matrixl.cpp

此源码为矩阵的求逆算法

// matrixl.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include "math.h"
#include <stdlib.h> 
//#include <math.h> 
//#include <stdio.h> 
int brinv(double a[4][4], int n) 
{ int *is,*js,i,j,k,l,u,v; 
double d,p; 
is=new int[n]; 
js=new int[n];
for (k=0; k<=n-1; k++) 
{ d=0.0; 
for (i=k; i<=n-1; i++) 
for (j=k; j<=n-1; j++) 
{ l=i*n+j; p=fabs(a[i][j]); 
if (p>d) { d=p; is[k]=i; js[k]=j;} 
} 
if (d+1.0==1.0) 
{ delete is; delete js; printf("err**not inv\n"); 
return(0); 
} 
if (is[k]!=k) 
for (j=0; j<=n-1; j++) 
{ u=k*n+j; v=is[k]*n+j; 
p=a[k][j]; a[k][j]=a[is[k]][j]; a[is[k]][j]=p; 
} 
if (js[k]!=k) 
for (i=0; i<=n-1; i++) 
{ u=i*n+k; v=i*n+js[k]; 
p=a[i][k]; a[i][k]=a[i][js[k]]; a[i][js[k]]=p; 
} 
l=k*n+k; 
a[k][k]=1.0/a[k][k]; 
for (j=0; j<=n-1; j++) 
if (j!=k) 
{ u=k*n+j; a[k][j]=a[k][j]*a[k][k];} 
for (i=0; i<=n-1; i++) 
if (i!=k) 
for (j=0; j<=n-1; j++) 
if (j!=k) 
{ u=i*n+j; 
a[i][j]=a[i][j]-a[i][k]*a[k][j]; 
} 
for (i=0; i<=n-1; i++) 
if (i!=k) 
{ u=i*n+k; a[i][k]=-a[i][k]*a[k][k];} 
} 
for (k=n-1; k>=0; k--) 
{ if (js[k]!=k) 
for (j=0; j<=n-1; j++) 
{ u=k*n+j; v=js[k]*n+j; 
p=a[k][j]; a[k][j]=a[js[k]][j]; a[js[k]][j]=p; 
} 
if (is[k]!=k) 
for (i=0; i<=n-1; i++) 
{ u=i*n+k; v=i*n+is[k]; 
p=a[i][k]; a[i][k]=a[i][is[k]]; a[i][is[k]]=p; 
} 
} 
delete is; delete js; 
return(1); 
} 

void brmul(double a[4][4], double b[4][4],int m,int n,int k,double c[4][4]) 
{ int i,j,l,u; 
for (i=0; i<=m-1; i++) 
for (j=0; j<=k-1; j++) 
{ u=i*k+j; c[i][j]=0.0; 
for (l=0; l<=n-1; l++) 
c[i][j]=c[i][j]+a[i][l]*b[l][j]; 
} 
return; 
} 

int main() 
{ int i,j; 
static double a[4][4]={ {0.2368,0.2471,0.2568,1.2671}, 
{1.1161,0.1254,0.1397,0.1490}, 
{0.1582,1.1675,0.1768,0.1871}, 
{0.1968,0.2071,1.2168,0.2271}}; 
static double b[4][4],c[4][4]; 
for (i=0; i<=3; i++) 
for (j=0; j<=3; j++) 
b[i][j]=a[i][j]; 
i=brinv(a,4); 
if (i!=0) 
{ printf("MAT A IS:\n"); 
for (i=0; i<=3; i++) 
{ for (j=0; j<=3; j++) 
printf("%13.7e ",b[i][j]); 
printf("\n"); 
} 
printf("\n"); 
printf("MAT A- IS:\n"); 
for (i=0; i<=3; i++) 
{ for (j=0; j<=3; j++) 
printf("%13.7e ",a[i][j]); 
printf("\n"); 
} 
printf("\n"); 
printf("MAT AA- IS:\n"); 
brmul(b,a,4,4,4,c); 
for (i=0; i<=3; i++) 
{ for (j=0; j<=3; j++) 
printf("%13.7e ",c[i][j]); 
printf("\n"); 
} 
} 
}