zlahqr.c

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

/* lapack/complex16/zlahqr.f -- translated by f2c (version 20050501).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
        or, if you install libf2c.a in a standard place, with -lf2c -lm
        -- in that order, at the end of the command line, as in
                cc *.o -lf2c -lm
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

                http://www.netlib.org/f2c/libf2c.zip
*/

#ifdef __cplusplus
extern "C" {
#endif
#include "v3p_netlib.h"

/* Table of constant values */

static integer c__2 = 2;
static integer c__1 = 1;

/*<    >*/
/* Subroutine */ int zlahqr_(logical *wantt, logical *wantz, integer *n, 
        integer *ilo, integer *ihi, doublecomplex *h__, integer *ldh, 
        doublecomplex *w, integer *iloz, integer *ihiz, doublecomplex *z__, 
        integer *ldz, integer *info)
{
    /* System generated locals */
    integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5;
    doublereal d__1, d__2, d__3, d__4, d__5, d__6;
    doublecomplex z__1, z__2, z__3, z__4;

    /* Builtin functions */
    double d_imag(doublecomplex *);
    void z_sqrt(doublecomplex *, doublecomplex *), d_cnjg(doublecomplex *, 
            doublecomplex *);
    double z_abs(doublecomplex *);

    /* Local variables */
    integer i__, j, k, l, m;
    doublereal s;
    doublecomplex t, u, v[2], x, y;
    integer i1=0, i2=0;
    doublecomplex t1;
    doublereal t2;
    doublecomplex v2;
    doublereal h10;
    doublecomplex h11;
    doublereal h21;
    doublecomplex h22;
    integer nh, nz;
    doublecomplex h11s;
    integer itn, its;
    doublereal ulp;
    doublecomplex sum;
    doublereal tst1;
    doublecomplex temp;
    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 
            doublecomplex *, integer *);
    doublereal rtemp, rwork[1];
    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 
            doublecomplex *, integer *);
    extern doublereal dlamch_(char *, ftnlen);
    extern /* Subroutine */ int zlarfg_(integer *, doublecomplex *, 
            doublecomplex *, integer *, doublecomplex *);
    extern /* Double Complex */ VOID zladiv_(doublecomplex *, doublecomplex *,
             doublecomplex *);
    extern doublereal zlanhs_(char *, integer *, doublecomplex *, integer *, 
            doublereal *, ftnlen);
    doublereal smlnum;


/*  -- LAPACK auxiliary routine (version 3.0) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/*     Courant Institute, Argonne National Lab, and Rice University */
/*     June 30, 1999 */

/*     .. Scalar Arguments .. */
/*<       LOGICAL            WANTT, WANTZ >*/
/*<       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N >*/
/*     .. */
/*     .. Array Arguments .. */
/*<       COMPLEX*16         H( LDH, * ), W( * ), Z( LDZ, * ) >*/
/*     .. */

/*  Purpose */
/*  ======= */

/*  ZLAHQR is an auxiliary routine called by ZHSEQR to update the */
/*  eigenvalues and Schur decomposition already computed by ZHSEQR, by */
/*  dealing with the Hessenberg submatrix in rows and columns ILO to IHI. */

/*  Arguments */
/*  ========= */

/*  WANTT   (input) LOGICAL */
/*          = .TRUE. : the full Schur form T is required; */
/*          = .FALSE.: only eigenvalues are required. */

/*  WANTZ   (input) LOGICAL */
/*          = .TRUE. : the matrix of Schur vectors Z is required; */
/*          = .FALSE.: Schur vectors are not required. */

/*  N       (input) INTEGER */
/*          The order of the matrix H.  N >= 0. */

/*  ILO     (input) INTEGER */
/*  IHI     (input) INTEGER */
/*          It is assumed that H is already upper triangular in rows and */
/*          columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1). */
/*          ZLAHQR works primarily with the Hessenberg submatrix in rows */
/*          and columns ILO to IHI, but applies transformations to all of */
/*          H if WANTT is .TRUE.. */
/*          1 <= ILO <= max(1,IHI); IHI <= N. */

