zhseqr.c

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

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

/*     Initialize Z, if necessary */

/*<    >*/
    if (initz) {
        zlaset_("Full", n, n, &c_b1, &c_b2, &z__[z_offset], ldz, (ftnlen)4);
    }

/*     Store the eigenvalues isolated by ZGEBAL. */

/*<       DO 10 I = 1, ILO - 1 >*/
    i__1 = *ilo - 1;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          W( I ) = H( I, I ) >*/
        i__2 = i__;
        i__3 = i__ + i__ * h_dim1;
        w[i__2].r = h__[i__3].r, w[i__2].i = h__[i__3].i;
/*<    10 CONTINUE >*/
/* L10: */
    }
/*<       DO 20 I = IHI + 1, N >*/
    i__1 = *n;
    for (i__ = *ihi + 1; i__ <= i__1; ++i__) {
/*<          W( I ) = H( I, I ) >*/
        i__2 = i__;
        i__3 = i__ + i__ * h_dim1;
        w[i__2].r = h__[i__3].r, w[i__2].i = h__[i__3].i;
/*<    20 CONTINUE >*/
/* L20: */
    }

/*     Quick return if possible. */

/*<    >*/
    if (*n == 0) {
        return 0;
    }
/*<       IF( ILO.EQ.IHI ) THEN >*/
    if (*ilo == *ihi) {
/*<          W( ILO ) = H( ILO, ILO ) >*/
        i__1 = *ilo;
        i__2 = *ilo + *ilo * h_dim1;
        w[i__1].r = h__[i__2].r, w[i__1].i = h__[i__2].i;
/*<          RETURN >*/
        return 0;
/*<       END IF >*/
    }

/*     Set rows and columns ILO to IHI to zero below the first */
/*     subdiagonal. */

/*<       DO 40 J = ILO, IHI - 2 >*/
    i__1 = *ihi - 2;
    for (j = *ilo; j <= i__1; ++j) {
/*<          DO 30 I = J + 2, N >*/
        i__2 = *n;
        for (i__ = j + 2; i__ <= i__2; ++i__) {
/*<             H( I, J ) = ZERO >*/
            i__3 = i__ + j * h_dim1;
            h__[i__3].r = 0., h__[i__3].i = 0.;
/*<    30    CONTINUE >*/
/* L30: */
        }
/*<    40 CONTINUE >*/
/* L40: */
    }
/*<       NH = IHI - ILO + 1 >*/
    nh = *ihi - *ilo + 1;

/*     I1 and I2 are the indices of the first row and last column of H */
/*     to which transformations must be applied. If eigenvalues only are */
/*     being computed, I1 and I2 are re-set inside the main loop. */

/*<       IF( WANTT ) THEN >*/
    if (wantt) {
/*<          I1 = 1 >*/
        i1 = 1;
/*<          I2 = N >*/
        i2 = *n;
/*<       ELSE >*/
    } else {
/*<          I1 = ILO >*/
        i1 = *ilo;
/*<          I2 = IHI >*/
        i2 = *ihi;
/*<       END IF >*/
    }

/*     Ensure that the subdiagonal elements are real. */

/*<       DO 50 I = ILO + 1, IHI >*/
    i__1 = *ihi;
    for (i__ = *ilo + 1; i__ <= i__1; ++i__) {
/*<          TEMP = H( I, I-1 ) >*/
        i__2 = i__ + (i__ - 1) * h_dim1;
        temp.r = h__[i__2].r, temp.i = h__[i__2].i;
/*<          IF( DIMAG( TEMP ).NE.RZERO ) THEN >*/
        if (d_imag(&temp) != 0.) {
/*<             RTEMP = DLAPY2( DBLE( TEMP ), DIMAG( TEMP ) ) >*/
            d__1 = temp.r;
            d__2 = d_imag(&temp);
            rtemp = dlapy2_(&d__1, &d__2);
/*<             H( I, I-1 ) = RTEMP >*/
            i__2 = i__ + (i__ - 1) * h_dim1;
            h__[i__2].r = rtemp, h__[i__2].i = 0.;
/*<             TEMP = TEMP / RTEMP >*/
            z__1.r = temp.r / rtemp, z__1.i = temp.i / rtemp;
            temp.r = z__1.r, temp.i = z__1.i;
/*<    >*/
            if (i2 > i__) {
                i__2 = i2 - i__;
                d_cnjg(&z__1, &temp);
                zscal_(&i__2, &z__1, &h__[i__ + (i__ + 1) * h_dim1], ldh);
            }
/*<             CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 ) >*/
            i__2 = i__ - i1;
            zscal_(&i__2, &temp, &h__[i1 + i__ * h_dim1], &c__1);
/*<    >*/
            if (i__ < *ihi) {
                i__2 = i__ + 1 + i__ * h_dim1;
                i__3 = i__ + 1 + i__ * h_dim1;
                z__1.r = temp.r * h__[i__3].r - temp.i * h__[i__3].i, z__1.i =
                         temp.r * h__[i__3].i + temp.i * h__[i__3].r;
                h__[i__2].r = z__1.r, h__[i__2].i = z__1.i;
            }
/*<    >*/
            if (wantz) {
                zscal_(&nh, &temp, &z__[*ilo + i__ * z_dim1], &c__1);
            }
/*<          END IF >*/
        }
/*<    50 CONTINUE >*/
/* L50: */
    }

/*     Determine the order of the multi-shift QR algorithm to be used. */

