zsvdc.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/* linpack/zsvdc.f -- translated by f2c (version 20050501).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
        or, if you install libf2c.a in a standard place, with -lf2c -lm
        -- in that order, at the end of the command line, as in
                cc *.o -lf2c -lm
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

                http://www.netlib.org/f2c/libf2c.zip
*/

#ifdef __cplusplus
extern "C" {
#endif
#include "v3p_netlib.h"

/* Table of constant values */

static integer c__1 = 1;
static doublecomplex c_b10 = {1.,0.};
static doublecomplex c_b60 = {-1.,0.};

/*<       subroutine zsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info) >*/
/* Subroutine */ int zsvdc_(doublecomplex *x, integer *ldx, integer *n, 
        integer *p, doublecomplex *s, doublecomplex *e, doublecomplex *u, 
        integer *ldu, doublecomplex *v, integer *ldv, doublecomplex *work, 
        integer *job, integer *info)
{
    /* System generated locals */
    integer x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2, 
            i__3, i__4, i__5;
    doublereal d__1, d__2, d__3, d__4;
    doublecomplex z__1, z__2, z__3;

    /* Builtin functions */
    double z_abs(doublecomplex *);
    void z_div(doublecomplex *, doublecomplex *, doublecomplex *), d_cnjg(
            doublecomplex *, doublecomplex *);
    double sqrt(doublereal);

    /* Local variables */
    doublereal b, c__, f, g;
    integer i__, j, k, l=0, m;
    doublecomplex r__, t;
    doublereal t1, el;
    integer kk;
    doublereal cs;
    integer ll, mm, ls=0;
    doublereal sl;
    integer lu;
    doublereal sm, sn;
    integer lm1, mm1, lp1, mp1, nct, ncu, lls, nrt;
    doublereal emm1, smm1;
    integer kase, jobu, iter;
    doublereal test;
    integer nctp1, nrtp1;
    doublereal scale;
    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 
            doublecomplex *, integer *);
    doublereal shift;
    extern /* Subroutine */ int drotg_(doublereal *, doublereal *, doublereal 
            *, doublereal *);
    integer maxit;
    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 
            doublecomplex *, integer *, doublecomplex *, integer *);
    logical wantu, wantv;
    extern /* Subroutine */ int zdrot_(integer *, doublecomplex *, integer *, 
            doublecomplex *, integer *, doublereal *, doublereal *), zswap_(
            integer *, doublecomplex *, integer *, doublecomplex *, integer *)
            ;
    doublereal ztest;
    extern /* Subroutine */ int zaxpy_(integer *, doublecomplex *, 
            doublecomplex *, integer *, doublecomplex *, integer *);
    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);

/*<       integer ldx,n,p,ldu,ldv,job,info >*/
/*<       complex*16 x(ldx,1),s(1),e(1),u(ldu,1),v(ldv,1),work(1) >*/


/*     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 */

/*<    >*/
/*<       complex*16 zdotc,t,r >*/
/*<    >*/
/*<       logical wantu,wantv >*/

/*<       complex*16 csign,zdum,zdum1,zdum2 >*/
/*<       double precision cabs1 >*/
/*<       double precision dreal,dimag >*/
/*<       complex*16 zdumr,zdumi >*/
/*<       dreal(zdumr) = zdumr >*/
/*<       dimag(zdumi) = (0.0d0,-1.0d0)*zdumi >*/
/*<       cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum)) >*/
/*<       csign(zdum1,zdum2) = cdabs(zdum1)*(zdum2/cdabs(zdum2)) >*/

/*     set the maximum number of iterations. */

/*<       maxit = 1000 >*/
    /* Parameter adjustments */
    x_dim1 = *ldx;
    x_offset = 1 + x_dim1;
    x -= x_offset;
    --s;
    --e;
    u_dim1 = *ldu;
    u_offset = 1 + u_dim1;
    u -= u_offset;
    v_dim1 = *ldv;
    v_offset = 1 + v_dim1;
    v -= v_offset;
    --work;

    /* Function Body */
    maxit = 1000;

/*     determine what is to be computed. */

/*<       wantu = .false. >*/
    wantu = FALSE_;
/*<       wantv = .false. >*/
    wantv = FALSE_;
/*<       jobu = mod(job,100)/10 >*/
    jobu = *job % 100 / 10;
/*<       ncu = n >*/
    ncu = *n;
/*<       if (jobu .gt. 1) ncu = min0(n,p) >*/
    if (jobu > 1) {
        ncu = min(*n,*p);
    }
/*<       if (jobu .ne. 0) wantu = .true. >*/
    if (jobu != 0) {
        wantu = TRUE_;
    }
/*<       if (mod(job,10) .ne. 0) wantv = .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 >*/
    *info = 0;
/*<       nct = min0(n-1,p) >*/
/* Computing MIN */
    i__1 = *n - 1;
    nct = min(i__1,*p);
/*<       nrt = max0(0,min0(p-2,n)) >*/
/* Computing MAX */
/* Computing MIN */
    i__3 = *p - 2;
    i__1 = 0, i__2 = min(i__3,*n);
    nrt = max(i__1,i__2);
/*<       lu = max0(nct,nrt) >*/
    lu = max(nct,nrt);
/*<       if (lu .lt. 1) go to 170 >*/
    if (lu < 1) {
        goto L170;
    }
