zhseqr.c

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

/* lapack/complex16/zhseqr.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 doublecomplex c_b1 = {0.,0.};
static doublecomplex c_b2 = {1.,0.};
static integer c__1 = 1;
static integer c__4 = 4;
static integer c_n1 = -1;
static integer c__2 = 2;
static integer c__8 = 8;
static integer c__15 = 15;
static logical c_false = FALSE_;

/*<    >*/
/* Subroutine */ int zhseqr_(char *job, char *compz, integer *n, integer *ilo,
         integer *ihi, doublecomplex *h__, integer *ldh, doublecomplex *w, 
        doublecomplex *z__, integer *ldz, doublecomplex *work, integer *lwork,
         integer *info, ftnlen job_len, ftnlen compz_len)
{
    /* System generated locals */
    address a__1[2];
    integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4[2], 
            i__5, i__6;
    doublereal d__1, d__2, d__3, d__4;
    doublecomplex z__1;
    char ch__1[2];

    /* Builtin functions */
    double d_imag(doublecomplex *);
    void d_cnjg(doublecomplex *, doublecomplex *);
    /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);

    /* Local variables */
    integer i__, j, k, l;
    doublecomplex s[225]        /* was [15][15] */, v[16];
    integer i1, i2, ii, nh, nr, ns, nv;
    doublecomplex vv[16];
    integer itn;
    doublecomplex tau;
    integer its;
    doublereal ulp, tst1;
    integer maxb, ierr;
    doublereal unfl;
    doublecomplex temp;
    doublereal ovfl;
    extern logical lsame_(char *, char *, ftnlen, ftnlen);
    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 
            doublecomplex *, integer *);
    integer itemp;
    doublereal rtemp;
    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 
            doublecomplex *, doublecomplex *, integer *, doublecomplex *, 
            integer *, doublecomplex *, doublecomplex *, integer *, ftnlen);
    logical initz, wantt, wantz;
    doublereal rwork[1];
    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 
            doublecomplex *, integer *);
    extern doublereal dlapy2_(doublereal *, doublereal *);
    extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);
    extern doublereal dlamch_(char *, ftnlen);
    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
            integer *, integer *, ftnlen, ftnlen);
    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 
            doublecomplex *, integer *), zlarfg_(integer *, doublecomplex *, 
            doublecomplex *, integer *, doublecomplex *);
    extern integer izamax_(integer *, doublecomplex *, integer *);
    extern doublereal zlanhs_(char *, integer *, doublecomplex *, integer *, 
            doublereal *, ftnlen);
    extern /* Subroutine */ int zlahqr_(logical *, logical *, integer *, 
            integer *, integer *, doublecomplex *, integer *, doublecomplex *,
             integer *, integer *, doublecomplex *, integer *, integer *), 
            zlacpy_(char *, integer *, integer *, doublecomplex *, integer *, 
            doublecomplex *, integer *, ftnlen), zlaset_(char *, integer *, 
            integer *, doublecomplex *, doublecomplex *, doublecomplex *, 
            integer *, ftnlen), zlarfx_(char *, integer *, integer *, 
            doublecomplex *, doublecomplex *, doublecomplex *, integer *, 
            doublecomplex *, ftnlen);
    doublereal smlnum;
    logical lquery;
    (void)job_len;
    (void)compz_len;

/*  -- LAPACK 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 .. */
/*<       CHARACTER          COMPZ, JOB >*/
/*<       INTEGER            IHI, ILO, INFO, LDH, LDZ, LWORK, N >*/
/*     .. */
/*     .. Array Arguments .. */
/*<       COMPLEX*16         H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * ) >*/
/*     .. */

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

/*  ZHSEQR computes the eigenvalues of a complex upper Hessenberg */
/*  matrix H, and, optionally, the matrices T and Z from the Schur */
/*  decomposition H = Z T Z**H, where T is an upper triangular matrix */
/*  (the Schur form), and Z is the unitary matrix of Schur vectors. */

/*  Optionally Z may be postmultiplied into an input unitary matrix Q, */
/*  so that this routine can give the Schur factorization of a matrix A */
/*  which has been reduced to the Hessenberg form H by the unitary */
/*  matrix Q:  A = Q*H*Q**H = (QZ)*T*(QZ)**H. */

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

/*  JOB     (input) CHARACTER*1 */
/*          = 'E': compute eigenvalues only; */
/*          = 'S': compute eigenvalues and the Schur form T. */

/*  COMPZ   (input) CHARACTER*1 */
/*          = 'N': no Schur vectors are computed; */
/*          = 'I': Z is initialized to the unit matrix and the matrix Z */
/*                 of Schur vectors of H is returned; */
/*          = 'V': Z must contain an unitary matrix Q on entry, and */
/*                 the product Q*Z is returned. */

/*  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 1:ILO-1 and IHI+1:N. ILO and IHI are normally */
/*          set by a previous call to ZGEBAL, and then passed to CGEHRD */
/*          when the matrix output by ZGEBAL is reduced to Hessenberg */
/*          form. Otherwise ILO and IHI should be set to 1 and N */
/*          respectively. */
/*          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. */

/*  H       (input/output) COMPLEX*16 array, dimension (LDH,N) */
/*          On entry, the upper Hessenberg matrix H. */
/*          On exit, if JOB = 'S', H contains the upper triangular matrix */
/*          T from the Schur decomposition (the Schur form). If */
/*          JOB = 'E', 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. If JOB = 'S', 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). */

