雅克比求逆

#includestdio.h#includestdlib.hdouble*MatrixOpp(double*A,intm,intn);/*矩阵求逆*/double*MatrixInver(double*A,intm,intn);/*矩阵转置*/doubleSurplus(doubleA[],intm,intn);/*求矩阵行列式*/double*MatrixOpp(doubleA[],intm,intn)/*矩阵求逆*/{inti,j,x,y,k;double*SP=NULL,*AB=NULL,*B=NULL,X,*C;SP=(double*)malloc(m*n*sizeof(double));AB=(double*)malloc(m*n*sizeof(double));B=(double*)malloc(m*n*sizeof(double));X=Surplus(A,m,n);X=1/X;for(i=0;im;i++)for(j=0;jn;j++){for(k=0;km*n;k++)B[k]=A[k];{for(x=0;xn;x++)B[i*n+x]=0;for(y=0;ym;y++)B[m*y+j]=0;B[i*n+j]=1;SP[i*n+j]=Surplus(B,m,n);AB[i*n+j]=X*SP[i*n+j];}}C=MatrixInver(AB,m,n);returnC;}double*MatrixInver(doubleA[],intm,intn)/*矩阵转置*/{inti,j;double*B=NULL;B=(double*)malloc(m*n*sizeof(double));for(i=0;in;i++)for(j=0;jm;j++)B[i*m+j]=A[j*n+i];returnB;}doubleSurplus(doubleA[],intm,intn)/*求矩阵行列式*/{inti,j,k,p,r;doubleX,temp=1,temp1=1,s=0,s1=0;if(n==2){for(i=0;im;i++)for(j=0;jn;j++)if((i+j)%2)temp1*=A[i*n+j];elsetemp*=A[i*n+j];X=temp-temp1;}else{for(k=0;kn;k++){for(i=0,j=k;im,jn;i++,j++)temp*=A[i*n+j];if(m-i){for(p=m-i,r=m-1;p0;p--,r--)temp*=A[r*n+p-1];}s+=temp;temp=1;}for(k=n-1;k=0;k--){for(i=0,j=k;im,j=0;i++,j--)temp1*=A[i*n+j];if(m-i){for(p=m-1,r=i;rm;p--,r++)temp1*=A[r*n+p];}s1+=temp1;temp1=1;}X=s-s1;}returnX;}/*Test*/intmain(){inti,j;doublearr[5][5],*result,*t=arr;for(i=0;i5;i++)for(j=0;j5;j++)scanf(%lf,&arr[i][j]);result=MatrixOpp((double*)arr,5,5);//求逆/*...其他操作，如显示结果*/printf(\n\nTheresultis:\n);for(i=0;i5*5;i++){printf(%lf\t,*(result+i));if(i%5==4)printf(\n);}free(result);system(PAUSE);return0;}自己写的矩阵类，这是头文件#ifndefXXXXX#defineXXXXX#defineMAX_SIZE4#includeiostream#includemath.husingnamespacestd;doubleabs(doublea);classCMatrix;classCVector;CVectoroperator+(CVector&a,CVector&b);//向量加法CVectoroperator-(CVector&a,CVector&b);//向量减法classCVector{public:CVector();CVector(int);CVector(int,double*b);int_size();voiddiapaly();CVector&operator=(CVector&a);friendCVectoroperator*(CMatrix&a,CVector&b);//矩阵乘以向量friendvoidEliminate(CMatrixm,CVector&b,CVector&solution);friendvoidJacobiIter(CMatrixm,CVectorb,CVector&solution);friendCVectoroperator+(CVector&a,CVector&b);//向量加法friendCVectoroperator-(CVector&a,CVector&b);//向量减法CVector&operator+=(constCVector&a);CVector&operator-=(constCVector&a);frienddoubleModule(CVector);//向量的模//~CVector();private:intn;//thesizedouble*value;};doubleModule(CVectora);//向量的模CMatrixoperator+(CMatrix&,CMatrix&);//重载矩阵加法CMatrixoperator-(CMatrix&a,CMatrix&b);//矩阵减法classCMatrix{public:CMatrix();CMatrix(int,int);CMatrix(int,int,doublea[][MAX_SIZE]);CMatrix(CMatrix&a);friendvoidEliminate(CMatrixm,CVector&b,CVector&solution);friendvoidJacobiIter(CMatrixm,CVectorb,CVector&solution);friendCMatrixoperator+(CMatrix&a,CMatrix&b);friendCMatrixInvMatrix(CMatrix);//矩阵求逆friendCMatrixoperator*(CMatrix&a,CMatrix&b);//矩阵乘法friendCMatrixoperator-(CMatrix&a,CMatrix&b);//矩阵减法friendCVectoroperator*(CMatrix&a,CVector&b);//矩阵乘以向量CMatrix&operator=(constCMatrix&);CMatrix&operator+=(constCMatrix&);CMatrix&operator-=(constCMatrix&);boolIsDiag();//判断是否为对角矩阵intRow();intCol();voiddispaly();intRank();//返回矩阵的秩doubleAbs();//求矩阵的行列式的值，并打印出一个上三角矩阵private:introw;//行数intcol;//列数doublevalue[MAX_SIZE][MAX_SIZE];//值};CMatrixInvMatrix(CMatrix);//矩阵求逆CMatrixoperator*(CMatrix&a,CMatrix&b);//矩阵乘法CVectoroperator*(CMatrix&a,CVector&b);//矩阵乘以向量CMatrixoperator-(CMatrix&a,CMatrix&b);//矩阵减法voidEliminate(CMatrixm,CVector&b,CVector&solution);//选主元消元法voidJacobiIter(CMatrixm,CVectorb,CVector&solution);//解方程组的Jacibi迭代法#endifCMatrixInvMatrix(CMatrixa){//逆if(a.Row()=0||a.Col()=0)throw矩阵输入错误!无法执行;if(a.Row()!=a.Col())throw非方阵,无法求逆!;CMatrixrhs(a.Row(),a.Col());inti,j,k,t;doubleax,temp,s;intflow;for(i=0;ia.Row();i++){//先单位化for(j=0;ja.Col();j++)rhs.value[i][j]=0;rhs.value[i][i]=1;}for(k=0;ka.Col();k++){//对于每一列s=0;for(j=k;ja.Row();j++){//选取最大的主元if(abs(a.value[j][k])=s){s=abs(a.value[j][k]);flow=j;}}for(t=0;ta.Col();t++){//交换k行和flow行temp=a.value[k][t];a.value[k][t]=a.value[flow][t];a.value[flow][t]=temp;temp=rhs.value[k][t];rhs.value[k][t]=rhs.value[flow][t];rhs.value[flow][t]=temp;}for(i=k+1;ia.Row();i++){//从k+1行开始消ax=a.value[i][k]/a.value[k][k];for(t=0;ta.Col();t++){a.value[i][t]=a.value[i][t]-ax*a.value[k][t];rhs.value[i][t]=rhs.value[i][t]-ax*rhs.value[k][t];}}}//********得到一个上三角阵，下面将之化为对角阵***********for(i=a.Row()-2;i=0;i--){for(j=i+1;ja.Col();j++){ax=a.value[i][j]/a.value[j][j];for(t=0;ta.Col();t++){a.value[i][t]=a.value[i][t]-ax*a.value[j][t];rhs.value[i][t]=rhs.value[i][t]-ax*rhs.value[j][t];}}}for(i=0;ia.Row();i++){for(j=0;ja.Col();j++){rhs.value[i][j]/=a.value[i][i];}for(j=0;ja.Col();j++){a.value[i][j]/=a.value[i][i];}}returnrhs;}