rg.c

来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 1,457 行 · 第 1/3 页

/*          j          if the limit of 30*n iterations is exhausted */
/*                     while the j-th eigenvalue is being sought. */

/*     calls cdiv for complex division. */

/*     questions and comments should be directed to burton s. garbow, */
/*     mathematics and computer science div, argonne national laboratory */

/*     this version dated august 1983. */

/*     ------------------------------------------------------------------ */

    *ierr = 0;
    norm = 0.;
    k = 0;
/*     .......... store roots isolated by balanc */
/*                and compute matrix norm .......... */
    for (i = 0; i < *n; ++i) {
        for (j = k; j < *n; ++j) {
            norm += abs(h[i + j * *nm]);
        }
        k = i;
        if (i+1 < *low || i >= *igh) {
            wr[i] = h[i + i * *nm];
            wi[i] = 0.;
        }
    }

    en = *igh - 1;
    t = 0.;
    itn = *n * 30;
/*     .......... search for next eigenvalues .......... */
L60:
    if (en+1 < *low) {
        goto L340;
    }
    its = 0;
    na = en - 1;
    enm2 = na - 1;
/*     .......... look for single small sub-diagonal element */
/*                for l=en step -1 until low do -- .......... */
L70:
    for (l = en; l+1 > *low; --l) {
        s = abs(h[l-1 + (l-1) * *nm]) + abs(h[l + l * *nm]);
        if (s == 0.) {
            s = norm;
        }
        tst1 = s;
        tst2 = tst1 + abs(h[l + (l-1) * *nm]);
        if (tst2 == tst1) {
            break;
        }
    }
/*     .......... form shift .......... */
    x = h[en + en * *nm];
    if (l == en) {
        goto L270;
    }
    y = h[na + na * *nm];
    w = h[en + na * *nm] * h[na + en * *nm];
    if (l == na) {
        goto L280;
    }
    if (itn == 0) {
        goto L1000;
    }
    if (its != 10 && its != 20) {
        goto L130;
    }
/*     .......... form exceptional shift .......... */
    t += x;

    for (i = *low - 1; i <= en; ++i) {
        h[i + i * *nm] -= x;
    }

    s = abs(h[en + na * *nm]) + abs(h[na + enm2 * *nm]);
    x = s * .75;
    y = x;
    w = s * -.4375 * s;
L130:
    ++its;
    --itn;
/*     .......... look for two consecutive small */
/*                sub-diagonal elements. */
/*                for m=en-2 step -1 until l do -- .......... */
    for (m = enm2; m >= l; --m) {
        zz = h[m + m * *nm];
        r = x - zz;
        s = y - zz;
        p = (r * s - w) / h[m+1 + m * *nm] + h[m + (m+1) * *nm];
        q = h[m+1 + (m+1) * *nm] - zz - r - s;
        r = h[m+2 + (m+1) * *nm];
        s = abs(p) + abs(q) + abs(r);
        p /= s;
        q /= s;
        r /= s;
        if (m == l) {
            goto L150;
        }
        tst1 = abs(p) * (abs(h[m-1 + (m-1) * *nm]) + abs(zz) + abs(h[m+1 + (m+1) * *nm]));
        tst2 = tst1 + abs(h[m + (m-1) * *nm]) * (abs(q) + abs(r));
        if (tst2 == tst1) {
            goto L150;
        }
    }

L150:
    for (i = m+2; i <= en; ++i) {
        h[i + (i-2) * *nm] = 0.;
        if (i != m+2) {
            h[i + (i-3) * *nm] = 0.;
        }
    }
