zgeev.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

    }
/*<       IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN >*/
    if (*lwork < minwrk && ! lquery) {
/*<          INFO = -12 >*/
        *info = -12;
/*<       END IF >*/
    }
/*<       IF( INFO.NE.0 ) THEN >*/
    if (*info != 0) {
/*<          CALL XERBLA( 'ZGEEV ', -INFO ) >*/
        i__1 = -(*info);
        xerbla_("ZGEEV ", &i__1, (ftnlen)6);
/*<          RETURN >*/
        return 0;
/*<       ELSE IF( LQUERY ) THEN >*/
    } else if (lquery) {
/*<          RETURN >*/
        return 0;
/*<       END IF >*/
    }

/*     Quick return if possible */

/*<    >*/
    if (*n == 0) {
        return 0;
    }

/*     Get machine constants */

/*<       EPS = DLAMCH( 'P' ) >*/
    eps = dlamch_("P", (ftnlen)1);
/*<       SMLNUM = DLAMCH( 'S' ) >*/
    smlnum = dlamch_("S", (ftnlen)1);
/*<       BIGNUM = ONE / SMLNUM >*/
    bignum = 1. / smlnum;
/*<       CALL DLABAD( SMLNUM, BIGNUM ) >*/
    dlabad_(&smlnum, &bignum);
/*<       SMLNUM = SQRT( SMLNUM ) / EPS >*/
    smlnum = sqrt(smlnum) / eps;
/*<       BIGNUM = ONE / SMLNUM >*/
    bignum = 1. / smlnum;

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

/*<       ANRM = ZLANGE( 'M', N, N, A, LDA, DUM ) >*/
    anrm = zlange_("M", n, n, &a[a_offset], lda, dum, (ftnlen)1);
/*<       SCALEA = .FALSE. >*/
    scalea = FALSE_;
/*<       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN >*/
    if (anrm > 0. && anrm < smlnum) {
/*<          SCALEA = .TRUE. >*/
        scalea = TRUE_;
/*<          CSCALE = SMLNUM >*/
        cscale = smlnum;
/*<       ELSE IF( ANRM.GT.BIGNUM ) THEN >*/
    } else if (anrm > bignum) {
/*<          SCALEA = .TRUE. >*/
        scalea = TRUE_;
/*<          CSCALE = BIGNUM >*/
        cscale = bignum;
/*<       END IF >*/
    }
/*<    >*/
    if (scalea) {
        zlascl_("G", &c__0, &c__0, &anrm, &cscale, n, n, &a[a_offset], lda, &
                ierr, (ftnlen)1);
    }

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

/*<       IBAL = 1 >*/
    ibal = 1;
/*<       CALL ZGEBAL( 'B', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR ) >*/
    zgebal_("B", n, &a[a_offset], lda, &ilo, &ihi, &rwork[ibal], &ierr, (
            ftnlen)1);

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

/*<       ITAU = 1 >*/
    itau = 1;
/*<       IWRK = ITAU + N >*/
    iwrk = itau + *n;
/*<    >*/
    i__1 = *lwork - iwrk + 1;
    zgehrd_(n, &ilo, &ihi, &a[a_offset], lda, &work[itau], &work[iwrk], &i__1,
             &ierr);

