contrib/contact/clib/ssvdc.c
Go to the documentation of this file.
1 /* ssvdc is a function to compute the singular value decomposition
2  of a matrix.
3 
4  Ver.1.00, Jun,3,1988.
5  Ver.1.01, Jul.14,1988. remove the side effect to the matrix X
6  Ver.1.02, Jul.15,1988. debug the matrix base shift routine */
7 
8 #include "arith.h"
9 #define MAXIT 30 /* the maximum number of iteratons */
10 
11 ssvdc(xarg,ldx,n,p,s,e,u,ldu,v,ldv,work,job)
12 int n,p,ldx,ldu,ldv,job;
13 MATRIX xarg,s,e,u,v,work;
14 {
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;
17  int ii,jj;
18 
19  int max0(),min0();
20 
21  REAL t,r;
22  REAL b,c,cs,el,emm1,f,g,scale,shift,sl,sm,sn,smm1,t1,test;
23  REAL ztest;
24  MATRIX x;
25 
26  REAL amax1(),sdot(),snrm2();
27 
28  LOGICAL wantu,wantv;
29 
30  /* shift the matrix base from 0 to 1 */
31 
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];
35 
36  /* set the maximum number of iterations */
37 
38  maxit=MAXIT;
39 
40  /* determine what is to be computed */
41 
42  wantu=FALSE;
43  wantv=FALSE;
44  jobu=(job%100)/10;
45  ncu=n;
46 
47  if(jobu > 1) ncu=min0(2,n,p);
48  if(jobu != 0) wantu=TRUE;
49  if((job%10) != 0) wantv=TRUE;
50 
51  /* reduce x to bidiagonal form,
52  storing the diagonal elements in s and
53  the super-diagonal elements in e */
54 
55  info=0;
56  nct=min0(2,n-1,p);
57  nrt=max0(2,0,min0(2,p-2,n));
58 
59  lu=max0(2,nct,nrt);
60 
61  if(lu >= 1){
62  for(l=1; l<=lu; l++){
63  lp1=l+1;
64  if(l <= nct){
65 
66  /* compute the transformation for the l-th column and
67  place the l-th diagonal in s[l][1] */
68 
69  s[l][1]=snrm2(n-l+1,x,l,l,1);
70 
71  if(s[l][1] != 0.0e0){
72  if(x[l][l] != 0.0e0)
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];
76  }
77  s[l][1]= -s[l][1];
78  }
79  if(p >= lp1){
80  for(j=lp1; j<=p; j++){
81  if(l <= nct)
82  if(s[l][1] != 0.0e0){
83 
84  /* apply the transformation */
85 
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);
88  }
89 
90  /* place the l-th row of x into e for the subsequent
91  calculation of the row transformation */
92 
93  e[j][1]=x[l][j];
94  }
95  }
96  if(wantu && (l <= nct))
97 
98  /* place the transformation in u for subsequent
99  back multiplication */
100 
101  for(i=l;i<=n;i++)
102  u[i][l]=x[i][l];
103  if(l <= nrt){
104 
105  /* compute the l-th row transformation and
106  place the l-th super-diagonal in e[l][1] */
107 
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];
114  }
115  e[l][1]= -e[l][1];
116  if((lp1 <= n) && (e[l][1] != 0.0e0)){
117 
118  /* apply the transformation */
119 
120  for(i=lp1; i<=n; i++)
121  work[i][1]=0.0e0;
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);
126  }
127  if(wantv)
128 
129  /* place the transformation in v for subsequent
130  back multiplication */
131 
132  for(i=lp1; i<=p; i++)
133  v[i][l]=e[i][1];
134  }
135  }
136  }
137 
138  /* set up the final bidiagonal matrix or order m. */
139 
140  m=min0(2,p,n+1);
141 
142  nctp1=nct+1;
143  nrtp1=nrt+1;
144  if(nct < p)
145  s[nctp1][1]=x[nctp1][nctp1];
146  if(n < m)
147  s[m][1]=0.0e0;
148  if(nrtp1 < m)
149  e[nrtp1][1]=x[nrtp1][m];
150  e[m][1]=0.0e0;
151 
152  /* if required, generate u. */
153 
154  if(wantu){
155  if(ncu >= nctp1){
156  for(j=nctp1; j<=ncu; j++){
157  for(i=1; i<=n; i++)
158  u[i][j]=0.0e0;
159  u[j][j]=1.0e0;
160  }
161  }
162  if(nct >= 1){
163  for(ll=1; ll<=nct; ll++){
164  l=nct-ll+1;
165  if(s[l][1] != 0.0e0){
166  lp1=l+1;
167  if(ncu >= lp1){
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);
171  }
172  }
173  sscal(n-l+1,-1.0e0,u,l,l,1);
174  u[l][l]=1.0e0+u[l][l];
175  lm1=l-1;
176  if(lm1 >= 1)
177  for(i=1; i<=lm1; i++)
178  u[i][l]=0.0e0;
179  }
180  else{
181  for(i=1; i<=n; i++)
182  u[i][l]=0.0e0;
183  u[l][l]=1.0e0;
184  }
185  }
186  }
187  }
188 
189  /* if it is required, generate v. */
190 
191  if(wantv){
192  for(ll=1; ll<=p; ll++){
193  l=p-ll+1;
194  lp1=l+1;
195  if(l <= nrt)
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);
200  }
201  }
202  for(i=1; i<=p; i++)
203  v[i][l]=0.0e0;
204  v[l][l]=1.0e0;
205  }
206  }
207 
208  /* main iteration loop for the singular values */
209 
210  mm=m;
211  iter=0;
212 
213  /* quit if all the singular values have been found. */
214 
215  while(m != 0){
216 
217  if(iter >= maxit){
218  info=m;
219  break;
220  }
221 
222  /* this section of the program inspects for
223  negligible elements in the s and e arrays. on
224  completion the variables kase and l are set as follows
225 
226  kase=1 if s(m) and e(l-1) are negligible and l<m
227  kase=2 if s(l) is negligible and l<m
228  kase=3 if e(l-1) is negligible, l<m, and
229  s(l), ..., s(m) are not negligible (QR step).
230  kase=4 if e(m-1) is negligible (convergence).
231  */
232 
233  for(ll=1; ll<=m; ll++){
234  l=m-ll;
235  if(l == 0)
236  break;
237  test=fabs(s[l][1])+fabs(s[l+1][1]);
238  ztest=test+fabs(e[l][1]);
239  if(ztest == test){
240  e[l][1]=0.0e0;
241  break;
242  }
243  }
244 
245  if(l == (m-1))
246  kase=4;
247  else{
248  lp1=l+1;
249  mp1=m+1;
250  for(lls=lp1; lls<=mp1; lls++){
251  ls=m-lls+lp1;
252  if(ls == l)
253  break;
254  test=0.0e0;
255  if(ls != m)
256  test += fabs(e[ls][1]);
257  if(ls != (l+1))
258  test += fabs(e[ls-1][1]);
259  ztest=test+fabs(s[ls][1]);
260  if(ztest == test){
261  s[ls][1]=0.0e0;
262  break;
263  }
264  }
265 
266  if(ls == l)
267  kase=3;
268  else if(ls == m)
269  kase=1;
270  else{
271  kase=2;
272  l=ls;
273  }
274  }
275  l=l+1;
276 
277  /* perform the task indicated by kase */
278 
279  switch(kase){
280 
281  /* deflate neglegible s[m] */
282 
283  case 1:
284  mm1=m-1;
285  f=e[m-1][1];
286  e[m-1][1]=0.0e0;
287  for(kk=l; kk<=mm1; kk++){
288  k=mm1-kk+l;
289  t1=s[k][1];
290  srotg(&t1,&f,&cs,&sn);
291  s[k][1]=t1;
292  if(k != l){
293  f= -sn*e[k-1][1];
294  e[k-1][1]=cs*e[k-1][1];
295  }
296  if(wantv)
297  srot(p,v,1,k,1,v,1,m,1,cs,sn);
298  }
299  break;
300 
301  /* split at negligible s[l] */
302 
303  case 2:
304  f=e[l-1][1];
305  e[l-1][1]=0.0e0;
306  for(k=l; k<=m; k++){
307  t1=s[k][1];
308  srotg(&t1,&f,&cs,&sn);
309  s[k][1]=t1;
310  f= -sn*e[k][1];
311  e[k][1]=cs*e[k][1];
312  if(wantu)
313  srot(n,u,1,k,1,u,1,l-1,1,cs,sn);
314  }
315  break;
316 
317  /* perform one QR step */
318 
319  case 3:
320 
321  /* calculate the shift */
322 
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]));
325  sm=s[m][1]/scale;
326  smm1=s[m-1][1]/scale;
327  emm1=e[m-1][1]/scale;
328  sl=s[l][1]/scale;
329  el=e[l][1]/scale;
330  b=((smm1+sm)*(smm1-sm)+emm1*emm1)/2.0e0;
331  c=(sm*emm1)*(sm*emm1);
332  shift=0.0e0;
333  if((b != 0.0e0) || (c != 0.0e0)){
334  shift=sqrt(b*b+c);
335  if(b < 0.0e0)
336  shift= -shift;
337  shift=c/(b+shift);
338  }
339  f=(sl+sm)*(sl-sm)-shift;
340  g=sl*el;
341 
342  /* chase zeros. */
343 
344  mm1=m-1;
345  for(k=l; k<=mm1; k++){
346  srotg(&f,&g,&cs,&sn);
347  if(k != l)
348  e[k-1][1]=f;
349  f=cs*s[k][1]+sn*e[k][1];
350  e[k][1]=cs*e[k][1]-sn*s[k][1];
351  g=sn*s[k+1][1];
352  s[k+1][1]=cs*s[k+1][1];
353  if(wantv)
354  srot(p,v,1,k,1,v,1,k+1,1,cs,sn);
355  srotg(&f,&g,&cs,&sn);
356  s[k][1]=f;
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];
359  g=sn*e[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);
363  }
364  }
365  e[m-1][1]=f;
366  iter=iter+1;
367  break;
368 
369  /* convergence */
370 
371  case 4:
372 
373  /* make the singular value positive */
374 
375  if(s[l][1] < 0.0e0){
376  s[l][1]= -s[l][1];
377  if(wantv)
378  sscal(p,-1.0e0,v,1,l,1);
379  }
380 
381  /* order the singular value */
382 
383  while(l != mm){
384  if(s[l][1] >= s[l+1][1])
385  break;
386 
387  t=s[l][1];
388  s[l][1]=s[l+1][1];
389  s[l+1][1]=t;
390  if(wantv && (l < p))
391  sswap(p,v,1,l,1,v,1,l+1,1);
392  if(wantu && (l < n))
393  sswap(n,u,1,l,1,u,1,l+1,1);
394  l=l+1;
395  }
396  iter=0;
397  m=m-1;
398  break;
399  }
400  }
401  /* shift the matrix base from 1 to 0 */
402  /* move the singular values from the 1st column to diagonal position */
403 
404  for(ii=0; ii<n; ii++)
405  for(jj=0; jj<n; jj++)
406  u[ii][jj]=u[ii+1][jj+1];
407 
408  for(ii=0; ii<n; ii++){
409  s[ii][ii]=s[ii+1][1];
410  for(jj=0; jj<p; jj++)
411  if(ii != jj)
412  s[ii][jj]=0.0e0;
413  }
414 
415  for(ii=0; ii<p; ii++)
416  for(jj=0; jj<p; jj++)
417  v[ii][jj]=v[ii+1][jj+1];
418 
419  return(info);
420 }
f
srotg(REAL *psa, REAL *psb, REAL *pc, REAL *ps)
Definition: srotg.c:6
sswap(int n, MATRIX sx, int isx, int jsx, int incx, MATRIX sy, int isy, int jsy, int incy)
Definition: sswap.c:7
GLfloat n[6][3]
Definition: cube.c:15
REAL amax1(va_alist)
Definition: amax1.c:7
#define MAXIT
ssvdc(MATRIX xarg, int ldx, int n, int p, MATRIX s, MATRIX e, MATRIX u, int ldu, MATRIX v, int ldv, MATRIX work, int job)
double REAL
Definition: arith.h:25
int max0(va_alist)
Definition: max0.c:7
#define TRUE
Definition: arith.h:19
REAL sdot(int n, MATRIX sx, int isx, int jsx, int incx, MATRIX sy, int isy, int jsy, int incy)
Definition: sdot.c:7
int LOGICAL
Definition: arith.h:31
long l
Definition: structsize.c:3
double sqrt()
sscal(int n, REAL sa, MATRIX sx, int isx, int jsx, int incx)
Definition: sscal.c:7
#define FALSE
Definition: arith.h:20
short s
Definition: structsize.c:2
int min0(va_alist)
Definition: min0.c:7
VECTOR MATRIX[MAX]
Definition: arith.h:27
GLfloat v[8][3]
Definition: cube.c:21
saxpy(int n, REAL sa, MATRIX sx, int isx, int jsx, int incx, MATRIX sy, int isy, int jsy, int incy)
Definition: saxpy.c:7
srot(int n, MATRIX sx, int isx, int jsx, int incx, MATRIX sy, int isy, int jsy, int incy, REAL c, REAL s)
Definition: srot.c:6
REAL snrm2(int n, MATRIX sx, int isx, int jsx, int incx)
Definition: snrm2.c:12
static int(* g)()
Definition: cfunc.c:9


euslisp
Author(s): Toshihiro Matsui
autogenerated on Mon Feb 28 2022 22:18:28