VC++编程中队矩阵的各种操作的矩阵类

////////////////////////////////////////////////////////////////////////////// // // 完全矩阵类库 V2.11 // 具体实现 Matrix.cpp // BaconChina // 2002.11.18 // ////////////////////////////////////////////////////////////////////////////// // // 2002.11.18 duckli 修订了SVDResolve, 解决了原来分解奇异值小于0的问题 // // 2002.11.17 duckli 改写了LinearFormula(CBCMATRIX& A)和SVDResolve, 使 // 之适合于解决非方阵的线性方程求解和SVD分解问题 // ////////////////////////////////////////////////////////////////////////////// #include "stdafx.h" #include "BCMatrix.h" #include "iostream.h" #include "stdio.h" #include "memory.h" #include "stdlib.h" #include "math.h" CBCMATRIX::CBCMATRIX(int n) { int i; height=n; width=n; p=new double*[height]; for (i=0;i<height;i++) { p[i]=new double[width]; } } CBCMATRIX::CBCMATRIX(int _height,int _width) { int i; height=_height; width=_width; p=new double*[height]; for (i=0;i<height;i++) { p[i]=new double[width]; } } CBCMATRIX::CBCMATRIX(const CBCMATRIX& SRC) { int i,j; height=SRC.height; width=SRC.width; p=new double*[height]; for (i=0;i<height;i++) { p[i]=new double[width]; } for (i=0;i<height;i++) { for (j=0;j<width;j++) { p[i][j]=SRC.p[i][j]; } } } CBCMATRIX::~CBCMATRIX(void) { int i; for (i=0;i<height;i++) delete p[i]; delete p; } void CBCMATRIX::Display(void) { int i,j; for (i=0;i<height;i++) { for (j=0;j<width;j++) { cout<<p[i][j]<<' '; } cout<<endl; } cout<<"---------------------------------------"<<endl; } void CBCMATRIX::SetZero(void) { int i,j; for (i=0;i<height;i++) { for (j=0;j<width;j++) { p[i][j]=0; } } } void CBCMATRIX::SetUnit(void) { int i,j; for (i=0;i<height;i++) { for (j=0;j<width;j++) { if (i==j) p[i][j]=1; else p[i][j]=0; } } } void CBCMATRIX::Resize(int n) //重置为n*n阶方阵 { int i; for (i=0;i<height;i++) delete p[i]; delete p; height=n; width=n; p=new double*[height]; for (i=0;i<height;i++) { p[i]=new double[width]; } } void CBCMATRIX::Resize(int _height,int _width) //重置为h*w阶矩阵 { int i; for (i=0;i<height;i++) delete p[i]; delete p; height=_height; width=_width; p=new double*[height]; for (i=0;i<height;i++) { p[i]=new double[width]; } } void CBCMATRIX::SwapRow(int a,int b) { int i; double temp; for (i=0;i<width;i++) { temp=p[a][i]; p[a][i]=p[b][i]; p[b][i]=temp; } } void CBCMATRIX::SwapCol(int a,int b) { int i; double temp; for (i=0;i<height;i++) { temp=p[i][a]; p[i][a]=p[i][b]; p[i][b]=temp; } } void CBCMATRIX::PlusRow(int a,int b,double r) { int i; for (i=0;i<width;i++) { p[b][i]+=p[a][i]*r; } } void CBCMATRIX::PlusCol(int a,int b,double r) { int i; for (i=0;i<height;i++) { p[i][b]+=p[i][a]*r; } } void CBCMATRIX::MultRow(int a,double r) { int i; for (i=0;i<width;i++) { p[a][i]*=r; } } void CBCMATRIX::MultCol(int a,double r) { int i; for (i=0;i<height;i++) { p[i][a]*=r; } } double CBCMATRIX::EnValue(void) { int i,j,flag; double result,ration; CBCMATRIX temp; temp=*this; result=1; for (i=0;i<temp.height;i++) { if (temp.p[i][i]==0) { flag=0; for (j=i+1;j<temp.height;j++) { if (temp.p[j][i]!=0) { temp.SwapRow(i,j); result*=-1; flag=1; break; } } } if (flag==0) return 0; //保证对角线元素不为零 for (j=i+1;j<temp.height;j++) { ration=-temp.p[j][i]/temp.p[i][i]; temp.PlusRow(i,j,ration); } result*=temp.p[i][i]; } return result; } CBCMATRIX& CBCMATRIX::operator =(CBCMATRIX& A) { int i,j; Resize(A.height,A.width); for (i=0;i<height;i++) { for (j=0;j<width;j++) { p[i][j]=A.p[i][j]; } } return *this; } CBCMATRIX operator ~(CBCMATRIX& A) { int i,j; CBCMATRIX temp(A.width,A.height); for (i=0;i<temp.height;i++) { for (j=0;j<temp.width;j++) { temp.p[i][j]=A.p[j][i]; } } return temp; } CBCMATRIX operator +(CBCMATRIX& A,CBCMATRIX& B) { int i,j; CBCMATRIX temp(A.height,A.width); //A.width==B.width A.height==B.height for (i=0;i<temp.height;i++) { for (j=0;j<temp.width;j++) { temp.p[i][j]=A.p[i][j]+B.p[i][j]; } } return temp; } CBCMATRIX operator *(CBCMATRIX& A,double r) { int i,j; CBCMATRIX temp(A.height,A.width); for (i=0;i<temp.height;i++) { for (j=0;j<temp.width;j++) { temp.p[i][j]=r*A.p[i][j]; } } return temp; } CBCMATRIX operator *(double r,CBCMATRIX& A) { int i,j; CBCMATRIX temp(A.height,A.width); for (i=0;i<temp.height;i++) { for (j=0;j<temp.width;j++) { temp.p[i][j]=r*A.p[i][j]; } } return temp; } CBCMATRIX operator *(CBCMATRIX& A,CBCMATRIX& B) { int i,j,k; CBCMATRIX temp(A.height,B.width); temp.SetZero(); //A.width==B.height for (i=0;i<temp.height;i++) { for (j=0;j<temp.width;j++) { for (k=0;k<A.width;k++) { temp.p[i][j]+=A.p[i][k]*B.p[k][j]; } } } return temp; } CBCMATRIX operator >(CBCMATRIX& A,CBCMATRIX& B) { int i,j; CBCMATRIX temp(A.height,A.width+B.width); //A.height==B.height for (i=0;i<temp.height;i++) { for (j=0;j<A.width;j++) { temp.p[i][j]=A.p[i][j]; } } for (i=0;i<temp.height;i++) { for (j=0;j<B.width;j++) { temp.p[i][A.width+j]=B.p[i][j]; } } return temp; } CBCMATRIX operator <(CBCMATRIX& A,CBCMATRIX& B) { int i,j; CBCMATRIX temp(A.height+B.height,A.width); //A.width==B.width for (i=0;i<A.height;i++) { for (j=0;j<temp.width;j++) { temp.p[i][j]=A.p[i][j]; } } for (i=0;i<B.height;i++) { for (j=0;j<temp.width;j++) { temp.p[i+A.height][j]=B.p[i][j]; } } return temp; } CBCMATRIX operator !(CBCMATRIX& A) { int i,j; double ration; CBCMATRIX temp(A.height),tempA; temp.SetUnit(); //A.height==A.width tempA=A; for (i=0;i<temp.height;i++) { if (tempA.p[i][i]==0) { for (j=i+1;j<temp.height;j++) { if (tempA.p[j][i]!=0) { tempA.SwapRow(i,j); temp.SwapRow(i,j); break; } } } //保证对角线元素不为零 ration=1/tempA.p[i][i]; tempA.MultRow(i,ration); temp.MultRow(i,ration); for (j=i+1;j<temp.height;j++) { ration=-tempA.p[j][i]; tempA.PlusRow(i,j,ration); temp.PlusRow(i,j,ration); } } //现已变成对角线元素为1的上三角阵 for (i=0;i<temp.height;i++) { for (j=0;j<i;j++) { ration=-tempA.p[j][i]; tempA.PlusRow(i,j,ration); temp.PlusRow(i,j,ration); } } return temp; } void CBCMATRIX::QRResolve(CBCMATRIX& Q,CBCMATRIX& R) { CBCMATRIX P; int i,j; double g,a,b,v,c,s; Q.Resize(height); Q.SetUnit(); P.Resize(height); //height==width R=*this; for (i=0;i<R.height;i++) { for (j=i+1;j<R.height;j++) { P.SetUnit(); g=__max(fabs(R.p[i][i]),fabs(R.p[j][i])); if (g==0) break; a=R.p[i][i]/g; b=R.p[j][i]/g; v=sqrt(a*a+b*b); c=a/v; s=b/v; P.p[i][i]=c; P.p[j][j]=c; P.p[j][i]=-s; P.p[i][j]=s; R=P*R; Q=Q*~P; } } } double CBCMATRIX::MaxBelowDiag(void) { int i,j; double maxbelowdiag; maxbelowdiag=0; for (i=0;i<height;i++) { for (j=0;j<i;j++) { if (fabs(p[i][j])>maxbelowdiag) { maxbelowdiag=fabs(p[i][j]); } } } return(maxbelowdiag); } void CBCMATRIX::EigenAbove(CBCMATRIX& A) { CBCMATRIX Q,R; A=*this; while (A.MaxBelowDiag()>MINIMUM) { A.QRResolve(Q,R); A=R*Q; } A.Rectify(); } void CBCMATRIX::Rectify(void) //浮点数误差修正 { int i,j; for (i=0;i<height;i++) { for (j=0;j<width;j++) { if (fabs(p[i][j])<MINIMUM) p[i][j]=0; } } } void CBCMATRIX::LUResolve(CBCMATRIX& L,CBCMATRIX& U) //LU分解,L