dsvdc.c

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

/* linpack/dsvdc.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
*/

/* This code expects correct IEEE rounding behaviour which is not
   always provided.  The source should be built with -ffloat-store.
   A note from the GCC man page:

   -ffloat-store
    Do  not  store floating point variables in registers.  This pre-
    vents undesirable excess precision on machines such as the 68000
    where  the floating registers (of the 68881) keep more precision
    than a double is supposed to have.

    For most programs, the excess precision does only  good,  but  a
    few  programs  rely  on  the precise definition of IEEE floating
    point.  Use `-ffloat-store' for such programs.  */

#ifdef __INTEL_COMPILER
# pragma optimize("", off)
#endif

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

/* Table of constant values */

static integer c__1 = 1;
static doublereal c_b44 = -1.;

/*<       subroutine dsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info) >*/
/* Subroutine */ int 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)
{
    /* System generated locals */
    integer x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2, 
            i__3;
    doublereal d__1, d__2, d__3, d__4, d__5, d__6, d__7;

    /* Builtin functions */
    double d_sign(doublereal *, doublereal *), sqrt(doublereal);

    /* Local variables */
    doublereal b, c__, f, g;
    integer i__, j, k, l=0, m;
    doublereal t, 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;
    extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
            integer *);
    integer jobu, iter;
    extern /* Subroutine */ int drot_(integer *, doublereal *, integer *, 
            doublereal *, integer *, doublereal *, doublereal *);
    doublereal test;
    extern doublereal dnrm2_(integer *, doublereal *, integer *);
    integer nctp1, nrtp1;
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
            integer *);
    doublereal scale, shift;
    extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *, 
            doublereal *, integer *), drotg_(doublereal *, doublereal *, 
            doublereal *, doublereal *);
    integer maxit;
    extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 
            integer *, doublereal *, integer *);
    logical wantu, wantv;
    doublereal ztest;

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


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

/*     internal variables */

/*<    >*/
/*<       double precision ddot,t,r >*/
/*<    >*/
/*<       logical wantu,wantv >*/


/*     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) = dnrm2(n-l+1,x(l,l),1) >*/
        i__2 = *n - l + 1;
        s[l] = dnrm2_(&i__2, &x[l + l * x_dim1], &c__1);
/*<             if (s(l) .eq. 0.0d0) go to 10 >*/
        if (s[l] == 0.) {
            goto L10;
        }
/*<                if (x(l,l) .ne. 0.0d0) s(l) = dsign(s(l),x(l,l)) >*/
        if (x[l + l * x_dim1] != 0.) {
            s[l] = d_sign(&s[l], &x[l + l * x_dim1]);
        }
/*<                call dscal(n-l+1,1.0d0/s(l),x(l,l),1) >*/
        i__2 = *n - l + 1;
        d__1 = 1. / s[l];
        dscal_(&i__2, &d__1, &x[l + l * x_dim1], &c__1);
/*<                x(l,l) = 1.0d0 + x(l,l) >*/
        x[l + l * x_dim1] += 1.;
/*<    10       continue >*/
L10:
/*<             s(l) = -s(l) >*/
        s[l] = -s[l];
/*<    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 (s(l) .eq. 0.0d0) go to 30 >*/
            if (s[l] == 0.) {
                goto L30;
            }

/*              apply the transformation. */

/*<                t = -ddot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l) >*/
            i__3 = *n - l + 1;
            t = -ddot_(&i__3, &x[l + l * x_dim1], &c__1, &x[l + j * x_dim1], &
                    c__1) / x[l + l * x_dim1];
/*<                call daxpy(n-l+1,t,x(l,l),1,x(l,j),1) >*/
            i__3 = *n - l + 1;
            daxpy_(&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) = x(l,j) >*/
            e[j] = x[l + j * x_dim1];
/*<    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) >*/
            u[i__ + l * u_dim1] = x[i__ + l * x_dim1];
/*<    60       continue >*/
/* L60: */
        }
/*<    70    continue >*/