ssvdc.c

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

L400:
/*<          if (l .ne. m - 1) go to 410 >*/
    if (l != m - 1) {
        goto L410;
    }
/*<             kase = 4 >*/
    kase = 4;
/*<          go to 480 >*/
    goto L480;
/*<   410    continue >*/
L410:
/*<             lp1 = l + 1 >*/
    lp1 = l + 1;
/*<             mp1 = m + 1 >*/
    mp1 = m + 1;
/*<             do 430 lls = lp1, mp1 >*/
    i__1 = mp1;
    for (lls = lp1; lls <= i__1; ++lls) {
/*<                ls = m - lls + lp1 >*/
        ls = m - lls + lp1;
/*           ...exit */
/*<                if (ls .eq. l) go to 440 >*/
        if (ls == l) {
            goto L440;
        }
/*<                test = 0.0e0 >*/
        test = (float)0.;
/*<                if (ls .ne. m) test = test + abs(e(ls)) >*/
        if (ls != m) {
            test += (r__1 = e[ls], dabs(r__1));
        }
/*<                if (ls .ne. l + 1) test = test + abs(e(ls-1)) >*/
        if (ls != l + 1) {
            test += (r__1 = e[ls - 1], dabs(r__1));
        }
/*<                ztest = test + abs(s(ls)) >*/
        ztest = test + (r__1 = s[ls], dabs(r__1));
/*<                if (ztest .ne. test) go to 420 >*/
        if (ztest != test) {
            goto L420;
        }
/*<                   s(ls) = 0.0e0 >*/
        s[ls] = (float)0.;
/*           ......exit */
/*<                   go to 440 >*/
        goto L440;
/*<   420          continue >*/
L420:
/*<   430       continue >*/
/* L430: */
        ;
    }
/*<   440       continue >*/
L440:
/*<             if (ls .ne. l) go to 450 >*/
    if (ls != l) {
        goto L450;
    }
/*<                kase = 3 >*/
    kase = 3;
/*<             go to 470 >*/
    goto L470;
/*<   450       continue >*/
L450:
/*<             if (ls .ne. m) go to 460 >*/
    if (ls != m) {
        goto L460;
    }
/*<                kase = 1 >*/
    kase = 1;
/*<             go to 470 >*/
    goto L470;
/*<   460       continue >*/
L460:
/*<                kase = 2 >*/
    kase = 2;
/*<                l = ls >*/
    l = ls;
/*<   470       continue >*/
L470:
/*<   480    continue >*/
L480:
/*<          l = l + 1 >*/
    ++l;

/*        perform the task indicated by kase. */

/*<          go to (490,520,540,570), kase >*/
    switch (kase) {
        case 1:  goto L490;
        case 2:  goto L520;
        case 3:  goto L540;
        case 4:  goto L570;
    }

/*        deflate negligible s(m). */

/*<   490    continue >*/
L490:
/*<             mm1 = m - 1 >*/
    mm1 = m - 1;
/*<             f = e(m-1) >*/
    f = e[m - 1];
/*<             e(m-1) = 0.0e0 >*/
    e[m - 1] = (float)0.;
/*<             do 510 kk = l, mm1 >*/
    i__1 = mm1;
    for (kk = l; kk <= i__1; ++kk) {
/*<                k = mm1 - kk + l >*/
        k = mm1 - kk + l;
/*<                t1 = s(k) >*/
        t1 = s[k];
/*<                call srotg(t1,f,cs,sn) >*/
        srotg_(&t1, &f, &cs, &sn);
/*<                s(k) = t1 >*/
        s[k] = t1;
/*<                if (k .eq. l) go to 500 >*/
        if (k == l) {
            goto L500;
        }
/*<                   f = -sn*e(k-1) >*/
        f = -sn * e[k - 1];
/*<                   e(k-1) = cs*e(k-1) >*/
        e[k - 1] = cs * e[k - 1];
/*<   500          continue >*/
L500:
/*<                if (wantv) call srot(p,v(1,k),1,v(1,m),1,cs,sn) >*/
        if (wantv) {
            srot_(p, &v[k * v_dim1 + 1], &c__1, &v[m * v_dim1 + 1], &c__1, &
                    cs, &sn);
        }
/*<   510       continue >*/
/* L510: */
    }