/*  Z       (input/output) COMPLEX*16 array, dimension (LDZ,N) */
/*          If COMPZ = 'N': Z is not referenced. */
/*          If COMPZ = 'I': on entry, Z need not be set, and on exit, Z */
/*          contains the unitary matrix Z of the Schur vectors of H. */
/*          If COMPZ = 'V': on entry Z must contain an N-by-N matrix Q, */
/*          which is assumed to be equal to the unit matrix except for */
/*          the submatrix Z(ILO:IHI,ILO:IHI); on exit Z contains Q*Z. */
/*          Normally Q is the unitary matrix generated by ZUNGHR after */
/*          the call to ZGEHRD which formed the Hessenberg matrix H. */

/*  LDZ     (input) INTEGER */
/*          The leading dimension of the array Z. */
/*          LDZ >= max(1,N) if COMPZ = 'I' or 'V'; LDZ >= 1 otherwise. */

/*  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK) */
/*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */

/*  LWORK   (input) INTEGER */
/*          The dimension of the array WORK.  LWORK >= max(1,N). */

/*          If LWORK = -1, then a workspace query is assumed; the routine */
/*          only calculates the optimal size of the WORK array, returns */
/*          this value as the first entry of the WORK array, and no error */
/*          message related to LWORK is issued by XERBLA. */

/*  INFO    (output) INTEGER */
/*          = 0:  successful exit */
/*          < 0:  if INFO = -i, the i-th argument had an illegal value */
/*          > 0:  if INFO = i, ZHSEQR failed to compute all the */
/*                eigenvalues in a total of 30*(IHI-ILO+1) iterations; */
/*                elements 1:ilo-1 and i+1:n of W contain those */
/*                eigenvalues which have been successfully computed. */

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

/*     .. Parameters .. */
/*<       COMPLEX*16         ZERO, ONE >*/
/*<    >*/
/*<       DOUBLE PRECISION   RZERO, RONE, CONST >*/
/*<    >*/
/*<       INTEGER            NSMAX, LDS >*/
/*<       PARAMETER          ( NSMAX = 15, LDS = NSMAX ) >*/
/*     .. */
/*     .. Local Scalars .. */
/*<       LOGICAL            INITZ, LQUERY, WANTT, WANTZ >*/
/*<    >*/
/*<       DOUBLE PRECISION   OVFL, RTEMP, SMLNUM, TST1, ULP, UNFL >*/
/*<       COMPLEX*16         CDUM, TAU, TEMP >*/
/*     .. */
/*     .. Local Arrays .. */
/*<       DOUBLE PRECISION   RWORK( 1 ) >*/
/*<       COMPLEX*16         S( LDS, NSMAX ), V( NSMAX+1 ), VV( NSMAX+1 ) >*/
/*     .. */
/*     .. External Functions .. */
/*<       LOGICAL            LSAME >*/
/*<       INTEGER            ILAENV, IZAMAX >*/
/*<       DOUBLE PRECISION   DLAMCH, DLAPY2, ZLANHS >*/
/*<       EXTERNAL           LSAME, ILAENV, IZAMAX, DLAMCH, DLAPY2, ZLANHS >*/
/*     .. */
/*     .. External Subroutines .. */
/*<    >*/
/*     .. */
/*     .. Intrinsic Functions .. */
/*<       INTRINSIC          ABS, DBLE, DCONJG, DIMAG, MAX, MIN >*/
/*     .. */
/*     .. Statement Functions .. */
/*<       DOUBLE PRECISION   CABS1 >*/
/*     .. */
/*     .. Statement Function definitions .. */
/*<       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) >*/
/*     .. */
/*     .. Executable Statements .. */

/*     Decode and test the input parameters */

/*<       WANTT = LSAME( JOB, 'S' ) >*/
    /* 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;
    --work;

    /* Function Body */
    wantt = lsame_(job, "S", (ftnlen)1, (ftnlen)1);
/*<       INITZ = LSAME( COMPZ, 'I' ) >*/
    initz = lsame_(compz, "I", (ftnlen)1, (ftnlen)1);
/*<       WANTZ = INITZ .OR. LSAME( COMPZ, 'V' ) >*/
    wantz = initz || lsame_(compz, "V", (ftnlen)1, (ftnlen)1);

/*<       INFO = 0 >*/
    *info = 0;
/*<       WORK( 1 ) = MAX( 1, N ) >*/
    i__1 = max(1,*n);
    work[1].r = (doublereal) i__1, work[1].i = 0.;
/*<       LQUERY = ( LWORK.EQ.-1 ) >*/
    lquery = *lwork == -1;
/*<       IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN >*/
    if (! lsame_(job, "E", (ftnlen)1, (ftnlen)1) && ! wantt) {
/*<          INFO = -1 >*/
        *info = -1;
/*<       ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN >*/
    } else if (! lsame_(compz, "N", (ftnlen)1, (ftnlen)1) && ! wantz) {
/*<          INFO = -2 >*/
        *info = -2;
/*<       ELSE IF( N.LT.0 ) THEN >*/
    } else if (*n < 0) {
/*<          INFO = -3 >*/
        *info = -3;
/*<       ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN >*/
    } else if (*ilo < 1 || *ilo > max(1,*n)) {
/*<          INFO = -4 >*/
        *info = -4;
/*<       ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN >*/
    } else if (*ihi < min(*ilo,*n) || *ihi > *n) {
/*<          INFO = -5 >*/
        *info = -5;
/*<       ELSE IF( LDH.LT.MAX( 1, N ) ) THEN >*/
    } else if (*ldh < max(1,*n)) {
/*<          INFO = -7 >*/
        *info = -7;
/*<       ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. LDZ.LT.MAX( 1, N ) ) THEN >*/
    } else if (*ldz < 1 || (wantz && *ldz < max(1,*n))) {