dsvdc.c

InsightToolkit-1.4.0(有大量的优化算法程序)

        s[nct] = x[nct + nct * *ldx];
    }
    if (*n < m+1) {
        s[m] = 0.;
    }
    if (nrt < m) {
        e[nrt] = x[nrt + m * *ldx];
    }
    e[m] = 0.;

/*     if required, generate u. */

    if (wantu)
    for (j = nct; j < ncu; ++j) {
        for (i = 0; i < *n; ++i) {
            u[i + j * *ldu] = 0.;
        }
        u[j + j * *ldu] = 1.;
    }
    if (wantu)
    for (l = nct-1; l >= 0; --l) {
        if (s[l] == 0.) {
            for (i = 0; i < *n; ++i) {
                u[i + l * *ldu] = 0.;
            }
            u[l + l * *ldu] = 1.;
            continue;
        }
        for (j = l+1; j < ncu; ++j) {
            i__1 = *n - l;
            t = -ddot_(&i__1, &u[l + l * *ldu], &c__1, &u[l + j * *ldu], &c__1) / u[l + l * *ldu];
            daxpy_(&i__1, &t, &u[l + l * *ldu], &c__1, &u[l + j * *ldu], &c__1);
        }
        i__1 = *n - l;
        dscal_(&i__1, &c_m1, &u[l + l * *ldu], &c__1);
        u[l + l * *ldu] += 1.;
        for (i = 0; i < l; ++i) {
            u[i + l * *ldu] = 0.;
        }
    }

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

    if (wantv)
    for (l = *p-1; l >= 0; --l) {
        lp1 = l+1;
        if (l < nrt && e[l] != 0.)
        for (j = lp1; j < *p; ++j) {
            i__1 = *p - lp1;
            t = -ddot_(&i__1, &v[lp1 + l * *ldv], &c__1, &v[lp1 + j * *ldv], &c__1) / v[lp1 + l * *ldv];
            daxpy_(&i__1, &t, &v[lp1 + l * *ldv], &c__1, &v[lp1 + j * *ldv], &c__1);
        }
        for (i = 0; i < *p; ++i) {
            v[i + l * *ldv] = 0.;
        }
        v[l + l * *ldv] = 1.;
    }

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

    mm = m;
    iter = 0;
L360:

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

    if (m < 0) {
        return;
    }

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

    if (iter >= maxit) {
        *info = m+1;
        return;
    }

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

    for (l = m-1; l >= 0; --l) {
        test = abs(s[l]) + abs(s[l+1]);
        ztest = test + abs(e[l]);
        if (fsm_ieee_doubles_equal(&ztest, &test)) {
/* WAS: if (ztest == test) { */
            e[l] = 0.;
            break;
        }
    }
    if (l == m-1) {
        kase = 4;
        goto L480;
    }
    for (ls = m; ls > l; --ls) {
        test = 0.;
        if (ls != m) {
            test += abs(e[ls]);
        }
        if (ls != l+1) {
            test += abs(e[ls-1]);
        }
        ztest = test + abs(s[ls]);
        if (fsm_ieee_doubles_equal(&ztest, &test)) {
/* WAS: if (ztest == test) { */
            s[ls] = 0.;
            break;
        }
    }
    if (ls == l) {
        kase = 3;
    }
    else if (ls == m) {
        kase = 1;
    }
    else {
        kase = 2;
        l = ls;
    }
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:
    f = e[m-1];
    e[m-1] = 0.;
    for (k = m-1; k >= l; --k) {
        t1 = s[k];
        drotg_(&t1, &f, &cs, &sn);
        s[k] = t1;
        if (k != l) {
            f = -sn * e[k-1];
            e[k-1] *= cs;
        }
        if (wantv) {
            drot_(p, &v[k * *ldv], &c__1, &v[m * *ldv], &c__1, &cs, &sn);
        }
    }
    goto L360;

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

L520:
    f = e[l-1];
    e[l-1] = 0.;
    for (k = l; k <= m; ++k) {
        t1 = s[k];
        drotg_(&t1, &f, &cs, &sn);
        s[k] = t1;
        f = -sn * e[k];
        e[k] *= cs;
        if (wantu) {
            drot_(n, &u[k * *ldu], &c__1, &u[(l-1) * *ldu], &c__1, &cs, &sn);
        }
    }
    goto L360;

/*        perform one qr step. */

L540:

/*           calculate the shift. */

    scale = max(max(max(max(abs(s[m]),abs(s[m-1])),abs(e[m-1])),abs(s[l])),abs(e[l]));
    sm = s[m] / scale;
    smm1 = s[m-1] / scale;
    emm1 = e[m-1] / scale;
    sl = s[l] / scale;
    el = e[l] / scale;
    b = ((smm1 + sm) * (smm1 - sm) + emm1 * emm1) / 2.;
    c = sm * emm1; c *= c;
    if (b == 0. && c == 0.) {
        shift = 0.;
    }
    else {
        shift = sqrt(b * b + c);
        if (b < 0.) {
            shift = -shift;
        }
        shift = c / (b + shift);
    }
    f = (sl + sm) * (sl - sm) + shift;
    g = sl * el;

/*           chase zeros. */

    for (k = l; k < m; ++k) {
        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;
        if (wantv) {
            drot_(p, &v[k * *ldv], &c__1, &v[(k+1) * *ldv], &c__1, &cs, &sn);
        }
        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;
        if (wantu && k+1 < *n) {
            drot_(n, &u[k * *ldu], &c__1, &u[(k+1) * *ldu], &c__1, &cs, &sn);
        }
    }
    e[m-1] = f;
    ++iter;
    goto L360;

/*        convergence. */

L570:

/*           make the singular value  positive. */

    if (s[l] < 0.) {
        s[l] = -s[l];
        if (wantv) {
            dscal_(p, &c_m1, &v[l * *ldv], &c__1);
        }
    }

/*           order the singular value. */

L590:
    if (l == mm) {
        goto L600;
    }
    if (s[l] >= s[l+1]) {
        goto L600;
    }
    t = s[l];
    s[l] = s[l+1];
    ++l;
    s[l] = t;
    if (wantv && l < *p) {
        dswap_(p, &v[(l-1) * *ldv], &c__1, &v[l * *ldv], &c__1);
    }
    if (wantu && l < *n) {
        dswap_(n, &u[(l-1) * *ldu], &c__1, &u[l * *ldu], &c__1);
    }
    goto L590;
L600:
    iter = 0;
    --m;
    goto L360;
} /* dsvdc_ */

int fsm_ieee_doubles_equal(const doublereal *x, const doublereal *y)
{
  return *x == *y;
}