/*<       do 160 l = 1, lu >*/
    i__1 = lu;
    for (l = 1; l <= i__1; ++l) {
/*<          lp1 = l + 1 >*/
        lp1 = l + 1;
/*<          if (l .gt. nct) go to 20 >*/
        if (l > nct) {
            goto L20;
        }

/*           compute the transformation for the l-th column and */
/*           place the l-th diagonal in s(l). */

/*<             s(l) = dcmplx(dznrm2(n-l+1,x(l,l),1),0.0d0) >*/
        i__2 = l;
        i__3 = *n - l + 1;
        d__1 = dznrm2_(&i__3, &x[l + l * x_dim1], &c__1);
        z__1.r = d__1, z__1.i = 0.;
        s[i__2].r = z__1.r, s[i__2].i = z__1.i;
/*<             if (cabs1(s(l)) .eq. 0.0d0) go to 10 >*/
        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].i * 0. + 
                s[i__3].r * -1.;
        if ((d__1 = s[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) == 0.) 
                {
            goto L10;
        }
/*<                if (cabs1(x(l,l)) .ne. 0.0d0) s(l) = csign(s(l),x(l,l)) >*/
        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].i * 0. + 
                x[i__3].r * -1.;
        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;
        }
/*<                call zscal(n-l+1,1.0d0/s(l),x(l,l),1) >*/
        i__2 = *n - l + 1;
        z_div(&z__1, &c_b10, &s[l]);
        zscal_(&i__2, &z__1, &x[l + l * x_dim1], &c__1);
/*<                x(l,l) = (1.0d0,0.0d0) + x(l,l) >*/
        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;
/*<    10       continue >*/
L10:
/*<             s(l) = -s(l) >*/
        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;
/*<    20    continue >*/
L20:
/*<          if (p .lt. lp1) go to 50 >*/
        if (*p < lp1) {
            goto L50;
        }
/*<          do 40 j = lp1, p >*/
        i__2 = *p;
        for (j = lp1; j <= i__2; ++j) {
/*<             if (l .gt. nct) go to 30 >*/
            if (l > nct) {
                goto L30;
            }
/*<             if (cabs1(s(l)) .eq. 0.0d0) go to 30 >*/
            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].i * 
                    0. + s[i__4].r * -1.;
            if ((d__1 = s[i__3].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) == 
                    0.) {
                goto L30;
            }

/*              apply the transformation. */

/*<                t = -zdotc(n-l+1,x(l,l),1,x(l,j),1)/x(l,l) >*/
            i__3 = *n - l + 1;
            zdotc_(&z__3, &i__3, &x[l + l * x_dim1], &c__1, &x[l + j * x_dim1]
                    , &c__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;
/*<                call zaxpy(n-l+1,t,x(l,l),1,x(l,j),1) >*/
            i__3 = *n - l + 1;
            zaxpy_(&i__3, &t, &x[l + l * x_dim1], &c__1, &x[l + j * x_dim1], &
                    c__1);
/*<    30       continue >*/
L30:

/*           place the l-th row of x into  e for the */
/*           subsequent calculation of the row transformation. */

/*<             e(j) = dconjg(x(l,j)) >*/
            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;
/*<    40    continue >*/
/* L40: */
        }
/*<    50    continue >*/
L50:
/*<          if (.not.wantu .or. l .gt. nct) go to 70 >*/
        if (! wantu || l > nct) {
            goto L70;
        }

/*           place the transformation in u for subsequent back */
/*           multiplication. */

/*<             do 60 i = l, n >*/
        i__2 = *n;
        for (i__ = l; i__ <= i__2; ++i__) {
/*<                u(i,l) = x(i,l) >*/
            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;
/*<    60       continue >*/
/* L60: */
        }
/*<    70    continue >*/
L70:
/*<          if (l .gt. nrt) go to 150 >*/
        if (l > nrt) {
            goto L150;
        }

/*           compute the l-th row transformation and place the */
/*           l-th super-diagonal in e(l). */

/*<             e(l) = dcmplx(dznrm2(p-l,e(lp1),1),0.0d0) >*/
        i__2 = l;
        i__3 = *p - l;
        d__1 = dznrm2_(&i__3, &e[lp1], &c__1);
        z__1.r = d__1, z__1.i = 0.;
        e[i__2].r = z__1.r, e[i__2].i = z__1.i;
/*<             if (cabs1(e(l)) .eq. 0.0d0) go to 80 >*/
        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].i * 0. + 
                e[i__3].r * -1.;
        if ((d__1 = e[i__2].r, abs(d__1)) + (d__2 = z__1.r, abs(d__2)) == 0.) 
                {
            goto L80;
        }
/*<                if (cabs1(e(lp1)) .ne. 0.0d0) e(l) = csign(e(l),e(lp1)) >*/
        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].i * 0. + 
                e[i__3].r * -1.;
        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;
        }
/*<                call zscal(p-l,1.0d0/e(l),e(lp1),1) >*/
        i__2 = *p - l;
        z_div(&z__1, &c_b10, &e[l]);
        zscal_(&i__2, &z__1, &e[lp1], &c__1);
/*<                e(lp1) = (1.0d0,0.0d0) + e(lp1) >*/
        i__2 = lp1;
        i__3 = lp1;
        z__1.r = e[i__3].r + 1., z__1.i = e[i__3].i + 0.;