/*<       NS = ILAENV( 4, 'ZHSEQR', JOB // COMPZ, N, ILO, IHI, -1 ) >*/
/* Writing concatenation */
    i__4[0] = 1, a__1[0] = job;
    i__4[1] = 1, a__1[1] = compz;
    s_cat(ch__1, a__1, i__4, &c__2, (ftnlen)2);
    ns = ilaenv_(&c__4, "ZHSEQR", ch__1, n, ilo, ihi, &c_n1, (ftnlen)6, (
            ftnlen)2);
/*<       MAXB = ILAENV( 8, 'ZHSEQR', JOB // COMPZ, N, ILO, IHI, -1 ) >*/
/* Writing concatenation */
    i__4[0] = 1, a__1[0] = job;
    i__4[1] = 1, a__1[1] = compz;
    s_cat(ch__1, a__1, i__4, &c__2, (ftnlen)2);
    maxb = ilaenv_(&c__8, "ZHSEQR", ch__1, n, ilo, ihi, &c_n1, (ftnlen)6, (
            ftnlen)2);
/*<       IF( NS.LE.1 .OR. NS.GT.NH .OR. MAXB.GE.NH ) THEN >*/
    if (ns <= 1 || ns > nh || maxb >= nh) {

/*        Use the standard double-shift algorithm */

/*<    >*/
        zlahqr_(&wantt, &wantz, n, ilo, ihi, &h__[h_offset], ldh, &w[1], ilo, 
                ihi, &z__[z_offset], ldz, info);
/*<          RETURN >*/
        return 0;
/*<       END IF >*/
    }
/*<       MAXB = MAX( 2, MAXB ) >*/
    maxb = max(2,maxb);
/*<       NS = MIN( NS, MAXB, NSMAX ) >*/
/* Computing MIN */
    i__1 = min(ns,maxb);
    ns = min(i__1,15);

/*     Now 1 < NS <= MAXB < NH. */

/*     Set machine-dependent constants for the stopping criterion. */
/*     If norm(H) <= sqrt(OVFL), overflow should not occur. */

/*<       UNFL = DLAMCH( 'Safe minimum' ) >*/
    unfl = dlamch_("Safe minimum", (ftnlen)12);
/*<       OVFL = RONE / UNFL >*/
    ovfl = 1. / unfl;
/*<       CALL DLABAD( UNFL, OVFL ) >*/
    dlabad_(&unfl, &ovfl);
/*<       ULP = DLAMCH( 'Precision' ) >*/
    ulp = dlamch_("Precision", (ftnlen)9);
/*<       SMLNUM = UNFL*( NH / ULP ) >*/
    smlnum = unfl * (nh / ulp);

/*     ITN is the total number of multiple-shift QR iterations allowed. */

/*<       ITN = 30*NH >*/
    itn = nh * 30;

/*     The main loop begins here. I is the loop index and decreases from */
/*     IHI to ILO in steps of at most MAXB. Each iteration of the loop */
/*     works with the active submatrix in rows and columns L to I. */
/*     Eigenvalues I+1 to IHI have already converged. Either L = ILO, or */
/*     H(L,L-1) is negligible so that the matrix splits. */

/*<       I = IHI >*/
    i__ = *ihi;
/*<    60 CONTINUE >*/
L60:
/*<    >*/
    if (i__ < *ilo) {
        goto L180;
    }

/*     Perform multiple-shift QR iterations on rows and columns ILO to I */
/*     until a submatrix of order at most MAXB splits off at the bottom */
/*     because a subdiagonal element has become negligible. */

/*<       L = ILO >*/
    l = *ilo;
/*<       DO 160 ITS = 0, ITN >*/
    i__1 = itn;
    for (its = 0; its <= i__1; ++its) {

/*        Look for a single small subdiagonal element. */

/*<          DO 70 K = I, L + 1, -1 >*/
        i__2 = l + 1;
        for (k = i__; k >= i__2; --k) {
/*<             TST1 = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) ) >*/
            i__3 = k - 1 + (k - 1) * h_dim1;
            i__5 = k + k * h_dim1;
            tst1 = (d__1 = h__[i__3].r, abs(d__1)) + (d__2 = d_imag(&h__[k - 
                    1 + (k - 1) * h_dim1]), abs(d__2)) + ((d__3 = h__[i__5].r,
                     abs(d__3)) + (d__4 = d_imag(&h__[k + k * h_dim1]), abs(
                    d__4)));
/*<    >*/
            if (tst1 == 0.) {
                i__3 = i__ - l + 1;
                tst1 = zlanhs_("1", &i__3, &h__[l + l * h_dim1], ldh, rwork, (
                        ftnlen)1);
            }
/*<    >*/
            i__3 = k + (k - 1) * h_dim1;
/* Computing MAX */
            d__2 = ulp * tst1;
            if ((d__1 = h__[i__3].r, abs(d__1)) <= max(d__2,smlnum)) {
                goto L80;
            }
/*<    70    CONTINUE >*/
/* L70: */
        }
/*<    80    CONTINUE >*/
L80:
/*<          L = K >*/
        l = k;
/*<          IF( L.GT.ILO ) THEN >*/
        if (l > *ilo) {

/*           H(L,L-1) is negligible. */

/*<             H( L, L-1 ) = ZERO >*/
            i__2 = l + (l - 1) * h_dim1;
            h__[i__2].r = 0., h__[i__2].i = 0.;
/*<          END IF >*/
        }

/*        Exit from loop if a submatrix of order <= MAXB has split off. */

/*<    >*/
        if (l >= i__ - maxb + 1) {
            goto L170;
        }

/*        Now the active submatrix is in rows and columns L to I. If */
/*        eigenvalues only are being computed, only the active submatrix */
/*        need be transformed. */

/*<          IF( .NOT.WANTT ) THEN >*/