gausspivot.cpp

利用高斯列主元消去法进行LU分解

#include<iostream>
#include<math.h>
#include<stdio.h>
#define N 10
using namespace std;

void input_file(double a[N][N],double b[N],int *n) 
{
     int i,j; 
     char infile[20]; 

     FILE *input; 
     cout<<"input infile name:"; 
     cin>>infile; 
     input=fopen(infile,"r"); 
     if(input==NULL)  {
        cout<<"Can't open the file."<<endl; 
        exit(0);
     }
     fscanf(input,"%d",n);   
     for(i=0;i<*n;i++) {
        for(j=0;j<*n;j++) 
           fscanf(input,"%lf",&a[i][j]);  
        fscanf(input,"%lf",&b[i]);
     }
     fclose(input); 
     
} 
void input_keyboard(double a[N][N],double b[N],int *n) 
{
     int i,j;
     cout<<"请输入维数：";
     cin>>*n;
      cout<<"请输入增广矩阵:"<<endl;
      for(i=0;i<*n;i++){
         for(j=0;j<*n;j++)
            cin>>a[i][j];
         cin>>b[i];
      }
}
int gauss(double a[N][N],int q[N],int n)
{
     int i,j,max_row=0,m,k,t111,t11;
     double c,max_num;
     
     for(i=0;i<n-1;i++){ 
         max_num=fabs(a[i][i]);
         max_row=i;        
         for(j=i+1;j<n;j++){
            if(fabs(a[j][i])>max_num){
               max_num=fabs(a[j][i]);
               max_row=j;
            }                   
         }
         q[i]=max_row;
         if(max_row!=i)
            for(k=0;k<n;k++){
               c=a[i][k];
               a[i][k]=a[max_row][k];
               a[max_row][k]=c;
            }
         for(m=i+1;m<n;m++) {/*进行消元计算*/ 
            c=a[m][i]/a[i][i]; 
            for(j=i;j<n;j++) 
               a[m][j]=a[m][j]-c*a[i][j]; 
            a[m][i]=c;
         } 
      }
      q[i]=i;
}
void solve(double b[N],double LU[N][N],int q[N],int n)
{
     int i,j,k;
     double L[N][N]={0},U[N][N]={0},y[N],c;
          //LUx=b,Ly=b,Ux=y
     for(i=0;i<n;i++)
        for(j=0;j<n;j++){
           if(j==i)
              L[i][j]=1,U[i][j]=LU[i][j];
           else if(j>i)
              U[i][j]=LU[i][j];
           else
              L[i][j]=LU[i][j];
        }         
     for(i=0;i<n;i++)
        if(q[i]!=i){
           c=b[i];
           b[i]=b[q[i]];
           b[q[i]]=c;
        }
     for(i=0;i<n;i++) {/*回代过程*/ 
       y[i]=b[i]; 
       for(j=0;j<i;j++) 
          y[i]=y[i]-L[i][j]*y[j]; 
       //y[i]=b[i]; 
    } 
     for(k=n-1;k>=0;k--) {/*回代过程*/ 
       b[k]=y[k]; 
       for(j=n-1;j>k;j--) 
          b[k]=b[k]-U[k][j]*b[j]; 
       b[k]=b[k]/U[k][k]; 
    } 
}
main()
{
      int n,i,j,max,q[N]={0},choice;
      
      double A[N][N]={0},B[N]={0};
      cout<<"请选择输入方式：1.文件输入。2.键盘输入。";
      cin>>choice;
      if(choice==1)
         input_file(A,B,&n);
      else
         input_keyboard(A,B,&n);    
      gauss(A,q,n);
      solve(B,A,q,n);
      
      for(i=0;i<n;i++) {
         printf("\nx[%d]=%f\n",i+1,B[i]);
         printf("\n"); 
      }
      system("pause"); 
      return 0;
}