C++矩阵类(4)

}

for(long i=0;i

return result; }

bool Matrix::isPtv() {

assert(width==height);//是方阵才可以计算 bool result=true; Matrix m;

for(long i=1;i<=height;i++){ m=this->subMatrix(i); if(m.Deter()<=0){ result=false; break; } }

return result; }

Matrix Matrix::T() {

double * t=new double[width*height]; for(long i=0;i

*(t+height*j+i)=*(p+width*i+j); } }

Matrix m(t,width,height); delete [] t; return m; }

Matrix Matrix::Eigen() {

Matrix matrix=*this;

assert(matrix.getHeight()==matrix.getWidth());

double **a,*wr,*wi;

long nn,m,l,k,j,its,i,mmin;

double z,y,x,w,v,u,t,s,r,q,p,anorm;

long n=matrix.getWidth(); a=new double*[n+1]; wr=new double[n+1]; wi=new double[n+1]; for(i=0;i

a[i]=new double[n+1];

for(i=0;i

a[i+1][j+1]=matrix[i][j]; } }

anorm=fabs(a[1][1]); for (i=2;i<=n;i++) for (j=(i-1);j<=n;j++) anorm += fabs(a[i][j]); nn=n; t=0.0;

while (nn >= 1) { its=0; do {

for (l=nn;l>=2;l--) {

s=fabs(a[l-1][l-1])+fabs(a[l][l]); if (s == 0.0) s=anorm;

if (fabs(a[l][l-1]) + s == s) break; }

x=a[nn][nn]; if (l == nn) { wr[nn]=x+t; wi[nn--]=0.0; } else {

y=a[nn-1][nn-1];

w=a[nn][nn-1]*a[nn-1][nn]; if (l == (nn-1)) { p=0.5*(y-x); q=p*p+w;

z=sqrt(fabs(q)); x += t;

if (q >= 0.0) { z=p+SIGN(z,p);

wr[nn-1]=wr[nn]=x+z;

if (z) wr[nn]=x-w/z; wi[nn-1]=wi[nn]=0.0; } else {

wr[nn-1]=wr[nn]=x+p; wi[nn-1]= -(wi[nn]=z); }

nn -= 2; } else {

if (its == 30) printf(\ if (its == 10 || its == 20) { t += x;

for (i=1;i<=nn;i++) a[i][i] -= x;

s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]); y=x=0.75*s; w = -0.4375*s*s; } ++its;

for (m=(nn-2);m>=l;m--) { z=a[m][m]; r=x-z; s=y-z;

p=(r*s-w)/a[m+1][m]+a[m][m+1]; q=a[m+1][m+1]-z-r-s; r=a[m+2][m+1];

s=fabs(p)+fabs(q)+fabs(r); p /= s; q /= s; r /= s;

if (m == l) break;

u=fabs(a[m][m-1])*(fabs(q)+fabs(r));

v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1])); if (u+v == v) break; }

for (i=m+2;i<=nn;i++) { a[i][i-2]=0.0;

if (i != (m+2)) a[i][i-3]=0.0; }

for (k=m;k<=nn-1;k++) { if (k != m) { p=a[k][k-1]; q=a[k+1][k-1]; r=0.0;

if (k != (nn-1)) r=a[k+2][k-1]; if (x=fabs(p)+fabs(q)+fabs(r)) {

p /= x; q /= x; r /= x; } }

if (s=SIGN(sqrt(p*p+q*q+r*r),p)) { if (k == m) { if (l != m)

a[k][k-1] = -a[k][k-1]; } else

a[k][k-1] = -s*x; p += s; x=p/s; y=q/s; z=r/s; q /= p; r /= p;

for (j=k;j<=nn;j++) { p=a[k][j]+q*a[k+1][j]; if (k != (nn-1)) { p += r*a[k+2][j]; a[k+2][j] -= p*z; }

a[k+1][j] -= p*y; a[k][j] -= p*x; }

mmin = nn

a[i][k+1] -= p*q; a[i][k] -= p; } } } } }

} while (l < nn-1); }

Matrix eigenvalue(1,n); for(i=0;i

eigenvalue[0][i]=wr[i+1]; }

return eigenvalue; }

Matrix Matrix::Adjoint() {

return this->Inv()*this->Deter(); }

void Matrix::Unity() {

double total;

total=sqrt(sum(multiply(*this,*this))); for(long i=0;iheight;i++){ for(long j=0;jwidth;j++){ (*this)[i][j]/=total; } }

return; }

void Matrix::Store(const char * filename){ ofstream fout(filename); int w=5;

for(long i=0;i

fout.width(w); fout<

fout.width(15);

fout<<(*this)[i][j]<

cout.width(width); cout<

cout.width(width);

cout<<(*this)[i][j]<

fout.close(); }

Matrix Matrix::Inv() {

assert(width==height && width>0); Matrix m=*this;