/*     .......... double qr step involving rows l to en and */
/*                columns m to en .......... */
    for (k = m; k <= na; ++k) {
        notlas = k != na;
        if (k == m) {
            goto L170;
        }
        p = h[k + (k-1) * *nm];
        q = h[k+1 + (k-1) * *nm];
        r = 0.;
        if (notlas) {
            r = h[k+2 + (k-1) * *nm];
        }
        x = abs(p) + abs(q) + abs(r);
        if (x == 0.) {
            continue;
        }
        p /= x;
        q /= x;
        r /= x;
L170:
        d__1 = sqrt(p * p + q * q + r * r);
        s = d_sign(&d__1, &p);
        if (k == m) {
            goto L180;
        }
        h[k + (k-1) * *nm] = -s * x;
        goto L190;
L180:
        if (l != m) {
            h[k + (k-1) * *nm] = -h[k + (k-1) * *nm];
        }
L190:
        p += s;
        x = p / s;
        y = q / s;
        zz = r / s;
        q /= p;
        r /= p;
        if (notlas) {
            goto L225;
        }
/*     .......... row modification .......... */
        for (j = k; j < *n; ++j) {
            p = h[k + j * *nm] + q * h[k+1 + j * *nm];
            h[k + j * *nm] -= p * x;
            h[k+1 + j * *nm] -= p * y;
        }

        j = min(en,k+3);
/*     .......... column modification .......... */
        for (i = 0; i <= j; ++i) {
            p = x * h[i + k * *nm] + y * h[i + (k+1) * *nm];
            h[i + k * *nm] -= p;
            h[i + (k+1) * *nm] -= p * q;
        }
/*     .......... accumulate transformations .......... */
        for (i = *low - 1; i < *igh; ++i) {
            p = x * z[i + k * *nm] + y * z[i + (k+1) * *nm];
            z[i + k * *nm] -= p;
            z[i + (k+1) * *nm] -= p * q;
        }
        continue;
L225:
/*     .......... row modification .......... */
        for (j = k; j < *n; ++j) {
            p = h[k + j * *nm] + q * h[k+1 + j * *nm] + r * h[k+2 + j * *nm];
            h[k + j * *nm] -= p * x;
            h[k+1 + j * *nm] -= p * y;
            h[k+2 + j * *nm] -= p * zz;
        }

        j = min(en,k+3);
/*     .......... column modification .......... */
        for (i = 0; i <= j; ++i) {
            p = x * h[i + k * *nm] + y * h[i + (k+1) * *nm] + zz * h[i + (k+2) * *nm];
            h[i + k * *nm] -= p;
            h[i + (k+1) * *nm] -= p * q;
            h[i + (k+2) * *nm] -= p * r;
        }
/*     .......... accumulate transformations .......... */
        for (i = *low - 1; i < *igh; ++i) {
            p = x * z[i + k * *nm] + y * z[i + (k+1) * *nm] + zz * z[i + (k+2) * *nm];
            z[i + k * *nm] -= p;
            z[i + (k+1) * *nm] -= p * q;
            z[i + (k+2) * *nm] -= p * r;
        }
    }

    goto L70;
/*     .......... one root found .......... */
L270:
    h[en + en * *nm] = x + t;
    wr[en] = h[en + en * *nm];
    wi[en] = 0.;
    en = na;
    goto L60;
/*     .......... two roots found .......... */
L280:
    p = (y - x) / 2.;
    q = p * p + w;
    zz = sqrt(abs(q));
    h[en + en * *nm] = x + t;
    x = h[en + en * *nm];
    h[na + na * *nm] = y + t;
    if (q < 0.) {
        goto L320;
    }
/*     .......... real pair .......... */
    zz = p + d_sign(&zz, &p);
    wr[na] = x + zz;
    wr[en] = wr[na];
    if (zz != 0.) {
        wr[en] = x - w / zz;
    }
    wi[na] = 0.;
    wi[en] = 0.;
    x = h[en + na * *nm];
    s = abs(x) + abs(zz);
    p = x / s;
    q = zz / s;
    r = sqrt(p * p + q * q);
    p /= r;
    q /= r;
/*     .......... row modification .......... */
    for (j = na; j < *n; ++j) {
        zz = h[na + j * *nm];
        h[na + j * *nm] = q * zz + p * h[en + j * *nm];
        h[en + j * *nm] = q * h[en + j * *nm] - p * zz;
    }
