#include "iostream.h"
#include"fstream.h"
#include "math.h"

#define N 100

void Jacobi()
{
  int i,j,k,g;
  double a[N][N],b[N],x[N],X[N];
  int n;
  cout<<"輸入行列式的階數:";
  cin>>n;
  ifstream istrm("data.txt");
  for(i=0;i<n;i++)
  {
   for(j=0;j<n;j++)
    istrm>>a[i][j];
    istrm>>b[i];
  }
 
  for(i=0;i<n;i++)
  {
   b[i]/=a[i][i];
   for(j=0;j<n;j++)
   {  if(j==i) continue;
    a[i][j]/=a[i][i];
   }
   a[i][i]=0;
  }
  cout<<"設置初始向量:";
  for(i=0;i<n;i++)
   cin>>x[i];
  int m;
  cout<<"輸入跌代的次數:";
   cin>>m;
  for(g=0;g<m;g++)
  { 
   for(i=0;i<n;i++)
   {
      double sum=0;
      for(j=0;j<n;j++)
   {
    if(j==i) continue;
    sum+=a[i][j]*x[j];
   }
   X[i]=b[i]-sum;
 
   }
   cout<<g<<" ";
   for(k=0;k<n;k++)
   {
    x[k]=X[k];
         cout<<X[k]<<" ";
   }
   cout<<endl;
 
  }
}
void GaussSeidel()
{
  int i,j,k;
  double a[N][N],b[N],x[N];
  int n;
  cout<<"輸入行列式的階數:";
  cin>>n;
  ifstream istrm("data.txt");
  for(i=0;i<n;i++)
  {
   for(j=0;j<n;j++)
    istrm>>a[i][j];
    istrm>>b[i];
  }
  for(i=0;i<n;i++)
  {
   b[i]/=a[i][i];
   for(j=0;j<n;j++)
   {  if(j==i) continue;
    a[i][j]/=a[i][i];
   }
   a[i][i]=0;
  }
  cout<<"設置初始向量";
  for(i=0;i<n;i++)
   cin>>x[i];
  int m;
  cout<<"跌代的次數:";
   cin>>m;
  for(int g=0;g<m;g++)
  {
     for(i=0;i<n;i++)
   {
      double sum=0;
      for(j=0;j<n;j++)
   {
    if(j==i) continue;
    else  sum+=a[i][j]*x[j];
   
   }
   x[i]=b[i]-sum;
 
   }
  cout<<g<<" ";
   for(k=0;k<n;k++)
   {
   
         cout<<x[k]<<" ";
   }
   cout<<endl;
 
  }
}
void SOR()
{
  int i,j,k,cnt=0;
  double a[N][N],b[N],x[N],X[N];
  int n;
  cout<<"輸入行列式的階數:";
  cin>>n;
  ifstream istrm("data.txt");
  for(i=0;i<n;i++)
  {
   for(j=0;j<n;j++)
    istrm>>a[i][j];
    istrm>>b[i];
  }
  for(i=0;i<n;i++)
  {
   b[i]/=a[i][i];
   for(j=0;j<n;j++)
   {  if(j==i) continue;
    a[i][j]/=a[i][i];
   }
   a[i][i]=0;
  }
  cout<<"設置初始向量";
  for(i=0;i<n;i++)
   cin>>x[i];
  double e,m=1;
  cout<<"輸入允許的誤差:";
  cin>>e;
  int w;
  cout<<"輸入松弛因子:";
  cin>>w;
 
  while(m>e)
  {
  cnt++;
     for(i=0;i<n;i++)
   {
         double sum=0;
         for(j=0;j<n;j++)
   {
     if(j==i) continue;
     if(j>i)  sum+=a[i][j]*x[j];
     if(j<i)  sum+=a[i][j]*X[j];
   
   }
      X[i]=x[i]+w*(b[i]-sum-x[i]);
 
  }cout<<cnt<<" ";
     for(k=0,m=0;k<n;k++)
   {
    m+=(X[k]-x[k])*(X[k]-x[k]);
    x[k]=X[k];
         cout<<x[k]<<" ";

   }
     m=sqrt(m);
     cout<<endl;
  }
  cout<<"跌代的次數為"<<cnt<<endl;
}
void Gauss()
{
  int i,j,k;
  float a[N][N],b[N],x[N];
  int n;
  cout<<"輸入行列式的階數:";
  cin>>n;
  ifstream istrm("data.txt");
  for(i=0;i<n;i++)
  {
   for(j=0;j<n;j++)
  istrm>>a[i][j];
   istrm>>b[i];
  }
  for(i=1;i<n;i++)
   for(j=i;j<n;j++)
   { float q=a[j][i-1];
     for(k=i-1;k<n;k++) a[j][k]+=-q*a[i-1][k]/a[i-1][i-1];
  b[j]+=-q*b[i-1]/a[i-1][i-1];
   }
   x[n-1]=b[n-1]/a[n-1][n-1];
   for(i=n-2;i>=0;i--)
   {
     float sum=0;
  for(j=n-1;j>i;j--) sum+=x[j]*a[i][j];
        x[i]=(b[i]-sum)/a[i][i];
   }
  for(i=0;i<n;i++)
    cout<<x[i]<<" ";
}
void SanJiao()
{
  int i,j,r;
  float a[N][N],b[N],l[N][N],u[N][N];
  int n;
  cout<<"輸入行列式的階數:";
  cin>>n;
  ifstream istrm("data.txt");
  for(i=0;i<n;i++)
  {
   for(j=0;j<n;j++)
  istrm>>a[i][j];
   istrm>>b[i];
  }
 for(i=0;i<n;i++)
  for(j=0;j<n;j++)
  {
    if(i<j) l[i][j]=0;
    if(i==j) l[i][j]=1;
    else u[i][j]=0;
  }
  u[0][0]=a[0][0];
   for(i=1;i<n;i++)
  {
    u[0][i]=a[0][i];
 l[i][0]=a[i][0]/u[0][0];
  }
  for(r=1;r<n;r++)
  {
    for(i=r;i<n;i++)
 {
   u[r][i]=a[r][i];
   for(j=0;j<=r-1;j++)
   {  
    u[r][i]-=l[r][j]*u[j][i];
   
   }

 }
 for(i=r+1;i<n;i++)
 {
   l[i][r]=a[i][r];
   for(j=0;j<=r-1;j++)
   {l[i][r]-=l[i][j]*u[j][r];}
   l[i][r]/=u[r][r];
 }
  }
  for(i=0;i<n;i++)
 {  cout<<"|";
  for(j=0;j<n;j++)
   {
    cout<<l[i][j]<<" ";
   }
  cout<<"|";
  cout<<" "<<"|";
  for(j=0;j<n;j++)
   {
     cout<<u[i][j]<<" ";
   }
    cout<<"|"<<endl;
 }
}
void Sqrt()
{
    int i,j,k;
  double temp;
  double a[10][10],l[10][10]={0};
  int n;
  cout<<"輸入行列式的階數:";
  cin>>n;
  ifstream istrm("data.txt");
  for(i=1;i<=n;i++)
    for(j=1;j<=n;j++)
  istrm>>a[i][j];
  for(j=1;j<=n;j++)
  {temp=0;
   for(k=1;k<j;k++) temp+=(l[j][k]*l[j][k]);
    l[j][j]=sqrt(a[j][j]-temp);
 for(i=j+1;i<=n;i++)
 {
  temp=0;
   for(k=1;k<j;k++) temp+=(l[j][k]*l[i][k]);
   l[i][j]=(a[i][j]-temp)/l[j][j];
  
 }
  }
  for(i=1;i<=n;i++)
 {  cout<<"|";
  for(j=1;j<=n;j++)
   {
    cout<<l[i][j]<<" ";
   }
  cout<<"|"<<endl;
  
 
  }
}
void main()
{
  int n,m;
  cout<<"***********提示菜單***************"<<endl;
  cout<<"*******1-線性方程的直接解法*******"<<endl;
  cout<<"*******2-線性方程的間接解法*******"<<endl;
  cout<<"*******0-退出程序*****************"<<endl;
  cout<<"**********************************"<<endl;
  while(1)
  {
    cout<<"輸入你要進行的操作:";
 cin>>n;
 if(n==0) break;
 switch(n)
 {
 case 1:cout<<"***********************"<<endl;
     cout<<"***1-高斯消去法********"<<endl;
     cout<<"***2-矩陣三角分解法****"<<endl;
     cout<<"***3-平方根分解法******"<<endl;
     cout<<"***0-返回上一層********"<<endl;
     cout<<"***********************"<<endl;
     cin>>m;
     switch(m)
     {
     case 1:Gauss();break;
     case 2:SanJiao();break;
     case 3:Sqrt();break;
     case 0:break;
     default:cout<<"輸入錯誤請重新輸入"<<endl;      
     }
     break;
 case 2:cout<<"***********************"<<endl;
     cout<<"***1-雅可比迭代法******"<<endl;
     cout<<"***2-高斯-塞德爾迭代法*"<<endl;
     cout<<"***3-松弛迭代法********"<<endl;
     cout<<"***0-返回上一層********"<<endl;
     cout<<"***********************"<<endl;
     cin>>m;
     switch(m)
     {
     case 1:Jacobi();break;
     case 2:GaussSeidel();break;
     case 3:SOR();break;
           case 0:break;
     default:cout<<"輸入錯誤請重新輸入"<<endl;
     }break;
 default:cout<<"輸入錯誤請重新輸入"<<endl;
 }
  }

}