#include <time.h>
#include <math.h>
#include <stdlib.h>
#include <iostream>
using namespace std;
void swap(double &a,double &b)
{
double temp=a;
a=b;
b=temp;
}
/**********************************************
*函數名：InverseMatrix
*函數介紹：求逆矩陣（高斯—約當法）
*輸入參數：矩陣首地址（二維數組）matrix，階數row
*輸出參數：matrix原矩陣的逆矩陣
*返回值：成功，0；失敗，1
*調用函數：swap(double &a,double &b)
*作者：
*完成時間：2009-10-04
**********************************************/
int InverseMatrix(double *matrix,const int &row)
{
double *m=new double[row*row];
double *ptemp,*pt=m;
int i,j;
ptemp=matrix;
for (i=0;i<row;i++)
{
for (j=0;j<row;j++)
{
*pt=*ptemp;
ptemp++;
pt++;
}
}
int k;
int *is=new int[row],*js=new int[row];
for (k=0;k<row;k++)
{
double max=0;
//全選主元
//尋找最大元素
for (i=k;i<row;i++)
{
for (j=k;j<row;j++)
{
if (fabs(*(m+i*row+j))>max)
{
max=*(m+i*row+j);
is[k]=i;
js[k]=j;
}
}
}
if (0 == max)
{
return 1;
}
//行交換
if (is[k]!=k)
{
for (i=0;i<row;i++)
{
swap(*(m+k*row+i),*(m+is[k]*row+i));
}
}
//列交換
if (js[k]!=k)
{
for (i=0;i<row;i++)
{
swap(*(m+i*row+k),*(m+i*row+js[k]));
}
}
*(m+k*row+k)=1/(*(m+k*row+k));
for (j=0;j<row;j++)
{
if (j!=k)
{
*(m+k*row+j)*=*((m+k*row+k));
}
}
for (i=0;i<row;i++)
{
if (i!=k)
{
for (j=0;j<row;j++)
{
if(j!=k)
{
*(m+i*row+j)-=*(m+i*row+k)**(m+k*row+j);
}
}
}
}
for (i=0;i<row;i++)
{
if(i!=k)
{
*(m+i*row+k)*=-(*(m+k*row+k));
}
}
}
int r;
//恢復
for (r=row-1;r>=0;r--)
{
if (js[r]!=r)
{
for (j=0;j<row;j++)
{
swap(*(m+r*row+j),*(m+js[r]*row+j));
}
}
if (is[r]!=r)
{
for (i=0;i<row;i++)
{
swap(*(m+i*row+r),*(m+i*row+is[r]));
}
}
}
ptemp=matrix;
pt=m;
for (i=0;i<row;i++)
{
for (j=0;j<row;j++)
{
*ptemp=*pt;
ptemp++;
pt++;
}
}
delete []is;
delete []js;
delete []m;
return 0;
}
int main()
{
double a[6*6];
int i;
for(i=0;i<36;i++)
{
a[i] = rand()%1000+1;
}
for(i=0;i<36;i++)
{
cout<<a[i]<<" ";
if ((i+1)%6 == 0)
{
cout<<endl;
}
}
cout<<endl;
InverseMatrix(a,6);
for(i=0;i<36;i++)
{
cout<<a[i]<<" ";
if ((i+1)%6 == 0)
{
cout<<endl;
}
}
cout<<endl;
system("pause");
return 0;
}