rg.c

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

        int_[m] = i+1;
        if (i == m) {
            goto L130;
        }
/*     .......... interchange rows and columns of a .......... */
        for (j = mm1; j < *n; ++j) {
            y = a[i + j * *nm];
            a[i + j * *nm] = a[m + j * *nm];
            a[m + j * *nm] = y;
        }

        for (j = 0; j < *igh; ++j) {
            y = a[j + i * *nm];
            a[j + i * *nm] = a[j + m * *nm];
            a[j + m * *nm] = y;
        }
/*     .......... end interchange .......... */
L130:
        if (x != 0.)
        for (i = m+1; i < *igh; ++i) {
            y = a[i + mm1 * *nm];
            if (y == 0.) {
                continue;
            }
            y /= x;
            a[i + mm1 * *nm] = y;

            for (j = m; j < *n; ++j) {
                a[i + j * *nm] -= y * a[m + j * *nm];
            }

            for (j = 0; j < *igh; ++j) {
                a[j + m * *nm] += y * a[j + i * *nm];
            }
        }
    }
} /* elmhes_ */

/* Subroutine */ void eltran_(nm, n, low, igh, a, int_, z)
integer *nm, *n, *low, *igh;
doublereal *a;
integer *int_;
doublereal *z;
{
    /* Local variables */
    static integer i, j, kl, mp;

/*     this subroutine is a translation of the algol procedure elmtrans, */
/*     num. math. 16, 181-204(1970) by peters and wilkinson. */
/*     handbook for auto. comp., vol.ii-linear algebra, 372-395(1971). */

/*     this subroutine accumulates the stabilized elementary */
/*     similarity transformations used in the reduction of a */
/*     real general matrix to upper hessenberg form by  elmhes. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        low and igh are integers determined by the balancing */
/*          subroutine  balanc.  if  balanc  has not been used, */
/*          set low=1, igh=n. */

/*        a contains the multipliers which were used in the */
/*          reduction by  elmhes  in its lower triangle */
/*          below the subdiagonal. */

/*        int contains information on the rows and columns */
/*          interchanged in the reduction by  elmhes. */
/*          only elements low through igh are used. */

/*     on output */

/*        z contains the transformation matrix produced in the */
/*          reduction by  elmhes. */

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

/*     this version dated august 1983. */

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

/*     .......... initialize z to identity matrix .......... */

    for (j = 0; j < *n; ++j) {
        for (i = 0; i < *n; ++i) {
            z[i + j * *nm] = 0.;
        }

        z[j + j * *nm] = 1.;
    }

    kl = *igh - *low - 1;
    if (kl < 1) {
        return;
    }
/*     .......... for mp=igh-1 step -1 until low+1 do -- .......... */
    for (mp = *igh - 2; mp > *igh -kl-2; --mp) {
        for (i = mp+1; i < *igh; ++i) {
            z[i + mp * *nm] = a[i + (mp - 1) * *nm];
        }

        i = int_[mp] - 1;
        if (i == mp) {
            continue;
        }

        for (j = mp; j < *igh; ++j) {
            z[mp + j * *nm] = z[i + j * *nm];
            z[i + j * *nm] = 0.;
        }

        z[i + mp * *nm] = 1.;
    }
} /* eltran_ */

/* Subroutine */ void hqr_(nm, n, low, igh, h, wr, wi, ierr)
integer *nm, *n, *low, *igh;
doublereal *h, *wr, *wi;
integer *ierr;
{
    /* System generated locals */
    doublereal d__1;

    /* Local variables */
    static doublereal norm;
    static integer i, j, k, l, m;
    static doublereal p, q, r, s, t, w, x, y;
    static integer na, en;
    static doublereal zz;
    static logical notlas;
    static integer itn, its, enm2;
    static doublereal tst1, tst2;

/*  RESTORED CORRECT INDICES OF LOOPS (200,210,230,240). (9/29/89 BSG) */

/*     this subroutine is a translation of the algol procedure hqr, */
/*     num. math. 14, 219-231(1970) by martin, peters, and wilkinson. */
/*     handbook for auto. comp., vol.ii-linear algebra, 359-371(1971). */

/*     this subroutine finds the eigenvalues of a real */
/*     upper hessenberg matrix by the qr method. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        low and igh are integers determined by the balancing */
/*          subroutine  balanc.  if  balanc  has not been used, */
/*          set low=1, igh=n. */

/*        h contains the upper hessenberg matrix.  information about */
/*          the transformations used in the reduction to hessenberg */
/*          form by  elmhes  or  orthes, if performed, is stored */
/*          in the remaining triangle under the hessenberg matrix. */

/*     on output */

/*        h has been destroyed.  therefore, it must be saved */
/*          before calling  hqr  if subsequent calculation and */
/*          back transformation of eigenvectors is to be performed. */

/*        wr and wi contain the real and imaginary parts, */
/*          respectively, of the eigenvalues.  the eigenvalues */
/*          are unordered except that complex conjugate pairs */
/*          of values appear consecutively with the eigenvalue */
/*          having the positive imaginary part first.  if an */
/*          error exit is made, the eigenvalues should be correct */
/*          for indices ierr+1,...,n. */

/*        ierr is set to */
/*          zero       for normal return, */
/*          j          if the limit of 30*n iterations is exhausted */
/*                     while the j-th eigenvalue is being sought. */

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

/*     this version dated september 1989. */

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

    *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 L1001;
    }
    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 <= en; ++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 = l; 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;
        }
        continue;