/*<          go to 610 >*/
    goto L610;

/*        split at negligible s(l). */

/*<   520    continue >*/
L520:
/*<             f = e(l-1) >*/
    f = e[l - 1];
/*<             e(l-1) = 0.0e0 >*/
    e[l - 1] = (float)0.;
/*<             do 530 k = l, m >*/
    i__1 = m;
    for (k = l; k <= i__1; ++k) {
/*<                t1 = s(k) >*/
        t1 = s[k];
/*<                call srotg(t1,f,cs,sn) >*/
        srotg_(&t1, &f, &cs, &sn);
/*<                s(k) = t1 >*/
        s[k] = t1;
/*<                f = -sn*e(k) >*/
        f = -sn * e[k];
/*<                e(k) = cs*e(k) >*/
        e[k] = cs * e[k];
/*<                if (wantu) call srot(n,u(1,k),1,u(1,l-1),1,cs,sn) >*/
        if (wantu) {
            srot_(n, &u[k * u_dim1 + 1], &c__1, &u[(l - 1) * u_dim1 + 1], &
                    c__1, &cs, &sn);
        }
/*<   530       continue >*/
/* L530: */
    }
/*<          go to 610 >*/
    goto L610;

/*        perform one qr step. */

/*<   540    continue >*/
L540:

/*           calculate the shift. */

/*<    >*/
/* Computing MAX */
    r__6 = (r__1 = s[m], dabs(r__1)), r__7 = (r__2 = s[m - 1], dabs(r__2)), 
            r__6 = max(r__6,r__7), r__7 = (r__3 = e[m - 1], dabs(r__3)), r__6 
            = max(r__6,r__7), r__7 = (r__4 = s[l], dabs(r__4)), r__6 = max(
            r__6,r__7), r__7 = (r__5 = e[l], dabs(r__5));
    scale = dmax(r__6,r__7);
/*<             sm = s(m)/scale >*/
    sm = s[m] / scale;
/*<             smm1 = s(m-1)/scale >*/
    smm1 = s[m - 1] / scale;
/*<             emm1 = e(m-1)/scale >*/
    emm1 = e[m - 1] / scale;
/*<             sl = s(l)/scale >*/
    sl = s[l] / scale;
/*<             el = e(l)/scale >*/
    el = e[l] / scale;
/*<             b = ((smm1 + sm)*(smm1 - sm) + emm1**2)/2.0e0 >*/
/* Computing 2nd power */
    r__1 = emm1;
    b = ((smm1 + sm) * (smm1 - sm) + r__1 * r__1) / (float)2.;
/*<             c = (sm*emm1)**2 >*/
/* Computing 2nd power */
    r__1 = sm * emm1;
    c__ = r__1 * r__1;
/*<             shift = 0.0e0 >*/
    shift = (float)0.;
/*<             if (b .eq. 0.0e0 .and. c .eq. 0.0e0) go to 550 >*/
    if (b == (float)0. && c__ == (float)0.) {
        goto L550;
    }
/*<                shift = sqrt(b**2+c) >*/
/* Computing 2nd power */
    r__1 = b;
    shift = sqrt(r__1 * r__1 + c__);
/*<                if (b .lt. 0.0e0) shift = -shift >*/
    if (b < (float)0.) {
        shift = -shift;
    }
/*<                shift = c/(b + shift) >*/
    shift = c__ / (b + shift);
/*<   550       continue >*/
L550:
/*<             f = (sl + sm)*(sl - sm) + shift >*/
    f = (sl + sm) * (sl - sm) + shift;
/*<             g = sl*el >*/
    g = sl * el;

/*           chase zeros. */

/*<             mm1 = m - 1 >*/
    mm1 = m - 1;
