用C++实现矩阵奇异值的求解

fmxsl

要求是,给出一个12*8的矩阵,然后进行求解.

fmxsl

#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][j]>scale))
         scale=fabs(a[k][j]);
     }
           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;



}