00001
00002
00003
00004
00005
00006
00007
00008 #include "arith.h"
00009 #define MAXIT 30
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
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
00037
00038 maxit=MAXIT;
00039
00040
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
00052
00053
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
00067
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
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
00091
00092
00093 e[j][1]=x[l][j];
00094 }
00095 }
00096 if(wantu && (l <= nct))
00097
00098
00099
00100
00101 for(i=l;i<=n;i++)
00102 u[i][l]=x[i][l];
00103 if(l <= nrt){
00104
00105
00106
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
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
00130
00131
00132 for(i=lp1; i<=p; i++)
00133 v[i][l]=e[i][1];
00134 }
00135 }
00136 }
00137
00138
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
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
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
00209
00210 mm=m;
00211 iter=0;
00212
00213
00214
00215 while(m != 0){
00216
00217 if(iter >= maxit){
00218 info=m;
00219 break;
00220 }
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
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
00278
00279 switch(kase){
00280
00281
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
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
00318
00319 case 3:
00320
00321
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
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
00370
00371 case 4:
00372
00373
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
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
00402
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 }