zhseqr.c

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

        if (! wantt) {
/*<             I1 = L >*/
            i1 = l;
/*<             I2 = I >*/
            i2 = i__;
/*<          END IF >*/
        }

/*<          IF( ITS.EQ.20 .OR. ITS.EQ.30 ) THEN >*/
        if (its == 20 || its == 30) {

/*           Exceptional shifts. */

/*<             DO 90 II = I - NS + 1, I >*/
            i__2 = i__;
            for (ii = i__ - ns + 1; ii <= i__2; ++ii) {
/*<    >*/
                i__3 = ii;
                i__5 = ii + (ii - 1) * h_dim1;
                i__6 = ii + ii * h_dim1;
                d__3 = ((d__1 = h__[i__5].r, abs(d__1)) + (d__2 = h__[i__6].r,
                         abs(d__2))) * 1.5;
                w[i__3].r = d__3, w[i__3].i = 0.;
/*<    90       CONTINUE >*/
/* L90: */
            }
/*<          ELSE >*/
        } else {

/*           Use eigenvalues of trailing submatrix of order NS as shifts. */

/*<    >*/
            zlacpy_("Full", &ns, &ns, &h__[i__ - ns + 1 + (i__ - ns + 1) * 
                    h_dim1], ldh, s, &c__15, (ftnlen)4);
/*<    >*/
            zlahqr_(&c_false, &c_false, &ns, &c__1, &ns, s, &c__15, &w[i__ - 
                    ns + 1], &c__1, &ns, &z__[z_offset], ldz, &ierr);
/*<             IF( IERR.GT.0 ) THEN >*/
            if (ierr > 0) {

/*              If ZLAHQR failed to compute all NS eigenvalues, use the */
/*              unconverged diagonal elements as the remaining shifts. */

/*<                DO 100 II = 1, IERR >*/
                i__2 = ierr;
                for (ii = 1; ii <= i__2; ++ii) {
/*<                   W( I-NS+II ) = S( II, II ) >*/
                    i__3 = i__ - ns + ii;
                    i__5 = ii + ii * 15 - 16;
                    w[i__3].r = s[i__5].r, w[i__3].i = s[i__5].i;
/*<   100          CONTINUE >*/
/* L100: */
                }
/*<             END IF >*/
            }
/*<          END IF >*/
        }

/*        Form the first column of (G-w(1)) (G-w(2)) . . . (G-w(ns)) */
/*        where G is the Hessenberg submatrix H(L:I,L:I) and w is */
/*        the vector of shifts (stored in W). The result is */
/*        stored in the local array V. */

/*<          V( 1 ) = ONE >*/
        v[0].r = 1., v[0].i = 0.;
/*<          DO 110 II = 2, NS + 1 >*/
        i__2 = ns + 1;
        for (ii = 2; ii <= i__2; ++ii) {
/*<             V( II ) = ZERO >*/
            i__3 = ii - 1;
            v[i__3].r = 0., v[i__3].i = 0.;
/*<   110    CONTINUE >*/
/* L110: */
        }
/*<          NV = 1 >*/
        nv = 1;
/*<          DO 130 J = I - NS + 1, I >*/
        i__2 = i__;
        for (j = i__ - ns + 1; j <= i__2; ++j) {
/*<             CALL ZCOPY( NV+1, V, 1, VV, 1 ) >*/
            i__3 = nv + 1;
            zcopy_(&i__3, v, &c__1, vv, &c__1);
/*<    >*/
            i__3 = nv + 1;
            i__5 = j;
            z__1.r = -w[i__5].r, z__1.i = -w[i__5].i;
            zgemv_("No transpose", &i__3, &nv, &c_b2, &h__[l + l * h_dim1], 
                    ldh, vv, &c__1, &z__1, v, &c__1, (ftnlen)12);
/*<             NV = NV + 1 >*/
            ++nv;

/*           Scale V(1:NV) so that max(abs(V(i))) = 1. If V is zero, */
/*           reset it to the unit vector. */

/*<             ITEMP = IZAMAX( NV, V, 1 ) >*/
            itemp = izamax_(&nv, v, &c__1);
/*<             RTEMP = CABS1( V( ITEMP ) ) >*/
            i__3 = itemp - 1;
            rtemp = (d__1 = v[i__3].r, abs(d__1)) + (d__2 = d_imag(&v[itemp - 
                    1]), abs(d__2));
/*<             IF( RTEMP.EQ.RZERO ) THEN >*/
            if (rtemp == 0.) {
/*<                V( 1 ) = ONE >*/
                v[0].r = 1., v[0].i = 0.;
/*<                DO 120 II = 2, NV >*/
                i__3 = nv;
                for (ii = 2; ii <= i__3; ++ii) {
/*<                   V( II ) = ZERO >*/
                    i__5 = ii - 1;
                    v[i__5].r = 0., v[i__5].i = 0.;
/*<   120          CONTINUE >*/
/* L120: */
                }
/*<             ELSE >*/
            } else {
/*<                RTEMP = MAX( RTEMP, SMLNUM ) >*/
                rtemp = max(rtemp,smlnum);
/*<                CALL ZDSCAL( NV, RONE / RTEMP, V, 1 ) >*/
                d__1 = 1. / rtemp;
                zdscal_(&nv, &d__1, v, &c__1);
/*<             END IF >*/
            }
/*<   130    CONTINUE >*/
/* L130: */
        }