/*<             do 560 k = l, mm1 >*/
    i__1 = mm1;
    for (k = l; k <= i__1; ++k) {
/*<                call srotg(f,g,cs,sn) >*/
        srotg_(&f, &g, &cs, &sn);
/*<                if (k .ne. l) e(k-1) = f >*/
        if (k != l) {
            e[k - 1] = f;
        }
/*<                f = cs*s(k) + sn*e(k) >*/
        f = cs * s[k] + sn * e[k];
/*<                e(k) = cs*e(k) - sn*s(k) >*/
        e[k] = cs * e[k] - sn * s[k];
/*<                g = sn*s(k+1) >*/
        g = sn * s[k + 1];
/*<                s(k+1) = cs*s(k+1) >*/
        s[k + 1] = cs * s[k + 1];
/*<                if (wantv) call srot(p,v(1,k),1,v(1,k+1),1,cs,sn) >*/
        if (wantv) {
            srot_(p, &v[k * v_dim1 + 1], &c__1, &v[(k + 1) * v_dim1 + 1], &
                    c__1, &cs, &sn);
        }
/*<                call srotg(f,g,cs,sn) >*/
        srotg_(&f, &g, &cs, &sn);
/*<                s(k) = f >*/
        s[k] = f;
/*<                f = cs*e(k) + sn*s(k+1) >*/
        f = cs * e[k] + sn * s[k + 1];
/*<                s(k+1) = -sn*e(k) + cs*s(k+1) >*/
        s[k + 1] = -sn * e[k] + cs * s[k + 1];
/*<                g = sn*e(k+1) >*/
        g = sn * e[k + 1];
/*<                e(k+1) = cs*e(k+1) >*/
        e[k + 1] = cs * e[k + 1];
/*<    >*/
        if (wantu && k < *n) {
            srot_(n, &u[k * u_dim1 + 1], &c__1, &u[(k + 1) * u_dim1 + 1], &
                    c__1, &cs, &sn);
        }
/*<   560       continue >*/
/* L560: */
    }
/*<             e(m-1) = f >*/
    e[m - 1] = f;
/*<             iter = iter + 1 >*/
    ++iter;
/*<          go to 610 >*/
    goto L610;

/*        convergence. */

/*<   570    continue >*/
L570:

/*           make the singular value  positive. */

/*<             if (s(l) .ge. 0.0e0) go to 580 >*/
    if (s[l] >= (float)0.) {
        goto L580;
    }
/*<                s(l) = -s(l) >*/
    s[l] = -s[l];
/*<                if (wantv) call sscal(p,-1.0e0,v(1,l),1) >*/
    if (wantv) {
        sscal_(p, &c_b44, &v[l * v_dim1 + 1], &c__1);
    }
/*<   580       continue >*/
L580:

/*           order the singular value. */

/*<   590       if (l .eq. mm) go to 600 >*/
L590:
    if (l == mm) {
        goto L600;
    }
/*           ...exit */
/*<                if (s(l) .ge. s(l+1)) go to 600 >*/
    if (s[l] >= s[l + 1]) {
        goto L600;
    }
/*<                t = s(l) >*/
    t = s[l];
/*<                s(l) = s(l+1) >*/
    s[l] = s[l + 1];
/*<                s(l+1) = t >*/
    s[l + 1] = t;
/*<    >*/
    if (wantv && l < *p) {
        sswap_(p, &v[l * v_dim1 + 1], &c__1, &v[(l + 1) * v_dim1 + 1], &c__1);
    }
/*<    >*/
    if (wantu && l < *n) {
        sswap_(n, &u[l * u_dim1 + 1], &c__1, &u[(l + 1) * u_dim1 + 1], &c__1);
    }
/*<                l = l + 1 >*/
    ++l;
/*<             go to 590 >*/
    goto L590;
/*<   600       continue >*/
L600:
/*<             iter = 0 >*/
    iter = 0;
/*<             m = m - 1 >*/
    --m;
/*<   610    continue >*/
L610:
/*<       go to 360 >*/
    goto L360;
/*<   620 continue >*/
L620:
/*<       return >*/
    return 0;
/*<       end >*/
} /* ssvdc_ */

#ifdef __cplusplus
        }
#endif