ssvdc.c

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

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

/*<             e(l) = snrm2(p-l,e(lp1),1) >*/
        i__2 = *p - l;
        e[l] = snrm2_(&i__2, &e[lp1], &c__1);
/*<             if (e(l) .eq. 0.0e0) go to 80 >*/
        if (e[l] == (float)0.) {
            goto L80;
        }
/*<                if (e(lp1) .ne. 0.0e0) e(l) = sign(e(l),e(lp1)) >*/
        if (e[lp1] != (float)0.) {
            e[l] = r_sign(&e[l], &e[lp1]);
        }
/*<                call sscal(p-l,1.0e0/e(l),e(lp1),1) >*/
        i__2 = *p - l;
        r__1 = (float)1. / e[l];
        sscal_(&i__2, &r__1, &e[lp1], &c__1);
/*<                e(lp1) = 1.0e0 + e(lp1) >*/
        e[lp1] += (float)1.;
/*<    80       continue >*/
L80:
/*<             e(l) = -e(l) >*/
        e[l] = -e[l];
/*<             if (lp1 .gt. n .or. e(l) .eq. 0.0e0) go to 120 >*/
        if (lp1 > *n || e[l] == (float)0.) {
            goto L120;
        }

/*              apply the transformation. */

/*<                do 90 i = lp1, n >*/
        i__2 = *n;
        for (i__ = lp1; i__ <= i__2; ++i__) {
/*<                   work(i) = 0.0e0 >*/
            work[i__] = (float)0.;
/*<    90          continue >*/
/* L90: */
        }
/*<                do 100 j = lp1, p >*/
        i__2 = *p;
        for (j = lp1; j <= i__2; ++j) {
/*<                   call saxpy(n-l,e(j),x(lp1,j),1,work(lp1),1) >*/
            i__3 = *n - l;
            saxpy_(&i__3, &e[j], &x[lp1 + j * x_dim1], &c__1, &work[lp1], &
                    c__1);
/*<   100          continue >*/
/* L100: */
        }
/*<                do 110 j = lp1, p >*/
        i__2 = *p;
        for (j = lp1; j <= i__2; ++j) {
/*<                   call saxpy(n-l,-e(j)/e(lp1),work(lp1),1,x(lp1,j),1) >*/
            i__3 = *n - l;
            r__1 = -e[j] / e[lp1];
            saxpy_(&i__3, &r__1, &work[lp1], &c__1, &x[lp1 + j * x_dim1], &
                    c__1);
/*<   110          continue >*/
/* L110: */
        }
/*<   120       continue >*/
L120:
/*<             if (.not.wantv) go to 140 >*/
        if (! wantv) {
            goto L140;
        }

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

/*<                do 130 i = lp1, p >*/
        i__2 = *p;
        for (i__ = lp1; i__ <= i__2; ++i__) {
/*<                   v(i,l) = e(i) >*/
            v[i__ + l * v_dim1] = e[i__];
/*<   130          continue >*/
/* L130: */
        }
/*<   140       continue >*/
L140:
/*<   150    continue >*/
L150:
/*<   160 continue >*/
/* L160: */
        ;
    }
/*<   170 continue >*/
L170:

/*     set up the final bidiagonal matrix or order m. */

/*<       m = min0(p,n+1) >*/
/* Computing MIN */
    i__1 = *p, i__2 = *n + 1;
    m = min(i__1,i__2);
/*<       nctp1 = nct + 1 >*/
    nctp1 = nct + 1;
/*<       nrtp1 = nrt + 1 >*/
    nrtp1 = nrt + 1;
/*<       if (nct .lt. p) s(nctp1) = x(nctp1,nctp1) >*/
    if (nct < *p) {
        s[nctp1] = x[nctp1 + nctp1 * x_dim1];
    }
/*<       if (n .lt. m) s(m) = 0.0e0 >*/
    if (*n < m) {
        s[m] = (float)0.;
    }
