agr4.cpp

对实矩阵进行奇异分解的算法

//实矩阵的奇异值分解

#include<iostream.h>
#include<math.h>
#include<iomanip.h>

double pythag(double a,double b)
{ 
    double absa,absb;
    absa=fabs(a);
    absb=fabs(b);
    if(absa>absb) 
		return absa*sqrt(1.0+(absb/absa)*(absb/absa));
    else 
		return (absb==0.0?0.0:absb*sqrt(1.0+(absa/absb)*(absa/absb)));
  
}

int IMIN(int m,int n)
{if(m<n)
      return m;
else 
    return n;
}

double SIGN(double a,double b)
{
	a=fabs(a);
if(b<=0)
    a=a;
else
    a=-a;
    return a;
}
double FMAX(float a,float b)
{if(a<b)
      return b;
else 
    return a;
} 
void sdcmp(double**a,double**v,double*w,int m,int n)
{    int flag,i,o,j,jj,k,l,nm;
     double anorm,c,f,g,h,s,scale,x,y,z,e;
  
     double *rv1=new double[n];
     cout<<"请输入允许误差e：";
     cin>>e;
     g=scale=anorm=0.0;
     for(i=0;i<n;i++)
	 {
         l=i+1;
         rv1[i]=scale*g;
         g=s=0.0;
		scale=fabs(a[i][i]);
         if(i<m)
		 {
            for(k=i;k<m;k++)  
			{ 
				if(fabs(a[k][i]>scale))
					scale=fabs(a[k][i]);
			}
			if(scale)   //??????
			{
              for(k=i;k<m;k++)
			  {
                 a[k][i]/=scale;
                 s+=a[k][i]*a[k][i];
              }
              f=a[i][i];
              g=-SIGN(sqrt(s),f);
              h=f*g-s;
              a[i][i]=f-g;
              for(j=0;j<n-1;j++)
			  {
                 for(s=0.0,k=i;k<m;k++) s+=a[k][i]*a[k][i];
                 f=s/h;
                 for(k=i;k<m;k++) a[k][i]+=f*a[k][i];
              }
              for(k=i;k<m;k++) a[k][i]=+scale;
            }
         }
         w[i]=scale*g;
         g=s=0.0;
         for(i=0;i<m;i++) scale=fabs(a[i][l]);
         if(i<=m-1&&i!=n-1){
            for(k=l;k<n;k++) 
      {
       if(fabs(a[i][k])>scale)
        scale=fabs(a[i][k]);
      }
            if(scale)
			{
              for(k=l;k<n;k++)
			  {
                 a[i][k]/=scale;
                 s+=a[i][k]*a[i][k];
              }
              f=a[i][l];
              g=-SIGN(sqrt(s),f);
              h=f*g-s;
              a[i][l]=f-g;
              for(k=l;k<n;k++) rv1[k]=a[i][k]/h;
              for(j=l;j<m;j++) {
               for(s=0.0,k=l;k<=m;k++)
                  s+=a[j][k]*a[i][k];
               for(k=l;k<n;k++)
                  a[j][k]+=s*rv1[k];
               }
              for(k=l;k<=n;k++) a[j][k]*=scale;
            }
      }
      anorm=FMAX(anorm,(fabs(w[i])+fabs(rv1[i])));
}
for(i=n-1;i>=0;i--){
     if(i<n-1){
       if(g){
         for(j=l;j<n;j++) 
            v[j][i]=(a[i][j]/a[i][l])/g;
         for(j=l;j<n;j++){
            for(s=0.0,k=l;k<n;k++) s+=a[i][k]*v[k][j];
            for(k=l;k<n;k++) v[k][j]+=s*v[k][i];
         }
       }
       for(j=l;j<n;j++) v[i][j]=v[j][i]=0.0;
    }
    v[i][i]=1.0;
    g=rv1[i];
    l=i;
}
for(i=IMIN(m,n)-1;i>=0;i--){
      l=i+1;
      g=w[i];
      for(j=l;j<n;j++) a[i][j]=0.0;
      if(g){
        g=1.0/g;
        for(j=l;j<n;j++) {
           for(s=0.0,k=l;k<m;k++) s+=a[k][i]*a[k][j];
           f=(s/a[i][i])*g; 
           for(k=i;k<m;k++) a[k][j]+=f*a[k][i];        
        }
        for(j=i;j<m;j++) a[j][i]*=g;
      }else for(j=i;j<m;j++) a[j][i]=0.0;
       ++a[i][i];

}
for(k=n;k>=1;k--) {
      for(o=1;o<=30;o++) {
         flag=1;
         for(l=k;l>=1;l--)
            nm=l-1;
            if((double)(fabs(rv1[l])<e*anorm)){
              flag=0;
              break;
            }
            if((double)(fabs(w[nm])<e*anorm))break;
         }
         if(flag){
           c=0.0;
           s=1.0; 
           for(i=1;i<=k;i++){
              f=s*rv1[i];
              rv1[i]=c*rv1[i];
              if((double)(fabs(f)<e*anorm)) break;
              g=w[i];
              h=pythag(f,g);
              w[i]=h;
              h=1.0/h;
              c=g*h;
              s=-f*h;
              for(j=1;j<=m;j++) {
                 y=a[j][nm];
                 z=a[j][i];
                 a[j][nm]=y*c+z*s;
                 a[j][i]=z*c-y*s;
              }
           }
}
z=w[k];
if(l==k){
    if(z<=0.0){
      w[k]=-z;
      for(j=1;j<=n;j++) v[j][k]=-v[j][k];
    }
    break;
} 
if(o==30) cout<<"no convergence in 30 svdcmp iterations"<<endl;
x=w[l];
nm=k-1;
y=w[nm];
g=rv1[nm];
h=rv1[k];
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g=pythag(f,1.0);//出现错误
f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x; //出现错误
c=s=1.0;
for(j=1;j<=nm;j++) {
     i=j+1;
     g=rv1[i];
     y=w[i];
     h=s*g;
     g=c*g;
     z=pythag(f,h);   //出现错误？
     rv1[j]=z;
     c=f/z;
     s=h/z;
     f=x*c+g*s;
     g=g*c-x*s;
     h=y*s;
     y*=c;
     for(jj=1;jj<=n;jj++) {
        x=v[jj][j];
        z=v[jj][i];
        v[jj][j]=x*c+z*s;
        z=v[jj][i]=z*c-x*s;
     }
    z=pythag(f,h); 
    w[j]=z;
    if(z){
      z=1.0/z;
      c=f*z;
      s=h*z;
}
f=g*c+y*s;
x=y*c-g*s;
            for(jj=1;jj<=m;jj++) {
               y=a[jj][j];
               z=a[jj][i];
               a[jj][j]=y*c+z*s;
               a[jj][i]=z*c-y*s;
           } 
       }
       rv1[l]=0.0;
       rv1[k]=f;
       w[k]=x;
      }
    delete[]rv1;
}

