11 ssvdc(xarg,ldx,
n,p,
s,e,u,ldu,
v,ldv,work,job)
12 int n,p,ldx,ldu,ldv,job;
15 int i,iter,j,jobu,k,kase,kk,
l,ll,lls,lm1,lp1,ls,lu,m,maxit,info;
16 int mm,mm1,mp1,nct,nctp1,ncu,nrt,nrtp1;
22 REAL b,c,cs,el,emm1,
f,
g,scale,shift,sl,sm,sn,smm1,t1,test;
32 for(ii=
n-1; ii>=0; ii--)
33 for(jj=p-1; jj>=0; jj--)
34 x[ii+1][jj+1]=xarg[ii][jj];
47 if(jobu > 1) ncu=
min0(2,
n,p);
48 if(jobu != 0) wantu=
TRUE;
49 if((job%10) != 0) wantv=
TRUE;
73 s[
l][1]=copysign(
s[l][1],x[l][l]);
74 sscal(
n-l+1,1.0e0/
s[l][1],x,l,l,1);
75 x[
l][
l]=1.0e0+x[
l][
l];
80 for(j=lp1; j<=p; j++){
86 t= -
sdot(
n-l+1,x,l,l,1,x,l,j,1)/x[
l][
l];
87 saxpy(
n-l+1,t,x,l,l,1,x,l,j,1);
96 if(wantu && (l <= nct))
108 e[
l][1]=
snrm2(p-l,e,lp1,1,1);
109 if(e[l][1] != 0.0e0){
110 if(e[lp1][1] != 0.0e0)
111 e[
l][1]=copysign(e[l][1],e[lp1][1]);
112 sscal(p-l,1.0e0/e[l][1],e,lp1,1,1);
113 e[lp1][1]=1.0e0+e[lp1][1];
116 if((lp1 <=
n) && (e[l][1] != 0.0e0)){
120 for(i=lp1; i<=
n; i++)
122 for(j=lp1; j<=p; j++)
123 saxpy(
n-l,e[j][1],x,lp1,j,1,work,lp1,1,1);
124 for(j=lp1; j<=p; j++)
125 saxpy(
n-l,-e[j][1]/e[lp1][1],work,lp1,1,1,x,lp1,j,1);
132 for(i=lp1; i<=p; i++)
145 s[nctp1][1]=x[nctp1][nctp1];
149 e[nrtp1][1]=x[nrtp1][m];
156 for(j=nctp1; j<=ncu; j++){
163 for(ll=1; ll<=nct; ll++){
165 if(
s[l][1] != 0.0e0){
168 for(j=lp1; j<=ncu; j++){
169 t= -
sdot(
n-l+1,u,l,l,1,u,l,j,1)/u[
l][
l];
170 saxpy(
n-l+1,t,u,l,l,1,u,l,j,1);
173 sscal(
n-l+1,-1.0e0,u,l,l,1);
174 u[
l][
l]=1.0e0+u[
l][
l];
177 for(i=1; i<=lm1; i++)
192 for(ll=1; ll<=p; ll++){
196 if(e[l][1] != 0.0e0){
197 for(j=lp1; j<=p; j++){
198 t= -
sdot(p-l,
v,lp1,l,1,
v,lp1,j,1)/
v[lp1][
l];
199 saxpy(p-l,t,
v,lp1,l,1,
v,lp1,j,1);
233 for(ll=1; ll<=m; ll++){
237 test=fabs(
s[l][1])+fabs(
s[l+1][1]);
238 ztest=test+fabs(e[l][1]);
250 for(lls=lp1; lls<=mp1; lls++){
256 test += fabs(e[ls][1]);
258 test += fabs(e[ls-1][1]);
259 ztest=test+fabs(
s[ls][1]);
287 for(kk=l; kk<=mm1; kk++){
290 srotg(&t1,&f,&cs,&sn);
294 e[k-1][1]=cs*e[k-1][1];
297 srot(p,
v,1,k,1,
v,1,m,1,cs,sn);
308 srotg(&t1,&f,&cs,&sn);
313 srot(
n,u,1,k,1,u,1,l-1,1,cs,sn);
323 scale=
amax1(5,fabs(
s[m][1]),fabs(
s[m-1][1]),fabs(e[m-1][1]),
324 fabs(
s[l][1]),fabs(e[l][1]));
326 smm1=
s[m-1][1]/scale;
327 emm1=e[m-1][1]/scale;
330 b=((smm1+sm)*(smm1-sm)+emm1*emm1)/2.0e0;
331 c=(sm*emm1)*(sm*emm1);
333 if((b != 0.0e0) || (c != 0.0e0)){
339 f=(sl+sm)*(sl-sm)-shift;
345 for(k=l; k<=mm1; k++){
346 srotg(&f,&g,&cs,&sn);
349 f=cs*
s[k][1]+sn*e[k][1];
350 e[k][1]=cs*e[k][1]-sn*
s[k][1];
352 s[k+1][1]=cs*s[k+1][1];
354 srot(p,
v,1,k,1,
v,1,k+1,1,cs,sn);
355 srotg(&f,&g,&cs,&sn);
357 f=cs*e[k][1]+sn*s[k+1][1];
358 s[k+1][1]= -sn*e[k][1]+cs*s[k+1][1];
360 e[k+1][1]=cs*e[k+1][1];
361 if(wantu && (k <
n)){
362 srot(
n,u,1,k,1,u,1,k+1,1,cs,sn);
384 if(
s[l][1] >=
s[l+1][1])
393 sswap(
n,u,1,l,1,u,1,l+1,1);
404 for(ii=0; ii<
n; ii++)
405 for(jj=0; jj<
n; jj++)
406 u[ii][jj]=u[ii+1][jj+1];
408 for(ii=0; ii<
n; ii++){
409 s[ii][ii]=
s[ii+1][1];
410 for(jj=0; jj<p; jj++)
415 for(ii=0; ii<p; ii++)
416 for(jj=0; jj<p; jj++)
417 v[ii][jj]=
v[ii+1][jj+1];
srotg(REAL *psa, REAL *psb, REAL *pc, REAL *ps)
sswap(int n, MATRIX sx, int isx, int jsx, int incx, MATRIX sy, int isy, int jsy, int incy)
REAL sdot(int n, MATRIX sx, int isx, int jsx, int incx, MATRIX sy, int isy, int jsy, int incy)
sscal(int n, REAL sa, MATRIX sx, int isx, int jsx, int incx)
saxpy(int n, REAL sa, MATRIX sx, int isx, int jsx, int incx, MATRIX sy, int isy, int jsy, int incy)
srot(int n, MATRIX sx, int isx, int jsx, int incx, MATRIX sy, int isy, int jsy, int incy, REAL c, REAL s)
REAL snrm2(int n, MATRIX sx, int isx, int jsx, int incx)