ssvdc.c
Go to the documentation of this file.
00001 /*  ssvdc is a function to compute the singular value decomposition
00002     of a matrix.
00003 
00004     Ver.1.00, Jun,3,1988.  
00005     Ver.1.01, Jul.14,1988.  remove the side effect to the matrix X   
00006     Ver.1.02, Jul.15,1988.  debug the matrix base shift routine     */
00007 
00008 #include "arith.h"
00009 #define MAXIT 30  /* the maximum number of iteratons */
00010 
00011 ssvdc(xarg,ldx,n,p,s,e,u,ldu,v,ldv,work,job)
00012 int n,p,ldx,ldu,ldv,job;
00013 MATRIX xarg,s,e,u,v,work;
00014 {
00015     int i,iter,j,jobu,k,kase,kk,l,ll,lls,lm1,lp1,ls,lu,m,maxit,info;
00016     int mm,mm1,mp1,nct,nctp1,ncu,nrt,nrtp1;
00017     int ii,jj;
00018 
00019     int max0(),min0();
00020 
00021     REAL t,r;
00022     REAL b,c,cs,el,emm1,f,g,scale,shift,sl,sm,sn,smm1,t1,test;
00023     REAL ztest;
00024     MATRIX x;
00025 
00026     REAL amax1(),sdot(),snrm2();
00027 
00028     LOGICAL wantu,wantv;
00029 
00030     /* shift the matrix base from 0 to 1 */
00031 
00032     for(ii=n-1; ii>=0; ii--)
00033         for(jj=p-1; jj>=0; jj--)
00034             x[ii+1][jj+1]=xarg[ii][jj];
00035 
00036     /*  set the maximum number of iterations */
00037 
00038     maxit=MAXIT;
00039 
00040     /* determine what is to be computed */
00041 
00042     wantu=FALSE;
00043     wantv=FALSE;
00044     jobu=(job%100)/10;
00045     ncu=n;
00046 
00047     if(jobu > 1) ncu=min0(2,n,p);
00048     if(jobu != 0) wantu=TRUE;
00049     if((job%10) != 0) wantv=TRUE;
00050 
00051     /* reduce x to bidiagonal form,
00052        storing the diagonal elements in s and
00053        the super-diagonal elements in e */
00054 
00055     info=0;
00056     nct=min0(2,n-1,p);
00057     nrt=max0(2,0,min0(2,p-2,n));
00058 
00059     lu=max0(2,nct,nrt);
00060 
00061     if(lu >= 1){
00062         for(l=1; l<=lu; l++){
00063             lp1=l+1;
00064             if(l <= nct){
00065 
00066                 /* compute the transformation for the l-th column and
00067                    place the l-th diagonal in s[l][1] */
00068 
00069                 s[l][1]=snrm2(n-l+1,x,l,l,1);
00070 
00071                 if(s[l][1] != 0.0e0){
00072                     if(x[l][l] != 0.0e0)
00073                         s[l][1]=copysign(s[l][1],x[l][l]);
00074                     sscal(n-l+1,1.0e0/s[l][1],x,l,l,1);
00075                     x[l][l]=1.0e0+x[l][l];
00076                 }
00077                 s[l][1]= -s[l][1];
00078             }
00079             if(p >= lp1){
00080                 for(j=lp1; j<=p; j++){
00081                     if(l <= nct)
00082                         if(s[l][1] != 0.0e0){
00083 
00084                             /* apply the transformation */
00085 
00086                             t= -sdot(n-l+1,x,l,l,1,x,l,j,1)/x[l][l];
00087                             saxpy(n-l+1,t,x,l,l,1,x,l,j,1);
00088                         }
00089 
00090                     /* place the l-th row of x into e for the subsequent
00091                        calculation of the row transformation */
00092 
00093                     e[j][1]=x[l][j];
00094                 }
00095             }
00096             if(wantu && (l <= nct))
00097 
00098                 /* place the transformation in u for subsequent 
00099                    back multiplication */
00100 
00101                 for(i=l;i<=n;i++)
00102                     u[i][l]=x[i][l];
00103             if(l <= nrt){
00104 
00105                 /* compute the l-th row transformation and
00106                    place the l-th super-diagonal in e[l][1] */
00107 
00108                 e[l][1]=snrm2(p-l,e,lp1,1,1);
00109                 if(e[l][1] != 0.0e0){
00110                     if(e[lp1][1] != 0.0e0)
00111                         e[l][1]=copysign(e[l][1],e[lp1][1]);
00112                     sscal(p-l,1.0e0/e[l][1],e,lp1,1,1);
00113                     e[lp1][1]=1.0e0+e[lp1][1];
00114                 }
00115                 e[l][1]= -e[l][1];
00116                 if((lp1 <= n) && (e[l][1] != 0.0e0)){
00117 
00118                     /* apply the transformation */
00119 
00120                     for(i=lp1; i<=n; i++)
00121                         work[i][1]=0.0e0;
00122                     for(j=lp1; j<=p; j++)
00123                         saxpy(n-l,e[j][1],x,lp1,j,1,work,lp1,1,1);
00124                     for(j=lp1; j<=p; j++)
00125                         saxpy(n-l,-e[j][1]/e[lp1][1],work,lp1,1,1,x,lp1,j,1);
00126                 }
00127                 if(wantv)
00128 
00129                     /* place the transformation in v for subsequent
00130                        back multiplication */
00131 
00132                     for(i=lp1; i<=p; i++)
00133                         v[i][l]=e[i][1];
00134             }
00135         }
00136     }
00137 
00138     /* set up the final bidiagonal matrix or order m. */
00139 
00140     m=min0(2,p,n+1);
00141 
00142     nctp1=nct+1;
00143     nrtp1=nrt+1;
00144     if(nct < p)
00145         s[nctp1][1]=x[nctp1][nctp1];
00146     if(n < m)
00147         s[m][1]=0.0e0;
00148     if(nrtp1 < m)
00149         e[nrtp1][1]=x[nrtp1][m];
00150     e[m][1]=0.0e0;
00151 
00152     /* if required, generate u. */
00153 
00154     if(wantu){
00155         if(ncu >= nctp1){
00156             for(j=nctp1; j<=ncu; j++){
00157                 for(i=1; i<=n; i++)
00158                     u[i][j]=0.0e0;
00159                 u[j][j]=1.0e0;
00160             }
00161         }
00162         if(nct >= 1){
00163             for(ll=1; ll<=nct; ll++){
00164                 l=nct-ll+1;
00165                 if(s[l][1] != 0.0e0){
00166                     lp1=l+1;
00167                     if(ncu >= lp1){
00168                         for(j=lp1; j<=ncu; j++){
00169                             t= -sdot(n-l+1,u,l,l,1,u,l,j,1)/u[l][l];
00170                             saxpy(n-l+1,t,u,l,l,1,u,l,j,1);
00171                         }
00172                     }
00173                     sscal(n-l+1,-1.0e0,u,l,l,1);
00174                     u[l][l]=1.0e0+u[l][l];
00175                     lm1=l-1;
00176                     if(lm1 >= 1)
00177                         for(i=1; i<=lm1; i++)
00178                             u[i][l]=0.0e0;
00179                 }
00180                 else{
00181                     for(i=1; i<=n; i++)
00182                         u[i][l]=0.0e0;
00183                     u[l][l]=1.0e0;
00184                 }
00185             }
00186         }
00187     }
00188 
00189     /* if it is required, generate v. */
00190 
00191     if(wantv){
00192         for(ll=1; ll<=p; ll++){
00193             l=p-ll+1;
00194             lp1=l+1;
00195             if(l <= nrt)
00196                 if(e[l][1] != 0.0e0){
00197                     for(j=lp1; j<=p; j++){
00198                         t= -sdot(p-l,v,lp1,l,1,v,lp1,j,1)/v[lp1][l];
00199                         saxpy(p-l,t,v,lp1,l,1,v,lp1,j,1);
00200                     }
00201                 }
00202             for(i=1; i<=p; i++)
00203                 v[i][l]=0.0e0;
00204             v[l][l]=1.0e0;
00205         }
00206     }
00207 
00208     /* main iteration loop for the singular values */
00209 
00210     mm=m;
00211     iter=0;
00212 
00213     /* quit if all the singular values have been found. */
00214 
00215     while(m != 0){
00216 
00217         if(iter >= maxit){
00218             info=m;
00219             break;
00220         }
00221 
00222         /* this section of the program inspects for 
00223            negligible elements in the s and e arrays. on
00224            completion the variables kase and l are set as follows
00225 
00226            kase=1       if s(m) and e(l-1) are negligible and l<m
00227            kase=2       if s(l) is negligible and l<m
00228            kase=3       if e(l-1) is negligible, l<m, and
00229                         s(l), ..., s(m) are not negligible (QR step).
00230            kase=4       if e(m-1) is negligible (convergence).
00231                                                                 */
00232 
00233         for(ll=1; ll<=m; ll++){
00234             l=m-ll;
00235             if(l == 0)
00236                 break;
00237             test=fabs(s[l][1])+fabs(s[l+1][1]);
00238             ztest=test+fabs(e[l][1]);
00239             if(ztest == test){
00240                 e[l][1]=0.0e0;
00241                 break;
00242             }
00243         }
00244 
00245         if(l == (m-1))
00246             kase=4;
00247         else{
00248             lp1=l+1;
00249             mp1=m+1;
00250             for(lls=lp1; lls<=mp1; lls++){
00251                 ls=m-lls+lp1;
00252                 if(ls == l)
00253                     break;
00254                 test=0.0e0;
00255                 if(ls != m)
00256                     test += fabs(e[ls][1]);
00257                 if(ls != (l+1))
00258                     test += fabs(e[ls-1][1]);
00259                 ztest=test+fabs(s[ls][1]);
00260                 if(ztest == test){
00261                     s[ls][1]=0.0e0;
00262                     break;
00263                 }
00264             }
00265             
00266             if(ls == l)
00267                 kase=3;
00268             else if(ls == m)
00269                 kase=1;
00270             else{
00271                 kase=2;
00272                 l=ls;
00273             }
00274         }
00275         l=l+1;
00276 
00277         /* perform the task indicated by kase */
00278 
00279         switch(kase){
00280 
00281             /* deflate neglegible s[m] */
00282 
00283           case 1:
00284             mm1=m-1;
00285             f=e[m-1][1];
00286             e[m-1][1]=0.0e0;
00287             for(kk=l; kk<=mm1; kk++){
00288                 k=mm1-kk+l;
00289                 t1=s[k][1];
00290                 srotg(&t1,&f,&cs,&sn);
00291                 s[k][1]=t1;
00292                 if(k != l){
00293                     f= -sn*e[k-1][1];
00294                     e[k-1][1]=cs*e[k-1][1];
00295                 }
00296                 if(wantv)
00297                     srot(p,v,1,k,1,v,1,m,1,cs,sn);
00298             }
00299             break;
00300 
00301             /* split at negligible s[l] */
00302 
00303           case 2:
00304             f=e[l-1][1];
00305             e[l-1][1]=0.0e0;
00306             for(k=l; k<=m; k++){
00307                 t1=s[k][1];
00308                 srotg(&t1,&f,&cs,&sn);
00309                 s[k][1]=t1;
00310                 f= -sn*e[k][1];
00311                 e[k][1]=cs*e[k][1];
00312                 if(wantu)
00313                     srot(n,u,1,k,1,u,1,l-1,1,cs,sn);
00314             }
00315             break;
00316 
00317             /* perform one QR step */
00318 
00319           case 3:
00320 
00321             /* calculate the shift */
00322 
00323             scale=amax1(5,fabs(s[m][1]),fabs(s[m-1][1]),fabs(e[m-1][1]),
00324                          fabs(s[l][1]),fabs(e[l][1]));
00325             sm=s[m][1]/scale;
00326             smm1=s[m-1][1]/scale;
00327             emm1=e[m-1][1]/scale;
00328             sl=s[l][1]/scale;
00329             el=e[l][1]/scale;
00330             b=((smm1+sm)*(smm1-sm)+emm1*emm1)/2.0e0;
00331             c=(sm*emm1)*(sm*emm1);
00332             shift=0.0e0;
00333             if((b != 0.0e0) || (c != 0.0e0)){
00334                 shift=sqrt(b*b+c);
00335                 if(b < 0.0e0)
00336                     shift= -shift;
00337                 shift=c/(b+shift);
00338             }
00339             f=(sl+sm)*(sl-sm)-shift;
00340             g=sl*el;
00341 
00342             /* chase zeros. */
00343 
00344             mm1=m-1;
00345             for(k=l; k<=mm1; k++){
00346                 srotg(&f,&g,&cs,&sn);
00347                 if(k != l)
00348                     e[k-1][1]=f;
00349                 f=cs*s[k][1]+sn*e[k][1];
00350                 e[k][1]=cs*e[k][1]-sn*s[k][1];
00351                 g=sn*s[k+1][1];
00352                 s[k+1][1]=cs*s[k+1][1];
00353                 if(wantv)
00354                     srot(p,v,1,k,1,v,1,k+1,1,cs,sn);
00355                 srotg(&f,&g,&cs,&sn);
00356                 s[k][1]=f;
00357                 f=cs*e[k][1]+sn*s[k+1][1];
00358                 s[k+1][1]= -sn*e[k][1]+cs*s[k+1][1];
00359                 g=sn*e[k+1][1];
00360                 e[k+1][1]=cs*e[k+1][1];
00361                 if(wantu && (k < n)){
00362                     srot(n,u,1,k,1,u,1,k+1,1,cs,sn);
00363                 }
00364             }
00365             e[m-1][1]=f;
00366             iter=iter+1;
00367             break;
00368 
00369             /* convergence */
00370 
00371           case 4:
00372 
00373             /* make the singular value positive */
00374 
00375             if(s[l][1] < 0.0e0){
00376                 s[l][1]= -s[l][1];
00377                 if(wantv)
00378                     sscal(p,-1.0e0,v,1,l,1);
00379             }
00380 
00381             /* order the singular value */
00382 
00383             while(l != mm){
00384                 if(s[l][1] >= s[l+1][1])
00385                     break;
00386 
00387                 t=s[l][1];
00388                 s[l][1]=s[l+1][1];
00389                 s[l+1][1]=t;
00390                 if(wantv && (l < p))
00391                     sswap(p,v,1,l,1,v,1,l+1,1);
00392                 if(wantu && (l < n))
00393                     sswap(n,u,1,l,1,u,1,l+1,1);
00394                 l=l+1;
00395             }
00396             iter=0;
00397             m=m-1;
00398             break;
00399         }
00400     }
00401     /* shift the matrix base from 1 to 0 */
00402     /* move the singular values from the 1st column to diagonal position */
00403 
00404     for(ii=0; ii<n; ii++)
00405         for(jj=0; jj<n; jj++)
00406             u[ii][jj]=u[ii+1][jj+1];
00407 
00408     for(ii=0; ii<n; ii++){
00409         s[ii][ii]=s[ii+1][1];
00410         for(jj=0; jj<p; jj++)
00411             if(ii != jj)
00412                 s[ii][jj]=0.0e0;
00413     }
00414 
00415     for(ii=0; ii<p; ii++)
00416         for(jj=0; jj<p; jj++)
00417             v[ii][jj]=v[ii+1][jj+1];
00418 
00419     return(info);
00420 }


euslisp
Author(s): Toshihiro Matsui
autogenerated on Thu Mar 9 2017 04:57:50