/*<       IF( WANTVL ) THEN >*/
    if (wantvl) {

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

/*<          SIDE = 'L' >*/
        *(unsigned char *)side = 'L';
/*<          CALL ZLACPY( 'L', N, N, A, LDA, VL, LDVL ) >*/
        zlacpy_("L", n, n, &a[a_offset], lda, &vl[vl_offset], ldvl, (ftnlen)1)
                ;

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

/*<    >*/
        i__1 = *lwork - iwrk + 1;
        zunghr_(n, &ilo, &ihi, &vl[vl_offset], 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 >*/
        iwrk = itau;
/*<    >*/
        i__1 = *lwork - iwrk + 1;
        zhseqr_("S", "V", n, &ilo, &ihi, &a[a_offset], lda, &w[1], &vl[
                vl_offset], ldvl, &work[iwrk], &i__1, info, (ftnlen)1, (
                ftnlen)1);

/*<          IF( WANTVR ) THEN >*/
        if (wantvr) {

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

/*<             SIDE = 'B' >*/
            *(unsigned char *)side = 'B';
/*<             CALL ZLACPY( 'F', N, N, VL, LDVL, VR, LDVR ) >*/
            zlacpy_("F", n, n, &vl[vl_offset], ldvl, &vr[vr_offset], ldvr, (
                    ftnlen)1);
/*<          END IF >*/
        }

/*<       ELSE IF( WANTVR ) THEN >*/
    } else if (wantvr) {

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

/*<          SIDE = 'R' >*/
        *(unsigned char *)side = 'R';
/*<          CALL ZLACPY( 'L', N, N, A, LDA, VR, LDVR ) >*/
        zlacpy_("L", n, n, &a[a_offset], lda, &vr[vr_offset], ldvr, (ftnlen)1)
                ;

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

/*<    >*/
        i__1 = *lwork - iwrk + 1;
        zunghr_(n, &ilo, &ihi, &vr[vr_offset], 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 >*/
        iwrk = itau;
/*<    >*/
        i__1 = *lwork - iwrk + 1;
        zhseqr_("S", "V", n, &ilo, &ihi, &a[a_offset], lda, &w[1], &vr[
                vr_offset], ldvr, &work[iwrk], &i__1, info, (ftnlen)1, (
                ftnlen)1);

/*<       ELSE >*/
    } else {

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

/*<          IWRK = ITAU >*/
        iwrk = itau;
/*<    >*/
        i__1 = *lwork - iwrk + 1;
        zhseqr_("E", "N", n, &ilo, &ihi, &a[a_offset], lda, &w[1], &vr[
                vr_offset], ldvr, &work[iwrk], &i__1, info, (ftnlen)1, (
                ftnlen)1);
/*<       END IF >*/
    }

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

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

/*<       IF( WANTVL .OR. WANTVR ) THEN >*/
    if (wantvl || wantvr) {

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

/*<          IRWORK = IBAL + N >*/
        irwork = ibal + *n;
/*<    >*/
        ztrevc_(side, "B", select, n, &a[a_offset], lda, &vl[vl_offset], ldvl,
                 &vr[vr_offset], ldvr, n, &nout, &work[iwrk], &rwork[irwork], 
                &ierr, (ftnlen)1, (ftnlen)1);
/*<       END IF >*/
    }

/*<       IF( WANTVL ) THEN >*/
    if (wantvl) {

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

/*<    >*/
        zgebak_("B", "L", n, &ilo, &ihi, &rwork[ibal], n, &vl[vl_offset], 
                ldvl, &ierr, (ftnlen)1, (ftnlen)1);

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

/*<          DO 20 I = 1, N >*/
        i__1 = *n;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<             SCL = ONE / DZNRM2( N, VL( 1, I ), 1 ) >*/
            scl = 1. / dznrm2_(n, &vl[i__ * vl_dim1 + 1], &c__1);
/*<             CALL ZDSCAL( N, SCL, VL( 1, I ), 1 ) >*/
            zdscal_(n, &scl, &vl[i__ * vl_dim1 + 1], &c__1);
/*<             DO 10 K = 1, N >*/
            i__2 = *n;
            for (k = 1; k <= i__2; ++k) {
/*<    >*/
                i__3 = k + i__ * vl_dim1;
/* Computing 2nd power */
                d__1 = vl[i__3].r;
/* Computing 2nd power */
                d__2 = d_imag(&vl[k + i__ * vl_dim1]);
                rwork[irwork + k - 1] = d__1 * d__1 + d__2 * d__2;
/*<    10       CONTINUE >*/
/* L10: */
            }
/*<             K = IDAMAX( N, RWORK( IRWORK ), 1 ) >*/
            k = idamax_(n, &rwork[irwork], &c__1);
/*<             TMP = DCONJG( VL( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) ) >*/
            d_cnjg(&z__2, &vl[k + i__ * vl_dim1]);
            d__1 = sqrt(rwork[irwork + k - 1]);
            z__1.r = z__2.r / d__1, z__1.i = z__2.i / d__1;
            tmp.r = z__1.r, tmp.i = z__1.i;
/*<             CALL ZSCAL( N, TMP, VL( 1, I ), 1 ) >*/
            zscal_(n, &tmp, &vl[i__ * vl_dim1 + 1], &c__1);
/*<             VL( K, I ) = DCMPLX( DBLE( VL( K, I ) ), ZERO ) >*/
            i__2 = k + i__ * vl_dim1;
            i__3 = k + i__ * vl_dim1;
            d__1 = vl[i__3].r;
            z__1.r = d__1, z__1.i = 0.;
            vl[i__2].r = z__1.r, vl[i__2].i = z__1.i;
/*<    20    CONTINUE >*/
/* L20: */
        }
/*<       END IF >*/
    }

/*<       IF( WANTVR ) THEN >*/
    if (wantvr) {

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

/*<    >*/
        zgebak_("B", "R", n, &ilo, &ihi, &rwork[ibal], n, &vr[vr_offset], 
                ldvr, &ierr, (ftnlen)1, (ftnlen)1);

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

/*<          DO 40 I = 1, N >*/
        i__1 = *n;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<             SCL = ONE / DZNRM2( N, VR( 1, I ), 1 ) >*/
            scl = 1. / dznrm2_(n, &vr[i__ * vr_dim1 + 1], &c__1);
/*<             CALL ZDSCAL( N, SCL, VR( 1, I ), 1 ) >*/
            zdscal_(n, &scl, &vr[i__ * vr_dim1 + 1], &c__1);
/*<             DO 30 K = 1, N >*/
            i__2 = *n;
            for (k = 1; k <= i__2; ++k) {
/*<    >*/
                i__3 = k + i__ * vr_dim1;
/* Computing 2nd power */
                d__1 = vr[i__3].r;
/* Computing 2nd power */
                d__2 = d_imag(&vr[k + i__ * vr_dim1]);
                rwork[irwork + k - 1] = d__1 * d__1 + d__2 * d__2;
/*<    30       CONTINUE >*/
/* L30: */
            }
/*<             K = IDAMAX( N, RWORK( IRWORK ), 1 ) >*/
            k = idamax_(n, &rwork[irwork], &c__1);
/*<             TMP = DCONJG( VR( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) ) >*/
            d_cnjg(&z__2, &vr[k + i__ * vr_dim1]);
            d__1 = sqrt(rwork[irwork + k - 1]);
            z__1.r = z__2.r / d__1, z__1.i = z__2.i / d__1;
            tmp.r = z__1.r, tmp.i = z__1.i;
/*<             CALL ZSCAL( N, TMP, VR( 1, I ), 1 ) >*/
            zscal_(n, &tmp, &vr[i__ * vr_dim1 + 1], &c__1);
/*<             VR( K, I ) = DCMPLX( DBLE( VR( K, I ) ), ZERO ) >*/
            i__2 = k + i__ * vr_dim1;
            i__3 = k + i__ * vr_dim1;
            d__1 = vr[i__3].r;
            z__1.r = d__1, z__1.i = 0.;
            vr[i__2].r = z__1.r, vr[i__2].i = z__1.i;
/*<    40    CONTINUE >*/
/* L40: */
        }
/*<       END IF >*/
    }

/*     Undo scaling if necessary */

/*<    50 CONTINUE >*/
L50:
/*<       IF( SCALEA ) THEN >*/
    if (scalea) {
/*<    >*/
        i__1 = *n - *info;
/* Computing MAX */
        i__3 = *n - *info;
        i__2 = max(i__3,1);
        zlascl_("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &w[*info + 1]
                , &i__2, &ierr, (ftnlen)1);
/*<          IF( INFO.GT.0 ) THEN >*/
        if (*info > 0) {
/*<             CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR ) >*/
            i__1 = ilo - 1;
            zlascl_("G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &w[1], n,
                     &ierr, (ftnlen)1);
/*<          END IF >*/
        }
/*<       END IF >*/
    }

/*<       WORK( 1 ) = MAXWRK >*/
    work[1].r = (doublereal) maxwrk, work[1].i = 0.;
/*<       RETURN >*/
    return 0;

/*     End of ZGEEV */

/*<       END >*/
} /* zgeev_ */

#ifdef __cplusplus
        }
#endif