/*        Multiple-shift QR step */

/*<          DO 150 K = L, I - 1 >*/
        i__2 = i__ - 1;
        for (k = l; k <= i__2; ++k) {

/*           The first iteration of this loop determines a reflection G */
/*           from the vector V and applies it from left and right to H, */
/*           thus creating a nonzero bulge below the subdiagonal. */

/*           Each subsequent iteration determines a reflection G to */
/*           restore the Hessenberg form in the (K-1)th column, and thus */
/*           chases the bulge one step toward the bottom of the active */
/*           submatrix. NR is the order of G. */

/*<             NR = MIN( NS+1, I-K+1 ) >*/
/* Computing MIN */
            i__3 = ns + 1, i__5 = i__ - k + 1;
            nr = min(i__3,i__5);
/*<    >*/
            if (k > l) {
                zcopy_(&nr, &h__[k + (k - 1) * h_dim1], &c__1, v, &c__1);
            }
/*<             CALL ZLARFG( NR, V( 1 ), V( 2 ), 1, TAU ) >*/
            zlarfg_(&nr, v, &v[1], &c__1, &tau);
/*<             IF( K.GT.L ) THEN >*/
            if (k > l) {
/*<                H( K, K-1 ) = V( 1 ) >*/
                i__3 = k + (k - 1) * h_dim1;
                h__[i__3].r = v[0].r, h__[i__3].i = v[0].i;
/*<                DO 140 II = K + 1, I >*/
                i__3 = i__;
                for (ii = k + 1; ii <= i__3; ++ii) {
/*<                   H( II, K-1 ) = ZERO >*/
                    i__5 = ii + (k - 1) * h_dim1;
                    h__[i__5].r = 0., h__[i__5].i = 0.;
/*<   140          CONTINUE >*/
/* L140: */
                }
/*<             END IF >*/
            }
/*<             V( 1 ) = ONE >*/
            v[0].r = 1., v[0].i = 0.;

/*           Apply G' from the left to transform the rows of the matrix */
/*           in columns K to I2. */

/*<    >*/
            i__3 = i2 - k + 1;
            d_cnjg(&z__1, &tau);
            zlarfx_("Left", &nr, &i__3, v, &z__1, &h__[k + k * h_dim1], ldh, &
                    work[1], (ftnlen)4);

/*           Apply G from the right to transform the columns of the */
/*           matrix in rows I1 to min(K+NR,I). */

/*<    >*/
/* Computing MIN */
            i__5 = k + nr;
            i__3 = min(i__5,i__) - i1 + 1;
            zlarfx_("Right", &i__3, &nr, v, &tau, &h__[i1 + k * h_dim1], ldh, 
                    &work[1], (ftnlen)5);

/*<             IF( WANTZ ) THEN >*/
            if (wantz) {

/*              Accumulate transformations in the matrix Z */

/*<    >*/
                zlarfx_("Right", &nh, &nr, v, &tau, &z__[*ilo + k * z_dim1], 
                        ldz, &work[1], (ftnlen)5);
/*<             END IF >*/
            }
/*<   150    CONTINUE >*/
/* L150: */
        }

/*        Ensure that H(I,I-1) is real. */

/*<          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( WANTZ ) THEN >*/
            if (wantz) {
/*<                CALL ZSCAL( NH, TEMP, Z( ILO, I ), 1 ) >*/
                zscal_(&nh, &temp, &z__[*ilo + i__ * z_dim1], &c__1);
/*<             END IF >*/
            }
/*<          END IF >*/
        }

/*<   160 CONTINUE >*/
/* L160: */
    }

/*     Failure to converge in remaining number of iterations */

/*<       INFO = I >*/
    *info = i__;
/*<       RETURN >*/
    return 0;

/*<   170 CONTINUE >*/
L170:

/*     A submatrix of order <= MAXB in rows and columns L to I has split */
/*     off. Use the double-shift QR algorithm to handle it. */

/*<    >*/
    zlahqr_(&wantt, &wantz, n, &l, &i__, &h__[h_offset], ldh, &w[1], ilo, ihi,
             &z__[z_offset], ldz, info);
/*<    >*/
    if (*info > 0) {
        return 0;
    }

/*     Decrement number of remaining iterations, and return to start of */
/*     the main loop with a new value of I. */

/*<       ITN = ITN - ITS >*/
    itn -= its;
/*<       I = L - 1 >*/
    i__ = l - 1;
/*<       GO TO 60 >*/
    goto L60;

/*<   180 CONTINUE >*/
L180:
/*<       WORK( 1 ) = MAX( 1, N ) >*/
    i__1 = max(1,*n);
    work[1].r = (doublereal) i__1, work[1].i = 0.;
/*<       RETURN >*/
    return 0;

/*     End of ZHSEQR */

/*<       END >*/
} /* zhseqr_ */

#ifdef __cplusplus
        }
#endif