/* Singular Value Decomposition from LINPACK */
#include "linalg.h"
VOID linpack_dsvdc P13C(double *, x,
int, ldx,
int, n,
int, p,
double *, s,
double *, e,
double *, u,
int, ldu,
double *, v,
int, ldv,
double *, work,
int, job,
int *, info)
{
/* System generated locals */
int x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2, i__3;
double d__1, d__2, d__3, d__4, d__5, d__6, d__7;
/* Local variables */
int kase;
int jobu, iter;
double test;
int nctp1;
double b, c;
int nrtp1;
double f, g;
int i, j, k, l, m;
double t, scale;
double shift;
int maxit;
int wantu, wantv;
double t1, ztest, el;
int kk;
double cs;
int ll, mm, ls;
double sl;
int lu;
double sm, sn;
int lm1, mm1, lp1, mp1, nct, ncu, lls, nrt;
double emm1, smm1;
/* dsvdc is a subroutine to reduce a double precision nxp matrix x */
/* by orthogonal transformations u and v to diagonal form. the */
/* diagonal elements s(i) are the singular values of x. the */
/* columns of u are the corresponding left singular vectors, */
/* and the columns of v the right singular vectors. */
/* on entry */
/* x double precision(ldx,p), where ldx.ge.n. */
/* x contains the matrix whose singular value */
/* decomposition is to be computed. x is */
/* destroyed by dsvdc. */
/* ldx integer. */
/* ldx is the leading dimension of the array x. */
/* n integer. */
/* n is the number of rows of the matrix x. */
/* p integer. */
/* p is the number of columns of the matrix x. */
/* ldu integer. */
/* ldu is the leading dimension of the array u. */
/* (see below). */
/* ldv integer. */
/* ldv is the leading dimension of the array v. */
/* (see below). */
/* work double precision(n). */
/* work is a scratch array. */
/* job integer. */
/* job controls the computation of the singular */
/* vectors. it has the decimal expansion ab */
/* with the following meaning */
/* a.eq.0 do not compute the left singular */
/* vectors. */
/* a.eq.1 return the n left singular vectors */
/* in u. */
/* a.ge.2 return the first min(n,p) singular */
/* vectors in u. */
/* b.eq.0 do not compute the right singular */
/* vectors. */
/* b.eq.1 return the right singular vectors */
/* in v. */
/* on return */
/* s double precision(mm), where mm=min(n+1,p). */
/* the first min(n,p) entries of s contain the */
/* singular values of x arranged in descending */
/* order of magnitude. */
/* e double precision(p), */
/* e ordinarily contains zeros. however see the */
/* discussion of info for exceptions. */
/* u double precision(ldu,k), where ldu.ge.n. if */
/* joba.eq.1 then k.eq.n, if joba.ge.2 */
/* then k.eq.min(n,p). */
/* u contains the matrix of left singular vectors. */
/* u is not referenced if joba.eq.0. if n.le.p */
/* or if joba.eq.2, then u may be identified with x */
/* in the subroutine call. */
/* v double precision(ldv,p), where ldv.ge.p. */
/* v contains the matrix of right singular vectors. */
/* v is not referenced if job.eq.0. if p.le.n, */
/* then v may be identified with x in the */
/* subroutine call. */
/* info integer. */
/* the singular values (and their corresponding */
/* singular vectors) s(info+1),s(info+2),...,s(m) */
/* are correct (here m=min(n,p)). thus if */
/* info.eq.0, all the singular values and their */
/* vectors are correct. in any event, the matrix */
/* b = trans(u)*x*v is the bidiagonal matrix */
/* with the elements of s on its diagonal and the */
/* elements of e on its super-diagonal (trans(u) */
/* is the transpose of u). thus the singular */
/* values of x and b are the same. */
/* linpack. this version dated 08/14/78 . */
/* correction made to shift 2/84. */
/* g.w. stewart, university of maryland, argonne national lab. */
/* dsvdc uses the following functions and subprograms. */
/* external drot */
/* blas daxpy,ddot,dscal,dswap,dnrm2,drotg */
/* fortran dabs,dmax1,max0,min0,mod,dsqrt */
/* Parameter adjustments */
x_dim1 = ldx;
x_offset = x_dim1 + 1;
x -= x_offset;
--s;
--e;
u_dim1 = ldu;
u_offset = u_dim1 + 1;
u -= u_offset;
v_dim1 = ldv;
v_offset = v_dim1 + 1;
v -= v_offset;
--work;
/* keep compiler happy */
l = ls = 0;
/* set the maximum number of iterations. */
maxit = 30;
/* determine what is to be computed. */
wantu = FALSE;
wantv = FALSE;
jobu = job % 100 / 10;
ncu = n;
if (jobu > 1) {
ncu = min(n,p);
}
if (jobu != 0) {
wantu = TRUE;
}
if (job % 10 != 0) {
wantv = TRUE;
}
/* reduce x to bidiagonal form, storing the diagonal elements */
/* in s and the super-diagonal elements in e. */
*info = 0;
/* Computing MIN */
i__1 = n - 1;
nct = min(i__1,p);
/* Computing MAX */
/* Computing MIN */
i__3 = p - 2;
i__1 = 0, i__2 = min(i__3,n);
nrt = max(i__1,i__2);
lu = max(nct,nrt);
if (lu < 1) {
goto L170;
}
i__1 = lu;
for (l = 1; l <= i__1; ++l) {
lp1 = l + 1;
if (l > nct) {
goto L20;
}
/* compute the transformation for the l-th column and */
/* place the l-th diagonal in s(l). */
i__2 = n - l + 1;
s[l] = blas_dnrm2(i__2, &x[l + l * x_dim1], 1);
if (s[l] == 0.) {
goto L10;
}
if (x[l + l * x_dim1] != 0.) {
s[l] = d_sign(&s[l], &x[l + l * x_dim1]);
}
i__2 = n - l + 1;
d__1 = 1. / s[l];
blas_dscal(i__2, d__1, &x[l + l * x_dim1], 1);
x[l + l * x_dim1] += 1.;
L10:
s[l] = -s[l];
L20:
if (p < lp1) {
goto L50;
}
i__2 = p;
for (j = lp1; j <= i__2; ++j) {
if (l > nct) {
goto L30;
}
if (s[l] == 0.) {
goto L30;
}
/* apply the transformation. */
i__3 = n - l + 1;
t = -blas_ddot(i__3, &x[l + l * x_dim1], 1, &x[l + j * x_dim1],
1) / x[l + l * x_dim1];
i__3 = n - l + 1;
blas_daxpy(i__3, t, &x[l + l * x_dim1], 1, &x[l + j * x_dim1], 1);
L30:
/* place the l-th row of x into e for the */
/* subsequent calculation of the row transformation. */
e[j] = x[l + j * x_dim1];
/* L40: */
}
L50:
if (! wantu || l > nct) {
goto L70;
}
/* place the transformation in u for subsequent back */
/* multiplication. */
i__2 = n;
for (i = l; i <= i__2; ++i) {
u[i + l * u_dim1] = x[i + l * x_dim1];
/* L60: */
}
L70:
if (l > nrt) {
goto L150;
}
/* compute the l-th row transformation and place the */
/* l-th super-diagonal in e(l). */
i__2 = p - l;
e[l] = blas_dnrm2(i__2, &e[lp1], 1);
if (e[l] == 0.) {
goto L80;
}
if (e[lp1] != 0.) {
e[l] = d_sign(&e[l], &e[lp1]);
}
i__2 = p - l;
d__1 = 1. / e[l];
blas_dscal(i__2, d__1, &e[lp1], 1);
e[lp1] += 1.;
L80:
e[l] = -e[l];
if (lp1 > n || e[l] == 0.) {
goto L120;
}
/* apply the transformation. */
i__2 = n;
for (i = lp1; i <= i__2; ++i) {
work[i] = 0.;
/* L90: */
}
i__2 = p;
for (j = lp1; j <= i__2; ++j) {
i__3 = n - l;
blas_daxpy(i__3, e[j], &x[lp1 + j * x_dim1], 1, &work[lp1], 1);
/* L100: */
}
i__2 = p;
for (j = lp1; j <= i__2; ++j) {
i__3 = n - l;
d__1 = -e[j] / e[lp1];
blas_daxpy(i__3, d__1, &work[lp1], 1, &x[lp1 + j * x_dim1], 1);
/* L110: */
}
L120:
if (! wantv) {
goto L140;
}
/* place the transformation in v for subsequent */
/* back multiplication. */
i__2 = p;
for (i = lp1; i <= i__2; ++i) {
v[i + l * v_dim1] = e[i];
/* L130: */
}
L140:
L150:
/* L160: */
;
}
L170:
/* set up the final bidiagonal matrix or order m. */
/* Computing MIN */
i__1 = p, i__2 = n + 1;
m = min(i__1,i__2);
nctp1 = nct + 1;
nrtp1 = nrt + 1;
if (nct < p) {
s[nctp1] = x[nctp1 + nctp1 * x_dim1];
}
if (n < m) {
s[m] = 0.;
}
if (nrtp1 < m) {
e[nrtp1] = x[nrtp1 + m * x_dim1];
}
e[m] = 0.;
/* if required, generate u. */
if (! wantu) {
goto L300;
}
if (ncu < nctp1) {
goto L200;
}
i__1 = ncu;
for (j = nctp1; j <= i__1; ++j) {
i__2 = n;
for (i = 1; i <= i__2; ++i) {
u[i + j * u_dim1] = 0.;
/* L180: */
}
u[j + j * u_dim1] = 1.;
/* L190: */
}
L200:
if (nct < 1) {
goto L290;
}
i__1 = nct;
for (ll = 1; ll <= i__1; ++ll) {
l = nct - ll + 1;
if (s[l] == 0.) {
goto L250;
}
lp1 = l + 1;
if (ncu < lp1) {
goto L220;
}
i__2 = ncu;
for (j = lp1; j <= i__2; ++j) {
i__3 = n - l + 1;
t = -blas_ddot(i__3, &u[l + l * u_dim1], 1, &u[l + j * u_dim1],
1) / u[l + l * u_dim1];
i__3 = n - l + 1;
blas_daxpy(i__3, t, &u[l + l * u_dim1], 1, &u[l + j * u_dim1], 1);
/* L210: */
}
L220:
i__2 = n - l + 1;
blas_dscal(i__2, -1.0, &u[l + l * u_dim1], 1);
u[l + l * u_dim1] += 1.;
lm1 = l - 1;
if (lm1 < 1) {
goto L240;
}
i__2 = lm1;
for (i = 1; i <= i__2; ++i) {
u[i + l * u_dim1] = 0.;
/* L230: */
}
L240:
goto L270;
L250:
i__2 = n;
for (i = 1; i <= i__2; ++i) {
u[i + l * u_dim1] = 0.;
/* L260: */
}
u[l + l * u_dim1] = 1.;
L270:
/* L280: */
;
}
L290:
L300:
/* if it is required, generate v. */
if (! wantv) {
goto L350;
}
i__1 = p;
for (ll = 1; ll <= i__1; ++ll) {
l = p - ll + 1;
lp1 = l + 1;
if (l > nrt) {
goto L320;
}
if (e[l] == 0.) {
goto L320;
}
i__2 = p;
for (j = lp1; j <= i__2; ++j) {
i__3 = p - l;
t = -blas_ddot(i__3, &v[lp1 + l * v_dim1], 1, &v[lp1 + j *
v_dim1], 1) / v[lp1 + l * v_dim1];
i__3 = p - l;
blas_daxpy(i__3, t, &v[lp1 + l * v_dim1], 1, &v[lp1 + j *
v_dim1], 1);
/* L310: */
}
L320:
i__2 = p;
for (i = 1; i <= i__2; ++i) {
v[i + l * v_dim1] = 0.;
/* L330: */
}
v[l + l * v_dim1] = 1.;
/* L340: */
}
L350:
/* main iteration loop for the singular values. */
mm = m;
iter = 0;
L360:
/* quit if all the singular values have been found. */
/* ...exit */
if (m == 0) {
goto L620;
}
/* if too many iterations have been performed, set */
/* flag and return. */
if (iter < maxit) {
goto L370;
}
*info = m;
/* ......exit */
goto L620;
L370:
/* this section of the program inspects for */
/* negligible elements in the s and e arrays. on */
/* completion the variables kase and l are set as follows. */
/* kase = 1 if s(m) and e(l-1) are negligible and l.lt.m */
/* kase = 2 if s(l) is negligible and l.lt.m */
/* kase = 3 if e(l-1) is negligible, l.lt.m, and */
/* s(l), ..., s(m) are not negligible (qr step). */
/* kase = 4 if e(m-1) is negligible (convergence). */
i__1 = m;
for (ll = 1; ll <= i__1; ++ll) {
l = m - ll;
/* ...exit */
if (l == 0) {
goto L400;
}
test = (d__1 = s[l], abs(d__1)) + (d__2 = s[l + 1], abs(d__2));
ztest = test + (d__1 = e[l], abs(d__1));
if (ztest != test) {
goto L380;
}
e[l] = 0.;
/* ......exit */
goto L400;
L380:
/* L390: */
;
}
L400:
if (l != m - 1) {
goto L410;
}
kase = 4;
goto L480;
L410:
lp1 = l + 1;
mp1 = m + 1;
i__1 = mp1;
for (lls = lp1; lls <= i__1; ++lls) {
ls = m - lls + lp1;
/* ...exit */
if (ls == l) {
goto L440;
}
test = 0.;
if (ls != m) {
test += (d__1 = e[ls], abs(d__1));
}
if (ls != l + 1) {
test += (d__1 = e[ls - 1], abs(d__1));
}
ztest = test + (d__1 = s[ls], abs(d__1));
if (ztest != test) {
goto L420;
}
s[ls] = 0.;
/* ......exit */
goto L440;
L420:
/* L430: */
;
}
L440:
if (ls != l) {
goto L450;
}
kase = 3;
goto L470;
L450:
if (ls != m) {
goto L460;
}
kase = 1;
goto L470;
L460:
kase = 2;
l = ls;
L470:
L480:
++l;
/* perform the task indicated by kase. */
switch ((int)kase) {
case 1: goto L490;
case 2: goto L520;
case 3: goto L540;
case 4: goto L570;
}
/* deflate negligible s(m). */
L490:
mm1 = m - 1;
f = e[m - 1];
e[m - 1] = 0.;
i__1 = mm1;
for (kk = l; kk <= i__1; ++kk) {
k = mm1 - kk + l;
t1 = s[k];
blas_drotg(&t1, &f, &cs, &sn);
s[k] = t1;
if (k == l) {
goto L500;
}
f = -sn * e[k - 1];
e[k - 1] = cs * e[k - 1];
L500:
if (wantv) {
blas_drot(p, &v[k * v_dim1 + 1], 1, &v[m * v_dim1 + 1], 1, cs, sn);
}
/* L510: */
}
goto L610;
/* split at negligible s(l). */
L520:
f = e[l - 1];
e[l - 1] = 0.;
i__1 = m;
for (k = l; k <= i__1; ++k) {
t1 = s[k];
blas_drotg(&t1, &f, &cs, &sn);
s[k] = t1;
f = -sn * e[k];
e[k] = cs * e[k];
if (wantu) {
blas_drot(n, &u[k * u_dim1 + 1], 1, &u[(l - 1) * u_dim1 + 1], 1, cs, sn);
}
/* L530: */
}
goto L610;
/* perform one qr step. */
L540:
/* calculate the shift. */
/* Computing MAX */
d__6 = (d__1 = s[m], abs(d__1)), d__7 = (d__2 = s[m - 1], abs(d__2)),
d__6 = max(d__6,d__7), d__7 = (d__3 = e[m - 1], abs(d__3)), d__6 =
max(d__6,d__7), d__7 = (d__4 = s[l], abs(d__4)), d__6 = max(d__6,
d__7), d__7 = (d__5 = e[l], abs(d__5));
scale = max(d__6,d__7);
sm = s[m] / scale;
smm1 = s[m - 1] / scale;
emm1 = e[m - 1] / scale;
sl = s[l] / scale;
el = e[l] / scale;
/* Computing 2nd power */
d__1 = emm1;
b = ((smm1 + sm) * (smm1 - sm) + d__1 * d__1) / 2.;
/* Computing 2nd power */
d__1 = sm * emm1;
c = d__1 * d__1;
shift = 0.;
if (b == 0. && c == 0.) {
goto L550;
}
/* Computing 2nd power */
d__1 = b;
shift = sqrt(d__1 * d__1 + c);
if (b < 0.) {
shift = -shift;
}
shift = c / (b + shift);
L550:
f = (sl + sm) * (sl - sm) + shift;
g = sl * el;
/* chase zeros. */
mm1 = m - 1;
i__1 = mm1;
for (k = l; k <= i__1; ++k) {
blas_drotg(&f, &g, &cs, &sn);
if (k != l) {
e[k - 1] = f;
}
f = cs * s[k] + sn * e[k];
e[k] = cs * e[k] - sn * s[k];
g = sn * s[k + 1];
s[k + 1] = cs * s[k + 1];
if (wantv) {
blas_drot(p, &v[k * v_dim1 + 1], 1, &v[(k + 1) * v_dim1 + 1], 1, cs, sn);
}
blas_drotg(&f, &g, &cs, &sn);
s[k] = f;
f = cs * e[k] + sn * s[k + 1];
s[k + 1] = -sn * e[k] + cs * s[k + 1];
g = sn * e[k + 1];
e[k + 1] = cs * e[k + 1];
if (wantu && k < n) {
blas_drot(n, &u[k * u_dim1 + 1], 1, &u[(k + 1) * u_dim1 + 1], 1, cs, sn);
}
/* L560: */
}
e[m - 1] = f;
++iter;
goto L610;
/* convergence. */
L570:
/* make the singular value positive. */
if (s[l] >= 0.) {
goto L580;
}
s[l] = -s[l];
if (wantv) {
blas_dscal(p, -1.0, &v[l * v_dim1 + 1], 1);
}
L580:
/* order the singular value. */
L590:
if (l == mm) {
goto L600;
}
/* ...exit */
if (s[l] >= s[l + 1]) {
goto L600;
}
t = s[l];
s[l] = s[l + 1];
s[l + 1] = t;
if (wantv && l < p) {
blas_dswap(p, &v[l * v_dim1 + 1], 1, &v[(l + 1) * v_dim1 + 1], 1);
}
if (wantu && l < n) {
blas_dswap(n, &u[l * u_dim1 + 1], 1, &u[(l + 1) * u_dim1 + 1], 1);
}
++l;
goto L590;
L600:
iter = 0;
--m;
L610:
goto L360;
L620:
return;
}
VOID linpack_zsvdc P13C(dcomplex *, x,
int, ldx,
int, n,
int, p,
dcomplex *, s,
dcomplex *, e,
dcomplex *, u,
int, ldu,
dcomplex *, v,
int, ldv,
dcomplex *, work,
int, job,
int *, info)
{
/* System generated locals */
int x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2,
i__3, i__4, i__5;
double r__1, r__2;
double d__1, d__2, d__3, d__4;
dcomplex z__1, z__2, z__3;
static dcomplex c_b10 = {1.,0.};
static dcomplex c_b60 = {-1.,0.};
/* Local variables */
int kase, jobu, iter;
double test;
int nctp1;
double b, c;
int nrtp1;
double f, g;
int i, j, k, l, m;
dcomplex r, t;
double scale;
double shift;
int maxit;
int wantu, wantv;
double t1, ztest;
double el;
int kk;
double cs;
int ll, mm, ls;
double sl;
int lu;
double sm, sn;
int lm1, mm1, lp1, mp1, nct, ncu, lls, nrt;
double emm1, smm1;
/* zsvdc is a subroutine to reduce a complex*16 nxp matrix x by */
/* unitary transformations u and v to diagonal form. the */
/* diagonal elements s(i) are the singular values of x. the */
/* columns of u are the corresponding left singular vectors, */
/* and the columns of v the right singular vectors. */
/* on entry */
/* x complex*16(ldx,p), where ldx.ge.n. */
/* x contains the matrix whose singular value */
/* decomposition is to be computed. x is */
/* destroyed by zsvdc. */
/* ldx integer. */
/* ldx is the leading dimension of the array x. */
/* n integer. */
/* n is the number of rows of the matrix x. */
/* p integer. */
/* p is the number of columns of the matrix x. */
/* ldu integer. */
/* ldu is the leading dimension of the array u */
/* (see below). */
/* ldv integer. */
/* ldv is the leading dimension of the array v */
/* (see below). */
/* work complex*16(n). */
/* work is a scratch array. */
/* job integer. */
/* job controls the computation of the singular */
/* vectors. it has the decimal expansion ab */
/* with the following meaning */
/* a.eq.0 do not compute the left singular */
/* vectors. */
/* a.eq.1 return the n left singular vectors */
/* in u. */
/* a.ge.2 returns the first min(n,p) */
/* left singular vectors in u. */
/* b.eq.0 do not compute the right singular */
/* vectors. */
/* b.eq.1 return the right singular vectors */
/* in v. */
/* on return */
/* s complex*16(mm), where mm=min(n+1,p). */
/* the first min(n,p) entries of s contain the */
/* singular values of x arranged in descending */
/* order of magnitude. */
/* e complex*16(p). */
/* e ordinarily contains zeros. however see the */
/* discussion of info for exceptions. */
/* u complex*16(ldu,k), where ldu.ge.n. if joba.eq.1 */
/* then k.eq.n, if joba.ge.2 then */
/* k.eq.min(n,p). */
/* u contains the matrix of left singular vectors. */
/* u is not referenced if joba.eq.0. if n.le.p */
/* or if joba.gt.2, then u may be identified with x */
/* in the subroutine call. */
/* v complex*16(ldv,p), where ldv.ge.p. */
/* v contains the matrix of right singular vectors. */
/* v is not referenced if jobb.eq.0. if p.le.n, */
/* then v may be identified whth x in the */
/* subroutine call. */
/* info integer. */
/* the singular values (and their corresponding */
/* singular vectors) s(info+1),s(info+2),...,s(m) */
/* are correct (here m=min(n,p)). thus if */
/* info.eq.0, all the singular values and their */
/* vectors are correct. in any event, the matrix */
/* b = ctrans(u)*x*v is the bidiagonal matrix */
/* with the elements of s on its diagonal and the */
/* elements of e on its super-diagonal (ctrans(u) */
/* is the conjugate-transpose of u). thus the */
/* singular values of x and b are the same. */
/* linpack. this version dated 03/19/79 . */
/* correction to shift calculation made 2/85. */
/* g.w. stewart, university of maryland, argonne national lab. */
/* zsvdc uses the following functions and subprograms. */
/* external zdrot */
/* blas zaxpy,zdotc,zscal,zswap,dznrm2,drotg */
/* fortran dabs,dmax1,cdabs,dcmplx */
/* fortran dconjg,max0,min0,mod,dsqrt */
/* internal variables */
/* set the maximum number of iterations. */
/* Parameter adjustments */
x_dim1 = ldx;
x_offset = x_dim1 + 1;
x -= x_offset;
--s;
--e;
u_dim1 = ldu;
u_offset = u_dim1 + 1;
u -= u_offset;
v_dim1 = ldv;
v_offset = v_dim1 + 1;
v -= v_offset;
--work;
/* keep compiler happy */
l = ls = 0;
/* Function Body */
maxit = 30;
/* determine what is to be computed. */
wantu = FALSE;
wantv = FALSE;
jobu = job % 100 / 10;
ncu = n;
if (jobu > 1) {
ncu = min(n,p);
}
if (jobu != 0) {
wantu = TRUE;
}
if (job % 10 != 0) {
wantv = TRUE;
}
/* reduce x to bidiagonal form, storing the diagonal elements */
/* in s and the super-diagonal elements in e. */
*info = 0;
/* Computing MIN */
i__1 = n - 1;
nct = min(i__1,p);
/* Computing MAX */
/* Computing MIN */
i__3 = p - 2;
i__1 = 0, i__2 = min(i__3,n);
nrt = max(i__1,i__2);
lu = max(nct,nrt);
if (lu < 1) {
goto L170;
}
i__1 = lu;
for (l = 1; l <= i__1; ++l) {
lp1 = l + 1;
if (l > nct) {
goto L20;
}
/* compute the transformation for the l-th column and */
/* place the l-th diagonal in s(l). */
i__2 = l;
i__3 = n - l + 1;
d__1 = blas_dznrm2(i__3, &x[l + l * x_dim1], 1);
z__1.r = d__1, z__1.i = 0.;
s[i__2].r = z__1.r, s[i__2].i = z__1.i;
i__2 = l;
i__3 = l;
z__1.r = s[i__3].r * 0. - s[i__3].i * -1., z__1.i = s[i__3].r * -1. +
s[i__3].i * 0.;
if ((d__1 = s[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) == 0.)
{
goto L10;
}
i__2 = l + l * x_dim1;
i__3 = l + l * x_dim1;
z__1.r = x[i__3].r * 0. - x[i__3].i * -1., z__1.i = x[i__3].r * -1. +
x[i__3].i * 0.;
if ((d__1 = x[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) != 0.)
{
i__4 = l;
d__3 = z_abs(&s[l]);
i__5 = l + l * x_dim1;
d__4 = z_abs(&x[l + l * x_dim1]);
z__3.r = x[i__5].r / d__4, z__3.i = x[i__5].i / d__4;
z__2.r = d__3 * z__3.r, z__2.i = d__3 * z__3.i;
s[i__4].r = z__2.r, s[i__4].i = z__2.i;
}
i__2 = n - l + 1;
z_div(&z__1, &c_b10, &s[l]);
blas_zscal(i__2, &z__1, &x[l + l * x_dim1], 1);
i__2 = l + l * x_dim1;
i__3 = l + l * x_dim1;
z__1.r = x[i__3].r + 1., z__1.i = x[i__3].i + 0.;
x[i__2].r = z__1.r, x[i__2].i = z__1.i;
L10:
i__2 = l;
i__3 = l;
z__1.r = -s[i__3].r, z__1.i = -s[i__3].i;
s[i__2].r = z__1.r, s[i__2].i = z__1.i;
L20:
if (p < lp1) {
goto L50;
}
i__2 = p;
for (j = lp1; j <= i__2; ++j) {
if (l > nct) {
goto L30;
}
i__3 = l;
i__4 = l;
z__1.r = s[i__4].r * 0. - s[i__4].i * -1., z__1.i = s[i__4].r *
-1. + s[i__4].i * 0.;
if ((d__1 = s[i__3].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) ==
0.) {
goto L30;
}
/* apply the transformation. */
i__3 = n - l + 1;
blas_zdotc(&z__3, i__3, &x[l + l * x_dim1], 1, &x[l + j * x_dim1]
, 1);
z__2.r = -z__3.r, z__2.i = -z__3.i;
z_div(&z__1, &z__2, &x[l + l * x_dim1]);
t.r = z__1.r, t.i = z__1.i;
i__3 = n - l + 1;
blas_zaxpy(i__3, &t, &x[l + l * x_dim1], 1, &x[l + j * x_dim1], 1);
L30:
/* place the l-th row of x into e for the */
/* subsequent calculation of the row transformation. */
i__3 = j;
d_cnjg(&z__1, &x[l + j * x_dim1]);
e[i__3].r = z__1.r, e[i__3].i = z__1.i;
/* L40: */
}
L50:
if (! wantu || l > nct) {
goto L70;
}
/* place the transformation in u for subsequent back */
/* multiplication. */
i__2 = n;
for (i = l; i <= i__2; ++i) {
i__3 = i + l * u_dim1;
i__4 = i + l * x_dim1;
u[i__3].r = x[i__4].r, u[i__3].i = x[i__4].i;
/* L60: */
}
L70:
if (l > nrt) {
goto L150;
}
/* compute the l-th row transformation and place the */
/* l-th super-diagonal in e(l). */
i__2 = l;
i__3 = p - l;
d__1 = blas_dznrm2(i__3, &e[lp1], 1);
z__1.r = d__1, z__1.i = 0.;
e[i__2].r = z__1.r, e[i__2].i = z__1.i;
i__2 = l;
i__3 = l;
z__1.r = e[i__3].r * 0. - e[i__3].i * -1., z__1.i = e[i__3].r * -1. +
e[i__3].i * 0.;
if ((d__1 = e[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) == 0.)
{
goto L80;
}
i__2 = lp1;
i__3 = lp1;
z__1.r = e[i__3].r * 0. - e[i__3].i * -1., z__1.i = e[i__3].r * -1. +
e[i__3].i * 0.;
if ((d__1 = e[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) != 0.)
{
i__4 = l;
d__3 = z_abs(&e[l]);
i__5 = lp1;
d__4 = z_abs(&e[lp1]);
z__3.r = e[i__5].r / d__4, z__3.i = e[i__5].i / d__4;
z__2.r = d__3 * z__3.r, z__2.i = d__3 * z__3.i;
e[i__4].r = z__2.r, e[i__4].i = z__2.i;
}
i__2 = p - l;
z_div(&z__1, &c_b10, &e[l]);
blas_zscal(i__2, &z__1, &e[lp1], 1);
i__2 = lp1;
i__3 = lp1;
z__1.r = e[i__3].r + 1., z__1.i = e[i__3].i + 0.;
e[i__2].r = z__1.r, e[i__2].i = z__1.i;
L80:
i__2 = l;
d_cnjg(&z__2, &e[l]);
z__1.r = -z__2.r, z__1.i = -z__2.i;
e[i__2].r = z__1.r, e[i__2].i = z__1.i;
i__2 = l;
i__3 = l;
z__1.r = e[i__3].r * 0. - e[i__3].i * -1., z__1.i = e[i__3].r * -1. +
e[i__3].i * 0.;
if (lp1 > n || (d__1 = e[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(
d__2)) == 0.) {
goto L120;
}
/* apply the transformation. */
i__2 = n;
for (i = lp1; i <= i__2; ++i) {
i__3 = i;
work[i__3].r = 0., work[i__3].i = 0.;
/* L90: */
}
i__2 = p;
for (j = lp1; j <= i__2; ++j) {
i__3 = n - l;
blas_zaxpy(i__3, &e[j], &x[lp1 + j * x_dim1], 1, &work[lp1], 1);
/* L100: */
}
i__2 = p;
for (j = lp1; j <= i__2; ++j) {
i__3 = n - l;
i__4 = j;
z__3.r = -e[i__4].r, z__3.i = -e[i__4].i;
z_div(&z__2, &z__3, &e[lp1]);
d_cnjg(&z__1, &z__2);
blas_zaxpy(i__3, &z__1, &work[lp1], 1, &x[lp1 + j * x_dim1], 1);
/* L110: */
}
L120:
if (! wantv) {
goto L140;
}
/* place the transformation in v for subsequent */
/* back multiplication. */
i__2 = p;
for (i = lp1; i <= i__2; ++i) {
i__3 = i + l * v_dim1;
i__4 = i;
v[i__3].r = e[i__4].r, v[i__3].i = e[i__4].i;
/* L130: */
}
L140:
L150:
/* L160: */
;
}
L170:
/* set up the final bidiagonal matrix or order m. */
/* Computing MIN */
i__1 = p, i__2 = n + 1;
m = min(i__1,i__2);
nctp1 = nct + 1;
nrtp1 = nrt + 1;
if (nct < p) {
i__1 = nctp1;
i__2 = nctp1 + nctp1 * x_dim1;
s[i__1].r = x[i__2].r, s[i__1].i = x[i__2].i;
}
if (n < m) {
i__1 = m;
s[i__1].r = 0., s[i__1].i = 0.;
}
if (nrtp1 < m) {
i__1 = nrtp1;
i__2 = nrtp1 + m * x_dim1;
e[i__1].r = x[i__2].r, e[i__1].i = x[i__2].i;
}
i__1 = m;
e[i__1].r = 0., e[i__1].i = 0.;
/* if required, generate u. */
if (! wantu) {
goto L300;
}
if (ncu < nctp1) {
goto L200;
}
i__1 = ncu;
for (j = nctp1; j <= i__1; ++j) {
i__2 = n;
for (i = 1; i <= i__2; ++i) {
i__3 = i + j * u_dim1;
u[i__3].r = 0., u[i__3].i = 0.;
/* L180: */
}
i__2 = j + j * u_dim1;
u[i__2].r = 1., u[i__2].i = 0.;
/* L190: */
}
L200:
if (nct < 1) {
goto L290;
}
i__1 = nct;
for (ll = 1; ll <= i__1; ++ll) {
l = nct - ll + 1;
i__2 = l;
i__3 = l;
z__1.r = s[i__3].r * 0. - s[i__3].i * -1., z__1.i = s[i__3].r * -1. +
s[i__3].i * 0.;
if ((d__1 = s[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) == 0.)
{
goto L250;
}
lp1 = l + 1;
if (ncu < lp1) {
goto L220;
}
i__2 = ncu;
for (j = lp1; j <= i__2; ++j) {
i__3 = n - l + 1;
blas_zdotc(&z__3, i__3, &u[l + l * u_dim1], 1, &u[l + j * u_dim1]
, 1);
z__2.r = -z__3.r, z__2.i = -z__3.i;
z_div(&z__1, &z__2, &u[l + l * u_dim1]);
t.r = z__1.r, t.i = z__1.i;
i__3 = n - l + 1;
blas_zaxpy(i__3, &t, &u[l + l * u_dim1], 1, &u[l + j * u_dim1], 1);
/* L210: */
}
L220:
i__2 = n - l + 1;
blas_zscal(i__2, &c_b60, &u[l + l * u_dim1], 1);
i__2 = l + l * u_dim1;
i__3 = l + l * u_dim1;
z__1.r = u[i__3].r + 1., z__1.i = u[i__3].i + 0.;
u[i__2].r = z__1.r, u[i__2].i = z__1.i;
lm1 = l - 1;
if (lm1 < 1) {
goto L240;
}
i__2 = lm1;
for (i = 1; i <= i__2; ++i) {
i__3 = i + l * u_dim1;
u[i__3].r = 0., u[i__3].i = 0.;
/* L230: */
}
L240:
goto L270;
L250:
i__2 = n;
for (i = 1; i <= i__2; ++i) {
i__3 = i + l * u_dim1;
u[i__3].r = 0., u[i__3].i = 0.;
/* L260: */
}
i__2 = l + l * u_dim1;
u[i__2].r = 1., u[i__2].i = 0.;
L270:
/* L280: */
;
}
L290:
L300:
/* if it is required, generate v. */
if (! wantv) {
goto L350;
}
i__1 = p;
for (ll = 1; ll <= i__1; ++ll) {
l = p - ll + 1;
lp1 = l + 1;
if (l > nrt) {
goto L320;
}
i__2 = l;
i__3 = l;
z__1.r = e[i__3].r * 0. - e[i__3].i * -1., z__1.i = e[i__3].r * -1. +
e[i__3].i * 0.;
if ((d__1 = e[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) == 0.)
{
goto L320;
}
i__2 = p;
for (j = lp1; j <= i__2; ++j) {
i__3 = p - l;
blas_zdotc(&z__3, i__3, &v[lp1 + l * v_dim1], 1, &v[lp1 + j *
v_dim1], 1);
z__2.r = -z__3.r, z__2.i = -z__3.i;
z_div(&z__1, &z__2, &v[lp1 + l * v_dim1]);
t.r = z__1.r, t.i = z__1.i;
i__3 = p - l;
blas_zaxpy(i__3, &t, &v[lp1 + l * v_dim1], 1, &v[lp1 + j *
v_dim1], 1);
/* L310: */
}
L320:
i__2 = p;
for (i = 1; i <= i__2; ++i) {
i__3 = i + l * v_dim1;
v[i__3].r = 0., v[i__3].i = 0.;
/* L330: */
}
i__2 = l + l * v_dim1;
v[i__2].r = 1., v[i__2].i = 0.;
/* L340: */
}
L350:
/* transform s and e so that they are double precision. */
i__1 = m;
for (i = 1; i <= i__1; ++i) {
i__2 = i;
i__3 = i;
z__1.r = s[i__3].r * 0. - s[i__3].i * -1., z__1.i = s[i__3].r * -1. +
s[i__3].i * 0.;
if ((d__1 = s[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) == 0.)
{
goto L360;
}
d__1 = (double) z_abs(&s[i]);
z__1.r = d__1, z__1.i = 0.;
t.r = z__1.r, t.i = z__1.i;
z_div(&z__1, &s[i], &t);
r.r = z__1.r, r.i = z__1.i;
i__2 = i;
s[i__2].r = t.r, s[i__2].i = t.i;
if (i < m) {
i__2 = i;
z_div(&z__1, &e[i], &r);
e[i__2].r = z__1.r, e[i__2].i = z__1.i;
}
if (wantu) {
blas_zscal(n, &r, &u[i * u_dim1 + 1], 1);
}
L360:
/* ...exit */
if (i == m) {
goto L390;
}
i__2 = i;
i__3 = i;
z__1.r = e[i__3].r * 0. - e[i__3].i * -1., z__1.i = e[i__3].r * -1. +
e[i__3].i * 0.;
if ((d__1 = e[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) == 0.)
{
goto L370;
}
d__1 = (double) z_abs(&e[i]);
z__1.r = d__1, z__1.i = 0.;
t.r = z__1.r, t.i = z__1.i;
z_div(&z__1, &t, &e[i]);
r.r = z__1.r, r.i = z__1.i;
i__2 = i;
e[i__2].r = t.r, e[i__2].i = t.i;
i__2 = i + 1;
i__3 = i + 1;
z__1.r = s[i__3].r * r.r - s[i__3].i * r.i, z__1.i = s[i__3].r * r.i
+ s[i__3].i * r.r;
s[i__2].r = z__1.r, s[i__2].i = z__1.i;
if (wantv) {
blas_zscal(p, &r, &v[(i + 1) * v_dim1 + 1], 1);
}
L370:
/* L380: */
;
}
L390:
/* main iteration loop for the singular values. */
mm = m;
iter = 0;
L400:
/* quit if all the singular values have been found. */
/* ...exit */
if (m == 0) {
goto L660;
}
/* if too many iterations have been performed, set */
/* flag and return. */
if (iter < maxit) {
goto L410;
}
*info = m;
/* ......exit */
goto L660;
L410:
/* this section of the program inspects for */
/* negligible elements in the s and e arrays. on */
/* completion the variables kase and l are set as follows. */
/* kase = 1 if s(m) and e(l-1) are negligible and l.lt.m */
/* kase = 2 if s(l) is negligible and l.lt.m */
/* kase = 3 if e(l-1) is negligible, l.lt.m, and */
/* s(l), ..., s(m) are not negligible (qr step). */
/* kase = 4 if e(m-1) is negligible (convergence). */
i__1 = m;
for (ll = 1; ll <= i__1; ++ll) {
l = m - ll;
/* ...exit */
if (l == 0) {
goto L440;
}
test = z_abs(&s[l]) + z_abs(&s[l + 1]);
ztest = test + z_abs(&e[l]);
if (ztest != test) {
goto L420;
}
i__2 = l;
e[i__2].r = 0., e[i__2].i = 0.;
/* ......exit */
goto L440;
L420:
/* L430: */
;
}
L440:
if (l != m - 1) {
goto L450;
}
kase = 4;
goto L520;
L450:
lp1 = l + 1;
mp1 = m + 1;
i__1 = mp1;
for (lls = lp1; lls <= i__1; ++lls) {
ls = m - lls + lp1;
/* ...exit */
if (ls == l) {
goto L480;
}
test = 0.;
if (ls != m) {
test += z_abs(&e[ls]);
}
if (ls != l + 1) {
test += z_abs(&e[ls - 1]);
}
ztest = test + z_abs(&s[ls]);
if (ztest != test) {
goto L460;
}
i__2 = ls;
s[i__2].r = 0., s[i__2].i = 0.;
/* ......exit */
goto L480;
L460:
/* L470: */
;
}
L480:
if (ls != l) {
goto L490;
}
kase = 3;
goto L510;
L490:
if (ls != m) {
goto L500;
}
kase = 1;
goto L510;
L500:
kase = 2;
l = ls;
L510:
L520:
++l;
/* perform the task indicated by kase. */
switch ((int)kase) {
case 1: goto L530;
case 2: goto L560;
case 3: goto L580;
case 4: goto L610;
}
/* deflate negligible s(m). */
L530:
mm1 = m - 1;
i__1 = m - 1;
f = e[i__1].r;
i__1 = m - 1;
e[i__1].r = 0., e[i__1].i = 0.;
i__1 = mm1;
for (kk = l; kk <= i__1; ++kk) {
k = mm1 - kk + l;
i__2 = k;
t1 = s[i__2].r;
blas_drotg(&t1, &f, &cs, &sn);
i__2 = k;
z__1.r = t1, z__1.i = 0.;
s[i__2].r = z__1.r, s[i__2].i = z__1.i;
if (k == l) {
goto L540;
}
i__2 = k - 1;
f = -sn * e[i__2].r;
i__2 = k - 1;
i__3 = k - 1;
z__1.r = cs * e[i__3].r, z__1.i = cs * e[i__3].i;
e[i__2].r = z__1.r, e[i__2].i = z__1.i;
L540:
if (wantv) {
blas_zdrot(p, &v[k * v_dim1 + 1], 1, &v[m * v_dim1 + 1], 1, cs, sn);
}
/* L550: */
}
goto L650;
/* split at negligible s(l). */
L560:
i__1 = l - 1;
f = e[i__1].r;
i__1 = l - 1;
e[i__1].r = 0., e[i__1].i = 0.;
i__1 = m;
for (k = l; k <= i__1; ++k) {
i__2 = k;
t1 = s[i__2].r;
blas_drotg(&t1, &f, &cs, &sn);
i__2 = k;
z__1.r = t1, z__1.i = 0.;
s[i__2].r = z__1.r, s[i__2].i = z__1.i;
i__2 = k;
f = -sn * e[i__2].r;
i__2 = k;
i__3 = k;
z__1.r = cs * e[i__3].r, z__1.i = cs * e[i__3].i;
e[i__2].r = z__1.r, e[i__2].i = z__1.i;
if (wantu) {
blas_zdrot(n, &u[k * u_dim1 + 1], 1, &u[(l - 1) * u_dim1 + 1], 1, cs, sn);
}
/* L570: */
}
goto L650;
/* perform one qr step. */
L580:
/* calculate the shift. */
/* Computing MAX */
r__1 = z_abs(&s[m]), r__2 = z_abs(&s[m - 1]), r__1 = max(r__1,r__2),
r__2 = z_abs(&e[m - 1]), r__1 = max(r__1,r__2), r__2 = z_abs(&s[
l]), r__1 = max(r__1,r__2), r__2 = z_abs(&e[l]);
scale = max(r__1,r__2);
i__1 = m;
sm = s[i__1].r / scale;
i__1 = m - 1;
smm1 = s[i__1].r / scale;
i__1 = m - 1;
emm1 = e[i__1].r / scale;
i__1 = l;
sl = s[i__1].r / scale;
i__1 = l;
el = e[i__1].r / scale;
/* Computing 2nd power */
d__1 = emm1;
b = ((smm1 + sm) * (smm1 - sm) + d__1 * d__1) / 2.;
/* Computing 2nd power */
d__1 = sm * emm1;
c = d__1 * d__1;
shift = 0.;
if (b == 0. && c == 0.) {
goto L590;
}
/* Computing 2nd power */
d__1 = b;
shift = sqrt(d__1 * d__1 + c);
if (b < 0.) {
shift = -shift;
}
shift = c / (b + shift);
L590:
f = (sl + sm) * (sl - sm) + shift;
g = sl * el;
/* chase zeros. */
mm1 = m - 1;
i__1 = mm1;
for (k = l; k <= i__1; ++k) {
blas_drotg(&f, &g, &cs, &sn);
if (k != l) {
i__2 = k - 1;
z__1.r = f, z__1.i = 0.;
e[i__2].r = z__1.r, e[i__2].i = z__1.i;
}
i__2 = k;
i__3 = k;
f = cs * s[i__2].r + sn * e[i__3].r;
i__2 = k;
i__3 = k;
z__2.r = cs * e[i__3].r, z__2.i = cs * e[i__3].i;
i__4 = k;
z__3.r = sn * s[i__4].r, z__3.i = sn * s[i__4].i;
z__1.r = z__2.r - z__3.r, z__1.i = z__2.i - z__3.i;
e[i__2].r = z__1.r, e[i__2].i = z__1.i;
i__2 = k + 1;
g = sn * s[i__2].r;
i__2 = k + 1;
i__3 = k + 1;
z__1.r = cs * s[i__3].r, z__1.i = cs * s[i__3].i;
s[i__2].r = z__1.r, s[i__2].i = z__1.i;
if (wantv) {
blas_zdrot(p, &v[k * v_dim1 + 1], 1, &v[(k + 1) * v_dim1 + 1], 1, cs, sn);
}
blas_drotg(&f, &g, &cs, &sn);
i__2 = k;
z__1.r = f, z__1.i = 0.;
s[i__2].r = z__1.r, s[i__2].i = z__1.i;
i__2 = k;
i__3 = k + 1;
f = cs * e[i__2].r + sn * s[i__3].r;
i__2 = k + 1;
d__1 = -sn;
i__3 = k;
z__2.r = d__1 * e[i__3].r, z__2.i = d__1 * e[i__3].i;
i__4 = k + 1;
z__3.r = cs * s[i__4].r, z__3.i = cs * s[i__4].i;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
s[i__2].r = z__1.r, s[i__2].i = z__1.i;
i__2 = k + 1;
g = sn * e[i__2].r;
i__2 = k + 1;
i__3 = k + 1;
z__1.r = cs * e[i__3].r, z__1.i = cs * e[i__3].i;
e[i__2].r = z__1.r, e[i__2].i = z__1.i;
if (wantu && k < n) {
blas_zdrot(n, &u[k * u_dim1 + 1], 1, &u[(k + 1) * u_dim1 + 1], 1, cs, sn);
}
/* L600: */
}
i__1 = m - 1;
z__1.r = f, z__1.i = 0.;
e[i__1].r = z__1.r, e[i__1].i = z__1.i;
++iter;
goto L650;
/* convergence. */
L610:
/* make the singular value positive */
i__1 = l;
if (s[i__1].r >= 0.) {
goto L620;
}
i__1 = l;
i__2 = l;
z__1.r = -s[i__2].r, z__1.i = -s[i__2].i;
s[i__1].r = z__1.r, s[i__1].i = z__1.i;
if (wantv) {
blas_zscal(p, &c_b60, &v[l * v_dim1 + 1], 1);
}
L620:
/* order the singular value. */
L630:
if (l == mm) {
goto L640;
}
/* ...exit */
i__1 = l;
i__2 = l + 1;
if (s[i__1].r >= s[i__2].r) {
goto L640;
}
i__1 = l;
t.r = s[i__1].r, t.i = s[i__1].i;
i__1 = l;
i__2 = l + 1;
s[i__1].r = s[i__2].r, s[i__1].i = s[i__2].i;
i__1 = l + 1;
s[i__1].r = t.r, s[i__1].i = t.i;
if (wantv && l < p) {
blas_zswap(p, &v[l * v_dim1 + 1], 1, &v[(l + 1) * v_dim1 + 1], 1);
}
if (wantu && l < n) {
blas_zswap(n, &u[l * u_dim1 + 1], 1, &u[(l + 1) * u_dim1 + 1], 1);
}
++l;
goto L630;
L640:
iter = 0;
--m;
L650:
goto L400;
L660:
return;
}
syntax highlighted by Code2HTML, v. 0.9.1