/*  H       (input/output) COMPLEX*16 array, dimension (LDH,N) */
/*          On entry, the upper Hessenberg matrix H. */
/*          On exit, if WANTT is .TRUE., H is upper triangular in rows */
/*          and columns ILO:IHI, with any 2-by-2 diagonal blocks in */
/*          standard form. If WANTT is .FALSE., the contents of H are */
/*          unspecified on exit. */

/*  LDH     (input) INTEGER */
/*          The leading dimension of the array H. LDH >= max(1,N). */

/*  W       (output) COMPLEX*16 array, dimension (N) */
/*          The computed eigenvalues ILO to IHI are stored in the */
/*          corresponding elements of W. If WANTT is .TRUE., the */
/*          eigenvalues are stored in the same order as on the diagonal */
/*          of the Schur form returned in H, with W(i) = H(i,i). */

/*  ILOZ    (input) INTEGER */
/*  IHIZ    (input) INTEGER */
/*          Specify the rows of Z to which transformations must be */
/*          applied if WANTZ is .TRUE.. */
/*          1 <= ILOZ <= ILO; IHI <= IHIZ <= N. */

/*  Z       (input/output) COMPLEX*16 array, dimension (LDZ,N) */
/*          If WANTZ is .TRUE., on entry Z must contain the current */
/*          matrix Z of transformations accumulated by ZHSEQR, and on */
/*          exit Z has been updated; transformations are applied only to */
/*          the submatrix Z(ILOZ:IHIZ,ILO:IHI). */
/*          If WANTZ is .FALSE., Z is not referenced. */

/*  LDZ     (input) INTEGER */
/*          The leading dimension of the array Z. LDZ >= max(1,N). */

/*  INFO    (output) INTEGER */
/*          = 0: successful exit */
/*          > 0: if INFO = i, ZLAHQR failed to compute all the */
/*               eigenvalues ILO to IHI in a total of 30*(IHI-ILO+1) */
/*               iterations; elements i+1:ihi of W contain those */
/*               eigenvalues which have been successfully computed. */

/*  ===================================================================== */

/*     .. Parameters .. */
/*<       COMPLEX*16         ZERO, ONE >*/
/*<    >*/
/*<       DOUBLE PRECISION   RZERO, HALF >*/
/*<       PARAMETER          ( RZERO = 0.0D+0, HALF = 0.5D+0 ) >*/
/*<       DOUBLE PRECISION   DAT1 >*/
/*<       PARAMETER          ( DAT1 = 0.75D+0 ) >*/
/*     .. */
/*     .. Local Scalars .. */
/*<       INTEGER            I, I1, I2, ITN, ITS, J, K, L, M, NH, NZ >*/
/*<       DOUBLE PRECISION   H10, H21, RTEMP, S, SMLNUM, T2, TST1, ULP >*/
/*<    >*/
/*     .. */
/*     .. Local Arrays .. */
/*<       DOUBLE PRECISION   RWORK( 1 ) >*/
/*<       COMPLEX*16         V( 2 ) >*/
/*     .. */
/*     .. External Functions .. */
/*<       DOUBLE PRECISION   DLAMCH, ZLANHS >*/
/*<       COMPLEX*16         ZLADIV >*/
/*<       EXTERNAL           DLAMCH, ZLANHS, ZLADIV >*/
/*     .. */
/*     .. External Subroutines .. */
/*<       EXTERNAL           ZCOPY, ZLARFG, ZSCAL >*/
/*     .. */
/*     .. Intrinsic Functions .. */
/*<       INTRINSIC          ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SQRT >*/
/*     .. */
/*     .. Statement Functions .. */
/*<       DOUBLE PRECISION   CABS1 >*/
/*     .. */
/*     .. Statement Function definitions .. */
/*<       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) >*/
/*     .. */
/*     .. Executable Statements .. */

