ztrevc.c

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

        *info = -11;
    }
    if (*info != 0) {
        i__1 = -(*info);
        xerbla_("ZTREVC", &i__1);
        return;
    }

/*     Quick return if possible. */

    if (*n == 0) {
        return;
    }

/*     Set the constants to control overflow. */

    unfl = dlamch_("Safe minimum");
    ovfl = 1. / unfl;
    dlabad_(&unfl, &ovfl);
    ulp = dlamch_("Precision");
    smlnum = unfl * (*n / ulp);

/*     Store the diagonal elements of T in working array WORK. */

    for (i = 0; i < *n; ++i) {
        i__1 = i + *n;
        i__2 = i + i * *ldt; /* index [i,i] */
        work[i__1].r = t[i__2].r, work[i__1].i = t[i__2].i;
    }

/*     Compute 1-norm of each column of strictly upper triangular */
/*     part of T to control overflow in triangular solver. */

    rwork[0] = 0.;
    for (j = 1; j < *n; ++j) {
        rwork[j] = dzasum_(&j, &t[j * *ldt], &c__1);
    }

    if (rightv) {

/*        Compute right eigenvectors. */

        is = *m - 1;
        for (ki = *n - 1; ki >= 0; --ki) {

            if (somev) {
                if (! select[ki]) {
                    continue; /* next ki */
                }
            }
            i__1 = ki + ki * *ldt; /* index [ki,ki] */
            smin = ulp * (abs(t[i__1].r) + abs(t[i__1].i));
            smin = max(smin, smlnum);

            work[0].r = 1., work[0].i = 0.;

/*           Form right-hand side. */

            for (k = 0; k < ki; ++k) {
                i__1 = k + ki * *ldt; /* index [k,ki] */
                work[k].r = -t[i__1].r, work[k].i = -t[i__1].i;
            }

/*           Solve the triangular system: */
/*              (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK. */

            for (k = 0; k < ki; ++k) {
                i__1 = k + k * *ldt; /* index [k,k] */
                i__2 = ki + ki * *ldt; /* index [ki,ki] */
                t[i__1].r -= t[i__2].r,
                t[i__1].i -= t[i__2].i;
                if (abs(t[i__1].r) + abs(t[i__1].i) < smin) {
                    t[i__1].r = smin, t[i__1].i = 0.;
                }
            }

            if (ki > 0) {
                zlatrs_("Upper", "No transpose", "Non-unit", "Y", &ki, t, ldt, work, &scale, rwork, info);
                work[ki].r = scale, work[ki].i = 0.;
            }

/*           Copy the vector x or Q*x to VR and normalize. */

            if (! over) {
                k = ki+1;
                zcopy_(&k, work, &c__1, &vr[is * *ldvr], &c__1);

                ii = izamax_(&k, &vr[is * *ldvr], &c__1);
                i__1 = ii-1 + is * *ldvr; /* index [ii-1,is] */
                remax = 1. / (abs(vr[i__1].r) + abs(vr[i__1].i));
                zdscal_(&k, &remax, &vr[is * *ldvr], &c__1);

                for (k = ki+1; k < *n; ++k) {
                    i__1 = k + is * *ldvr; /* index [k,is] */
                    vr[i__1].r = 0., vr[i__1].i = 0.;
                }
            } else {
                if (ki > 0) {
                    z__1.r = scale, z__1.i = 0.;
                    zgemv_("N", n, &ki, &c_b17, vr, ldvr, work, &c__1, &z__1, &vr[ki * *ldvr], &c__1);
                }

                ii = izamax_(n, &vr[ki * *ldvr], &c__1);
                i__1 = ii-1 + ki * *ldvr; /* index [ii-1,ki] */
                remax = 1. / (abs(vr[i__1].r) + abs(vr[i__1].i));
                zdscal_(n, &remax, &vr[ki * *ldvr], &c__1);
            }

/*           Set back the original diagonal elements of T. */

            for (k = 0; k < ki; ++k) {
                i__1 = k + k * *ldt; /* index [k,k] */
                i__2 = k + *n;
                t[i__1].r = work[i__2].r, t[i__1].i = work[i__2].i;
            }

            --is;
        }
    }

    if (leftv) {

/*        Compute left eigenvectors. */

        is = 0;
        for (ki = 0; ki < *n; ++ki) {
            if (somev) {
                if (! select[ki]) {
                    continue; /* next ki */
                }
            }
            i__1 = ki + ki * *ldt; /* index [ki,ki] */
            smin = ulp * (abs(t[i__1].r) + abs(t[i__1].i));
            smin = max(smin, smlnum);
            work[*n - 1].r = 1., work[*n - 1].i = 0.;

/*           Form right-hand side. */

            for (k = ki+1; k < *n; ++k) {
                work[k].r = -t[ki + k * *ldt].r, work[k].i = t[ki + k * *ldt].i;
            }

/*           Solve the triangular system: */
/*              (T(KI+1:N,KI+1:N) - T(KI,KI))'*X = SCALE*WORK. */

            for (k = ki+1; k < *n; ++k) {
                i__1 = k + k * *ldt; /* index [k,k] */
                i__2 = ki + ki * *ldt; /* index [ki,ki] */
                t[i__1].r -= t[i__2].r,
                t[i__1].i -= t[i__2].i;
                if (abs(t[i__1].r) + abs(t[i__1].i) < smin) {
                    t[i__1].r = smin, t[i__1].i = 0.;
                }
            }

            k = ki + 1;
            if (k < *n) {
                i__1 = *n - k;
                zlatrs_("Upper", "Conjugate transpose", "Non-unit", "Y",
                        &i__1, &t[k + k * *ldt], ldt, &work[k], &scale, rwork, info);
                work[ki].r = scale, work[ki].i = 0.;
            }

/*           Copy the vector x or Q*x to VL and normalize. */

            if (! over) {
                i__1 = *n - ki;
                zcopy_(&i__1, &work[ki], &c__1, &vl[ki + is * *ldvl], &c__1);
                ii = izamax_(&i__1, &vl[ki + is * *ldvl], &c__1) + ki;
                i__2 = ii-1 + is * *ldvl; /* index [ii-1,is] */
                remax = 1. / (abs(vl[i__2].r) + abs(vl[i__2].i));
                zdscal_(&i__1, &remax, &vl[ki + is * *ldvl], &c__1);

                for (k = 0; k < ki; ++k) {
                    i__1 = k + is * *ldvl; /* index [k,is] */
                    vl[i__1].r = 0., vl[i__1].i = 0.;
                }
            } else {
                k = ki + 1;
                if (k < *n) {
                    i__1 = *n - k;
                    z__1.r = scale, z__1.i = 0.;
                    zgemv_("N", n, &i__1, &c_b17, &vl[k * *ldvl],
                           ldvl, &work[k], &c__1, &z__1, &vl[ki * *ldvl], &c__1);
                }

                ii = izamax_(n, &vl[ki * *ldvl], &c__1);
                i__1 = ii-1 + ki * *ldvl; /* index [ii-1,ki] */
                remax = 1. / (abs(vl[i__1].r) + abs(vl[i__1].i));
                zdscal_(n, &remax, &vl[ki * *ldvl], &c__1);
            }

/*           Set back the original diagonal elements of T. */

            for (k = ki+1; k < *n; ++k) {
                i__1 = k + k * *ldt; /* index [k,k] */
                i__2 = k + *n;
                t[i__1].r = work[i__2].r, t[i__1].i = work[i__2].i;
            }
            ++is;
        }
    }
} /* ztrevc_ */