/*<       if (nrtp1 .lt. m) e(nrtp1) = x(nrtp1,m) >*/
    if (nrtp1 < m) {
        e[nrtp1] = x[nrtp1 + m * x_dim1];
    }
/*<       e(m) = 0.0e0 >*/
    e[m] = (float)0.;

/*     if required, generate u. */

/*<       if (.not.wantu) go to 300 >*/
    if (! wantu) {
        goto L300;
    }
/*<          if (ncu .lt. nctp1) go to 200 >*/
    if (ncu < nctp1) {
        goto L200;
    }
/*<          do 190 j = nctp1, ncu >*/
    i__1 = ncu;
    for (j = nctp1; j <= i__1; ++j) {
/*<             do 180 i = 1, n >*/
        i__2 = *n;
        for (i__ = 1; i__ <= i__2; ++i__) {
/*<                u(i,j) = 0.0e0 >*/
            u[i__ + j * u_dim1] = (float)0.;
/*<   180       continue >*/
/* L180: */
        }
/*<             u(j,j) = 1.0e0 >*/
        u[j + j * u_dim1] = (float)1.;
/*<   190    continue >*/
/* L190: */
    }
/*<   200    continue >*/
L200:
/*<          if (nct .lt. 1) go to 290 >*/
    if (nct < 1) {
        goto L290;
    }
/*<          do 280 ll = 1, nct >*/
    i__1 = nct;
    for (ll = 1; ll <= i__1; ++ll) {
/*<             l = nct - ll + 1 >*/
        l = nct - ll + 1;
/*<             if (s(l) .eq. 0.0e0) go to 250 >*/
        if (s[l] == (float)0.) {
            goto L250;
        }
/*<                lp1 = l + 1 >*/
        lp1 = l + 1;
/*<                if (ncu .lt. lp1) go to 220 >*/
        if (ncu < lp1) {
            goto L220;
        }
/*<                do 210 j = lp1, ncu >*/
        i__2 = ncu;
        for (j = lp1; j <= i__2; ++j) {
/*<                   t = -sdot(n-l+1,u(l,l),1,u(l,j),1)/u(l,l) >*/
            i__3 = *n - l + 1;
            t = -sdot_(&i__3, &u[l + l * u_dim1], &c__1, &u[l + j * u_dim1], &
                    c__1) / u[l + l * u_dim1];
/*<                   call saxpy(n-l+1,t,u(l,l),1,u(l,j),1) >*/
            i__3 = *n - l + 1;
            saxpy_(&i__3, &t, &u[l + l * u_dim1], &c__1, &u[l + j * u_dim1], &
                    c__1);
/*<   210          continue >*/
/* L210: */
        }
/*<   220          continue >*/
L220:
/*<                call sscal(n-l+1,-1.0e0,u(l,l),1) >*/
        i__2 = *n - l + 1;
        sscal_(&i__2, &c_b44, &u[l + l * u_dim1], &c__1);
/*<                u(l,l) = 1.0e0 + u(l,l) >*/
        u[l + l * u_dim1] += (float)1.;
/*<                lm1 = l - 1 >*/
        lm1 = l - 1;
/*<                if (lm1 .lt. 1) go to 240 >*/
        if (lm1 < 1) {
            goto L240;
        }
/*<                do 230 i = 1, lm1 >*/
        i__2 = lm1;
        for (i__ = 1; i__ <= i__2; ++i__) {
/*<                   u(i,l) = 0.0e0 >*/
            u[i__ + l * u_dim1] = (float)0.;
/*<   230          continue >*/
/* L230: */
        }
/*<   240          continue >*/
L240:
/*<             go to 270 >*/
        goto L270;
/*<   250       continue >*/
L250:
/*<                do 260 i = 1, n >*/
        i__2 = *n;
        for (i__ = 1; i__ <= i__2; ++i__) {
/*<                   u(i,l) = 0.0e0 >*/
            u[i__ + l * u_dim1] = (float)0.;
/*<   260          continue >*/
/* L260: */
        }
/*<                u(l,l) = 1.0e0 >*/
        u[l + l * u_dim1] = (float)1.;
/*<   270       continue >*/
L270:
/*<   280    continue >*/
/* L280: */
        ;
    }
