#include<math.h>

#include<stdio.h>//////////////////////////////////////////////////////////////////////////////////////////

//矩陣結(jié)構(gòu)體

structMatrix

{

intm,n;//m為行數(shù),n為列數(shù)

double**pm;//指向矩陣二維數(shù)組的指針

};

//初始化矩陣mt，并置矩陣的行為m，列為n

voidInitMatrix(structMatrix*mt,intm,intn)

{

inti;

//指定矩陣的行和列

mt->m=m;

mt->n=n;

//為矩陣分配內(nèi)存

mt->pm=newdouble*[m];

for(i=0;i<m;i++)

{

mt->pm[i]=newdouble[n];

}

}

//用0初始化矩陣mt，并置矩陣的行為m，列為n

voidInitMatrix0(structMatrix*mt,intm,intn)

{

inti,j;

InitMatrix(mt,m,n);

for(i=0;i<m;i++)

for(j=0;j<n;j++)mt->pm[i][j]=0.0;

}

//銷毀矩陣mt

voidDestroyMatrix(structMatrix*mt)

{

inti;

//釋放矩陣內(nèi)存

for(i=0;i<mt->m;i++)

{

delete[]mt->pm[i];

}

delete[]mt->pm;

}

//矩陣相乘mt3=mt1*mt2

//成功返回1，失敗返回0

intMatrixMul(Matrix*mt1,Matrix*mt2,Matrix*mt3)

{

inti,j,k;

doubles;

//檢查行列號(hào)是否匹配

if(mt1->n!=mt2->m||mt1->m!=mt3->m||mt2->n!=mt3->n)return0;

for(i=0;i<mt1->m;i++)

for(j=0;j<mt2->n;j++)

{

s=0.0;

for(k=0;k<mt1->n;k++)s=s+mt1->pm[i][k]*mt2->pm[k][j];

mt3->pm[i][j]=s;

}

return1;

}

//矩陣轉(zhuǎn)置mt2=T(mt1)

//成功返回1，失敗返回0

intMatrixTran(Matrix*mt1,Matrix*mt2)

{

inti,j;

//檢查行列號(hào)是否匹配

if(mt1->m!=mt2->n||mt1->n!=mt2->m)return0;

for(i=0;i<mt1->m;i++)

for(j=0;j<mt1->n;j++)mt2->pm[j][i]=mt1->pm[i][j];

return1;

}

//控制臺(tái)顯示矩陣mt

voidConsoleShowMatrix(Matrix*mt)

{

inti,j;

for(i=0;i<mt->m;i++)

{

printf("\n");

for(j=0;j<mt->n;j++)printf("%2.4f",mt->pm[i][j]);

}

printf("\n");

}

//////////////////////////////////////////////////////////////////////////////////////////

//Guass列主元素消元法求解方程Ax=b,a=(A,b)

intGuassl(double**a,double*x,intn)

{

inti,j,k,numl,*h,t;

double*l,max;

l=newdouble[n];

h=newint[n];

for(i=0;i<n;i++)h[i]=i;//行標(biāo)

for(i=1;i<n;i++)

{

max=fabs(a[h[i-1]][i-1]);

numl=i-1;

//列元的最大值

for(j=i;j<n;j++)

{

if(fabs(a[h[j]][i-1])>max)

{

numl=h[j];

max=fabs(a[h[j]][i-1]);

}

}

if(max<0.00000000001)return0;

//交換行號(hào)

if(numl>i-1)

{

t=h[i];

h[i]=h[numl];

h[numl]=t;

}

for(j=i;j<n;j++)l[j]=a[h[j]][i-1]/a[h[i-1]][i-1];

for(j=i;j<n;j++)

for(k=i;k<n+1;k++)a[h[j]][k]=a[h[j]][k]-l[j]*a[h[i-1]][k];

}

for(i=n-1;i>=0;i--)

{

x[i]=a[h[i]][n];

for(j=i+1;j<n;j++)x[i]=x[i]-a[h[i]][j]*x[j];

x[i]=x[i]/a[h[i]][i];

}

//清除臨時(shí)數(shù)組內(nèi)存

delete[]l;

delete[]h;

return1;

}///////////////////////////////////////////////////////////////////////////////////////////

//最小二乘法求解矩陣Ax=b

intMinMul(MatrixA,Matrixb,double*x)

{

inti,j;

if(b.n!=1)return0;

if(A.m!=b.m)return0;

MatrixTranA;//定義A的轉(zhuǎn)置

InitMatrix0(&TranA,A.n,A.m);

MatrixTran(&A,&TranA);

MatrixTranA_A;//定義A的轉(zhuǎn)置與A的乘積矩陣

InitMatrix0(&TranA_A,A.n,A.n);

MatrixMul(&TranA,&A,&TranA_A);//A的轉(zhuǎn)置與A的乘積

MatrixTranA_b;//定義A的轉(zhuǎn)置與b的乘積矩陣

InitMatrix0(&TranA_b,A.n,1);

MatrixMul(&TranA,&b,&TranA_b);//A的轉(zhuǎn)置與b的乘積

DestroyMatrix(&TranA);//釋放A的轉(zhuǎn)置的內(nèi)存

MatrixTranA_A_b;//定義增廣矩陣

InitMatrix0(&TranA_A_b,TranA_A.m,TranA_A.m+1);

//增廣矩陣賦值

for(i=0;i<TranA_A_b.m;i++)

{

for(j=0;j<TranA_A_b.m;j++)TranA_A_b.pm[i][j]=TranA_A.pm[i][j];

TranA_A_b.pm[i][TranA_A_b.m]=TranA_b.pm[i][0];

}

DestroyMatrix(&TranA_A);

DestroyMatrix(&TranA_b);

//Guass列主消元法求解

Guassl(TranA_A_b.pm,x,TranA_A_b.m);

DestroyMatrix(&TranA_A_b);

return1;

}intMinMul(double**A,double*b,intm,intn,double*x)

{

intr,i;

MatrixAl,bl;

Al.pm=newdouble*