zgeev.c

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

    }

/*     Quick return if possible */

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

/*     Get machine constants */

    eps = dlamch_("P");
    smlnum = dlamch_("S");
    bignum = 1. / smlnum;
    dlabad_(&smlnum, &bignum);
    smlnum = sqrt(smlnum) / eps;
    bignum = 1. / smlnum;

/*     Scale A if max element outside range [SMLNUM,BIGNUM] */

    anrm = zlange_("M", n, n, a, lda, dum);
    scalea = FALSE_;
    if (anrm > 0. && anrm < smlnum) {
        scalea = TRUE_;
        cscale = smlnum;
    } else if (anrm > bignum) {
        scalea = TRUE_;
        cscale = bignum;
    }
    if (scalea) {
        zlascl_("G", &c__0, &c__0, &anrm, &cscale, n, n, a, lda, &ierr);
    }

/*     Balance the matrix */
/*     (CWorkspace: none) */
/*     (RWorkspace: need N) */

    ibal = 0;
    zgebal_("B", n, a, lda, &ilo, &ihi, &rwork[ibal], &ierr);

/*     Reduce to upper Hessenberg form */
/*     (CWorkspace: need 2*N, prefer N+N*NB) */
/*     (RWorkspace: none) */

    itau = 0;
    iwrk = itau + *n;
    i__1 = *lwork - iwrk;
    zgehrd_(n, &ilo, &ihi, a, lda, &work[itau], &work[iwrk], &i__1, &ierr);

    if (wantvl) {

/*        Want left eigenvectors */
/*        Copy Householder vectors to VL */

        *side = 'L';
        zlacpy_("L", n, n, a, lda, vl, ldvl);

/*        Generate unitary matrix in VL */
/*        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB) */
/*        (RWorkspace: none) */

        i__1 = *lwork - iwrk;
        zunghr_(n, &ilo, &ihi, vl, ldvl, &work[itau], &work[iwrk], &i__1, &ierr);

/*        Perform QR iteration, accumulating Schur vectors in VL */
/*        (CWorkspace: need 1, prefer HSWORK (see comments) ) */
/*        (RWorkspace: none) */

        iwrk = itau;
        i__1 = *lwork - iwrk;
        zhseqr_("S", "V", n, &ilo, &ihi, a, lda, w, vl, ldvl, &work[iwrk], &i__1, info);

        if (wantvr) {

/*           Want left and right eigenvectors */
/*           Copy Schur vectors to VR */

            *side = 'B';
            zlacpy_("F", n, n, vl, ldvl, vr, ldvr);
        }

    } else if (wantvr) {

/*        Want right eigenvectors */
/*        Copy Householder vectors to VR */

        *side = 'R';
        zlacpy_("L", n, n, a, lda, vr, ldvr);

/*        Generate unitary matrix in VR */
/*        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB) */
/*        (RWorkspace: none) */

        i__1 = *lwork - iwrk;
        zunghr_(n, &ilo, &ihi, vr, ldvr, &work[itau], &work[iwrk], &i__1, &ierr);

/*        Perform QR iteration, accumulating Schur vectors in VR */
/*        (CWorkspace: need 1, prefer HSWORK (see comments) ) */
/*        (RWorkspace: none) */

        iwrk = itau;
        i__1 = *lwork - iwrk;
        zhseqr_("S", "V", n, &ilo, &ihi, a, lda, w, vr, ldvr, &work[iwrk], &i__1, info);

    } else {

/*        Compute eigenvalues only */
/*        (CWorkspace: need 1, prefer HSWORK (see comments) ) */
/*        (RWorkspace: none) */

        iwrk = itau;
        i__1 = *lwork - iwrk;
        zhseqr_("E", "N", n, &ilo, &ihi, a, lda, w, vr, ldvr, &work[iwrk], &i__1, info);
    }

/*     If INFO > 0 from ZHSEQR, then quit */

    if (*info > 0) {
        goto L50;
    }

    if (wantvl || wantvr) {

/*        Compute left and/or right eigenvectors */
/*        (CWorkspace: need 2*N) */
/*        (RWorkspace: need 2*N) */

        irwork = ibal + *n;
        ztrevc_(side, "B", select, n, a, lda, vl, ldvl, vr, ldvr, n, &nout, &work[iwrk], &rwork[irwork], &ierr);
    }

    if (wantvl) {

/*        Undo balancing of left eigenvectors */
/*        (CWorkspace: none) */
/*        (RWorkspace: need N) */

        zgebak_("B", "L", n, &ilo, &ihi, &rwork[ibal], n, vl, ldvl, &ierr);

/*        Normalize left eigenvectors and make largest component real */

        for (i = 0; i < *n; ++i) {
            scl = 1. / dznrm2_(n, &vl[i * *ldvl], &c__1);
            zdscal_(n, &scl, &vl[i * *ldvl], &c__1);
            for (k = 0; k < *n; ++k) {
                i__1 = k + i * *ldvl; /* index [k,i] */
                rwork[irwork + k] = vl[i__1].r * vl[i__1].r + vl[i__1].i * vl[i__1].i;
            }
            k = idamax_(n, &rwork[irwork], &c__1) - 1;
            d_cnjg(&tmp, &vl[k + i * *ldvl]);
            d__1 = sqrt(rwork[irwork + k]);
            tmp.r /= d__1, tmp.i /= d__1;
            zscal_(n, &tmp, &vl[i * *ldvl], &c__1);
            vl[k + i * *ldvl].i = 0.;
        }
    }

    if (wantvr) {

/*        Undo balancing of right eigenvectors */
/*        (CWorkspace: none) */
/*        (RWorkspace: need N) */

        zgebak_("B", "R", n, &ilo, &ihi, &rwork[ibal], n, vr, ldvr, &ierr);

/*        Normalize right eigenvectors and make largest component real */

        for (i = 0; i < *n; ++i) {
            scl = 1. / dznrm2_(n, &vr[i * *ldvr], &c__1);
            zdscal_(n, &scl, &vr[i * *ldvr], &c__1);
            for (k = 0; k < *n; ++k) {
                i__1 = k + i * *ldvr; /* index [k,i] */
                rwork[irwork + k] = vr[i__1].r * vr[i__1].r + vr[i__1].i * vr[i__1].i;
            }
            k = idamax_(n, &rwork[irwork], &c__1) - 1;
            d_cnjg(&tmp, &vr[k + i * *ldvr]);
            d__1 = sqrt(rwork[irwork + k]);
            tmp.r /= d__1, tmp.i /= d__1;
            zscal_(n, &tmp, &vr[i * *ldvr], &c__1);
            vr[k + i * *ldvr].i = 0.;
        }
    }

/*     Undo scaling if necessary */

L50:
    if (scalea) {
        i__1 = *n - *info;
        i__2 = max(i__1, 1);
        zlascl_("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &w[*info], &i__2, &ierr);
        if (*info > 0) {
            i__1 = ilo - 1;
            zlascl_("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, w, n, &ierr);
        }
    }

    work[0].r = (doublereal) maxwrk, work[0].i = 0.;

} /* zgeev_ */