C++求矩阵的逆程序代码.doc

求解逆矩阵源代码：
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define TINY 1.0e-20
void inverse(double**,int);
void ludcmp(double**, int, int*, double*);
void lubksb(double**, int, int*, double*);
double **matrix(int,int,int,int);
double *vector(int,int);
void free_matrix(double**,int,int,int,int);
void free_vector(double*,int,int);
void inverse(double **mat, int dim)
{
int i,j,*indx;
double **y,d,*col;
y = matrix(0,dim-1,0,dim-1);
indx = (int *)malloc((unsigned)(dim*sizeof(int)));
col = vector(0,dim-1);
ludcmp(mat,dim,indx,&d);
for (j=0;j<dim;j++)
{
for (i=0;i<dim;i++) col[i] = 0.0;
col[j] = 1.0;
lubksb(mat,dim,indx,col);
for (i=0;i<dim;i++) y[i][j] = col[i];
}
for (i=0;i<dim;i++)
for (j=0;j<dim;j++)
mat[i][j] = y[i][j];
free_matrix(y,0,dim-1,0,dim-1);
free_vector(col,0,dim-1);
free(indx);
}
void ludcmp(double **a, int n, int *indx, double *d)
{
int i,imax,j,k;
double big,dum,sum,temp;
double *vv;
vv = (double*)malloc((unsigned)(n*sizeof(double)));
if (!vv)
{
fprintf(stderr,"Error Allocating Vector Memory\n");
exit(1);
}
*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)
{
fprintf(stderr,"Singular Matrix in Routine LUDCMP\n");
for (j=0;j<n;j++) printf(" %f ",a[i][j]); printf("/n");
exit(1);
}
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-1)
{
dum = 1.0 / a[j][j];
for (i=j+1;i<n;i++) a[i][j] *= dum;
}
}
free(vv);
}
void lubksb(double **a, int n, int *indx, double *b)
{
int i,ip,j,ii=-1;
double sum;

for (i=0;i<n;i++)
{
ip = indx[i];
sum = b[ip];
b[ip] = b[i];
if (ii>=0)
for (j=ii;j<i;j++) sum -= a[i][j] * b[j];
else if (sum) ii = i;
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];
}
}