/*     .......... column modification .......... */
    for (i = 0; i <= en; ++i) {
        zz = h[i + na * *nm];
        h[i + na * *nm] = q * zz + p * h[i + en * *nm];
        h[i + en * *nm] = q * h[i + en * *nm] - p * zz;
    }
/*     .......... accumulate transformations .......... */
    for (i = *low - 1; i < *igh; ++i) {
        zz = z[i + na * *nm];
        z[i + na * *nm] = q * zz + p * z[i + en * *nm];
        z[i + en * *nm] = q * z[i + en * *nm] - p * zz;
    }

    goto L330;
/*     .......... complex pair .......... */
L320:
    wr[na] = x + p;
    wr[en] = x + p;
    wi[na] = zz;
    wi[en] = -zz;
L330:
    en = enm2;
    goto L60;
/*     .......... all roots found.  backsubstitute to find */
/*                vectors of upper triangular form .......... */
L340:
    if (norm == 0.) {
        goto L1001;
    }
/*     .......... for en=n step -1 until 1 do -- .......... */
    for (en = *n - 1; en >= 0; --en) {
        p = wr[en]; q = wi[en];
        na = en - 1;
        if (q < 0.) {
            goto L710;
        } else if (q != 0) {
            continue ;
        }
/*     .......... real vector .......... */
        m = en;
        h[en + en * *nm] = 1.;
        if (en == 0) {
            continue ;
        }
/*     .......... for i=en-1 step -1 until 1 do -- .......... */
        for (i = na; i >= 0; --i) {
            w = h[i + i * *nm] - p;
            r = 0.;

            for (j = m; j <= en; ++j) {
                r += h[i + j * *nm] * h[j + en * *nm];
            }

            if (wi[i] >= 0.) {
                goto L630;
            }
            zz = w;
            s = r;
            continue;
L630:
            m = i;
            if (wi[i] != 0.) {
                goto L640;
            }
            t = w;
            if (t != 0.) {
                goto L635;
            }
            tst1 = norm;
            t = tst1;
L632:
            t *= .01;
            tst2 = norm + t;
            if (tst2 > tst1) {
                goto L632;
            }
L635:
            h[i + en * *nm] = -r / t;
            goto L680;
/*     .......... solve real equations .......... */
L640:
            x = h[i + (i+1) * *nm];
            y = h[i+1 + i * *nm];
            q = (wr[i] - p) * (wr[i] - p) + wi[i] * wi[i];
            t = (x * s - zz * r) / q;
            h[i + en * *nm] = t;
            if (abs(x) <= abs(zz)) {
                goto L650;
            }
            h[i+1 + en * *nm] = (-r - w * t) / x;
            goto L680;
L650:
            h[i+1 + en * *nm] = (-s - y * t) / zz;

/*     .......... overflow control .......... */
L680:
            t = abs(h[i + en * *nm]);
            if (t == 0.) {
                continue;
            }
            tst1 = t;
            tst2 = tst1 + 1. / tst1;
            if (tst2 > tst1) {
                continue;
            }
            for (j = i; j <= en; ++j) {
                h[j + en * *nm] /= t;
            }
        }
/*     .......... end real vector .......... */
        continue;
/*     .......... complex vector .......... */
L710:
        m = na;
/*     .......... last vector component chosen imaginary so that */
/*                eigenvector matrix is triangular .......... */
        if (abs(h[en + na * *nm]) <= abs(h[na + en * *nm])) {
            goto L720;
        }
        h[na + na * *nm] = q / h[en + na * *nm];
        h[na + en * *nm] = -(h[en + en * *nm] - p) / h[en + na * *nm];
        goto L730;
L720:
        d__1 = -h[na + en * *nm];
        d__2 = h[na + na * *nm] - p;
        cdiv_(&c_b130, &d__1, &d__2, &q, &h[na + na * *nm], &h[na + en * *nm]);
L730:
        h[en + na * *nm] = 0.;
        h[en + en * *nm] = 1.;
        enm2 = na - 1;
/*     .......... for i=en-2 step -1 until 1 do -- .......... */
        if (na != 0)
        for (i = enm2; i >= 0; --i) {
            w = h[i + i * *nm] - p;
            ra = 0.;
            sa = 0.;

            for (j = m; j <= en; ++j) {
                ra += h[i + j * *nm] * h[j + na * *nm];
                sa += h[i + j * *nm] * h[j + en * *nm];
            }

            if (wi[i] >= 0.) {
                goto L770;
            }
            zz = w;
            r = ra;
            s = sa;
            continue;
L770:
            m = i;
            if (wi[i] != 0.) {
                goto L780;
            }
            d__1 = -ra;
            d__2 = -sa;
            cdiv_(&d__1, &d__2, &w, &q, &h[i + na * *nm], &h[i + en * *nm]);
            goto L790;
/*     .......... solve complex equations .......... */
L780:
            x = h[i + (i+1) * *nm];
            y = h[i+1 + i * *nm];
            vr = (wr[i] - p) * (wr[i] - p) + wi[i] * wi[i] - q * q;
            vi = (wr[i] - p) * 2. * q;
            if (vr != 0. || vi != 0.) {
                goto L784;
            }
            tst1 = norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(zz));
            vr = tst1;
L783:
            vr *= .01;
            tst2 = tst1 + vr;
            if (tst2 > tst1) {
                goto L783;
            }
L784:
            d__1 = x * r - zz * ra + q * sa;
            d__2 = x * s - zz * sa - q * ra;
            cdiv_(&d__1, &d__2, &vr, &vi, &h[i + na * *nm], &h[i + en * *nm]);
            if (abs(x) <= abs(zz) + abs(q)) {
                goto L785;
            }
            h[i+1 + na * *nm] = (-ra - w * h[i + na * *nm] + q * h[i + en * *nm]) / x;
            h[i+1 + en * *nm] = (-sa - w * h[i + en * *nm] - q * h[i + na * *nm]) / x;
            goto L790;
L785:
            d__1 = -r - y * h[i + na * *nm];
            d__2 = -s - y * h[i + en * *nm];
            cdiv_(&d__1, &d__2, &zz, &q, &h[i+1 + na * *nm], &h[i+1 + en * *nm]);

/*     .......... overflow control .......... */
L790:
            d__1 = abs(h[i + na * *nm]), d__2 = abs(h[i + en * *nm]);
            t = max(d__1,d__2);
            if (t == 0.) {
                continue;
            }
            tst1 = t;
            tst2 = tst1 + 1. / tst1;
            if (tst2 <= tst1)
            for (j = i; j <= en; ++j) {
                h[j + na * *nm] /= t;
                h[j + en * *nm] /= t;
            }
        }
/*     .......... end complex vector .......... */
    }
/*     .......... end back substitution. */
/*                vectors of isolated roots .......... */
    for (i = 0; i < *n; ++i) {
        if (i+1 < *low || i >= *igh)
        for (j = i; j < *n; ++j) {
            z[i + j * *nm] = h[i + j * *nm];
        }
    }
/*     .......... multiply by transformation matrix to give */
/*                vectors of original full matrix. */
/*                for j=n step -1 until low do -- .......... */
    for (j = *n - 1; j+1 >= *low; --j) {
        m = min(j,*igh-1);

        for (i = *low - 1; i < *igh; ++i) {
            zz = 0.;

            for (k = *low - 1; k <= m; ++k) {
                zz += z[i + k * *nm] * h[k + j * *nm];
            }

            z[i + j * *nm] = zz;
        }
    }

    goto L1001;
/*     .......... set error -- all eigenvalues have not */
/*                converged after 30*n iterations .......... */
L1000:
    *ierr = en-1;
L1001:
    return;
} /* hqr2_ */