L225:
/*     .......... row modification .......... */
        for (j = k; j <= en; ++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 = l; 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;
        }
    }

    goto L70;
/*     .......... one root found .......... */
L270:
    wr[en] = x + t;
    wi[en] = 0.;
    en = na;
    goto L60;
/*     .......... two roots found .......... */
L280:
    p = (y - x) / 2.;
    q = p * p + w;
    zz = sqrt(abs(q));
    x += 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.;
    goto L330;
/*     .......... complex pair .......... */
L320:
    wr[na] = x + p;
    wr[en] = x + p;
    wi[na] = zz;
    wi[en] = -zz;
L330:
    en = enm2;
    goto L60;
/*     .......... set error -- all eigenvalues have not */
/*                converged after 30*n iterations .......... */
L1000:
    *ierr = en-1;
L1001:
    return;
} /* hqr_ */

/* Subroutine */ void hqr2_(nm, n, low, igh, h, wr, wi, z, ierr)
integer *nm, *n, *low, *igh;
doublereal *h, *wr, *wi, *z;
integer *ierr;
{
    /* System generated locals */
    doublereal d__1, d__2;

    /* Local variables */
    static doublereal norm;
    static integer i, j, k, l, m;
    static doublereal p, q, r, s, t, w, x, y;
    static integer na, en;
    static doublereal ra, sa;
    static doublereal vi, vr, zz;
    static logical notlas;
    static integer itn, its, enm2;
    static doublereal tst1, tst2;


/*     this subroutine is a translation of the algol procedure hqr2, */
/*     num. math. 16, 181-204(1970) by peters and wilkinson. */
/*     handbook for auto. comp., vol.ii-linear algebra, 372-395(1971). */

/*     this subroutine finds the eigenvalues and eigenvectors */
/*     of a real upper hessenberg matrix by the qr method.  the */
/*     eigenvectors of a real general matrix can also be found */
/*     if  elmhes  and  eltran  or  orthes  and  ortran  have */
/*     been used to reduce this general matrix to hessenberg form */
/*     and to accumulate the similarity transformations. */

/*     on input */

/*        nm must be set to the row dimension of two-dimensional */
/*          array parameters as declared in the calling program */
/*          dimension statement. */

/*        n is the order of the matrix. */

/*        low and igh are integers determined by the balancing */
/*          subroutine  balanc.  if  balanc  has not been used, */
/*          set low=1, igh=n. */

/*        h contains the upper hessenberg matrix. */

/*        z contains the transformation matrix produced by  eltran */
/*          after the reduction by  elmhes, or by  ortran  after the */
/*          reduction by  orthes, if performed.  if the eigenvectors */
/*          of the hessenberg matrix are desired, z must contain the */
/*          identity matrix. */

/*     on output */

/*        h has been destroyed. */

/*        wr and wi contain the real and imaginary parts, */
/*          respectively, of the eigenvalues.  the eigenvalues */
/*          are unordered except that complex conjugate pairs */
/*          of values appear consecutively with the eigenvalue */
/*          having the positive imaginary part first.  if an */
/*          error exit is made, the eigenvalues should be correct */
/*          for indices ierr+1,...,n. */

/*        z contains the real and imaginary parts of the eigenvectors. */
/*          if the i-th eigenvalue is real, the i-th column of z */
/*          contains its eigenvector.  if the i-th eigenvalue is complex */
/*          with positive imaginary part, the i-th and (i+1)-th */
/*          columns of z contain the real and imaginary parts of its */
/*          eigenvector.  the eigenvectors are unnormalized.  if an */
/*          error exit is made, none of the eigenvectors has been found.  */

/*        ierr is set to */
/*          zero       for normal return, */