divh.h

简单的矩阵求逆的c语言算法

#include <conio.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 4
double A[N][N]={{1,3,2,13},{7,2,1,-2},{9,15,3,-2},{-2,-2,11,5}};
inline void swap(double &a,double &b){double c=a;a=b;b=c;};
int DinV(double A[N][N],int n)
{
 int i,j,k;
 double d;
 int JS[N],IS[N];
 for (k=0;k<n;k++)
 {
  d=0;
  for (i=k;i<n;i++)
   for (j=k;j<n;j++)
  {
   if (fabs(A[i][j])>d)
   {
    d=fabs(A[i][j]);
    IS[k]=i;
    JS[k]=j;
   }
  }
  if (d+1.0==1.0) return 0;
  if (IS[k]!=k)
   for (j=0;j<n;j++)
    swap(A[k][j],A[IS[k]][j]);
  if (JS[k]!=k)
   for (i=0;i<n;i++)
    swap(A[i][k],A[i][JS[k]]);
  A[k][k]=1/A[k][k];
  for (j=0;j<n;j++)
   if (j!=k) A[k][j]=A[k][j]*A[k][k];
  for (i=0;i<n;i++)
   if (i!=k)
    for (j=0;j<n;j++)
     if (j!=k) A[i][j]=A[i][j]-A[i][k]*A[k][j];
  for (i=0;i<n;i++)
   if (i!=k) A[i][k]=-A[i][k]*A[k][k];
 }
 for (k=n-1;k>=0;k--)
 {
  for (j=0;j<n;j++)
   if (JS[k]!=k) swap(A[k][j],A[JS[k]][j]);
  for (i=0;i<n;i++)
   if (IS[k]!=k) swap(A[i][k],A[i][IS[k]]);
 }
 for (i=0;i<n;i++)
 {
  for (j=0;j<n;j++)
   printf("  %1.4f",A[i][j]);
  puts("");
 }
}

void gauss(double a[][N+1],double uk[N+1],double n) 
{ 
	int i,j,k; 
	double temp,s,l; 
	double yk[N+1]={0},
	for (i=1;i<N+1;i++)
	{
		yk[i]=1;
	}

	for(k=1;k<n+1;k++) 
	{ 
		//选列主元 
		i=k; 
		for(j=k+1;j<=n;j++) 
		{ if(fabs(a[j][k])>fabs(a[i][k])) 
		i=j; 
		} 
		
		//换行 
		if(i!=k) 
		{	for(j=k;j<=n;j++) 
		{ 
			temp=a[k][j]; 
			a[k][j]=a[i][j]; 
			a[i][j]=temp; 
			
		} 
		temp=yk[k];
		yk[k]=yk[i];
		yk[i]=temp;
		}
		//消元 
		for(i=k+1;i<=n;i++) 
		{ 
			l=1.0*(a[i][k])/(a[k][k]); 
			for(j=1;j<=n;j++) 
			{a[i][j]=a[i][j]-(a[k][j])*l;}
			yk[i]=yk[i]-yk[k]*l;
		} 
        
		
	} 
	//回代 
	uk[n]=yk[n]/a[n][n]; 
	for(k=n-1;k>=1;k--) 
	{ 
		s=0; 
		for(j=k;j<=n;j++) 
		{ 
			if(j==k) 
                continue;
			s+=a[k][j]*uk[j]; 
			
		} 
		uk[k]=(yk[k]-s)/a[k][k]; 
	} 
	
}