connvert.c

对矩阵进行LU分解并求逆

//对N维方阵进行LU分解,并求逆 

# include <stdio.h>
# include <math.h>

# define TINY  1.0e-20
# define N   ....


void ludcmp(double a[N][N],int N, int indx[N], double d[N])
 {
   int i, imax, j, k;
   double  dum,  big, sum, temp;
   double vv[N];

   *d=1.0;
   
   for (i=0; i<N; i++){
	   big=0.0;
	   for (j=0; j<n; j++)		   
		   if ((temp = fabs (a[i][j])) > big) big = temp;        
//   if (big ==0.0) error("singular matrix in routine ludcmp");
       vv[i] = 1.0/big;
   }


   for (j=0; j<N; j++){
	   for (i=0; i<j; i++ ){
		   sum = a[i][j];
		   for (k=0; k<i; k++) sum -= a[i][k]*a[k][j];
		   a[i][j] = sum;
	   }

	   big=0.0;
//
	   for (i=j; i<n; i++){
		   sum=a[i][j];
//
		   for (k=0; k<j; k++)
			   sum -= a[i][k]*a[k][j];
		   a[i][j]=sum;
		   if ((dum=vv[i]*fabs(sum)) >= big) {
		   big=dum; 
		   imax=i;
		   }
	   }

	   if (j != imax) {
		   for (k=0; k<n; k++){
			   dum=a[imax][k];
			   a[imax][k]=a[j][k];
			   a[j][k]=dum;
		   }
		   *d = -(*d);
		   vv[imax]=vv[j];
	   }

	   indx[j]=imax;
	   if (a[j][j] == 0.0) a[j][j]= TINY;
	   if (j != n) {
		   dum=1.0/(a[j][j]);
		   for (i=j+1; i<N; i++) a[i][j] *= dum;
	   }
   }
 }
//
void lubksb (double a[N][N], int N, int indx[N], double b[N])
 {
	 int i, ip, j;
	 double sum;
//   
	 for ( i=0; i<N; i++){
		 ip=indx[i];
		 sum=b[ip];
		 b[ip]=b[i];
             for (j=0; j<=i-1; j++) sum -= a[i][j]*b[j];
		 b[i]=sum;
	 }
	 for ( i=N-1; i>=0; i--){
		 sum=b[i];
		 for ( j=i+1; j<N; j++) sum -= a[i][j]*b[j];
		 b[i]=sum/(a[i][i]);
	 }
 }

 
void main()
{ 
	double a[N][N], y[N][N], d[N], col[N];
	int i, j, indx[N];

	    ......
	ludcmp (a,N,indx,d);
	for (i=0; i<N; i++){
		printf("\n");
		for (j=0; j<N; j++){
			printf("%10.6f",a[i][j]);
		}
	}
printf("\n");
	for (j=0;j<N; j++){
		for (i=0; i<N; i++) col[i]=0.0;
		col[j]=1.0;
		lubksb(a,N,indx,col);
		for (i=0; i<N; i++) y[i][j]=col[i];
	}
	for (i=0; i<N; i++){
		printf("\n");
		for (j=0; j<N; j++){
			printf("%10.6f",y[i][j]);
		}
	}
	printf("\n");
}