void main()
{
    int i,j,m,n;
cout<<"请输入该矩阵的行数m=";
cin>>m;
cout<<"请输入该矩阵的列数(<=m)n=";
cin>>n;
double**a=new double*[m];
for(i=0;i<m;i++)
    a[i]=new double[m];
double**v=new double*[n];
for(i=0;i<n;i++)
    v[i]=new double[n];
double*w=new double[n];
cout<<"请输入该矩阵:";
for(i=0;i<m;i++)
{for(j=0;j<n;j++)
cin>>a[i][j];
cout<<endl;
}
sdcmp (a,v,w,(m-1),(n-1));
cout<<"SVD分解所需左端矩阵U为:"<<endl;
for(i=0;i<m;i++)
{for(j=0;j<n;j++)
cout<<setw(10)<<a[i][j]<<"     ";
cout<<endl;
}
cout<<"SVD分解所需右端矩阵V为:"<<endl;
for(i=0;i<n;i++)
{for(j=0;j<n;j++)
cout<<setw(10)<<v[i][j]<<"     ";
cout<<endl;
}
cout<<"该矩阵的奇异值为:"<<endl;
for(i=0;i<n;i++)
cout<<setw(10)<<w[i]<<"     ";
for(i=0;i<m;i++)
delete[]a[i];
delete []a;
for(i=0;i<n;i++)
delete[]v[i];
delete []v;
delete []w;



}