/*<   290    continue >*/
L290:
/*<   300 continue >*/
L300:

/*     if it is required, generate v. */

/*<       if (.not.wantv) go to 350 >*/
    if (! wantv) {
        goto L350;
    }
/*<          do 340 ll = 1, p >*/
    i__1 = *p;
    for (ll = 1; ll <= i__1; ++ll) {
/*<             l = p - ll + 1 >*/
        l = *p - ll + 1;
/*<             lp1 = l + 1 >*/
        lp1 = l + 1;
/*<             if (l .gt. nrt) go to 320 >*/
        if (l > nrt) {
            goto L320;
        }
/*<             if (e(l) .eq. 0.0e0) go to 320 >*/
        if (e[l] == (float)0.) {
            goto L320;
        }
/*<                do 310 j = lp1, p >*/
        i__2 = *p;
        for (j = lp1; j <= i__2; ++j) {
/*<                   t = -sdot(p-l,v(lp1,l),1,v(lp1,j),1)/v(lp1,l) >*/
            i__3 = *p - l;
            t = -sdot_(&i__3, &v[lp1 + l * v_dim1], &c__1, &v[lp1 + j * 
                    v_dim1], &c__1) / v[lp1 + l * v_dim1];
/*<                   call saxpy(p-l,t,v(lp1,l),1,v(lp1,j),1) >*/
            i__3 = *p - l;
            saxpy_(&i__3, &t, &v[lp1 + l * v_dim1], &c__1, &v[lp1 + j * 
                    v_dim1], &c__1);
/*<   310          continue >*/
/* L310: */
        }
/*<   320       continue >*/
L320:
/*<             do 330 i = 1, p >*/
        i__2 = *p;
        for (i__ = 1; i__ <= i__2; ++i__) {
/*<                v(i,l) = 0.0e0 >*/
            v[i__ + l * v_dim1] = (float)0.;
/*<   330       continue >*/
/* L330: */
        }
/*<             v(l,l) = 1.0e0 >*/
        v[l + l * v_dim1] = (float)1.;
/*<   340    continue >*/
/* L340: */
    }
/*<   350 continue >*/
L350:

/*     main iteration loop for the singular values. */

/*<       mm = m >*/
    mm = m;
/*<       iter = 0 >*/
    iter = 0;
/*<   360 continue >*/
L360:

/*        quit if all the singular values have been found. */

/*     ...exit */
/*<          if (m .eq. 0) go to 620 >*/
    if (m == 0) {
        goto L620;
    }

/*        if too many iterations have been performed, set */
/*        flag and return. */

/*<          if (iter .lt. maxit) go to 370 >*/
    if (iter < maxit) {
        goto L370;
    }
/*<             info = m >*/
    *info = m;
/*     ......exit */
/*<             go to 620 >*/
    goto L620;
/*<   370    continue >*/
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). */

/*<          do 390 ll = 1, m >*/
    i__1 = m;
    for (ll = 1; ll <= i__1; ++ll) {
/*<             l = m - ll >*/
        l = m - ll;
/*        ...exit */
/*<             if (l .eq. 0) go to 400 >*/
        if (l == 0) {
            goto L400;
        }
/*<             test = abs(s(l)) + abs(s(l+1)) >*/
        test = (r__1 = s[l], dabs(r__1)) + (r__2 = s[l + 1], dabs(r__2));
/*<             ztest = test + abs(e(l)) >*/
        ztest = test + (r__1 = e[l], dabs(r__1));
/*<             if (ztest .ne. test) go to 380 >*/
        if (ztest != test) {
            goto L380;
        }
/*<                e(l) = 0.0e0 >*/
        e[l] = (float)0.;
/*        ......exit */
/*<                go to 400 >*/
        goto L400;
/*<   380       continue >*/
L380:
/*<   390    continue >*/
/* L390: */
        ;
    }
/*<   400    continue >*/