/*<       INFO = 0 >*/
    /* Parameter adjustments */
    h_dim1 = *ldh;
    h_offset = 1 + h_dim1;
    h__ -= h_offset;
    --w;
    z_dim1 = *ldz;
    z_offset = 1 + z_dim1;
    z__ -= z_offset;

    /* Function Body */
    *info = 0;

/*     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 >*/
    }

/*<       NH = IHI - ILO + 1 >*/
    nh = *ihi - *ilo + 1;
/*<       NZ = IHIZ - ILOZ + 1 >*/
    nz = *ihiz - *iloz + 1;

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

/*<       ULP = DLAMCH( 'Precision' ) >*/
    ulp = dlamch_("Precision", (ftnlen)9);
/*<       SMLNUM = DLAMCH( 'Safe minimum' ) / ULP >*/
    smlnum = dlamch_("Safe minimum", (ftnlen)12) / ulp;

/*     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 set inside the main loop. */

/*<       IF( WANTT ) THEN >*/
    if (*wantt) {
/*<          I1 = 1 >*/
        i1 = 1;
/*<          I2 = N >*/
        i2 = *n;
/*<       END IF >*/
    }

/*     ITN is the total number of 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 1. 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;
/*<    10 CONTINUE >*/
L10:
/*<    >*/
    if (i__ < *ilo) {
        goto L130;
    }

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

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

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

/*<          DO 20 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__4 = 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__4].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 L30;
            }
/*<    20    CONTINUE >*/
/* L20: */
        }
/*<    30    CONTINUE >*/
L30:
/*<          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 1 has split off. */

/*<    >*/
        if (l >= i__) {
            goto L120;
        }

/*        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 >*/
        if (! (*wantt)) {
/*<             I1 = L >*/
            i1 = l;
/*<             I2 = I >*/
            i2 = i__;
/*<          END IF >*/
        }

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

/*           Exceptional shift. */

/*<             S = DAT1*ABS( DBLE( H( I, I-1 ) ) ) >*/
            i__2 = i__ + (i__ - 1) * h_dim1;
            s = (d__1 = h__[i__2].r, abs(d__1)) * .75;
/*<             T = S + H( I, I ) >*/
            i__2 = i__ + i__ * h_dim1;
            z__1.r = s + h__[i__2].r, z__1.i = h__[i__2].i;
            t.r = z__1.r, t.i = z__1.i;
/*<          ELSE >*/
        } else {

/*           Wilkinson's shift. */

/*<             T = H( I, I ) >*/
            i__2 = i__ + i__ * h_dim1;
            t.r = h__[i__2].r, t.i = h__[i__2].i;
/*<             U = H( I-1, I )*DBLE( H( I, I-1 ) ) >*/
            i__2 = i__ - 1 + i__ * h_dim1;
            i__3 = i__ + (i__ - 1) * h_dim1;
            d__1 = h__[i__3].r;
            z__1.r = d__1 * h__[i__2].r, z__1.i = d__1 * h__[i__2].i;
            u.r = z__1.r, u.i = z__1.i;
/*<             IF( U.NE.ZERO ) THEN >*/
            if (u.r != 0. || u.i != 0.) {
/*<                X = HALF*( H( I-1, I-1 )-T ) >*/
                i__2 = i__ - 1 + (i__ - 1) * h_dim1;
                z__2.r = h__[i__2].r - t.r, z__2.i = h__[i__2].i - t.i;
                z__1.r = z__2.r * .5, z__1.i = z__2.i * .5;
                x.r = z__1.r, x.i = z__1.i;
/*<                Y = SQRT( X*X+U ) >*/
                z__3.r = x.r * x.r - x.i * x.i, z__3.i = x.r * x.i + x.i * 
                        x.r;
                z__2.r = z__3.r + u.r, z__2.i = z__3.i + u.i;