1 /* This file contains code from two sources not covered by the GNU GPL that the
2 DynaMechs library is distributed under: LINPACK's ssvdc routine
3 that was translated to C using f2c, and portions of the f2c header file that
4 are need to compile the converted code. The copyrights and licenses are
5 found in comments in this file.
6 */
7
8 #if defined(WIN32)
9 #include <windows.h>
10 #endif
11
12 #include <stdio.h>
13 #include <iostream>
14 #include <iomanip>
15 using namespace std;
16
17 #ifdef __cplusplus
18 extern "C" {
19 #endif
20
21 /* -------------------------------- f2c.h -------------------------------
22
23 The following section was extracted from the f2c package (obtained from
24 http://www.netlib.org/f2c), is not covered by the GNU GPL, and comes with
25 the following permissions statement:
26
27 Copyright 1990 - 1997 by AT&T, Lucent Technologies and Bellcore.
28
29 Permission to use, copy, modify, and distribute this software
30 and its documentation for any purpose and without fee is hereby
31 granted, provided that the above copyright notice appear in all
32 copies and that both that the copyright notice and this
33 permission notice and warranty disclaimer appear in supporting
34 documentation, and that the names of AT&T, Bell Laboratories,
35 Lucent or Bellcore or any of their entities not be used in
36 advertising or publicity pertaining to distribution of the
37 software without specific, written prior permission.
38
39 AT&T, Lucent and Bellcore disclaim all warranties with regard to
40 this software, including all implied warranties of
41 merchantability and fitness. In no event shall AT&T, Lucent or
42 Bellcore be liable for any special, indirect or consequential
43 damages or any damages whatsoever resulting from loss of use,
44 data or profits, whether in an action of contract, negligence or
45 other tortious action, arising out of or in connection with the
46 use or performance of this software.
47 */
48
49 /* f2c.h -- Standard Fortran to C header file */
50
51 /** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed."
52
53 - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */
54
55 #ifndef F2C_INCLUDE
56 #define F2C_INCLUDE
57
58 typedef long int integer;
59 typedef unsigned long uinteger;
60 typedef char *address;
61 typedef short int shortint;
62 typedef float real;
63 typedef double doublereal;
64 typedef struct { real r, i; } complex;
65 typedef struct { doublereal r, i; } doublecomplex;
66 typedef long int logical;
67 typedef short int shortlogical;
68 typedef char logical1;
69 typedef char integer1;
70 #if 0 /* Adjust for integer*8. */
71 typedef long long longint; /* system-dependent */
72 typedef unsigned long long ulongint; /* system-dependent */
73 #define qbit_clear(a,b) ((a) & ~((ulongint)1 << (b)))
74 #define qbit_set(a,b) ((a) | ((ulongint)1 << (b)))
75 #endif
76
77 #define TRUE_ (1)
78 #define FALSE_ (0)
79
80 /* Extern is for use with -E */
81 #ifndef Extern
82 #define Extern extern
83 #endif
84
85 /* I/O stuff */
86
87 #ifdef f2c_i2
88 /* for -i2 */
89 typedef short flag;
90 typedef short ftnlen;
91 typedef short ftnint;
92 #else
93 typedef long int flag;
94 typedef long int ftnlen;
95 typedef long int ftnint;
96 #endif
97
98 /*external read, write*/
99 typedef struct
100 { flag cierr;
101 ftnint ciunit;
102 flag ciend;
103 char *cifmt;
104 ftnint cirec;
105 } cilist;
106
107 /*internal read, write*/
108 typedef struct
109 { flag icierr;
110 char *iciunit;
111 flag iciend;
112 char *icifmt;
113 ftnint icirlen;
114 ftnint icirnum;
115 } icilist;
116
117 /*open*/
118 typedef struct
119 { flag oerr;
120 ftnint ounit;
121 char *ofnm;
122 ftnlen ofnmlen;
123 char *osta;
124 char *oacc;
125 char *ofm;
126 ftnint orl;
127 char *oblnk;
128 } olist;
129
130 /*close*/
131 typedef struct
132 { flag cerr;
133 ftnint cunit;
134 char *csta;
135 } cllist;
136
137 /*rewind, backspace, endfile*/
138 typedef struct
139 { flag aerr;
140 ftnint aunit;
141 } alist;
142
143 /* inquire */
144 typedef struct
145 { flag inerr;
146 ftnint inunit;
147 char *infile;
148 ftnlen infilen;
149 ftnint *inex; /*parameters in standard's order*/
150 ftnint *inopen;
151 ftnint *innum;
152 ftnint *innamed;
153 char *inname;
154 ftnlen innamlen;
155 char *inacc;
156 ftnlen inacclen;
157 char *inseq;
158 ftnlen inseqlen;
159 char *indir;
160 ftnlen indirlen;
161 char *infmt;
162 ftnlen infmtlen;
163 char *inform;
164 ftnint informlen;
165 char *inunf;
166 ftnlen inunflen;
167 ftnint *inrecl;
168 ftnint *innrec;
169 char *inblank;
170 ftnlen inblanklen;
171 } inlist;
172
173 #define VOID void
174
175 union Multitype { /* for multiple entry points */
176 integer1 g;
177 shortint h;
178 integer i;
179 /* longint j; */
180 real r;
181 doublereal d;
182 complex c;
183 doublecomplex z;
184 };
185
186 typedef union Multitype Multitype;
187
188 /*typedef long int Long;*/ /* No longer used; formerly in Namelist */
189
190 struct Vardesc { /* for Namelist */
191 char *name;
192 char *addr;
193 ftnlen *dims;
194 int type;
195 };
196 typedef struct Vardesc Vardesc;
197
198 struct Namelist {
199 char *name;
200 Vardesc **vars;
201 int nvars;
202 };
203 typedef struct Namelist Namelist;
204
205 #define abs(x) ((x) >= 0 ? (x) : -(x))
206 #define dabs(x) (doublereal)abs(x)
207 #define omin(a,b) ((a) <= (b) ? (a) : (b))
208 #define omax(a,b) ((a) >= (b) ? (a) : (b))
209 #define dmin(a,b) (doublereal)omin(a,b)
210 #define dmax(a,b) (doublereal)omax(a,b)
211 #define bit_test(a,b) ((a) >> (b) & 1)
212 #define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
213 #define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
214
215 /* procedure parameter types for -A and -C++ */
216
217 #define F2C_proc_par_types 1
218 #ifdef __cplusplus
219 typedef int /* Unknown procedure type */ (*U_fp)(...);
220 typedef shortint (*J_fp)(...);
221 typedef integer (*I_fp)(...);
222 typedef real (*R_fp)(...);
223 typedef doublereal (*D_fp)(...), (*E_fp)(...);
224 typedef /* Complex */ VOID (*C_fp)(...);
225 typedef /* Double Complex */ VOID (*Z_fp)(...);
226 typedef logical (*L_fp)(...);
227 typedef shortlogical (*K_fp)(...);
228 typedef /* Character */ VOID (*H_fp)(...);
229 typedef /* Subroutine */ int (*S_fp)(...);
230 #else
231 typedef int /* Unknown procedure type */ (*U_fp)();
232 typedef shortint (*J_fp)();
233 typedef integer (*I_fp)();
234 typedef real (*R_fp)();
235 typedef doublereal (*D_fp)(), (*E_fp)();
236 typedef /* Complex */ VOID (*C_fp)();
237 typedef /* Double Complex */ VOID (*Z_fp)();
238 typedef logical (*L_fp)();
239 typedef shortlogical (*K_fp)();
240 typedef /* Character */ VOID (*H_fp)();
241 typedef /* Subroutine */ int (*S_fp)();
242 #endif
243 /* E_fp is for real functions when -R is not specified */
244 typedef VOID C_f; /* complex function */
245 typedef VOID H_f; /* character function */
246 typedef VOID Z_f; /* double complex function */
247 typedef doublereal E_f; /* real function with -R not specified */
248
249 /* undef any lower-case symbols that your C compiler predefines, e.g.: */
250
251 #ifndef Skip_f2c_Undefs
252 #undef cray
253 #undef gcos
254 #undef mc68010
255 #undef mc68020
256 #undef mips
257 #undef pdp11
258 #undef sgi
259 #undef sparc
260 #undef sun
261 #undef sun2
262 #undef sun3
263 #undef sun4
264 #undef u370
265 #undef u3b
266 #undef u3b2
267 #undef u3b5
268 #undef unix
269 #undef vax
270 #endif
271 #endif
272
273
274 /* ----------------------------------------------------------------------
275 ssvdc - single precision SVD routine from LINPACK
276 dsvdc - double precision SVD routine from LINPACK
277
278 Obtained from http:://www.netlib.org/linpack (search for ssvdc and dsvdc)
279
280 LINPACK is a collection of Fortran subroutines that analyze and
281 solve linear equations and linear least-squares probles. The
282 package solves linear systems whose matrices are general, banded,
283 symmetric indefinite, symmetric positive definite, triangular,
284 and tridiagonal square. In addition, the package computes
285 the QR and singular value decompositions of rectangular matrices
286 and applies them to least-squares problems. LINPACK uses
287 column-oriented algorithms to increase efficiency by preserving
288 locality of reference.
289
290 LINPACK was designed for supercomputers in use in the 1970s and
291 early 1980s. LINPACK has been largely superceded by LAPACK
292 which has been designed to run efficiently on shared-memory, vector
293 supercomputers.
294
295 Developed by Jack Dongarra, Jim Bunch, Cleve Moler and Pete Stewart.
296 1 Feb 84
297
298 If you are interested in acquiring the entire LINPACK, it may
299 make more sense to talk with NAG. NAG distribute the software
300 on a mag tape for a nominal charge.
301 NAG
302 1400 Opus Place, Suite 200
303 Downers Grove, IL 60515-5702
304 708-971-2337, FAX 971-2706
305
306 ----------------------------------------------------------------------------*/
307
308 /* ssvdc.f -- translated by f2c (version 19990311).
309 You must link the resulting object file with the libraries:
310 -lf2c -lm (in that order)
311 */
312
313 /* Duane Marhefka Note: It is not necessary to link with f2c. Only one
314 function (r_sign) in the f2c library was required, so I appended that
315 function from the libf77.c file to the end of this file.
316
317 This file was created with the following commands:
318 > f2cx -C++ -R ssvdc.f
319 > move ssvdc.c ssvdc.cpp
320
321 And then the r_sign function was appended from libf77.c.
322 */
323
324 /* Scott McMillan Notes:
325
326 This wouldn't build on SGI because of a redifinition of
327 the sqrt function, so I "fixed" that. But this code will not converge on
328 SGI systems....I still don't know why.
329
330 On Windoze systems the min and max macros were redefinitions to I changed
331 them to omin and omax.
332 */
333
334 /* Table of constant values */
335
336 static integer c__1 = 1;
337 static real c_b44 = (float)-1.;
338 static real c_b122 = (float)1.;
339
ssvdc_(real * x,integer * ldx,integer * n,integer * p,real * s,real * e,real * u,integer * ldu,real * v,integer * ldv,real * work,integer * job,integer * info)340 /* Subroutine */ int ssvdc_(real *x, integer *ldx, integer *n, integer *p,
341 real *s, real *e, real *u, integer *ldu, real *v, integer *ldv, real *
342 work, integer *job, integer *info)
343 {
344 #if 0
345 cout << "\nBi_star (in routine);" << endl;
346 for (int rr = 0; rr < *n; rr++)
347 {
348 for (int cc = 0; cc < *n; cc++)
349 cout << setprecision(4) << setw(12) << x[rr*(*n) + cc];
350 cout << endl;
351 }
352 #endif
353
354 /* System generated locals */
355 integer x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2,
356 i__3;
357 real r__1, r__2, r__3, r__4, r__5, r__6, r__7;
358
359 /* Builtin functions */
360 real r_sign(real *, real *);
361 double sqrt(doublereal);
362
363 /* Local variables */
364 static integer kase, jobu, iter;
365 extern real sdot_(integer *, real *, integer *, real *, integer *);
366 static real test;
367 extern /* Subroutine */ int srot_(integer *, real *, integer *, real *,
368 integer *, real *, real *);
369 static integer nctp1;
370 extern real snrm2_(integer *, real *, integer *);
371 static real b, c__;
372 static integer nrtp1;
373 static real f, g;
374 static integer i__, j, k, l, m;
375 static real t, scale;
376 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
377 static real shift;
378 static integer maxit;
379 extern /* Subroutine */ int sswap_(integer *, real *, integer *, real *,
380 integer *);
381 static logical wantu, wantv;
382 extern /* Subroutine */ int srotg_(real *, real *, real *, real *),
383 saxpy_(integer *, real *, real *, integer *, real *, integer *);
384 static real t1, ztest, el;
385 static integer kk;
386 static real cs;
387 static integer ll, mm, ls;
388 static real sl;
389 static integer lu;
390 static real sm, sn;
391 static integer lm1, mm1, lp1, mp1, nct, ncu, lls, nrt;
392 static real emm1, smm1;
393
394
395
396 /* ssvdc is a subroutine to reduce a real nxp matrix x by */
397 /* orthogonal transformations u and v to diagonal form. the */
398 /* diagonal elements s(i) are the singular values of x. the */
399 /* columns of u are the corresponding left singular vectors, */
400 /* and the columns of v the right singular vectors. */
401
402 /* on entry */
403
404 /* x real(ldx,p), where ldx.ge.n. */
405 /* x contains the matrix whose singular value */
406 /* decomposition is to be computed. x is */
407 /* destroyed by ssvdc. */
408
409 /* ldx integer. */
410 /* ldx is the leading dimension of the array x. */
411
412 /* n integer. */
413 /* n is the number of rows of the matrix x. */
414
415 /* p integer. */
416 /* p is the number of columns of the matrix x. */
417
418 /* ldu integer. */
419 /* ldu is the leading dimension of the array u. */
420 /* (see below). */
421
422 /* ldv integer. */
423 /* ldv is the leading dimension of the array v. */
424 /* (see below). */
425
426 /* work real(n). */
427 /* work is a scratch array. */
428
429 /* job integer. */
430 /* job controls the computation of the singular */
431 /* vectors. it has the decimal expansion ab */
432 /* with the following meaning */
433
434 /* a.eq.0 do not compute the left singular */
435 /* vectors. */
436 /* a.eq.1 return the n left singular vectors */
437 /* in u. */
438 /* a.ge.2 return the first min(n,p) singular */
439 /* vectors in u. */
440 /* b.eq.0 do not compute the right singular */
441 /* vectors. */
442 /* b.eq.1 return the right singular vectors */
443 /* in v. */
444
445 /* on return */
446
447 /* s real(mm), where mm=min(n+1,p). */
448 /* the first min(n,p) entries of s contain the */
449 /* singular values of x arranged in descending */
450 /* order of magnitude. */
451
452 /* e real(p). */
453 /* e ordinarily contains zeros. however see the */
454 /* discussion of info for exceptions. */
455
456 /* u real(ldu,k), where ldu.ge.n. if joba.eq.1 then */
457 /* k.eq.n, if joba.ge.2 then */
458 /* k.eq.min(n,p). */
459 /* u contains the matrix of left singular vectors. */
460 /* u is not referenced if joba.eq.0. if n.le.p */
461 /* or if joba.eq.2, then u may be identified with x */
462 /* in the subroutine call. */
463
464 /* v real(ldv,p), where ldv.ge.p. */
465 /* v contains the matrix of right singular vectors. */
466 /* v is not referenced if job.eq.0. if p.le.n, */
467 /* then v may be identified with x in the */
468 /* subroutine call. */
469
470 /* info integer. */
471 /* the singular values (and their corresponding */
472 /* singular vectors) s(info+1),s(info+2),...,s(m) */
473 /* are correct (here m=min(n,p)). thus if */
474 /* info.eq.0, all the singular values and their */
475 /* vectors are correct. in any event, the matrix */
476 /* b = trans(u)*x*v is the bidiagonal matrix */
477 /* with the elements of s on its diagonal and the */
478 /* elements of e on its super-diagonal (trans(u) */
479 /* is the transpose of u). thus the singular */
480 /* values of x and b are the same. */
481
482 /* linpack. this version dated 03/19/79 . */
483 /* correction to shift calculation made 2/85. */
484 /* g.w. stewart, university of maryland, argonne national lab. */
485
486 /* ***** uses the following functions and subprograms. */
487
488 /* external srot */
489 /* blas saxpy,sdot,sscal,sswap,snrm2,srotg */
490 /* fortran abs,amax1,max0,min0,mod,sqrt */
491
492 /* internal variables */
493
494
495
496 /* set the maximum number of iterations. */
497
498 /* Parameter adjustments */
499 x_dim1 = *ldx;
500 x_offset = 1 + x_dim1 * 1;
501 x -= x_offset;
502 --s;
503 --e;
504 u_dim1 = *ldu;
505 u_offset = 1 + u_dim1 * 1;
506 u -= u_offset;
507 v_dim1 = *ldv;
508 v_offset = 1 + v_dim1 * 1;
509 v -= v_offset;
510 --work;
511
512 /* Function Body */
513 maxit = 30;
514
515 /* determine what is to be computed. */
516
517 wantu = FALSE_;
518 wantv = FALSE_;
519 jobu = *job % 100 / 10;
520 ncu = *n;
521 if (jobu > 1) {
522 ncu = omin(*n,*p);
523 }
524 if (jobu != 0) {
525 wantu = TRUE_;
526 }
527 if (*job % 10 != 0) {
528 wantv = TRUE_;
529 }
530
531 /* reduce x to bidiagonal form, storing the diagonal elements */
532 /* in s and the super-diagonal elements in e. */
533
534 *info = 0;
535 /* Computing MIN */
536 i__1 = *n - 1;
537 nct = omin(i__1,*p);
538 /* Computing MAX */
539 /* Computing MIN */
540 i__3 = *p - 2;
541 i__1 = 0, i__2 = omin(i__3,*n);
542 nrt = omax(i__1,i__2);
543 lu = omax(nct,nrt);
544 if (lu < 1) {
545 goto L170;
546 }
547 i__1 = lu;
548 for (l = 1; l <= i__1; ++l) {
549 lp1 = l + 1;
550 if (l > nct) {
551 goto L20;
552 }
553
554 /* compute the transformation for the l-th column and */
555 /* place the l-th diagonal in s(l). */
556
557 i__2 = *n - l + 1;
558 s[l] = snrm2_(&i__2, &x[l + l * x_dim1], &c__1);
559 if (s[l] == (float)0.) {
560 goto L10;
561 }
562 if (x[l + l * x_dim1] != (float)0.) {
563 s[l] = r_sign(&s[l], &x[l + l * x_dim1]);
564 }
565 i__2 = *n - l + 1;
566 r__1 = (float)1. / s[l];
567 sscal_(&i__2, &r__1, &x[l + l * x_dim1], &c__1);
568 x[l + l * x_dim1] += (float)1.;
569 L10:
570 s[l] = -s[l];
571 L20:
572 if (*p < lp1) {
573 goto L50;
574 }
575 i__2 = *p;
576 for (j = lp1; j <= i__2; ++j) {
577 if (l > nct) {
578 goto L30;
579 }
580 if (s[l] == (float)0.) {
581 goto L30;
582 }
583
584 /* apply the transformation. */
585
586 i__3 = *n - l + 1;
587 t = -sdot_(&i__3, &x[l + l * x_dim1], &c__1, &x[l + j * x_dim1], &
588 c__1) / x[l + l * x_dim1];
589 i__3 = *n - l + 1;
590 saxpy_(&i__3, &t, &x[l + l * x_dim1], &c__1, &x[l + j * x_dim1], &
591 c__1);
592 L30:
593
594 /* place the l-th row of x into e for the */
595 /* subsequent calculation of the row transformation. */
596
597 e[j] = x[l + j * x_dim1];
598 /* L40: */
599 }
600 L50:
601 if (! wantu || l > nct) {
602 goto L70;
603 }
604
605 /* place the transformation in u for subsequent back */
606 /* multiplication. */
607
608 i__2 = *n;
609 for (i__ = l; i__ <= i__2; ++i__) {
610 u[i__ + l * u_dim1] = x[i__ + l * x_dim1];
611 /* L60: */
612 }
613 L70:
614 if (l > nrt) {
615 goto L150;
616 }
617
618 /* compute the l-th row transformation and place the */
619 /* l-th super-diagonal in e(l). */
620
621 i__2 = *p - l;
622 e[l] = snrm2_(&i__2, &e[lp1], &c__1);
623 if (e[l] == (float)0.) {
624 goto L80;
625 }
626 if (e[lp1] != (float)0.) {
627 e[l] = r_sign(&e[l], &e[lp1]);
628 }
629 i__2 = *p - l;
630 r__1 = (float)1. / e[l];
631 sscal_(&i__2, &r__1, &e[lp1], &c__1);
632 e[lp1] += (float)1.;
633 L80:
634 e[l] = -e[l];
635 if (lp1 > *n || e[l] == (float)0.) {
636 goto L120;
637 }
638
639 /* apply the transformation. */
640
641 i__2 = *n;
642 for (i__ = lp1; i__ <= i__2; ++i__) {
643 work[i__] = (float)0.;
644 /* L90: */
645 }
646 i__2 = *p;
647 for (j = lp1; j <= i__2; ++j) {
648 i__3 = *n - l;
649 saxpy_(&i__3, &e[j], &x[lp1 + j * x_dim1], &c__1, &work[lp1], &
650 c__1);
651 /* L100: */
652 }
653 i__2 = *p;
654 for (j = lp1; j <= i__2; ++j) {
655 i__3 = *n - l;
656 r__1 = -e[j] / e[lp1];
657 saxpy_(&i__3, &r__1, &work[lp1], &c__1, &x[lp1 + j * x_dim1], &
658 c__1);
659 /* L110: */
660 }
661 L120:
662 if (! wantv) {
663 goto L140;
664 }
665
666 /* place the transformation in v for subsequent */
667 /* back multiplication. */
668
669 i__2 = *p;
670 for (i__ = lp1; i__ <= i__2; ++i__) {
671 v[i__ + l * v_dim1] = e[i__];
672 /* L130: */
673 }
674 L140:
675 L150:
676 /* L160: */
677 ;
678 }
679 L170:
680
681 /* set up the final bidiagonal matrix or order m. */
682
683 /* Computing MIN */
684 i__1 = *p, i__2 = *n + 1;
685 m = omin(i__1,i__2);
686 nctp1 = nct + 1;
687 nrtp1 = nrt + 1;
688 if (nct < *p) {
689 s[nctp1] = x[nctp1 + nctp1 * x_dim1];
690 }
691 if (*n < m) {
692 s[m] = (float)0.;
693 }
694 if (nrtp1 < m) {
695 e[nrtp1] = x[nrtp1 + m * x_dim1];
696 }
697 e[m] = (float)0.;
698
699 /* if required, generate u. */
700
701 if (! wantu) {
702 goto L300;
703 }
704 if (ncu < nctp1) {
705 goto L200;
706 }
707 i__1 = ncu;
708 for (j = nctp1; j <= i__1; ++j) {
709 i__2 = *n;
710 for (i__ = 1; i__ <= i__2; ++i__) {
711 u[i__ + j * u_dim1] = (float)0.;
712 /* L180: */
713 }
714 u[j + j * u_dim1] = (float)1.;
715 /* L190: */
716 }
717 L200:
718 if (nct < 1) {
719 goto L290;
720 }
721 i__1 = nct;
722 for (ll = 1; ll <= i__1; ++ll) {
723 l = nct - ll + 1;
724 if (s[l] == (float)0.) {
725 goto L250;
726 }
727 lp1 = l + 1;
728 if (ncu < lp1) {
729 goto L220;
730 }
731 i__2 = ncu;
732 for (j = lp1; j <= i__2; ++j) {
733 i__3 = *n - l + 1;
734 t = -sdot_(&i__3, &u[l + l * u_dim1], &c__1, &u[l + j * u_dim1], &
735 c__1) / u[l + l * u_dim1];
736 i__3 = *n - l + 1;
737 saxpy_(&i__3, &t, &u[l + l * u_dim1], &c__1, &u[l + j * u_dim1], &
738 c__1);
739 /* L210: */
740 }
741 L220:
742 i__2 = *n - l + 1;
743 sscal_(&i__2, &c_b44, &u[l + l * u_dim1], &c__1);
744 u[l + l * u_dim1] += (float)1.;
745 lm1 = l - 1;
746 if (lm1 < 1) {
747 goto L240;
748 }
749 i__2 = lm1;
750 for (i__ = 1; i__ <= i__2; ++i__) {
751 u[i__ + l * u_dim1] = (float)0.;
752 /* L230: */
753 }
754 L240:
755 goto L270;
756 L250:
757 i__2 = *n;
758 for (i__ = 1; i__ <= i__2; ++i__) {
759 u[i__ + l * u_dim1] = (float)0.;
760 /* L260: */
761 }
762 u[l + l * u_dim1] = (float)1.;
763 L270:
764 /* L280: */
765 ;
766 }
767 L290:
768 L300:
769
770 /* if it is required, generate v. */
771
772 if (! wantv) {
773 goto L350;
774 }
775 i__1 = *p;
776 for (ll = 1; ll <= i__1; ++ll) {
777 l = *p - ll + 1;
778 lp1 = l + 1;
779 if (l > nrt) {
780 goto L320;
781 }
782 if (e[l] == (float)0.) {
783 goto L320;
784 }
785 i__2 = *p;
786 for (j = lp1; j <= i__2; ++j) {
787 i__3 = *p - l;
788 t = -sdot_(&i__3, &v[lp1 + l * v_dim1], &c__1, &v[lp1 + j *
789 v_dim1], &c__1) / v[lp1 + l * v_dim1];
790 i__3 = *p - l;
791 saxpy_(&i__3, &t, &v[lp1 + l * v_dim1], &c__1, &v[lp1 + j *
792 v_dim1], &c__1);
793 /* L310: */
794 }
795 L320:
796 i__2 = *p;
797 for (i__ = 1; i__ <= i__2; ++i__) {
798 v[i__ + l * v_dim1] = (float)0.;
799 /* L330: */
800 }
801 v[l + l * v_dim1] = (float)1.;
802 /* L340: */
803 }
804 L350:
805
806 /* main iteration loop for the singular values. */
807
808 mm = m;
809 iter = 0;
810 L360:
811
812 /* quit if all the singular values have been found. */
813
814 /* ...exit */
815 if (m == 0) {
816 goto L620;
817 }
818
819 /* if too many iterations have been performed, set */
820 /* flag and return. */
821
822 if (iter < maxit) {
823 // cout << "Iteration = " << iter << endl;
824 goto L370;
825 }
826 *info = m;
827 /* ......exit */
828 goto L620;
829 L370:
830
831 /* this section of the program inspects for */
832 /* negligible elements in the s and e arrays. on */
833 /* completion the variables kase and l are set as follows. */
834
835 /* kase = 1 if s(m) and e(l-1) are negligible and l.lt.m */
836 /* kase = 2 if s(l) is negligible and l.lt.m */
837 /* kase = 3 if e(l-1) is negligible, l.lt.m, and */
838 /* s(l), ..., s(m) are not negligible (qr step). */
839 /* kase = 4 if e(m-1) is negligible (convergence). */
840
841 i__1 = m;
842 for (ll = 1; ll <= i__1; ++ll) {
843 l = m - ll;
844 /* ...exit */
845 if (l == 0) {
846 goto L400;
847 }
848 test = (r__1 = s[l], abs(r__1)) + (r__2 = s[l + 1], abs(r__2));
849 ztest = test + (r__1 = e[l], abs(r__1));
850 if (ztest != test) {
851 goto L380;
852 }
853 e[l] = (float)0.;
854 /* ......exit */
855 goto L400;
856 L380:
857 /* L390: */
858 ;
859 }
860 L400:
861 if (l != m - 1) {
862 goto L410;
863 }
864 kase = 4;
865 goto L480;
866 L410:
867 lp1 = l + 1;
868 mp1 = m + 1;
869 i__1 = mp1;
870 for (lls = lp1; lls <= i__1; ++lls) {
871 ls = m - lls + lp1;
872 /* ...exit */
873 if (ls == l) {
874 goto L440;
875 }
876 test = (float)0.;
877 if (ls != m) {
878 test += (r__1 = e[ls], abs(r__1));
879 }
880 if (ls != l + 1) {
881 test += (r__1 = e[ls - 1], abs(r__1));
882 }
883 ztest = test + (r__1 = s[ls], abs(r__1));
884 if (ztest != test) {
885 goto L420;
886 }
887 s[ls] = (float)0.;
888 /* ......exit */
889 goto L440;
890 L420:
891 /* L430: */
892 ;
893 }
894 L440:
895 if (ls != l) {
896 goto L450;
897 }
898 kase = 3;
899 goto L470;
900 L450:
901 if (ls != m) {
902 goto L460;
903 }
904 kase = 1;
905 goto L470;
906 L460:
907 kase = 2;
908 l = ls;
909 L470:
910 L480:
911 ++l;
912
913 /* perform the task indicated by kase. */
914
915 switch (kase) {
916 case 1: goto L490;
917 case 2: goto L520;
918 case 3: goto L540;
919 case 4: goto L570;
920 }
921
922 /* deflate negligible s(m). */
923
924 L490:
925 mm1 = m - 1;
926 f = e[m - 1];
927 e[m - 1] = (float)0.;
928 i__1 = mm1;
929 for (kk = l; kk <= i__1; ++kk) {
930 k = mm1 - kk + l;
931 t1 = s[k];
932 srotg_(&t1, &f, &cs, &sn);
933 s[k] = t1;
934 if (k == l) {
935 goto L500;
936 }
937 f = -sn * e[k - 1];
938 e[k - 1] = cs * e[k - 1];
939 L500:
940 if (wantv) {
941 srot_(p, &v[k * v_dim1 + 1], &c__1, &v[m * v_dim1 + 1], &c__1, &
942 cs, &sn);
943 }
944 /* L510: */
945 }
946 goto L610;
947
948 /* split at negligible s(l). */
949
950 L520:
951 f = e[l - 1];
952 e[l - 1] = (float)0.;
953 i__1 = m;
954 for (k = l; k <= i__1; ++k) {
955 t1 = s[k];
956 srotg_(&t1, &f, &cs, &sn);
957 s[k] = t1;
958 f = -sn * e[k];
959 e[k] = cs * e[k];
960 if (wantu) {
961 srot_(n, &u[k * u_dim1 + 1], &c__1, &u[(l - 1) * u_dim1 + 1], &
962 c__1, &cs, &sn);
963 }
964 /* L530: */
965 }
966 goto L610;
967
968 /* perform one qr step. */
969
970 L540:
971
972 /* calculate the shift. */
973
974 /* Computing MAX */
975 r__6 = (r__1 = s[m], abs(r__1)), r__7 = (r__2 = s[m - 1], abs(r__2)),
976 r__6 = omax(r__6,r__7), r__7 = (r__3 = e[m - 1], abs(r__3)), r__6 =
977 omax(r__6,r__7), r__7 = (r__4 = s[l], abs(r__4)), r__6 = omax(r__6,
978 r__7), r__7 = (r__5 = e[l], abs(r__5));
979 scale = omax(r__6,r__7);
980 sm = s[m] / scale;
981 smm1 = s[m - 1] / scale;
982 emm1 = e[m - 1] / scale;
983 sl = s[l] / scale;
984 el = e[l] / scale;
985 /* Computing 2nd power */
986 r__1 = emm1;
987 b = ((smm1 + sm) * (smm1 - sm) + r__1 * r__1) / (float)2.;
988 /* Computing 2nd power */
989 r__1 = sm * emm1;
990 c__ = r__1 * r__1;
991 shift = (float)0.;
992 if (b == (float)0. && c__ == (float)0.) {
993 goto L550;
994 }
995 /* Computing 2nd power */
996 r__1 = b;
997 if ( (r__1 * r__1 + c__) < 0 )
998 cout << "ERROR: sqrt(NEG)" << endl;
999 shift = (float)sqrt(r__1 * r__1 + c__);
1000 if (b < (float)0.) {
1001 shift = -shift;
1002 }
1003 shift = c__ / (b + shift);
1004 L550:
1005 f = (sl + sm) * (sl - sm) + shift;
1006 g = sl * el;
1007
1008 /* chase zeros. */
1009
1010 mm1 = m - 1;
1011 i__1 = mm1;
1012 for (k = l; k <= i__1; ++k) {
1013 srotg_(&f, &g, &cs, &sn);
1014 if (k != l) {
1015 e[k - 1] = f;
1016 }
1017 f = cs * s[k] + sn * e[k];
1018 e[k] = cs * e[k] - sn * s[k];
1019 g = sn * s[k + 1];
1020 s[k + 1] = cs * s[k + 1];
1021 if (wantv) {
1022 srot_(p, &v[k * v_dim1 + 1], &c__1, &v[(k + 1) * v_dim1 + 1], &
1023 c__1, &cs, &sn);
1024 }
1025 srotg_(&f, &g, &cs, &sn);
1026 s[k] = f;
1027 f = cs * e[k] + sn * s[k + 1];
1028 s[k + 1] = -sn * e[k] + cs * s[k + 1];
1029 g = sn * e[k + 1];
1030 e[k + 1] = cs * e[k + 1];
1031 if (wantu && k < *n) {
1032 srot_(n, &u[k * u_dim1 + 1], &c__1, &u[(k + 1) * u_dim1 + 1], &
1033 c__1, &cs, &sn);
1034 }
1035 /* L560: */
1036 }
1037 e[m - 1] = f;
1038 ++iter;
1039 goto L610;
1040
1041 /* convergence. */
1042
1043 L570:
1044
1045 /* make the singular value positive. */
1046
1047 if (s[l] >= (float)0.) {
1048 goto L580;
1049 }
1050 s[l] = -s[l];
1051 if (wantv) {
1052 sscal_(p, &c_b44, &v[l * v_dim1 + 1], &c__1);
1053 }
1054 L580:
1055
1056 /* order the singular value. */
1057
1058 L590:
1059 if (l == mm) {
1060 goto L600;
1061 }
1062 /* ...exit */
1063 if (s[l] >= s[l + 1]) {
1064 goto L600;
1065 }
1066 t = s[l];
1067 s[l] = s[l + 1];
1068 s[l + 1] = t;
1069 if (wantv && l < *p) {
1070 sswap_(p, &v[l * v_dim1 + 1], &c__1, &v[(l + 1) * v_dim1 + 1], &c__1);
1071 }
1072 if (wantu && l < *n) {
1073 sswap_(n, &u[l * u_dim1 + 1], &c__1, &u[(l + 1) * u_dim1 + 1], &c__1);
1074 }
1075 ++l;
1076 goto L590;
1077 L600:
1078 iter = 0;
1079 --m;
1080 L610:
1081 goto L360;
1082 L620:
1083 return 0;
1084 } /* ssvdc_ */
1085
saxpy_(integer * n,real * sa,real * sx,integer * incx,real * sy,integer * incy)1086 /* Subroutine */ int saxpy_(integer *n, real *sa, real *sx, integer *incx,
1087 real *sy, integer *incy)
1088 {
1089 /* System generated locals */
1090 integer i__1;
1091
1092 /* Local variables */
1093 static integer i__, m, ix, iy, mp1;
1094
1095
1096 /* constant times a vector plus a vector. */
1097 /* uses unrolled loop for increments equal to one. */
1098 /* jack dongarra, linpack, 3/11/78. */
1099 /* modified 12/3/93, array(1) declarations changed to array(*) */
1100
1101
1102 /* Parameter adjustments */
1103 --sy;
1104 --sx;
1105
1106 /* Function Body */
1107 if (*n <= 0) {
1108 return 0;
1109 }
1110 if (*sa == (float)0.) {
1111 return 0;
1112 }
1113 if (*incx == 1 && *incy == 1) {
1114 goto L20;
1115 }
1116
1117 /* code for unequal increments or equal increments */
1118 /* not equal to 1 */
1119
1120 ix = 1;
1121 iy = 1;
1122 if (*incx < 0) {
1123 ix = (-(*n) + 1) * *incx + 1;
1124 }
1125 if (*incy < 0) {
1126 iy = (-(*n) + 1) * *incy + 1;
1127 }
1128 i__1 = *n;
1129 for (i__ = 1; i__ <= i__1; ++i__) {
1130 sy[iy] += *sa * sx[ix];
1131 ix += *incx;
1132 iy += *incy;
1133 /* L10: */
1134 }
1135 return 0;
1136
1137 /* code for both increments equal to 1 */
1138
1139
1140 /* clean-up loop */
1141
1142 L20:
1143 m = *n % 4;
1144 if (m == 0) {
1145 goto L40;
1146 }
1147 i__1 = m;
1148 for (i__ = 1; i__ <= i__1; ++i__) {
1149 sy[i__] += *sa * sx[i__];
1150 /* L30: */
1151 }
1152 if (*n < 4) {
1153 return 0;
1154 }
1155 L40:
1156 mp1 = m + 1;
1157 i__1 = *n;
1158 for (i__ = mp1; i__ <= i__1; i__ += 4) {
1159 sy[i__] += *sa * sx[i__];
1160 sy[i__ + 1] += *sa * sx[i__ + 1];
1161 sy[i__ + 2] += *sa * sx[i__ + 2];
1162 sy[i__ + 3] += *sa * sx[i__ + 3];
1163 /* L50: */
1164 }
1165 return 0;
1166 } /* saxpy_ */
1167
sdot_(integer * n,real * sx,integer * incx,real * sy,integer * incy)1168 real sdot_(integer *n, real *sx, integer *incx, real *sy, integer *incy)
1169 {
1170 /* System generated locals */
1171 integer i__1;
1172 real ret_val;
1173
1174 /* Local variables */
1175 static integer i__, m;
1176 static real stemp;
1177 static integer ix, iy, mp1;
1178
1179
1180 /* forms the dot product of two vectors. */
1181 /* uses unrolled loops for increments equal to one. */
1182 /* jack dongarra, linpack, 3/11/78. */
1183 /* modified 12/3/93, array(1) declarations changed to array(*) */
1184
1185
1186 /* Parameter adjustments */
1187 --sy;
1188 --sx;
1189
1190 /* Function Body */
1191 stemp = (float)0.;
1192 ret_val = (float)0.;
1193 if (*n <= 0) {
1194 return ret_val;
1195 }
1196 if (*incx == 1 && *incy == 1) {
1197 goto L20;
1198 }
1199
1200 /* code for unequal increments or equal increments */
1201 /* not equal to 1 */
1202
1203 ix = 1;
1204 iy = 1;
1205 if (*incx < 0) {
1206 ix = (-(*n) + 1) * *incx + 1;
1207 }
1208 if (*incy < 0) {
1209 iy = (-(*n) + 1) * *incy + 1;
1210 }
1211 i__1 = *n;
1212 for (i__ = 1; i__ <= i__1; ++i__) {
1213 stemp += sx[ix] * sy[iy];
1214 ix += *incx;
1215 iy += *incy;
1216 /* L10: */
1217 }
1218 ret_val = stemp;
1219 return ret_val;
1220
1221 /* code for both increments equal to 1 */
1222
1223
1224 /* clean-up loop */
1225
1226 L20:
1227 m = *n % 5;
1228 if (m == 0) {
1229 goto L40;
1230 }
1231 i__1 = m;
1232 for (i__ = 1; i__ <= i__1; ++i__) {
1233 stemp += sx[i__] * sy[i__];
1234 /* L30: */
1235 }
1236 if (*n < 5) {
1237 goto L60;
1238 }
1239 L40:
1240 mp1 = m + 1;
1241 i__1 = *n;
1242 for (i__ = mp1; i__ <= i__1; i__ += 5) {
1243 stemp = stemp + sx[i__] * sy[i__] + sx[i__ + 1] * sy[i__ + 1] + sx[
1244 i__ + 2] * sy[i__ + 2] + sx[i__ + 3] * sy[i__ + 3] + sx[i__ +
1245 4] * sy[i__ + 4];
1246 /* L50: */
1247 }
1248 L60:
1249 ret_val = stemp;
1250 return ret_val;
1251 } /* sdot_ */
1252
snrm2_(integer * n,real * x,integer * incx)1253 real snrm2_(integer *n, real *x, integer *incx)
1254 {
1255 /* System generated locals */
1256 integer i__1, i__2;
1257 real ret_val, r__1;
1258
1259 /* Builtin functions */
1260 double sqrt(doublereal);
1261
1262 /* Local variables */
1263 static real norm, scale, absxi;
1264 static integer ix;
1265 static real ssq;
1266
1267 /* .. Scalar Arguments .. */
1268 /* .. Array Arguments .. */
1269 /* .. */
1270
1271 /* SNRM2 returns the euclidean norm of a vector via the function */
1272 /* name, so that */
1273
1274 /* SNRM2 := sqrt( x'*x ) */
1275
1276
1277
1278 /* -- This version written on 25-October-1982. */
1279 /* Modified on 14-October-1993 to inline the call to SLASSQ. */
1280 /* Sven Hammarling, Nag Ltd. */
1281
1282
1283 /* .. Parameters .. */
1284 /* .. Local Scalars .. */
1285 /* .. Intrinsic Functions .. */
1286 /* .. */
1287 /* .. Executable Statements .. */
1288 /* Parameter adjustments */
1289 --x;
1290
1291 /* Function Body */
1292 if (*n < 1 || *incx < 1) {
1293 norm = (float)0.;
1294 } else if (*n == 1) {
1295 norm = abs(x[1]);
1296 } else {
1297 scale = (float)0.;
1298 ssq = (float)1.;
1299 /* The following loop is equivalent to this call to the LAPACK */
1300 /* auxiliary routine: */
1301 /* CALL SLASSQ( N, X, INCX, SCALE, SSQ ) */
1302
1303 i__1 = (*n - 1) * *incx + 1;
1304 i__2 = *incx;
1305 for (ix = 1; i__2 < 0 ? ix >= i__1 : ix <= i__1; ix += i__2) {
1306 if (x[ix] != (float)0.) {
1307 absxi = (r__1 = x[ix], abs(r__1));
1308 if (scale < absxi) {
1309 /* Computing 2nd power */
1310 r__1 = scale / absxi;
1311 ssq = ssq * (r__1 * r__1) + (float)1.;
1312 scale = absxi;
1313 } else {
1314 /* Computing 2nd power */
1315 r__1 = absxi / scale;
1316 ssq += r__1 * r__1;
1317 }
1318 }
1319 /* L10: */
1320 }
1321 if ( ssq < 0 )
1322 cout << "ERROR: sqrt(NEG)" << endl;
1323 norm = float(scale * sqrt(ssq));
1324 }
1325
1326 ret_val = norm;
1327 return ret_val;
1328
1329 /* End of SNRM2. */
1330
1331 } /* snrm2_ */
1332
srot_(integer * n,real * sx,integer * incx,real * sy,integer * incy,real * c__,real * s)1333 /* Subroutine */ int srot_(integer *n, real *sx, integer *incx, real *sy,
1334 integer *incy, real *c__, real *s)
1335 {
1336 /* System generated locals */
1337 integer i__1;
1338
1339 /* Local variables */
1340 static integer i__;
1341 static real stemp;
1342 static integer ix, iy;
1343
1344
1345 /* applies a plane rotation. */
1346 /* jack dongarra, linpack, 3/11/78. */
1347 /* modified 12/3/93, array(1) declarations changed to array(*) */
1348
1349
1350 /* Parameter adjustments */
1351 --sy;
1352 --sx;
1353
1354 /* Function Body */
1355 if (*n <= 0) {
1356 return 0;
1357 }
1358 if (*incx == 1 && *incy == 1) {
1359 goto L20;
1360 }
1361
1362 /* code for unequal increments or equal increments not equal */
1363 /* to 1 */
1364
1365 ix = 1;
1366 iy = 1;
1367 if (*incx < 0) {
1368 ix = (-(*n) + 1) * *incx + 1;
1369 }
1370 if (*incy < 0) {
1371 iy = (-(*n) + 1) * *incy + 1;
1372 }
1373 i__1 = *n;
1374 for (i__ = 1; i__ <= i__1; ++i__) {
1375 stemp = *c__ * sx[ix] + *s * sy[iy];
1376 sy[iy] = *c__ * sy[iy] - *s * sx[ix];
1377 sx[ix] = stemp;
1378 ix += *incx;
1379 iy += *incy;
1380 /* L10: */
1381 }
1382 return 0;
1383
1384 /* code for both increments equal to 1 */
1385
1386 L20:
1387 i__1 = *n;
1388 for (i__ = 1; i__ <= i__1; ++i__) {
1389 stemp = *c__ * sx[i__] + *s * sy[i__];
1390 sy[i__] = *c__ * sy[i__] - *s * sx[i__];
1391 sx[i__] = stemp;
1392 /* L30: */
1393 }
1394 return 0;
1395 } /* srot_ */
1396
srotg_(real * sa,real * sb,real * c__,real * s)1397 /* Subroutine */ int srotg_(real *sa, real *sb, real *c__, real *s)
1398 {
1399 /* System generated locals */
1400 real r__1, r__2;
1401
1402 /* Builtin functions */
1403 double sqrt(doublereal);
1404 real r_sign(real *, real *);
1405
1406 /* Local variables */
1407 static real r__, scale, z__, roe;
1408
1409
1410 /* construct givens plane rotation. */
1411 /* jack dongarra, linpack, 3/11/78. */
1412
1413
1414 roe = *sb;
1415 if (abs(*sa) > abs(*sb)) {
1416 roe = *sa;
1417 }
1418 scale = abs(*sa) + abs(*sb);
1419 if (scale != (float)0.) {
1420 goto L10;
1421 }
1422 *c__ = (float)1.;
1423 *s = (float)0.;
1424 r__ = (float)0.;
1425 z__ = (float)0.;
1426 goto L20;
1427 L10:
1428 /* Computing 2nd power */
1429 r__1 = *sa / scale;
1430 /* Computing 2nd power */
1431 r__2 = *sb / scale;
1432 if ( (r__1 * r__1 + r__2 * r__2) < 0 )
1433 cout << "ERROR: sqrt(NEG)" << endl;
1434 r__ = float(scale * sqrt(r__1 * r__1 + r__2 * r__2));
1435 r__ = r_sign(&c_b122, &roe) * r__;
1436 *c__ = *sa / r__;
1437 *s = *sb / r__;
1438 z__ = (float)1.;
1439 if (abs(*sa) > abs(*sb)) {
1440 z__ = *s;
1441 }
1442 if (abs(*sb) >= abs(*sa) && *c__ != (float)0.) {
1443 z__ = (float)1. / *c__;
1444 }
1445 L20:
1446 *sa = r__;
1447 *sb = z__;
1448 return 0;
1449 } /* srotg_ */
1450
sscal_(integer * n,real * sa,real * sx,integer * incx)1451 /* Subroutine */ int sscal_(integer *n, real *sa, real *sx, integer *incx)
1452 {
1453 /* System generated locals */
1454 integer i__1, i__2;
1455
1456 /* Local variables */
1457 static integer i__, m, nincx, mp1;
1458
1459
1460 /* scales a vector by a constant. */
1461 /* uses unrolled loops for increment equal to 1. */
1462 /* jack dongarra, linpack, 3/11/78. */
1463 /* modified 3/93 to return if incx .le. 0. */
1464 /* modified 12/3/93, array(1) declarations changed to array(*) */
1465
1466
1467 /* Parameter adjustments */
1468 --sx;
1469
1470 /* Function Body */
1471 if (*n <= 0 || *incx <= 0) {
1472 return 0;
1473 }
1474 if (*incx == 1) {
1475 goto L20;
1476 }
1477
1478 /* code for increment not equal to 1 */
1479
1480 nincx = *n * *incx;
1481 i__1 = nincx;
1482 i__2 = *incx;
1483 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
1484 sx[i__] = *sa * sx[i__];
1485 /* L10: */
1486 }
1487 return 0;
1488
1489 /* code for increment equal to 1 */
1490
1491
1492 /* clean-up loop */
1493
1494 L20:
1495 m = *n % 5;
1496 if (m == 0) {
1497 goto L40;
1498 }
1499 i__2 = m;
1500 for (i__ = 1; i__ <= i__2; ++i__) {
1501 sx[i__] = *sa * sx[i__];
1502 /* L30: */
1503 }
1504 if (*n < 5) {
1505 return 0;
1506 }
1507 L40:
1508 mp1 = m + 1;
1509 i__2 = *n;
1510 for (i__ = mp1; i__ <= i__2; i__ += 5) {
1511 sx[i__] = *sa * sx[i__];
1512 sx[i__ + 1] = *sa * sx[i__ + 1];
1513 sx[i__ + 2] = *sa * sx[i__ + 2];
1514 sx[i__ + 3] = *sa * sx[i__ + 3];
1515 sx[i__ + 4] = *sa * sx[i__ + 4];
1516 /* L50: */
1517 }
1518 return 0;
1519 } /* sscal_ */
1520
sswap_(integer * n,real * sx,integer * incx,real * sy,integer * incy)1521 /* Subroutine */ int sswap_(integer *n, real *sx, integer *incx, real *sy,
1522 integer *incy)
1523 {
1524 /* System generated locals */
1525 integer i__1;
1526
1527 /* Local variables */
1528 static integer i__, m;
1529 static real stemp;
1530 static integer ix, iy, mp1;
1531
1532
1533 /* interchanges two vectors. */
1534 /* uses unrolled loops for increments equal to 1. */
1535 /* jack dongarra, linpack, 3/11/78. */
1536 /* modified 12/3/93, array(1) declarations changed to array(*) */
1537
1538
1539 /* Parameter adjustments */
1540 --sy;
1541 --sx;
1542
1543 /* Function Body */
1544 if (*n <= 0) {
1545 return 0;
1546 }
1547 if (*incx == 1 && *incy == 1) {
1548 goto L20;
1549 }
1550
1551 /* code for unequal increments or equal increments not equal */
1552 /* to 1 */
1553
1554 ix = 1;
1555 iy = 1;
1556 if (*incx < 0) {
1557 ix = (-(*n) + 1) * *incx + 1;
1558 }
1559 if (*incy < 0) {
1560 iy = (-(*n) + 1) * *incy + 1;
1561 }
1562 i__1 = *n;
1563 for (i__ = 1; i__ <= i__1; ++i__) {
1564 stemp = sx[ix];
1565 sx[ix] = sy[iy];
1566 sy[iy] = stemp;
1567 ix += *incx;
1568 iy += *incy;
1569 /* L10: */
1570 }
1571 return 0;
1572
1573 /* code for both increments equal to 1 */
1574
1575
1576 /* clean-up loop */
1577
1578 L20:
1579 m = *n % 3;
1580 if (m == 0) {
1581 goto L40;
1582 }
1583 i__1 = m;
1584 for (i__ = 1; i__ <= i__1; ++i__) {
1585 stemp = sx[i__];
1586 sx[i__] = sy[i__];
1587 sy[i__] = stemp;
1588 /* L30: */
1589 }
1590 if (*n < 3) {
1591 return 0;
1592 }
1593 L40:
1594 mp1 = m + 1;
1595 i__1 = *n;
1596 for (i__ = mp1; i__ <= i__1; i__ += 3) {
1597 stemp = sx[i__];
1598 sx[i__] = sy[i__];
1599 sy[i__] = stemp;
1600 stemp = sx[i__ + 1];
1601 sx[i__ + 1] = sy[i__ + 1];
1602 sy[i__ + 1] = stemp;
1603 stemp = sx[i__ + 2];
1604 sx[i__ + 2] = sy[i__ + 2];
1605 sy[i__ + 2] = stemp;
1606 /* L50: */
1607 }
1608 return 0;
1609 } /* sswap_ */
1610
1611 #ifdef KR_headers
1612 real r_sign(a,b) real *a, *b;
1613 #else
r_sign(real * a,real * b)1614 real r_sign(real *a, real *b)
1615 #endif
1616 {
1617 real x;
1618 x = (*a >= 0 ? *a : - *a);
1619 return( *b >= 0 ? x : -x);
1620 }
1621
1622 /* ------------------------------------------------------------------------ */
1623
1624 /* dsvdc.f -- translated by f2c (version 19990311).
1625 You must link the resulting object file with the libraries:
1626 -lf2c -lm (in that order)
1627 */
1628
1629 /* Duane Marhefka Note: It is not necessary to link with f2c. Only one
1630 function (d_sign) in the f2c library was required, so I appended that
1631 function from the libf77.c file to the end of this file.
1632
1633 This file was created with the following commands:
1634 > f2cx -C++ -R dsvdc.f
1635 > move dsvdc.c dsvdc.cpp
1636
1637 And then the d_sign function was appended from libf77.c.
1638 */
1639
1640 /* Table of constant values */
1641
1642 static doublereal c_b23 = 1.;
1643 static integer c__1d = 1;
1644 static doublereal c_b79 = -1.;
1645
daxpy_(integer * n,doublereal * da,doublereal * dx,integer * incx,doublereal * dy,integer * incy)1646 /* Subroutine */ int daxpy_(integer *n, doublereal *da, doublereal *dx,
1647 integer *incx, doublereal *dy, integer *incy)
1648 {
1649 /* System generated locals */
1650 integer i__1;
1651
1652 /* Local variables */
1653 static integer i__, m, ix, iy, mp1;
1654
1655
1656 /* constant times a vector plus a vector. */
1657 /* uses unrolled loops for increments equal to one. */
1658 /* jack dongarra, linpack, 3/11/78. */
1659 /* modified 12/3/93, array(1) declarations changed to array(*) */
1660
1661
1662 /* Parameter adjustments */
1663 --dy;
1664 --dx;
1665
1666 /* Function Body */
1667 if (*n <= 0) {
1668 return 0;
1669 }
1670 if (*da == 0.) {
1671 return 0;
1672 }
1673 if (*incx == 1 && *incy == 1) {
1674 goto L20;
1675 }
1676
1677 /* code for unequal increments or equal increments */
1678 /* not equal to 1 */
1679
1680 ix = 1;
1681 iy = 1;
1682 if (*incx < 0) {
1683 ix = (-(*n) + 1) * *incx + 1;
1684 }
1685 if (*incy < 0) {
1686 iy = (-(*n) + 1) * *incy + 1;
1687 }
1688 i__1 = *n;
1689 for (i__ = 1; i__ <= i__1; ++i__) {
1690 dy[iy] += *da * dx[ix];
1691 ix += *incx;
1692 iy += *incy;
1693 /* L10: */
1694 }
1695 return 0;
1696
1697 /* code for both increments equal to 1 */
1698
1699
1700 /* clean-up loop */
1701
1702 L20:
1703 m = *n % 4;
1704 if (m == 0) {
1705 goto L40;
1706 }
1707 i__1 = m;
1708 for (i__ = 1; i__ <= i__1; ++i__) {
1709 dy[i__] += *da * dx[i__];
1710 /* L30: */
1711 }
1712 if (*n < 4) {
1713 return 0;
1714 }
1715 L40:
1716 mp1 = m + 1;
1717 i__1 = *n;
1718 for (i__ = mp1; i__ <= i__1; i__ += 4) {
1719 dy[i__] += *da * dx[i__];
1720 dy[i__ + 1] += *da * dx[i__ + 1];
1721 dy[i__ + 2] += *da * dx[i__ + 2];
1722 dy[i__ + 3] += *da * dx[i__ + 3];
1723 /* L50: */
1724 }
1725 return 0;
1726 } /* daxpy_ */
1727
ddot_(integer * n,doublereal * dx,integer * incx,doublereal * dy,integer * incy)1728 doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy,
1729 integer *incy)
1730 {
1731 /* System generated locals */
1732 integer i__1;
1733 doublereal ret_val;
1734
1735 /* Local variables */
1736 static integer i__, m;
1737 static doublereal dtemp;
1738 static integer ix, iy, mp1;
1739
1740
1741 /* forms the dot product of two vectors. */
1742 /* uses unrolled loops for increments equal to one. */
1743 /* jack dongarra, linpack, 3/11/78. */
1744 /* modified 12/3/93, array(1) declarations changed to array(*) */
1745
1746
1747 /* Parameter adjustments */
1748 --dy;
1749 --dx;
1750
1751 /* Function Body */
1752 ret_val = 0.;
1753 dtemp = 0.;
1754 if (*n <= 0) {
1755 return ret_val;
1756 }
1757 if (*incx == 1 && *incy == 1) {
1758 goto L20;
1759 }
1760
1761 /* code for unequal increments or equal increments */
1762 /* not equal to 1 */
1763
1764 ix = 1;
1765 iy = 1;
1766 if (*incx < 0) {
1767 ix = (-(*n) + 1) * *incx + 1;
1768 }
1769 if (*incy < 0) {
1770 iy = (-(*n) + 1) * *incy + 1;
1771 }
1772 i__1 = *n;
1773 for (i__ = 1; i__ <= i__1; ++i__) {
1774 dtemp += dx[ix] * dy[iy];
1775 ix += *incx;
1776 iy += *incy;
1777 /* L10: */
1778 }
1779 ret_val = dtemp;
1780 return ret_val;
1781
1782 /* code for both increments equal to 1 */
1783
1784
1785 /* clean-up loop */
1786
1787 L20:
1788 m = *n % 5;
1789 if (m == 0) {
1790 goto L40;
1791 }
1792 i__1 = m;
1793 for (i__ = 1; i__ <= i__1; ++i__) {
1794 dtemp += dx[i__] * dy[i__];
1795 /* L30: */
1796 }
1797 if (*n < 5) {
1798 goto L60;
1799 }
1800 L40:
1801 mp1 = m + 1;
1802 i__1 = *n;
1803 for (i__ = mp1; i__ <= i__1; i__ += 5) {
1804 dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[
1805 i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ +
1806 4] * dy[i__ + 4];
1807 /* L50: */
1808 }
1809 L60:
1810 ret_val = dtemp;
1811 return ret_val;
1812 } /* ddot_ */
1813
dnrm2_(integer * n,doublereal * x,integer * incx)1814 doublereal dnrm2_(integer *n, doublereal *x, integer *incx)
1815 {
1816 /* System generated locals */
1817 integer i__1, i__2;
1818 doublereal ret_val, d__1;
1819
1820 /* Builtin functions */
1821 double sqrt(doublereal);
1822
1823 /* Local variables */
1824 static doublereal norm, scale, absxi;
1825 static integer ix;
1826 static doublereal ssq;
1827
1828 /* .. Scalar Arguments .. */
1829 /* .. Array Arguments .. */
1830 /* .. */
1831
1832 /* DNRM2 returns the euclidean norm of a vector via the function */
1833 /* name, so that */
1834
1835 /* DNRM2 := sqrt( x'*x ) */
1836
1837
1838
1839 /* -- This version written on 25-October-1982. */
1840 /* Modified on 14-October-1993 to inline the call to DLASSQ. */
1841 /* Sven Hammarling, Nag Ltd. */
1842
1843
1844 /* .. Parameters .. */
1845 /* .. Local Scalars .. */
1846 /* .. Intrinsic Functions .. */
1847 /* .. */
1848 /* .. Executable Statements .. */
1849 /* Parameter adjustments */
1850 --x;
1851
1852 /* Function Body */
1853 if (*n < 1 || *incx < 1) {
1854 norm = 0.;
1855 } else if (*n == 1) {
1856 norm = abs(x[1]);
1857 } else {
1858 scale = 0.;
1859 ssq = 1.;
1860 /* The following loop is equivalent to this call to the LAPACK */
1861 /* auxiliary routine: */
1862 /* CALL DLASSQ( N, X, INCX, SCALE, SSQ ) */
1863
1864 i__1 = (*n - 1) * *incx + 1;
1865 i__2 = *incx;
1866 for (ix = 1; i__2 < 0 ? ix >= i__1 : ix <= i__1; ix += i__2) {
1867 if (x[ix] != 0.) {
1868 absxi = (d__1 = x[ix], abs(d__1));
1869 if (scale < absxi) {
1870 /* Computing 2nd power */
1871 d__1 = scale / absxi;
1872 ssq = ssq * (d__1 * d__1) + 1.;
1873 scale = absxi;
1874 } else {
1875 /* Computing 2nd power */
1876 d__1 = absxi / scale;
1877 ssq += d__1 * d__1;
1878 }
1879 }
1880 /* L10: */
1881 }
1882 norm = scale * sqrt(ssq);
1883 }
1884
1885 ret_val = norm;
1886 return ret_val;
1887
1888 /* End of DNRM2. */
1889
1890 } /* dnrm2_ */
1891
drot_(integer * n,doublereal * dx,integer * incx,doublereal * dy,integer * incy,doublereal * c__,doublereal * s)1892 /* Subroutine */ int drot_(integer *n, doublereal *dx, integer *incx,
1893 doublereal *dy, integer *incy, doublereal *c__, doublereal *s)
1894 {
1895 /* System generated locals */
1896 integer i__1;
1897
1898 /* Local variables */
1899 static integer i__;
1900 static doublereal dtemp;
1901 static integer ix, iy;
1902
1903
1904 /* applies a plane rotation. */
1905 /* jack dongarra, linpack, 3/11/78. */
1906 /* modified 12/3/93, array(1) declarations changed to array(*) */
1907
1908
1909 /* Parameter adjustments */
1910 --dy;
1911 --dx;
1912
1913 /* Function Body */
1914 if (*n <= 0) {
1915 return 0;
1916 }
1917 if (*incx == 1 && *incy == 1) {
1918 goto L20;
1919 }
1920
1921 /* code for unequal increments or equal increments not equal */
1922 /* to 1 */
1923
1924 ix = 1;
1925 iy = 1;
1926 if (*incx < 0) {
1927 ix = (-(*n) + 1) * *incx + 1;
1928 }
1929 if (*incy < 0) {
1930 iy = (-(*n) + 1) * *incy + 1;
1931 }
1932 i__1 = *n;
1933 for (i__ = 1; i__ <= i__1; ++i__) {
1934 dtemp = *c__ * dx[ix] + *s * dy[iy];
1935 dy[iy] = *c__ * dy[iy] - *s * dx[ix];
1936 dx[ix] = dtemp;
1937 ix += *incx;
1938 iy += *incy;
1939 /* L10: */
1940 }
1941 return 0;
1942
1943 /* code for both increments equal to 1 */
1944
1945 L20:
1946 i__1 = *n;
1947 for (i__ = 1; i__ <= i__1; ++i__) {
1948 dtemp = *c__ * dx[i__] + *s * dy[i__];
1949 dy[i__] = *c__ * dy[i__] - *s * dx[i__];
1950 dx[i__] = dtemp;
1951 /* L30: */
1952 }
1953 return 0;
1954 } /* drot_ */
1955
drotg_(doublereal * da,doublereal * db,doublereal * c__,doublereal * s)1956 /* Subroutine */ int drotg_(doublereal *da, doublereal *db, doublereal *c__,
1957 doublereal *s)
1958 {
1959 /* System generated locals */
1960 doublereal d__1, d__2;
1961
1962 /* Builtin functions */
1963 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
1964
1965 /* Local variables */
1966 static doublereal r__, scale, z__, roe;
1967
1968
1969 /* construct givens plane rotation. */
1970 /* jack dongarra, linpack, 3/11/78. */
1971
1972
1973 roe = *db;
1974 if (abs(*da) > abs(*db)) {
1975 roe = *da;
1976 }
1977 scale = abs(*da) + abs(*db);
1978 if (scale != 0.) {
1979 goto L10;
1980 }
1981 *c__ = 1.;
1982 *s = 0.;
1983 r__ = 0.;
1984 z__ = 0.;
1985 goto L20;
1986 L10:
1987 /* Computing 2nd power */
1988 d__1 = *da / scale;
1989 /* Computing 2nd power */
1990 d__2 = *db / scale;
1991 r__ = scale * sqrt(d__1 * d__1 + d__2 * d__2);
1992 r__ = d_sign(&c_b23, &roe) * r__;
1993 *c__ = *da / r__;
1994 *s = *db / r__;
1995 z__ = 1.;
1996 if (abs(*da) > abs(*db)) {
1997 z__ = *s;
1998 }
1999 if (abs(*db) >= abs(*da) && *c__ != 0.) {
2000 z__ = 1. / *c__;
2001 }
2002 L20:
2003 *da = r__;
2004 *db = z__;
2005 return 0;
2006 } /* drotg_ */
2007
dscal_(integer * n,doublereal * da,doublereal * dx,integer * incx)2008 /* Subroutine */ int dscal_(integer *n, doublereal *da, doublereal *dx,
2009 integer *incx)
2010 {
2011 /* System generated locals */
2012 integer i__1, i__2;
2013
2014 /* Local variables */
2015 static integer i__, m, nincx, mp1;
2016
2017
2018 /* scales a vector by a constant. */
2019 /* uses unrolled loops for increment equal to one. */
2020 /* jack dongarra, linpack, 3/11/78. */
2021 /* modified 3/93 to return if incx .le. 0. */
2022 /* modified 12/3/93, array(1) declarations changed to array(*) */
2023
2024
2025 /* Parameter adjustments */
2026 --dx;
2027
2028 /* Function Body */
2029 if (*n <= 0 || *incx <= 0) {
2030 return 0;
2031 }
2032 if (*incx == 1) {
2033 goto L20;
2034 }
2035
2036 /* code for increment not equal to 1 */
2037
2038 nincx = *n * *incx;
2039 i__1 = nincx;
2040 i__2 = *incx;
2041 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
2042 dx[i__] = *da * dx[i__];
2043 /* L10: */
2044 }
2045 return 0;
2046
2047 /* code for increment equal to 1 */
2048
2049
2050 /* clean-up loop */
2051
2052 L20:
2053 m = *n % 5;
2054 if (m == 0) {
2055 goto L40;
2056 }
2057 i__2 = m;
2058 for (i__ = 1; i__ <= i__2; ++i__) {
2059 dx[i__] = *da * dx[i__];
2060 /* L30: */
2061 }
2062 if (*n < 5) {
2063 return 0;
2064 }
2065 L40:
2066 mp1 = m + 1;
2067 i__2 = *n;
2068 for (i__ = mp1; i__ <= i__2; i__ += 5) {
2069 dx[i__] = *da * dx[i__];
2070 dx[i__ + 1] = *da * dx[i__ + 1];
2071 dx[i__ + 2] = *da * dx[i__ + 2];
2072 dx[i__ + 3] = *da * dx[i__ + 3];
2073 dx[i__ + 4] = *da * dx[i__ + 4];
2074 /* L50: */
2075 }
2076 return 0;
2077 } /* dscal_ */
2078
dswap_(integer * n,doublereal * dx,integer * incx,doublereal * dy,integer * incy)2079 /* Subroutine */ int dswap_(integer *n, doublereal *dx, integer *incx,
2080 doublereal *dy, integer *incy)
2081 {
2082 /* System generated locals */
2083 integer i__1;
2084
2085 /* Local variables */
2086 static integer i__, m;
2087 static doublereal dtemp;
2088 static integer ix, iy, mp1;
2089
2090
2091 /* interchanges two vectors. */
2092 /* uses unrolled loops for increments equal one. */
2093 /* jack dongarra, linpack, 3/11/78. */
2094 /* modified 12/3/93, array(1) declarations changed to array(*) */
2095
2096
2097 /* Parameter adjustments */
2098 --dy;
2099 --dx;
2100
2101 /* Function Body */
2102 if (*n <= 0) {
2103 return 0;
2104 }
2105 if (*incx == 1 && *incy == 1) {
2106 goto L20;
2107 }
2108
2109 /* code for unequal increments or equal increments not equal */
2110 /* to 1 */
2111
2112 ix = 1;
2113 iy = 1;
2114 if (*incx < 0) {
2115 ix = (-(*n) + 1) * *incx + 1;
2116 }
2117 if (*incy < 0) {
2118 iy = (-(*n) + 1) * *incy + 1;
2119 }
2120 i__1 = *n;
2121 for (i__ = 1; i__ <= i__1; ++i__) {
2122 dtemp = dx[ix];
2123 dx[ix] = dy[iy];
2124 dy[iy] = dtemp;
2125 ix += *incx;
2126 iy += *incy;
2127 /* L10: */
2128 }
2129 return 0;
2130
2131 /* code for both increments equal to 1 */
2132
2133
2134 /* clean-up loop */
2135
2136 L20:
2137 m = *n % 3;
2138 if (m == 0) {
2139 goto L40;
2140 }
2141 i__1 = m;
2142 for (i__ = 1; i__ <= i__1; ++i__) {
2143 dtemp = dx[i__];
2144 dx[i__] = dy[i__];
2145 dy[i__] = dtemp;
2146 /* L30: */
2147 }
2148 if (*n < 3) {
2149 return 0;
2150 }
2151 L40:
2152 mp1 = m + 1;
2153 i__1 = *n;
2154 for (i__ = mp1; i__ <= i__1; i__ += 3) {
2155 dtemp = dx[i__];
2156 dx[i__] = dy[i__];
2157 dy[i__] = dtemp;
2158 dtemp = dx[i__ + 1];
2159 dx[i__ + 1] = dy[i__ + 1];
2160 dy[i__ + 1] = dtemp;
2161 dtemp = dx[i__ + 2];
2162 dx[i__ + 2] = dy[i__ + 2];
2163 dy[i__ + 2] = dtemp;
2164 /* L50: */
2165 }
2166 return 0;
2167 } /* dswap_ */
2168
dsvdc_(doublereal * x,integer * ldx,integer * n,integer * p,doublereal * s,doublereal * e,doublereal * u,integer * ldu,doublereal * v,integer * ldv,doublereal * work,integer * job,integer * info)2169 /* Subroutine */ int dsvdc_(doublereal *x, integer *ldx, integer *n, integer *
2170 p, doublereal *s, doublereal *e, doublereal *u, integer *ldu,
2171 doublereal *v, integer *ldv, doublereal *work, integer *job, integer *
2172 info)
2173 {
2174 /* System generated locals */
2175 integer x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2,
2176 i__3;
2177 doublereal d__1, d__2, d__3, d__4, d__5, d__6, d__7;
2178
2179 /* Builtin functions */
2180 double d_sign(doublereal *, doublereal *), sqrt(doublereal);
2181
2182 /* Local variables */
2183 static integer kase;
2184 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
2185 integer *);
2186 static integer jobu, iter;
2187 extern /* Subroutine */ int drot_(integer *, doublereal *, integer *,
2188 doublereal *, integer *, doublereal *, doublereal *);
2189 static doublereal test;
2190 extern doublereal dnrm2_(integer *, doublereal *, integer *);
2191 static integer nctp1;
2192 static doublereal b, c__;
2193 static integer nrtp1;
2194 static doublereal f, g;
2195 static integer i__, j, k, l, m;
2196 static doublereal t, scale;
2197 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
2198 integer *);
2199 static doublereal shift;
2200 extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *,
2201 doublereal *, integer *), drotg_(doublereal *, doublereal *,
2202 doublereal *, doublereal *);
2203 static integer maxit;
2204 extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
2205 integer *, doublereal *, integer *);
2206 static logical wantu, wantv;
2207 static doublereal t1, ztest, el;
2208 static integer kk;
2209 static doublereal cs;
2210 static integer ll, mm, ls;
2211 static doublereal sl;
2212 static integer lu;
2213 static doublereal sm, sn;
2214 static integer lm1, mm1, lp1, mp1, nct, ncu, lls, nrt;
2215 static doublereal emm1, smm1;
2216
2217
2218
2219 /* dsvdc is a subroutine to reduce a double precision nxp matrix x */
2220 /* by orthogonal transformations u and v to diagonal form. the */
2221 /* diagonal elements s(i) are the singular values of x. the */
2222 /* columns of u are the corresponding left singular vectors, */
2223 /* and the columns of v the right singular vectors. */
2224
2225 /* on entry */
2226
2227 /* x double precision(ldx,p), where ldx.ge.n. */
2228 /* x contains the matrix whose singular value */
2229 /* decomposition is to be computed. x is */
2230 /* destroyed by dsvdc. */
2231
2232 /* ldx integer. */
2233 /* ldx is the leading dimension of the array x. */
2234
2235 /* n integer. */
2236 /* n is the number of rows of the matrix x. */
2237
2238 /* p integer. */
2239 /* p is the number of columns of the matrix x. */
2240
2241 /* ldu integer. */
2242 /* ldu is the leading dimension of the array u. */
2243 /* (see below). */
2244
2245 /* ldv integer. */
2246 /* ldv is the leading dimension of the array v. */
2247 /* (see below). */
2248
2249 /* work double precision(n). */
2250 /* work is a scratch array. */
2251
2252 /* job integer. */
2253 /* job controls the computation of the singular */
2254 /* vectors. it has the decimal expansion ab */
2255 /* with the following meaning */
2256
2257 /* a.eq.0 do not compute the left singular */
2258 /* vectors. */
2259 /* a.eq.1 return the n left singular vectors */
2260 /* in u. */
2261 /* a.ge.2 return the first omin(n,p) singular */
2262 /* vectors in u. */
2263 /* b.eq.0 do not compute the right singular */
2264 /* vectors. */
2265 /* b.eq.1 return the right singular vectors */
2266 /* in v. */
2267
2268 /* on return */
2269
2270 /* s double precision(mm), where mm=omin(n+1,p). */
2271 /* the first omin(n,p) entries of s contain the */
2272 /* singular values of x arranged in descending */
2273 /* order of magnitude. */
2274
2275 /* e double precision(p), */
2276 /* e ordinarily contains zeros. however see the */
2277 /* discussion of info for exceptions. */
2278
2279 /* u double precision(ldu,k), where ldu.ge.n. if */
2280 /* joba.eq.1 then k.eq.n, if joba.ge.2 */
2281 /* then k.eq.omin(n,p). */
2282 /* u contains the matrix of left singular vectors. */
2283 /* u is not referenced if joba.eq.0. if n.le.p */
2284 /* or if joba.eq.2, then u may be identified with x */
2285 /* in the subroutine call. */
2286
2287 /* v double precision(ldv,p), where ldv.ge.p. */
2288 /* v contains the matrix of right singular vectors. */
2289 /* v is not referenced if job.eq.0. if p.le.n, */
2290 /* then v may be identified with x in the */
2291 /* subroutine call. */
2292
2293 /* info integer. */
2294 /* the singular values (and their corresponding */
2295 /* singular vectors) s(info+1),s(info+2),...,s(m) */
2296 /* are correct (here m=omin(n,p)). thus if */
2297 /* info.eq.0, all the singular values and their */
2298 /* vectors are correct. in any event, the matrix */
2299 /* b = trans(u)*x*v is the bidiagonal matrix */
2300 /* with the elements of s on its diagonal and the */
2301 /* elements of e on its super-diagonal (trans(u) */
2302 /* is the transpose of u). thus the singular */
2303 /* values of x and b are the same. */
2304
2305 /* linpack. this version dated 08/14/78 . */
2306 /* correction made to shift 2/84. */
2307 /* g.w. stewart, university of maryland, argonne national lab. */
2308
2309 /* dsvdc uses the following functions and subprograms. */
2310
2311 /* external drot */
2312 /* blas daxpy,ddot,dscal,dswap,dnrm2,drotg */
2313 /* fortran dabs,dmax1,max0,min0,mod,dsqrt */
2314
2315 /* internal variables */
2316
2317
2318
2319 /* set the maximum number of iterations. */
2320
2321 /* Parameter adjustments */
2322 x_dim1 = *ldx;
2323 x_offset = 1 + x_dim1 * 1;
2324 x -= x_offset;
2325 --s;
2326 --e;
2327 u_dim1 = *ldu;
2328 u_offset = 1 + u_dim1 * 1;
2329 u -= u_offset;
2330 v_dim1 = *ldv;
2331 v_offset = 1 + v_dim1 * 1;
2332 v -= v_offset;
2333 --work;
2334
2335 /* Function Body */
2336 maxit = 30;
2337
2338 /* determine what is to be computed. */
2339
2340 wantu = FALSE_;
2341 wantv = FALSE_;
2342 jobu = *job % 100 / 10;
2343 ncu = *n;
2344 if (jobu > 1) {
2345 ncu = omin(*n,*p);
2346 }
2347 if (jobu != 0) {
2348 wantu = TRUE_;
2349 }
2350 if (*job % 10 != 0) {
2351 wantv = TRUE_;
2352 }
2353
2354 /* reduce x to bidiagonal form, storing the diagonal elements */
2355 /* in s and the super-diagonal elements in e. */
2356
2357 *info = 0;
2358 /* Computing MIN */
2359 i__1 = *n - 1;
2360 nct = omin(i__1,*p);
2361 /* Computing MAX */
2362 /* Computing MIN */
2363 i__3 = *p - 2;
2364 i__1 = 0, i__2 = omin(i__3,*n);
2365 nrt = omax(i__1,i__2);
2366 lu = omax(nct,nrt);
2367 if (lu < 1) {
2368 goto L170;
2369 }
2370 i__1 = lu;
2371 for (l = 1; l <= i__1; ++l) {
2372 lp1 = l + 1;
2373 if (l > nct) {
2374 goto L20;
2375 }
2376
2377 /* compute the transformation for the l-th column and */
2378 /* place the l-th diagonal in s(l). */
2379
2380 i__2 = *n - l + 1;
2381 s[l] = dnrm2_(&i__2, &x[l + l * x_dim1], &c__1d);
2382 if (s[l] == 0.) {
2383 goto L10;
2384 }
2385 if (x[l + l * x_dim1] != 0.) {
2386 s[l] = d_sign(&s[l], &x[l + l * x_dim1]);
2387 }
2388 i__2 = *n - l + 1;
2389 d__1 = 1. / s[l];
2390 dscal_(&i__2, &d__1, &x[l + l * x_dim1], &c__1d);
2391 x[l + l * x_dim1] += 1.;
2392 L10:
2393 s[l] = -s[l];
2394 L20:
2395 if (*p < lp1) {
2396 goto L50;
2397 }
2398 i__2 = *p;
2399 for (j = lp1; j <= i__2; ++j) {
2400 if (l > nct) {
2401 goto L30;
2402 }
2403 if (s[l] == 0.) {
2404 goto L30;
2405 }
2406
2407 /* apply the transformation. */
2408
2409 i__3 = *n - l + 1;
2410 t = -ddot_(&i__3, &x[l + l * x_dim1], &c__1d, &x[l + j * x_dim1], &
2411 c__1d) / x[l + l * x_dim1];
2412 i__3 = *n - l + 1;
2413 daxpy_(&i__3, &t, &x[l + l * x_dim1], &c__1d, &x[l + j * x_dim1], &
2414 c__1d);
2415 L30:
2416
2417 /* place the l-th row of x into e for the */
2418 /* subsequent calculation of the row transformation. */
2419
2420 e[j] = x[l + j * x_dim1];
2421 /* L40: */
2422 }
2423 L50:
2424 if (! wantu || l > nct) {
2425 goto L70;
2426 }
2427
2428 /* place the transformation in u for subsequent back */
2429 /* multiplication. */
2430
2431 i__2 = *n;
2432 for (i__ = l; i__ <= i__2; ++i__) {
2433 u[i__ + l * u_dim1] = x[i__ + l * x_dim1];
2434 /* L60: */
2435 }
2436 L70:
2437 if (l > nrt) {
2438 goto L150;
2439 }
2440
2441 /* compute the l-th row transformation and place the */
2442 /* l-th super-diagonal in e(l). */
2443
2444 i__2 = *p - l;
2445 e[l] = dnrm2_(&i__2, &e[lp1], &c__1d);
2446 if (e[l] == 0.) {
2447 goto L80;
2448 }
2449 if (e[lp1] != 0.) {
2450 e[l] = d_sign(&e[l], &e[lp1]);
2451 }
2452 i__2 = *p - l;
2453 d__1 = 1. / e[l];
2454 dscal_(&i__2, &d__1, &e[lp1], &c__1d);
2455 e[lp1] += 1.;
2456 L80:
2457 e[l] = -e[l];
2458 if (lp1 > *n || e[l] == 0.) {
2459 goto L120;
2460 }
2461
2462 /* apply the transformation. */
2463
2464 i__2 = *n;
2465 for (i__ = lp1; i__ <= i__2; ++i__) {
2466 work[i__] = 0.;
2467 /* L90: */
2468 }
2469 i__2 = *p;
2470 for (j = lp1; j <= i__2; ++j) {
2471 i__3 = *n - l;
2472 daxpy_(&i__3, &e[j], &x[lp1 + j * x_dim1], &c__1d, &work[lp1], &
2473 c__1d);
2474 /* L100: */
2475 }
2476 i__2 = *p;
2477 for (j = lp1; j <= i__2; ++j) {
2478 i__3 = *n - l;
2479 d__1 = -e[j] / e[lp1];
2480 daxpy_(&i__3, &d__1, &work[lp1], &c__1d, &x[lp1 + j * x_dim1], &
2481 c__1d);
2482 /* L110: */
2483 }
2484 L120:
2485 if (! wantv) {
2486 goto L140;
2487 }
2488
2489 /* place the transformation in v for subsequent */
2490 /* back multiplication. */
2491
2492 i__2 = *p;
2493 for (i__ = lp1; i__ <= i__2; ++i__) {
2494 v[i__ + l * v_dim1] = e[i__];
2495 /* L130: */
2496 }
2497 L140:
2498 L150:
2499 /* L160: */
2500 ;
2501 }
2502 L170:
2503
2504 /* set up the final bidiagonal matrix or order m. */
2505
2506 /* Computing MIN */
2507 i__1 = *p, i__2 = *n + 1;
2508 m = omin(i__1,i__2);
2509 nctp1 = nct + 1;
2510 nrtp1 = nrt + 1;
2511 if (nct < *p) {
2512 s[nctp1] = x[nctp1 + nctp1 * x_dim1];
2513 }
2514 if (*n < m) {
2515 s[m] = 0.;
2516 }
2517 if (nrtp1 < m) {
2518 e[nrtp1] = x[nrtp1 + m * x_dim1];
2519 }
2520 e[m] = 0.;
2521
2522 /* if required, generate u. */
2523
2524 if (! wantu) {
2525 goto L300;
2526 }
2527 if (ncu < nctp1) {
2528 goto L200;
2529 }
2530 i__1 = ncu;
2531 for (j = nctp1; j <= i__1; ++j) {
2532 i__2 = *n;
2533 for (i__ = 1; i__ <= i__2; ++i__) {
2534 u[i__ + j * u_dim1] = 0.;
2535 /* L180: */
2536 }
2537 u[j + j * u_dim1] = 1.;
2538 /* L190: */
2539 }
2540 L200:
2541 if (nct < 1) {
2542 goto L290;
2543 }
2544 i__1 = nct;
2545 for (ll = 1; ll <= i__1; ++ll) {
2546 l = nct - ll + 1;
2547 if (s[l] == 0.) {
2548 goto L250;
2549 }
2550 lp1 = l + 1;
2551 if (ncu < lp1) {
2552 goto L220;
2553 }
2554 i__2 = ncu;
2555 for (j = lp1; j <= i__2; ++j) {
2556 i__3 = *n - l + 1;
2557 t = -ddot_(&i__3, &u[l + l * u_dim1], &c__1d, &u[l + j * u_dim1], &
2558 c__1d) / u[l + l * u_dim1];
2559 i__3 = *n - l + 1;
2560 daxpy_(&i__3, &t, &u[l + l * u_dim1], &c__1d, &u[l + j * u_dim1], &
2561 c__1d);
2562 /* L210: */
2563 }
2564 L220:
2565 i__2 = *n - l + 1;
2566 dscal_(&i__2, &c_b79, &u[l + l * u_dim1], &c__1d);
2567 u[l + l * u_dim1] += 1.;
2568 lm1 = l - 1;
2569 if (lm1 < 1) {
2570 goto L240;
2571 }
2572 i__2 = lm1;
2573 for (i__ = 1; i__ <= i__2; ++i__) {
2574 u[i__ + l * u_dim1] = 0.;
2575 /* L230: */
2576 }
2577 L240:
2578 goto L270;
2579 L250:
2580 i__2 = *n;
2581 for (i__ = 1; i__ <= i__2; ++i__) {
2582 u[i__ + l * u_dim1] = 0.;
2583 /* L260: */
2584 }
2585 u[l + l * u_dim1] = 1.;
2586 L270:
2587 /* L280: */
2588 ;
2589 }
2590 L290:
2591 L300:
2592
2593 /* if it is required, generate v. */
2594
2595 if (! wantv) {
2596 goto L350;
2597 }
2598 i__1 = *p;
2599 for (ll = 1; ll <= i__1; ++ll) {
2600 l = *p - ll + 1;
2601 lp1 = l + 1;
2602 if (l > nrt) {
2603 goto L320;
2604 }
2605 if (e[l] == 0.) {
2606 goto L320;
2607 }
2608 i__2 = *p;
2609 for (j = lp1; j <= i__2; ++j) {
2610 i__3 = *p - l;
2611 t = -ddot_(&i__3, &v[lp1 + l * v_dim1], &c__1d, &v[lp1 + j *
2612 v_dim1], &c__1d) / v[lp1 + l * v_dim1];
2613 i__3 = *p - l;
2614 daxpy_(&i__3, &t, &v[lp1 + l * v_dim1], &c__1d, &v[lp1 + j *
2615 v_dim1], &c__1d);
2616 /* L310: */
2617 }
2618 L320:
2619 i__2 = *p;
2620 for (i__ = 1; i__ <= i__2; ++i__) {
2621 v[i__ + l * v_dim1] = 0.;
2622 /* L330: */
2623 }
2624 v[l + l * v_dim1] = 1.;
2625 /* L340: */
2626 }
2627 L350:
2628
2629 /* main iteration loop for the singular values. */
2630
2631 mm = m;
2632 iter = 0;
2633 L360:
2634
2635 /* quit if all the singular values have been found. */
2636
2637 /* ...exit */
2638 if (m == 0) {
2639 goto L620;
2640 }
2641
2642 /* if too many iterations have been performed, set */
2643 /* flag and return. */
2644
2645 if (iter < maxit) {
2646 goto L370;
2647 }
2648 *info = m;
2649 /* ......exit */
2650 goto L620;
2651 L370:
2652
2653 /* this section of the program inspects for */
2654 /* negligible elements in the s and e arrays. on */
2655 /* completion the variables kase and l are set as follows. */
2656
2657 /* kase = 1 if s(m) and e(l-1) are negligible and l.lt.m */
2658 /* kase = 2 if s(l) is negligible and l.lt.m */
2659 /* kase = 3 if e(l-1) is negligible, l.lt.m, and */
2660 /* s(l), ..., s(m) are not negligible (qr step). */
2661 /* kase = 4 if e(m-1) is negligible (convergence). */
2662
2663 i__1 = m;
2664 for (ll = 1; ll <= i__1; ++ll) {
2665 l = m - ll;
2666 /* ...exit */
2667 if (l == 0) {
2668 goto L400;
2669 }
2670 test = (d__1 = s[l], abs(d__1)) + (d__2 = s[l + 1], abs(d__2));
2671 ztest = test + (d__1 = e[l], abs(d__1));
2672 if (ztest != test) {
2673 goto L380;
2674 }
2675 e[l] = 0.;
2676 /* ......exit */
2677 goto L400;
2678 L380:
2679 /* L390: */
2680 ;
2681 }
2682 L400:
2683 if (l != m - 1) {
2684 goto L410;
2685 }
2686 kase = 4;
2687 goto L480;
2688 L410:
2689 lp1 = l + 1;
2690 mp1 = m + 1;
2691 i__1 = mp1;
2692 for (lls = lp1; lls <= i__1; ++lls) {
2693 ls = m - lls + lp1;
2694 /* ...exit */
2695 if (ls == l) {
2696 goto L440;
2697 }
2698 test = 0.;
2699 if (ls != m) {
2700 test += (d__1 = e[ls], abs(d__1));
2701 }
2702 if (ls != l + 1) {
2703 test += (d__1 = e[ls - 1], abs(d__1));
2704 }
2705 ztest = test + (d__1 = s[ls], abs(d__1));
2706 if (ztest != test) {
2707 goto L420;
2708 }
2709 s[ls] = 0.;
2710 /* ......exit */
2711 goto L440;
2712 L420:
2713 /* L430: */
2714 ;
2715 }
2716 L440:
2717 if (ls != l) {
2718 goto L450;
2719 }
2720 kase = 3;
2721 goto L470;
2722 L450:
2723 if (ls != m) {
2724 goto L460;
2725 }
2726 kase = 1;
2727 goto L470;
2728 L460:
2729 kase = 2;
2730 l = ls;
2731 L470:
2732 L480:
2733 ++l;
2734
2735 /* perform the task indicated by kase. */
2736
2737 switch (kase) {
2738 case 1: goto L490;
2739 case 2: goto L520;
2740 case 3: goto L540;
2741 case 4: goto L570;
2742 }
2743
2744 /* deflate negligible s(m). */
2745
2746 L490:
2747 mm1 = m - 1;
2748 f = e[m - 1];
2749 e[m - 1] = 0.;
2750 i__1 = mm1;
2751 for (kk = l; kk <= i__1; ++kk) {
2752 k = mm1 - kk + l;
2753 t1 = s[k];
2754 drotg_(&t1, &f, &cs, &sn);
2755 s[k] = t1;
2756 if (k == l) {
2757 goto L500;
2758 }
2759 f = -sn * e[k - 1];
2760 e[k - 1] = cs * e[k - 1];
2761 L500:
2762 if (wantv) {
2763 drot_(p, &v[k * v_dim1 + 1], &c__1d, &v[m * v_dim1 + 1], &c__1d, &
2764 cs, &sn);
2765 }
2766 /* L510: */
2767 }
2768 goto L610;
2769
2770 /* split at negligible s(l). */
2771
2772 L520:
2773 f = e[l - 1];
2774 e[l - 1] = 0.;
2775 i__1 = m;
2776 for (k = l; k <= i__1; ++k) {
2777 t1 = s[k];
2778 drotg_(&t1, &f, &cs, &sn);
2779 s[k] = t1;
2780 f = -sn * e[k];
2781 e[k] = cs * e[k];
2782 if (wantu) {
2783 drot_(n, &u[k * u_dim1 + 1], &c__1d, &u[(l - 1) * u_dim1 + 1], &
2784 c__1d, &cs, &sn);
2785 }
2786 /* L530: */
2787 }
2788 goto L610;
2789
2790 /* perform one qr step. */
2791
2792 L540:
2793
2794 /* calculate the shift. */
2795
2796 /* Computing MAX */
2797 d__6 = (d__1 = s[m], abs(d__1)), d__7 = (d__2 = s[m - 1], abs(d__2)),
2798 d__6 = omax(d__6,d__7), d__7 = (d__3 = e[m - 1], abs(d__3)), d__6 =
2799 omax(d__6,d__7), d__7 = (d__4 = s[l], abs(d__4)), d__6 = omax(d__6,
2800 d__7), d__7 = (d__5 = e[l], abs(d__5));
2801 scale = omax(d__6,d__7);
2802 sm = s[m] / scale;
2803 smm1 = s[m - 1] / scale;
2804 emm1 = e[m - 1] / scale;
2805 sl = s[l] / scale;
2806 el = e[l] / scale;
2807 /* Computing 2nd power */
2808 d__1 = emm1;
2809 b = ((smm1 + sm) * (smm1 - sm) + d__1 * d__1) / 2.;
2810 /* Computing 2nd power */
2811 d__1 = sm * emm1;
2812 c__ = d__1 * d__1;
2813 shift = 0.;
2814 if (b == 0. && c__ == 0.) {
2815 goto L550;
2816 }
2817 /* Computing 2nd power */
2818 d__1 = b;
2819 shift = sqrt(d__1 * d__1 + c__);
2820 if (b < 0.) {
2821 shift = -shift;
2822 }
2823 shift = c__ / (b + shift);
2824 L550:
2825 f = (sl + sm) * (sl - sm) + shift;
2826 g = sl * el;
2827
2828 /* chase zeros. */
2829
2830 mm1 = m - 1;
2831 i__1 = mm1;
2832 for (k = l; k <= i__1; ++k) {
2833 drotg_(&f, &g, &cs, &sn);
2834 if (k != l) {
2835 e[k - 1] = f;
2836 }
2837 f = cs * s[k] + sn * e[k];
2838 e[k] = cs * e[k] - sn * s[k];
2839 g = sn * s[k + 1];
2840 s[k + 1] = cs * s[k + 1];
2841 if (wantv) {
2842 drot_(p, &v[k * v_dim1 + 1], &c__1d, &v[(k + 1) * v_dim1 + 1], &
2843 c__1d, &cs, &sn);
2844 }
2845 drotg_(&f, &g, &cs, &sn);
2846 s[k] = f;
2847 f = cs * e[k] + sn * s[k + 1];
2848 s[k + 1] = -sn * e[k] + cs * s[k + 1];
2849 g = sn * e[k + 1];
2850 e[k + 1] = cs * e[k + 1];
2851 if (wantu && k < *n) {
2852 drot_(n, &u[k * u_dim1 + 1], &c__1d, &u[(k + 1) * u_dim1 + 1], &
2853 c__1d, &cs, &sn);
2854 }
2855 /* L560: */
2856 }
2857 e[m - 1] = f;
2858 ++iter;
2859 goto L610;
2860
2861 /* convergence. */
2862
2863 L570:
2864
2865 /* make the singular value positive. */
2866
2867 if (s[l] >= 0.) {
2868 goto L580;
2869 }
2870 s[l] = -s[l];
2871 if (wantv) {
2872 dscal_(p, &c_b79, &v[l * v_dim1 + 1], &c__1d);
2873 }
2874 L580:
2875
2876 /* order the singular value. */
2877
2878 L590:
2879 if (l == mm) {
2880 goto L600;
2881 }
2882 /* ...exit */
2883 if (s[l] >= s[l + 1]) {
2884 goto L600;
2885 }
2886 t = s[l];
2887 s[l] = s[l + 1];
2888 s[l + 1] = t;
2889 if (wantv && l < *p) {
2890 dswap_(p, &v[l * v_dim1 + 1], &c__1d, &v[(l + 1) * v_dim1 + 1], &c__1d);
2891 }
2892 if (wantu && l < *n) {
2893 dswap_(n, &u[l * u_dim1 + 1], &c__1d, &u[(l + 1) * u_dim1 + 1], &c__1d);
2894 }
2895 ++l;
2896 goto L590;
2897 L600:
2898 iter = 0;
2899 --m;
2900 L610:
2901 goto L360;
2902 L620:
2903 return 0;
2904 } /* dsvdc_ */
2905
2906
2907 #ifdef KR_headers
2908 double d_sign(a,b) doublereal *a, *b;
2909 #else
d_sign(doublereal * a,doublereal * b)2910 double d_sign(doublereal *a, doublereal *b)
2911 #endif
2912 {
2913 double x;
2914 x = (*a >= 0 ? *a : - *a);
2915 return( *b >= 0 ? x : -x);
2916 }
2917
2918 #ifdef __cplusplus
2919 }
2920 #endif
2921