自学教程：C++ zlaset_函数代码示例

//.........这里部分代码省略........./*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. Scalars in Common .. *//*     .. *//*     .. Common blocks .. *//*     .. *//*     .. Executable Statements .. *//*     Quick return if possible */    /* Parameter adjustments */    r_dim1 = *lda;    r_offset = 1 + r_dim1;    r__ -= r_offset;    q_dim1 = *lda;    q_offset = 1 + q_dim1;    q -= q_offset;    af_dim1 = *lda;    af_offset = 1 + af_dim1;    af -= af_offset;    a_dim1 = *lda;    a_offset = 1 + a_dim1;    a -= a_offset;    --tau;    --work;    --rwork;    --result;    /* Function Body */    if (*m == 0 || *n == 0 || *k == 0) {	result[1] = 0.;	result[2] = 0.;	return 0;    }    eps = dlamch_("Epsilon");/*     Copy the last k rows of the factorization to the array Q */    zlaset_("Full", m, n, &c_b1, &c_b1, &q[q_offset], lda);    if (*k < *n) {	i__1 = *n - *k;	zlacpy_("Full", k, &i__1, &af[*m - *k + 1 + af_dim1], lda, &q[*m - *k 		+ 1 + q_dim1], lda);    }    if (*k > 1) {	i__1 = *k - 1;	i__2 = *k - 1;	zlacpy_("Lower", &i__1, &i__2, &af[*m - *k + 2 + (*n - *k + 1) * 		af_dim1], lda, &q[*m - *k + 2 + (*n - *k + 1) * q_dim1], lda);    }/*     Generate the last n rows of the matrix Q */    s_copy(srnamc_1.srnamt, "ZUNGRQ", (ftnlen)32, (ftnlen)6);    zungrq_(m, n, k, &q[q_offset], lda, &tau[*m - *k + 1], &work[1], lwork, &	    info);/*     Copy R(m-k+1:m,n-m+1:n) */    zlaset_("Full", k, m, &c_b9, &c_b9, &r__[*m - *k + 1 + (*n - *m + 1) * 	    r_dim1], lda);    zlacpy_("Upper", k, k, &af[*m - *k + 1 + (*n - *k + 1) * af_dim1], lda, &	    r__[*m - *k + 1 + (*n - *k + 1) * r_dim1], lda);/*     Compute R(m-k+1:m,n-m+1:n) - A(m-k+1:m,1:n) * Q(n-m+1:n,1:n)' */    zgemm_("No transpose", "Conjugate transpose", k, m, n, &c_b14, &a[*m - *k 	    + 1 + a_dim1], lda, &q[q_offset], lda, &c_b15, &r__[*m - *k + 1 + 	    (*n - *m + 1) * r_dim1], lda);/*     Compute norm( R - A*Q' ) / ( N * norm(A) * EPS ) . */    anorm = zlange_("1", k, n, &a[*m - *k + 1 + a_dim1], lda, &rwork[1]);    resid = zlange_("1", k, m, &r__[*m - *k + 1 + (*n - *m + 1) * r_dim1], 	    lda, &rwork[1]);    if (anorm > 0.) {	result[1] = resid / (doublereal) max(1,*n) / anorm / eps;    } else {	result[1] = 0.;    }/*     Compute I - Q*Q' */    zlaset_("Full", m, m, &c_b9, &c_b15, &r__[r_offset], lda);    zherk_("Upper", "No transpose", m, n, &c_b23, &q[q_offset], lda, &c_b24, &	    r__[r_offset], lda);/*     Compute norm( I - Q*Q' ) / ( N * EPS ) . */    resid = zlansy_("1", "Upper", m, &r__[r_offset], lda, &rwork[1]);    result[2] = resid / (doublereal) max(1,*n) / eps;    return 0;/*     End of ZRQT02 */} /* zrqt02_ */

//.........这里部分代码省略.........	     doublecomplex *, integer *);    static char xtype[1];    extern /* Subroutine */ int zlatb4_(char *, integer *, integer *, integer 	    *, char *, integer *, integer *, doublereal *, integer *, 	    doublereal *, char *), aladhd_(integer *, 	    char *);    static integer kd, nb, in, kl;    extern /* Subroutine */ int alaerh_(char *, char *, integer *, integer *, 	    char *, integer *, integer *, integer *, integer *, integer *, 	    integer *, integer *, integer *, integer *);    static logical prefac;    static integer iw, ku, nt;    static doublereal rcondc;    static logical nofact;    static char packit[1];    static integer iequed;    extern doublereal zlanhb_(char *, char *, integer *, integer *, 	    doublecomplex *, integer *, doublereal *), 	    zlange_(char *, integer *, integer *, doublecomplex *, integer *, 	    doublereal *);    extern /* Subroutine */ int zlaqhb_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublereal *, doublereal *, 	    doublereal *, char *), alasvm_(char *, integer *, 	    integer *, integer *, integer *);    static doublereal cndnum;    extern /* Subroutine */ int zlaipd_(integer *, doublecomplex *, integer *,	     integer *);    static doublereal ainvnm;    extern /* Subroutine */ int xlaenv_(integer *, integer *), zlacpy_(char *,	     integer *, integer *, doublecomplex *, integer *, doublecomplex *	    , integer *), zlarhs_(char *, char *, char *, char *, 	    integer *, integer *, integer *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, integer *, integer *), zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *), zpbequ_(char *, integer *, integer *, doublecomplex *, 	    integer *, doublereal *, doublereal *, doublereal *, integer *), zpbtrf_(char *, integer *, integer *, doublecomplex *, 	    integer *, integer *), zlatms_(integer *, integer *, char 	    *, integer *, char *, doublereal *, integer *, doublereal *, 	    doublereal *, integer *, integer *, char *, doublecomplex *, 	    integer *, doublecomplex *, integer *);    static doublereal result[6];    extern /* Subroutine */ int zpbtrs_(char *, integer *, integer *, integer 	    *, doublecomplex *, integer *, doublecomplex *, integer *, 	    integer *), zpbsvx_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, integer *, doublecomplex *, integer *,	     char *, doublereal *, doublecomplex *, integer *, doublecomplex *	    , integer *, doublereal *, doublereal *, doublereal *, 	    doublecomplex *, doublereal *, integer *),	     zerrvx_(char *, integer *);    static integer lda, ikd, nkd;    /* Fortran I/O blocks */    static cilist io___57 = { 0, 0, 0, fmt_9999, 0 };    static cilist io___60 = { 0, 0, 0, fmt_9997, 0 };    static cilist io___61 = { 0, 0, 0, fmt_9998, 0 };/*  -- LAPACK test routine (version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          June 30, 1999       Purpose       =======   

/* Subroutine */ int zdrvpt_(logical *dotype, integer *nn, integer *nval, 	integer *nrhs, doublereal *thresh, logical *tsterr, doublecomplex *a, 	doublereal *d__, doublecomplex *e, doublecomplex *b, doublecomplex *x, 	 doublecomplex *xact, doublecomplex *work, doublereal *rwork, integer 	*nout){    /* Initialized data */    static integer iseedy[4] = { 0,0,0,1 };    /* Format strings */    static char fmt_9999[] = "(1x,a6,/002, N =/002,i5,/002, type /002,i2,"	    "/002, test /002,i2,/002, ratio = /002,g12.5)";    static char fmt_9998[] = "(1x,a6,/002, FACT='/002,a1,/002', N =/002,i5"	    ",/002, type /002,i2,/002, test /002,i2,/002, ratio = /002,g12.5)";    /* System generated locals */    integer i__1, i__2, i__3, i__4, i__5;    doublereal d__1, d__2;    /* Builtin functions */    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);    double z_abs(doublecomplex *);    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);    /* Local variables */    integer i__, j, k, n;    doublereal z__[3];    integer k1, ia, in, kl, ku, ix, nt, lda;    char fact[1];    doublereal cond;    integer mode;    doublereal dmax__;    integer imat, info;    char path[3], dist[1], type__[1];    integer nrun, ifact;    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 	    integer *);    integer nfail, iseed[4];    extern doublereal dget06_(doublereal *, doublereal *);    doublereal rcond;    integer nimat;    doublereal anorm;    extern /* Subroutine */ int zget04_(integer *, integer *, doublecomplex *, 	     integer *, doublecomplex *, integer *, doublereal *, doublereal *), dcopy_(integer *, doublereal *, integer *, doublereal *, 	    integer *);    integer izero, nerrs;    extern /* Subroutine */ int zptt01_(integer *, doublereal *, 	    doublecomplex *, doublereal *, doublecomplex *, doublecomplex *, 	    doublereal *);    logical zerot;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zptt02_(char *, integer *, integer *, 	     doublereal *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *), zptt05_(	    integer *, integer *, doublereal *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublereal *, 	    doublereal *), zptsv_(integer *, integer *, doublereal *, 	    doublecomplex *, doublecomplex *, integer *, integer *), zlatb4_(	    char *, integer *, integer *, integer *, char *, integer *, 	    integer *, doublereal *, integer *, doublereal *, char *), aladhd_(integer *, char *), alaerh_(char 	    *, char *, integer *, integer *, char *, integer *, integer *, 	    integer *, integer *, integer *, integer *, integer *, integer *, 	    integer *);    extern integer idamax_(integer *, doublereal *, integer *);    doublereal rcondc;    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *), alasvm_(char *, integer *, integer *, 	     integer *, integer *), dlarnv_(integer *, integer *, 	    integer *, doublereal *);    doublereal ainvnm;    extern doublereal zlanht_(char *, integer *, doublereal *, doublecomplex *);    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern doublereal dzasum_(integer *, doublecomplex *, integer *);    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *), zlaptm_(char *, integer *, integer *, doublereal *, 	    doublereal *, doublecomplex *, doublecomplex *, integer *, 	    doublereal *, doublecomplex *, integer *), zlatms_(	    integer *, integer *, char *, integer *, char *, doublereal *, 	    integer *, doublereal *, doublereal *, integer *, integer *, char 	    *, doublecomplex *, integer *, doublecomplex *, integer *), zlarnv_(integer *, integer *, integer *, 	    doublecomplex *);    doublereal result[6];    extern /* Subroutine */ int zpttrf_(integer *, doublereal *, 	    doublecomplex *, integer *), zerrvx_(char *, integer *), 	    zpttrs_(char *, integer *, integer *, doublereal *, doublecomplex 	    *, doublecomplex *, integer *, integer *), zptsvx_(char *, 	     integer *, integer *, doublereal *, doublecomplex *, doublereal *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublereal *, doublereal *, doublereal *, 	    doublecomplex *, doublereal *, integer *);    /* Fortran I/O blocks */    static cilist io___35 = { 0, 0, 0, fmt_9999, 0 };    static cilist io___38 = { 0, 0, 0, fmt_9998, 0 };//.........这里部分代码省略.........

//.........这里部分代码省略........., doublecomplex *, integer *, doublecomplex *, doublereal *, 	    doublereal *), zgges_(char *, char *, char *, L_fp, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, integer *, 	     doublecomplex *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublereal *, logical *, integer *);    integer nmats, jsize;    extern /* Subroutine */ int zget54_(integer *, doublecomplex *, integer *, 	     doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublereal *);    integer nerrs, jtype, ntest, isort;    extern /* Subroutine */ int dlabad_(doublereal *, doublereal *), zlatm4_(	    integer *, integer *, integer *, integer *, logical *, doublereal 	    *, doublereal *, doublereal *, integer *, integer *, 	    doublecomplex *, integer *);    logical ilabad;    extern doublereal dlamch_(char *);    extern /* Subroutine */ int zunm2r_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    doublereal safmin, safmax;    integer knteig, ioldsd[4];    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 	    integer *, integer *);    extern /* Subroutine */ int alasvm_(char *, integer *, integer *, integer 	    *, integer *), xerbla_(char *, integer *), 	    zlarfg_(integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *);    extern /* Double Complex */ void zlarnd_(doublecomplex *, integer *, 	    integer *);    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *), 	    zlaset_(char *, integer *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *, integer *);    extern logical zlctes_(doublecomplex *, doublecomplex *);    integer minwrk, maxwrk;    doublereal ulpinv;    integer mtypes, ntestt;    /* Fortran I/O blocks */    static cilist io___41 = { 0, 0, 0, fmt_9999, 0 };    static cilist io___47 = { 0, 0, 0, fmt_9999, 0 };    static cilist io___51 = { 0, 0, 0, fmt_9998, 0 };    static cilist io___53 = { 0, 0, 0, fmt_9997, 0 };    static cilist io___54 = { 0, 0, 0, fmt_9996, 0 };    static cilist io___55 = { 0, 0, 0, fmt_9995, 0 };    static cilist io___56 = { 0, 0, 0, fmt_9994, 0 };    static cilist io___57 = { 0, 0, 0, fmt_9993, 0 };    static cilist io___58 = { 0, 0, 0, fmt_9992, 0 };    static cilist io___59 = { 0, 0, 0, fmt_9991, 0 };/*  -- LAPACK test routine (version 3.1.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     February 2007 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= */

//.........这里部分代码省略.........    static integer c__1 = 1;    static integer c_n1 = -1;    static integer c__0 = 0;        /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;    doublereal d__1;    /* Local variables */    static doublereal anrm, bnrm;    static integer brow;    static logical tpsd;    static integer i__, j, iascl, ibscl;    extern logical lsame_(char *, char *);    static integer wsize;    static doublereal rwork[1];    extern /* Subroutine */ int ztrsm_(char *, char *, char *, char *, 	    integer *, integer *, doublecomplex *, doublecomplex *, integer *,	     doublecomplex *, integer *), 	    dlabad_(doublereal *, doublereal *);    static integer nb;    extern doublereal dlamch_(char *);    static integer mn;    extern /* Subroutine */ int xerbla_(char *, integer *);    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 	    integer *, integer *, ftnlen, ftnlen);    static integer scllen;    static doublereal bignum;    extern doublereal zlange_(char *, integer *, integer *, doublecomplex *, 	    integer *, doublereal *);    extern /* Subroutine */ int zgelqf_(integer *, integer *, doublecomplex *,	     integer *, doublecomplex *, doublecomplex *, integer *, integer *	    ), zlascl_(char *, integer *, integer *, doublereal *, doublereal 	    *, integer *, integer *, doublecomplex *, integer *, integer *), zgeqrf_(integer *, integer *, doublecomplex *, integer *,	     doublecomplex *, doublecomplex *, integer *, integer *), zlaset_(	    char *, integer *, integer *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *);    static doublereal smlnum;    static logical lquery;    extern /* Subroutine */ int zunmlq_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, integer *), zunmqr_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, integer *);#define b_subscr(a_1,a_2) (a_2)*b_dim1 + a_1#define b_ref(a_1,a_2) b[b_subscr(a_1,a_2)]    a_dim1 = *lda;    a_offset = 1 + a_dim1 * 1;    a -= a_offset;    b_dim1 = *ldb;    b_offset = 1 + b_dim1 * 1;    b -= b_offset;    --work;    /* Function Body */    *info = 0;    mn = min(*m,*n);    lquery = *lwork == -1;    if (! (lsame_(trans, "N") || lsame_(trans, "C"))) {	*info = -1;    } else if (*m < 0) {	*info = -2;    } else if (*n < 0) {	*info = -3;    } else if (*nrhs < 0) {

/* Subroutine */ int zget22_(char *transa, char *transe, char *transw, 	integer *n, doublecomplex *a, integer *lda, doublecomplex *e, integer 	*lde, doublecomplex *w, doublecomplex *work, doublereal *rwork, 	doublereal *result){    /* System generated locals */    integer a_dim1, a_offset, e_dim1, e_offset, i__1, i__2, i__3, i__4;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1, z__2;    /* Builtin functions */    double d_imag(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    integer j;    doublereal ulp;    integer joff, jcol, jvec;    doublereal unfl;    integer jrow;    doublereal temp1;    extern logical lsame_(char *, char *);    char norma[1];    doublereal anorm;    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *);    char norme[1];    doublereal enorm;    doublecomplex wtemp;    extern doublereal dlamch_(char *), zlange_(char *, integer *, 	    integer *, doublecomplex *, integer *, doublereal *);    doublereal enrmin, enrmax;    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);    integer itrnse;    doublereal errnrm;    integer itrnsw;/*  -- LAPACK test routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZGET22 does an eigenvector check. *//*  The basic test is: *//*     RESULT(1) = | A E  -  E W | / ( |A| |E| ulp ) *//*  using the 1-norm.  It also tests the normalization of E: *//*     RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp ) *//*                  j *//*  where E(j) is the j-th eigenvector, and m-norm is the max-norm of a *//*  vector.  The max-norm of a complex n-vector x in this case is the *//*  maximum of |re(x(i)| + |im(x(i)| over i = 1, ..., n. *//*  Arguments *//*  ========== *//*  TRANSA  (input) CHARACTER*1 *//*          Specifies whether or not A is transposed. *//*          = 'N':  No transpose *//*          = 'T':  Transpose *//*          = 'C':  Conjugate transpose *//*  TRANSE  (input) CHARACTER*1 *//*          Specifies whether or not E is transposed. *//*          = 'N':  No transpose, eigenvectors are in columns of E *//*          = 'T':  Transpose, eigenvectors are in rows of E *//*          = 'C':  Conjugate transpose, eigenvectors are in rows of E *//*  TRANSW  (input) CHARACTER*1 *//*          Specifies whether or not W is transposed. *//*          = 'N':  No transpose *//*          = 'T':  Transpose, same as TRANSW = 'N' *//*          = 'C':  Conjugate transpose, use -WI(j) instead of WI(j) *//*  N       (input) INTEGER *//*          The order of the matrix A.  N >= 0. *//*  A       (input) COMPLEX*16 array, dimension (LDA,N) *//*          The matrix whose eigenvectors are in E. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A.  LDA >= max(1,N). *//*  E       (input) COMPLEX*16 array, dimension (LDE,N) *//*          The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors *///.........这里部分代码省略.........

//.........这里部分代码省略.........    wantc = *ncc > 0;    klu1 = *kl + *ku + 1;    *info = 0;    if (! wantq && ! wantpt && ! lsame_(vect, "N")) {	*info = -1;    } else if (*m < 0) {	*info = -2;    } else if (*n < 0) {	*info = -3;    } else if (*ncc < 0) {	*info = -4;    } else if (*kl < 0) {	*info = -5;    } else if (*ku < 0) {	*info = -6;    } else if (*ldab < klu1) {	*info = -8;    } else if (*ldq < 1 || wantq && *ldq < max(1,*m)) {	*info = -12;    } else if (*ldpt < 1 || wantpt && *ldpt < max(1,*n)) {	*info = -14;    } else if (*ldc < 1 || wantc && *ldc < max(1,*m)) {	*info = -16;    }    if (*info != 0) {	i__1 = -(*info);	xerbla_("ZGBBRD", &i__1);	return 0;    }/*     Initialize Q and P' to the unit matrix, if needed */    if (wantq) {	zlaset_("Full", m, m, &c_b1, &c_b2, &q[q_offset], ldq);    }    if (wantpt) {	zlaset_("Full", n, n, &c_b1, &c_b2, &pt[pt_offset], ldpt);    }/*     Quick return if possible. */    if (*m == 0 || *n == 0) {	return 0;    }    minmn = min(*m,*n);    if (*kl + *ku > 1) {/*        Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce *//*        first to lower bidiagonal form and then transform to upper *//*        bidiagonal */	if (*ku > 0) {	    ml0 = 1;	    mu0 = 2;	} else {	    ml0 = 2;	    mu0 = 1;	}/*        Wherever possible, plane rotations are generated and applied in *//*        vector operations of length NR over the index set J1:J2:KLU1. *//*        The complex sines of the plane rotations are stored in WORK, *//*        and the real cosines in RWORK. */

/* Subroutine */ int zlarrv_(integer *n, doublereal *vl, doublereal *vu, 	doublereal *d__, doublereal *l, doublereal *pivmin, integer *isplit, 	integer *m, integer *dol, integer *dou, doublereal *minrgp, 	doublereal *rtol1, doublereal *rtol2, doublereal *w, doublereal *werr, 	 doublereal *wgap, integer *iblock, integer *indexw, doublereal *gers, 	 doublecomplex *z__, integer *ldz, integer *isuppz, doublereal *work, 	integer *iwork, integer *info){    /* System generated locals */    integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5, i__6;    doublereal d__1, d__2;    doublecomplex z__1;    logical L__1;    /* Builtin functions */    double log(doublereal);    /* Local variables */    integer minwsize, i__, j, k, p, q, miniwsize, ii;    doublereal gl;    integer im, in;    doublereal gu, gap, eps, tau, tol, tmp;    integer zto;    doublereal ztz;    integer iend, jblk;    doublereal lgap;    integer done;    doublereal rgap, left;    integer wend, iter;    doublereal bstw;    integer itmp1, indld;    doublereal fudge;    integer idone;    doublereal sigma;    integer iinfo, iindr;    doublereal resid;    logical eskip;    doublereal right;    extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 	    doublereal *, integer *);    integer nclus, zfrom;    doublereal rqtol;    integer iindc1, iindc2, indin1, indin2;    logical stp2ii;    extern /* Subroutine */ int zlar1v_(integer *, integer *, integer *, 	    doublereal *, doublereal *, doublereal *, doublereal *, 	    doublereal *, doublereal *, doublereal *, doublecomplex *, 	    logical *, integer *, doublereal *, doublereal *, integer *, 	    integer *, doublereal *, doublereal *, doublereal *, doublereal *)	    ;    doublereal lambda;    extern doublereal dlamch_(char *);    integer ibegin, indeig;    logical needbs;    integer indlld;    doublereal sgndef, mingma;    extern /* Subroutine */ int dlarrb_(integer *, doublereal *, doublereal *, 	     integer *, integer *, doublereal *, doublereal *, integer *, 	    doublereal *, doublereal *, doublereal *, doublereal *, integer *, 	     doublereal *, doublereal *, integer *, integer *);    integer oldien, oldncl, wbegin;    doublereal spdiam;    integer negcnt;    extern /* Subroutine */ int dlarrf_(integer *, doublereal *, doublereal *, 	     doublereal *, integer *, integer *, doublereal *, doublereal *, 	    doublereal *, doublereal *, doublereal *, doublereal *, 	    doublereal *, doublereal *, doublereal *, doublereal *, 	    doublereal *, integer *);    integer oldcls;    doublereal savgap;    integer ndepth;    doublereal ssigma;    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *);    logical usedbs;    integer iindwk, offset;    doublereal gaptol;    integer newcls, oldfst, indwrk, windex, oldlst;    logical usedrq;    integer newfst, newftt, parity, windmn, windpl, isupmn, newlst, zusedl;    doublereal bstres;    integer newsiz, zusedu, zusedw;    doublereal nrminv;    logical tryrqc;    integer isupmx;    doublereal rqcorr;    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *///.........这里部分代码省略.........

/*<    >*//* Subroutine */ int zgghrd_(char *compq, char *compz, integer *n, integer *        ilo, integer *ihi, doublecomplex *a, integer *lda, doublecomplex *b,        integer *ldb, doublecomplex *q, integer *ldq, doublecomplex *z__,        integer *ldz, integer *info, ftnlen compq_len, ftnlen compz_len){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1,            z_offset, i__1, i__2, i__3;    doublecomplex z__1;    /* Builtin functions */    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    doublereal c__;    doublecomplex s;    logical ilq, ilz;    integer jcol, jrow;    extern /* Subroutine */ int zrot_(integer *, doublecomplex *, integer *,            doublecomplex *, integer *, doublereal *, doublecomplex *);    extern logical lsame_(const char *, const char *, ftnlen, ftnlen);    doublecomplex ctemp;    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);    integer icompq, icompz;    extern /* Subroutine */ int zlaset_(char *, integer *, integer *,            doublecomplex *, doublecomplex *, doublecomplex *, integer *,            ftnlen), zlartg_(doublecomplex *, doublecomplex *, doublereal *,            doublecomplex *, doublecomplex *);    (void)compq_len;    (void)compz_len;/*  -- LAPACK routine (version 3.2) -- *//*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- *//*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- *//*     November 2006 *//*     .. Scalar Arguments .. *//*<       CHARACTER          COMPQ, COMPZ >*//*<       INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N >*//*     .. *//*     .. Array Arguments .. *//*<    >*//*     .. *//*  Purpose *//*  ======= *//*  ZGGHRD reduces a pair of complex matrices (A,B) to generalized upper *//*  Hessenberg form using unitary transformations, where A is a *//*  general matrix and B is upper triangular.  The form of the *//*  generalized eigenvalue problem is *//*     A*x = lambda*B*x, *//*  and B is typically made upper triangular by computing its QR *//*  factorization and moving the unitary matrix Q to the left side *//*  of the equation. *//*  This subroutine simultaneously reduces A to a Hessenberg matrix H: *//*     Q**H*A*Z = H *//*  and transforms B to another upper triangular matrix T: *//*     Q**H*B*Z = T *//*  in order to reduce the problem to its standard form *//*     H*y = lambda*T*y *//*  where y = Z**H*x. *//*  The unitary matrices Q and Z are determined as products of Givens *//*  rotations.  They may either be formed explicitly, or they may be *//*  postmultiplied into input matrices Q1 and Z1, so that *//*       Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H *//*       Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H *//*  If Q1 is the unitary matrix from the QR factorization of B in the *//*  original equation A*x = lambda*B*x, then ZGGHRD reduces the original *//*  problem to generalized Hessenberg form. *//*  Arguments *//*  ========= *//*  COMPQ   (input) CHARACTER*1 *//*          = 'N': do not compute Q; *//*          = 'I': Q is initialized to the unit matrix, and the *//*                 unitary matrix Q is returned; *//*          = 'V': Q must contain a unitary matrix Q1 on entry, *//*                 and the product Q1*Q is returned. *//*  COMPZ   (input) CHARACTER*1 *//*          = 'N': do not compute Q; *//*          = 'I': Q is initialized to the unit matrix, and the *//*                 unitary matrix Q is returned; *//*          = 'V': Q must contain a unitary matrix Q1 on entry, *//*                 and the product Q1*Q is returned. *//*  N       (input) INTEGER *//*          The order of the matrices A and B.  N >= 0. *//*  ILO     (input) INTEGER *//*  IHI     (input) INTEGER *//*          ILO and IHI mark the rows and columns of A which are to be *//*          reduced.  It is assumed that A 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 ZGGBAL; otherwise they *///.........这里部分代码省略.........

doublereal ztzt01_(integer *m, integer *n, doublecomplex *a, doublecomplex *	af, integer *lda, doublecomplex *tau, doublecomplex *work, integer *	lwork){    /* System generated locals */    integer a_dim1, a_offset, af_dim1, af_offset, i__1, i__2, i__3, i__4;    doublereal ret_val;    /* Local variables */    static integer i__, j;    static doublereal norma, rwork[1];    extern /* Subroutine */ int zaxpy_(integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern doublereal dlamch_(char *);    extern /* Subroutine */ int xerbla_(char *, integer *);    extern doublereal zlange_(char *, integer *, integer *, doublecomplex *, 	    integer *, doublereal *);    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *), zlatzm_(char *, integer *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, doublecomplex *, 	    integer *, doublecomplex *);#define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1#define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]#define af_subscr(a_1,a_2) (a_2)*af_dim1 + a_1#define af_ref(a_1,a_2) af[af_subscr(a_1,a_2)]/*  -- LAPACK test routine (version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          September 30, 1994       Purpose       =======       ZTZT01 returns            || A - R*Q || / ( M * eps * ||A|| )       for an upper trapezoidal A that was factored with ZTZRQF.       Arguments       =========       M       (input) INTEGER               The number of rows of the matrices A and AF.       N       (input) INTEGER               The number of columns of the matrices A and AF.       A       (input) COMPLEX*16 array, dimension (LDA,N)               The original upper trapezoidal M by N matrix A.       AF      (input) COMPLEX*16 array, dimension (LDA,N)               The output of ZTZRQF for input matrix A.               The lower triangle is not referenced.       LDA     (input) INTEGER               The leading dimension of the arrays A and AF.       TAU     (input) COMPLEX*16 array, dimension (M)               Details of the  Householder transformations as returned by               ZTZRQF.       WORK    (workspace) COMPLEX*16 array, dimension (LWORK)       LWORK   (input) INTEGER               The length of the array WORK.  LWORK >= m*n + m.       =====================================================================          Parameter adjustments */    af_dim1 = *lda;    af_offset = 1 + af_dim1 * 1;    af -= af_offset;    a_dim1 = *lda;    a_offset = 1 + a_dim1 * 1;    a -= a_offset;    --tau;    --work;    /* Function Body */    ret_val = 0.;    if (*lwork < *m * *n + *m) {	xerbla_("ZTZT01", &c__8);	return ret_val;    }/*     Quick return if possible */    if (*m <= 0 || *n <= 0) {	return ret_val;    }    norma = zlange_("One-norm", m, n, &a[a_offset], lda, rwork);/*     Copy upper triangle R *///.........这里部分代码省略.........

/* Subroutine */ int zlalsd_(char *uplo, integer *smlsiz, integer *n, integer 	*nrhs, doublereal *d__, doublereal *e, doublecomplex *b, integer *ldb,	 doublereal *rcond, integer *rank, doublecomplex *work, doublereal *	rwork, integer *iwork, integer *info){    /* System generated locals */    integer b_dim1, b_offset, i__1, i__2, i__3, i__4, i__5, i__6;    doublereal d__1;    doublecomplex z__1;    /* Builtin functions */    double d_imag(doublecomplex *), log(doublereal), d_sign(doublereal *, 	    doublereal *);    /* Local variables */    static integer difl, difr, jcol, irwb, perm, nsub, nlvl, sqre, bxst, jrow,	     irwu, c__, i__, j, k;    static doublereal r__;    static integer s, u, jimag;    extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *, 	    integer *, doublereal *, doublereal *, integer *, doublereal *, 	    integer *, doublereal *, doublereal *, integer *);    static integer z__, jreal, irwib, poles, sizei, irwrb, nsize;    extern /* Subroutine */ int zdrot_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublereal *), zcopy_(	    integer *, doublecomplex *, integer *, doublecomplex *, integer *)	    ;    static integer irwvt, icmpq1, icmpq2;    static doublereal cs;    extern doublereal dlamch_(char *);    extern /* Subroutine */ int dlasda_(integer *, integer *, integer *, 	    integer *, doublereal *, doublereal *, doublereal *, integer *, 	    doublereal *, integer *, doublereal *, doublereal *, doublereal *,	     doublereal *, integer *, integer *, integer *, integer *, 	    doublereal *, doublereal *, doublereal *, doublereal *, integer *,	     integer *);    static integer bx;    static doublereal sn;    extern /* Subroutine */ int dlascl_(char *, integer *, integer *, 	    doublereal *, doublereal *, integer *, integer *, doublereal *, 	    integer *, integer *);    extern integer idamax_(integer *, doublereal *, integer *);    static integer st;    extern /* Subroutine */ int dlasdq_(char *, integer *, integer *, integer 	    *, integer *, integer *, doublereal *, doublereal *, doublereal *,	     integer *, doublereal *, integer *, doublereal *, integer *, 	    doublereal *, integer *);    static integer vt;    extern /* Subroutine */ int dlaset_(char *, integer *, integer *, 	    doublereal *, doublereal *, doublereal *, integer *), 	    dlartg_(doublereal *, doublereal *, doublereal *, doublereal *, 	    doublereal *), xerbla_(char *, integer *);    static integer givcol;    extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *);    extern /* Subroutine */ int zlalsa_(integer *, integer *, integer *, 	    integer *, doublecomplex *, integer *, doublecomplex *, integer *,	     doublereal *, integer *, doublereal *, integer *, doublereal *, 	    doublereal *, doublereal *, doublereal *, integer *, integer *, 	    integer *, integer *, doublereal *, doublereal *, doublereal *, 	    doublereal *, integer *, integer *), zlascl_(char *, integer *, 	    integer *, doublereal *, doublereal *, integer *, integer *, 	    doublecomplex *, integer *, integer *), dlasrt_(char *, 	    integer *, doublereal *, integer *), zlacpy_(char *, 	    integer *, integer *, doublecomplex *, integer *, doublecomplex *,	     integer *), zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);    static doublereal orgnrm;    static integer givnum, givptr, nm1, nrwork, irwwrk, smlszp, st1;    static doublereal eps;    static integer iwk;    static doublereal tol;#define b_subscr(a_1,a_2) (a_2)*b_dim1 + a_1#define b_ref(a_1,a_2) b[b_subscr(a_1,a_2)]/*  -- LAPACK routine (version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          October 31, 1999       Purpose       =======       ZLALSD uses the singular value decomposition of A to solve the least       squares problem of finding X to minimize the Euclidean norm of each       column of A*X-B, where A is N-by-N upper bidiagonal, and X and B       are N-by-NRHS. The solution X overwrites B.       The singular values of A smaller than RCOND times the largest       singular value are treated as zero in solving the least squares       problem; in this case a minimum norm solution is returned.       The actual singular values are returned in D in ascending order.       This code makes very mild assumptions about floating point       arithmetic. It will work on machines with a guard digit in       add/subtract, or on those binary machines without guard digits       which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.   //.........这里部分代码省略.........

/* Subroutine */ int zggev_(char *jobvl, char *jobvr, integer *n, 	doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, 	doublecomplex *alpha, doublecomplex *beta, doublecomplex *vl, integer 	*ldvl, doublecomplex *vr, integer *ldvr, doublecomplex *work, integer 	*lwork, doublereal *rwork, integer *info){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1, 	    vr_offset, i__1, i__2, i__3, i__4;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1;    /* Builtin functions */    double sqrt(doublereal), d_imag(doublecomplex *);    /* Local variables */    integer jc, in, jr, ihi, ilo;    doublereal eps;    logical ilv;    doublereal anrm, bnrm;    integer ierr, itau;    doublereal temp;    logical ilvl, ilvr;    integer iwrk;    extern logical lsame_(char *, char *);    integer ileft, icols, irwrk, irows;    extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);    extern doublereal dlamch_(char *);    extern /* Subroutine */ int zggbak_(char *, char *, integer *, integer *, 	    integer *, doublereal *, doublereal *, integer *, doublecomplex *, 	     integer *, integer *), zggbal_(char *, integer *, 	     doublecomplex *, integer *, doublecomplex *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, integer *);    logical ilascl, ilbscl;    extern /* Subroutine */ int xerbla_(char *, integer *);    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 	    integer *, integer *);    logical ldumma[1];    char chtemp[1];    doublereal bignum;    extern doublereal zlange_(char *, integer *, integer *, doublecomplex *, 	    integer *, doublereal *);    integer ijobvl, iright;    extern /* Subroutine */ int zgghrd_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, integer *, doublecomplex *, integer *, 	     doublecomplex *, integer *, doublecomplex *, integer *, integer *), zlascl_(char *, integer *, integer *, 	    doublereal *, doublereal *, integer *, integer *, doublecomplex *, 	     integer *, integer *);    integer ijobvr;    extern /* Subroutine */ int zgeqrf_(integer *, integer *, doublecomplex *, 	     integer *, doublecomplex *, doublecomplex *, integer *, integer *);    doublereal anrmto;    integer lwkmin;    doublereal bnrmto;    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *), 	    zlaset_(char *, integer *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *, integer *), ztgevc_(	    char *, char *, logical *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, integer *, integer *, doublecomplex *, 	     doublereal *, integer *), zhgeqz_(char *, char *, 	     char *, integer *, integer *, integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *, doublereal *, integer *);    doublereal smlnum;    integer lwkopt;    logical lquery;    extern /* Subroutine */ int zungqr_(integer *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *, integer *), zunmqr_(char *, char *, integer *, integer 	    *, integer *, doublecomplex *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, integer *);/*  -- LAPACK driver routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZGGEV computes for a pair of N-by-N complex nonsymmetric matrices *//*  (A,B), the generalized eigenvalues, and optionally, the left and/or *//*  right generalized eigenvectors. *//*  A generalized eigenvalue for a pair of matrices (A,B) is a scalar *//*  lambda or a ratio alpha/beta = lambda, such that A - lambda*B is *//*  singular. It is usually represented as the pair (alpha,beta), as *//*  there is a reasonable interpretation for beta=0, and even for both *//*  being zero. *///.........这里部分代码省略.........

 int zgges_(char *jobvsl, char *jobvsr, char *sort, L_fp 	selctg, int *n, doublecomplex *a, int *lda, doublecomplex *b, 	int *ldb, int *sdim, doublecomplex *alpha, doublecomplex *	beta, doublecomplex *vsl, int *ldvsl, doublecomplex *vsr, int 	*ldvsr, doublecomplex *work, int *lwork, double *rwork, 	int *bwork, int *info){    /* System generated locals */    int a_dim1, a_offset, b_dim1, b_offset, vsl_dim1, vsl_offset, 	    vsr_dim1, vsr_offset, i__1, i__2;    /* Builtin functions */    double sqrt(double);    /* Local variables */    int i__;    double dif[2];    int ihi, ilo;    double eps, anrm, bnrm;    int idum[1], ierr, itau, iwrk;    double pvsl, pvsr;    extern int lsame_(char *, char *);    int ileft, icols;    int cursl, ilvsl, ilvsr;    int irwrk, irows;    extern  int dlabad_(double *, double *);    extern double dlamch_(char *);    extern  int zggbak_(char *, char *, int *, int *, 	    int *, double *, double *, int *, doublecomplex *, 	     int *, int *), zggbal_(char *, int *, 	     doublecomplex *, int *, doublecomplex *, int *, int *, int *, double *, double *, double *, int *);    int ilascl, ilbscl;    extern  int xerbla_(char *, int *);    extern int ilaenv_(int *, char *, char *, int *, int *, 	    int *, int *);    extern double zlange_(char *, int *, int *, doublecomplex *, 	    int *, double *);    double bignum;    int ijobvl, iright;    extern  int zgghrd_(char *, char *, int *, int *, 	    int *, doublecomplex *, int *, doublecomplex *, int *, 	     doublecomplex *, int *, doublecomplex *, int *, int *), zlascl_(char *, int *, int *, 	    double *, double *, int *, int *, doublecomplex *, 	     int *, int *);    int ijobvr;    extern  int zgeqrf_(int *, int *, doublecomplex *, 	     int *, doublecomplex *, doublecomplex *, int *, int *);    double anrmto;    int lwkmin;    int lastsl;    double bnrmto;    extern  int zlacpy_(char *, int *, int *, 	    doublecomplex *, int *, doublecomplex *, int *), 	    zlaset_(char *, int *, int *, doublecomplex *, 	    doublecomplex *, doublecomplex *, int *), zhgeqz_(	    char *, char *, char *, int *, int *, int *, 	    doublecomplex *, int *, doublecomplex *, int *, 	    doublecomplex *, doublecomplex *, doublecomplex *, int *, 	    doublecomplex *, int *, doublecomplex *, int *, 	    double *, int *), ztgsen_(int 	    *, int *, int *, int *, int *, doublecomplex *, 	    int *, doublecomplex *, int *, doublecomplex *, 	    doublecomplex *, doublecomplex *, int *, doublecomplex *, 	    int *, int *, double *, double *, double *, 	    doublecomplex *, int *, int *, int *, int *);    double smlnum;    int wantst, lquery;    int lwkopt;    extern  int zungqr_(int *, int *, int *, 	    doublecomplex *, int *, doublecomplex *, doublecomplex *, 	    int *, int *), zunmqr_(char *, char *, int *, int 	    *, int *, doublecomplex *, int *, doublecomplex *, 	    doublecomplex *, int *, doublecomplex *, int *, int *);/*  -- LAPACK driver routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*     .. Function Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZGGES computes for a pair of N-by-N complex nonsymmetric matrices *//*  (A,B), the generalized eigenvalues, the generalized complex Schur *//*  form (S, T), and optionally left and/or right Schur vectors (VSL *//*  and VSR). This gives the generalized Schur factorization *//*          (A,B) = ( (VSL)*S*(VSR)**H, (VSL)*T*(VSR)**H ) *//*  where (VSR)**H is the conjugate-transpose of VSR. *///.........这里部分代码省略.........

/* Subroutine */ int zchkhs_(integer *nsizes, integer *nn, integer *ntypes, 	logical *dotype, integer *iseed, doublereal *thresh, integer *nounit, 	doublecomplex *a, integer *lda, doublecomplex *h__, doublecomplex *t1, 	 doublecomplex *t2, doublecomplex *u, integer *ldu, doublecomplex *	z__, doublecomplex *uz, doublecomplex *w1, doublecomplex *w3, 	doublecomplex *evectl, doublecomplex *evectr, doublecomplex *evecty, 	doublecomplex *evectx, doublecomplex *uu, doublecomplex *tau, 	doublecomplex *work, integer *nwork, doublereal *rwork, integer *	iwork, logical *select, doublereal *result, integer *info){    /* Initialized data */    static integer ktype[21] = { 1,2,3,4,4,4,4,4,6,6,6,6,6,6,6,6,6,6,9,9,9 };    static integer kmagn[21] = { 1,1,1,1,1,1,2,3,1,1,1,1,1,1,1,1,2,3,1,2,3 };    static integer kmode[21] = { 0,0,0,4,3,1,4,4,4,3,1,5,4,3,1,5,5,5,4,3,1 };    static integer kconds[21] = { 0,0,0,0,0,0,0,0,1,1,1,1,2,2,2,2,2,2,0,0,0 };    /* Format strings */    static char fmt_9999[] = "(/002 ZCHKHS: /002,a,/002 returned INFO=/002,i"	    "6,/002./002,/9x,/002N=/002,i6,/002, JTYPE=/002,i6,/002, ISEED="	    "(/002,3(i5,/002,/002),i5,/002)/002)";    static char fmt_9998[] = "(/002 ZCHKHS: /002,a,/002 Eigenvectors from"	    " /002,a,/002 incorrectly /002,/002normalized./002,//002 Bits of "	    "error=/002,0p,g10.3,/002,/002,9x,/002N=/002,i6,/002, JTYPE=/002,"	    "i6,/002, ISEED=(/002,3(i5,/002,/002),i5,/002)/002)";    static char fmt_9997[] = "(/002 ZCHKHS: Selected /002,a,/002 Eigenvector"	    "s from /002,a,/002 do not match other eigenvectors /002,9x,/002N="	    "/002,i6,/002, JTYPE=/002,i6,/002, ISEED=(/002,3(i5,/002,/002),i5,"	    "/002)/002)";    /* System generated locals */    integer a_dim1, a_offset, evectl_dim1, evectl_offset, evectr_dim1, 	    evectr_offset, evectx_dim1, evectx_offset, evecty_dim1, 	    evecty_offset, h_dim1, h_offset, t1_dim1, t1_offset, t2_dim1, 	    t2_offset, u_dim1, u_offset, uu_dim1, uu_offset, uz_dim1, 	    uz_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5, i__6;    doublereal d__1, d__2;    doublecomplex z__1;    /* Builtin functions */    double sqrt(doublereal);    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);    double z_abs(doublecomplex *);    /* Local variables */    integer i__, j, k, n, n1, jj, in, ihi, ilo;    doublereal ulp, cond;    integer jcol, nmax;    doublereal unfl, ovfl, temp1, temp2;    logical badnn, match;    integer imode;    doublereal dumma[4];    integer iinfo;    doublereal conds;    extern /* Subroutine */ int zget10_(integer *, integer *, doublecomplex *, 	     integer *, doublecomplex *, integer *, doublecomplex *, 	    doublereal *, doublereal *);    doublereal aninv, anorm;    extern /* Subroutine */ int zget22_(char *, char *, char *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, doublecomplex *, doublereal *, doublereal *), zgemm_(char *, char *, integer *, 	    integer *, integer *, doublecomplex *, doublecomplex *, integer *, 	     doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *);    integer nmats, jsize, nerrs, itype, jtype, ntest;    extern /* Subroutine */ int zhst01_(integer *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublereal *, doublereal *), zcopy_(integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *);    doublereal rtulp;    extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);    extern doublereal dlamch_(char *);    doublecomplex cdumma[4];    integer idumma[1];    extern /* Subroutine */ int dlafts_(char *, integer *, integer *, integer 	    *, integer *, doublereal *, integer *, doublereal *, integer *, 	    integer *);    integer ioldsd[4];    extern /* Subroutine */ int xerbla_(char *, integer *), zgehrd_(	    integer *, integer *, integer *, doublecomplex *, integer *, 	    doublecomplex *, doublecomplex *, integer *, integer *), dlasum_(	    char *, integer *, integer *, integer *), zlatme_(integer 	    *, char *, integer *, doublecomplex *, integer *, doublereal *, 	    doublecomplex *, char *, char *, char *, char *, doublereal *, 	    integer *, doublereal *, integer *, integer *, doublereal *, 	    doublecomplex *, integer *, doublecomplex *, integer *), zhsein_(char *, char *, char *, 	    logical *, integer *, doublecomplex *, integer *, doublecomplex *, 	     doublecomplex *, integer *, doublecomplex *, integer *, integer *, integer *, doublecomplex *, doublereal *, integer *, integer *, 	    integer *), zlacpy_(char *, integer *, 	    integer *, doublecomplex *, integer *, doublecomplex *, integer *), zlaset_(char *, integer *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *, integer *), zlatmr_(	    integer *, integer *, char *, integer *, char *, doublecomplex *, 	    integer *, doublereal *, doublecomplex *, char *, char *, 	    doublecomplex *, integer *, doublereal *, doublecomplex *, 	    integer *, doublereal *, char *, integer *, integer *, integer *, 	    doublereal *, doublereal *, char *, doublecomplex *, integer *, 	    integer *, integer *);    doublereal rtunfl, rtovfl, rtulpi, ulpinv;//.........这里部分代码省略.........

/* Subroutine */ int zchktz_(logical *dotype, integer *nm, integer *mval, 	integer *nn, integer *nval, doublereal *thresh, logical *tsterr, 	doublecomplex *a, doublecomplex *copya, doublereal *s, doublereal *	copys, doublecomplex *tau, doublecomplex *work, doublereal *rwork, 	integer *nout){    /* Initialized data */    static integer iseedy[4] = { 1988,1989,1990,1991 };    /* Format strings */    static char fmt_9999[] = "(/002 M =/002,i5,/002, N =/002,i5,/002, type"	    " /002,i2,/002, test /002,i2,/002, ratio =/002,g12.5)";    /* System generated locals */    integer i__1, i__2, i__3, i__4;    doublereal d__1;    /* Builtin functions          Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);    /* Local variables */    static integer mode, info;    static char path[3];    static integer nrun, i__;    extern /* Subroutine */ int alahd_(integer *, char *);    static integer k, m, n, nfail, iseed[4], imode, mnmin, nerrs, lwork;    extern doublereal zqrt12_(integer *, integer *, doublecomplex *, integer *	    , doublereal *, doublecomplex *, integer *, doublereal *), 	    zrzt01_(integer *, integer *, doublecomplex *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *), zrzt02_(	    integer *, integer *, doublecomplex *, integer *, doublecomplex *,	     doublecomplex *, integer *), ztzt01_(integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    doublecomplex *, integer *), ztzt02_(integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *);    extern /* Subroutine */ int zgeqr2_(integer *, integer *, doublecomplex *,	     integer *, doublecomplex *, doublecomplex *, integer *);    static integer im, in;    extern doublereal dlamch_(char *);    extern /* Subroutine */ int dlaord_(char *, integer *, doublereal *, 	    integer *), alasum_(char *, integer *, integer *, integer 	    *, integer *), zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *), 	    zlaset_(char *, integer *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *, integer *), zlatms_(	    integer *, integer *, char *, integer *, char *, doublereal *, 	    integer *, doublereal *, doublereal *, integer *, integer *, char 	    *, doublecomplex *, integer *, doublecomplex *, integer *);    static doublereal result[6];    extern /* Subroutine */ int zerrtz_(char *, integer *), ztzrqf_(	    integer *, integer *, doublecomplex *, integer *, doublecomplex *,	     integer *), ztzrzf_(integer *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *, integer *)	    ;    static integer lda;    static doublereal eps;    /* Fortran I/O blocks */    static cilist io___21 = { 0, 0, 0, fmt_9999, 0 };/*  -- LAPACK test routine (version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          June 30, 1999       Purpose       =======       ZCHKTZ tests ZTZRQF and ZTZRZF.       Arguments       =========       DOTYPE  (input) LOGICAL array, dimension (NTYPES)               The matrix types to be used for testing.  Matrices of type j               (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =               .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.       NM      (input) INTEGER               The number of values of M contained in the vector MVAL.       MVAL    (input) INTEGER array, dimension (NM)               The values of the matrix row dimension M.       NN      (input) INTEGER               The number of values of N contained in the vector NVAL.       NVAL    (input) INTEGER array, dimension (NN)               The values of the matrix column dimension N.       THRESH  (input) DOUBLE PRECISION               The threshold value for the test ratios.  A result is               included in the output file if RESULT >= THRESH.  To have               every test ratio printed, use THRESH = 0.   //.........这里部分代码省略.........

//.........这里部分代码省略........./*                       pass the threshold. */			    i__4 = nt;			    for (k = 1; k <= i__4; ++k) {				if (result[k - 1] >= *thresh) {				    if (nfail == 0 && nerrs == 0) {					aladhd_(nout, path);				    }				    io___48.ciunit = *nout;				    s_wsfe(&io___48);				    do_fio(&c__1, "ZPOSV ", (ftnlen)6);				    do_fio(&c__1, uplo, (ftnlen)1);				    do_fio(&c__1, (char *)&n, (ftnlen)sizeof(					    integer));				    do_fio(&c__1, (char *)&imat, (ftnlen)					    sizeof(integer));				    do_fio(&c__1, (char *)&k, (ftnlen)sizeof(					    integer));				    do_fio(&c__1, (char *)&result[k - 1], (					    ftnlen)sizeof(doublereal));				    e_wsfe();				    ++nfail;				}/* L60: */			    }			    nrun += nt;L70:			    ;			}/*                    --- Test ZPOSVX --- */			if (! prefac) {			    zlaset_(uplo, &n, &n, &c_b51, &c_b51, &afac[1], &				    lda);			}			zlaset_("Full", &n, nrhs, &c_b51, &c_b51, &x[1], &lda);			if (iequed > 1 && n > 0) {/*                       Equilibrate the matrix if FACT='F' and *//*                       EQUED='Y'. */			    zlaqhe_(uplo, &n, &a[1], &lda, &s[1], &scond, &				    amax, equed);			}/*                    Solve the system and compute the condition number *//*                    and error bounds using ZPOSVX. */			s_copy(srnamc_1.srnamt, "ZPOSVX", (ftnlen)32, (ftnlen)				6);			zposvx_(fact, uplo, &n, nrhs, &a[1], &lda, &afac[1], &				lda, equed, &s[1], &b[1], &lda, &x[1], &lda, &				rcond, &rwork[1], &rwork[*nrhs + 1], &work[1], 				 &rwork[(*nrhs << 1) + 1], &info);/*                    Check the error code from ZPOSVX. */			if (info != izero) {/* Writing concatenation */			    i__5[0] = 1, a__1[0] = fact;			    i__5[1] = 1, a__1[1] = uplo;			    s_cat(ch__1, a__1, i__5, &c__2, (ftnlen)2);			    alaerh_(path, "ZPOSVX", &info, &izero, ch__1, &n, 				    &n, &c_n1, &c_n1, nrhs, &imat, &nfail, &				    nerrs, nout);

/* Subroutine */ int zgghrd_(char *compq, char *compz, integer *n, integer *	ilo, integer *ihi, doublecomplex *a, integer *lda, doublecomplex *b, 	integer *ldb, doublecomplex *q, integer *ldq, doublecomplex *z__, 	integer *ldz, integer *info){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, 	    z_offset, i__1, i__2, i__3;    doublecomplex z__1;    /* Builtin functions */    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    doublereal c__;    doublecomplex s;    logical ilq, ilz;    integer jcol, jrow;    extern /* Subroutine */ int zrot_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublecomplex *);    extern logical lsame_(char *, char *);    doublecomplex ctemp;    extern /* Subroutine */ int xerbla_(char *, integer *);    integer icompq, icompz;    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *), zlartg_(doublecomplex *, doublecomplex *, doublereal *, 	    doublecomplex *, doublecomplex *);/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZGGHRD reduces a pair of complex matrices (A,B) to generalized upper *//*  Hessenberg form using unitary transformations, where A is a *//*  general matrix and B is upper triangular.  The form of the *//*  generalized eigenvalue problem is *//*     A*x = lambda*B*x, *//*  and B is typically made upper triangular by computing its QR *//*  factorization and moving the unitary matrix Q to the left side *//*  of the equation. *//*  This subroutine simultaneously reduces A to a Hessenberg matrix H: *//*     Q**H*A*Z = H *//*  and transforms B to another upper triangular matrix T: *//*     Q**H*B*Z = T *//*  in order to reduce the problem to its standard form *//*     H*y = lambda*T*y *//*  where y = Z**H*x. *//*  The unitary matrices Q and Z are determined as products of Givens *//*  rotations.  They may either be formed explicitly, or they may be *//*  postmultiplied into input matrices Q1 and Z1, so that *//*       Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H *//*       Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H *//*  If Q1 is the unitary matrix from the QR factorization of B in the *//*  original equation A*x = lambda*B*x, then ZGGHRD reduces the original *//*  problem to generalized Hessenberg form. *//*  Arguments *//*  ========= *//*  COMPQ   (input) CHARACTER*1 *//*          = 'N': do not compute Q; *//*          = 'I': Q is initialized to the unit matrix, and the *//*                 unitary matrix Q is returned; *//*          = 'V': Q must contain a unitary matrix Q1 on entry, *//*                 and the product Q1*Q is returned. *//*  COMPZ   (input) CHARACTER*1 *//*          = 'N': do not compute Q; *//*          = 'I': Q is initialized to the unit matrix, and the *//*                 unitary matrix Q is returned; *//*          = 'V': Q must contain a unitary matrix Q1 on entry, *//*                 and the product Q1*Q is returned. *//*  N       (input) INTEGER *//*          The order of the matrices A and B.  N >= 0. *//*  ILO     (input) INTEGER *//*  IHI     (input) INTEGER *//*          ILO and IHI mark the rows and columns of A which are to be *//*          reduced.  It is assumed that A 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 ZGGBAL; otherwise they *//*          should be set to 1 and N respectively. *//*          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA, N) *//*          On entry, the N-by-N general matrix to be reduced. *//*          On exit, the upper triangle and the first subdiagonal of A *//*          are overwritten with the upper Hessenberg matrix H, and the *///.........这里部分代码省略.........

//.........这里部分代码省略.........	rwork[1] = *bbnrm;	dlascl_("G", &c__0, &c__0, &bnrmto, &bnrm, &c__1, &c__1, &rwork[1], &		c__1, &ierr);	*bbnrm = rwork[1];    }/*     Reduce B to triangular form (QR decomposition of B) *//*     (Complex Workspace: need N, prefer N*NB ) */    irows = *ihi + 1 - *ilo;    if (ilv || ! wantsn) {	icols = *n + 1 - *ilo;    } else {	icols = irows;    }    itau = 1;    iwrk = itau + irows;    i__1 = *lwork + 1 - iwrk;    zgeqrf_(&irows, &icols, &b[*ilo + *ilo * b_dim1], ldb, &work[itau], &work[	    iwrk], &i__1, &ierr);/*     Apply the unitary transformation to A *//*     (Complex Workspace: need N, prefer N*NB) */    i__1 = *lwork + 1 - iwrk;    zunmqr_("L", "C", &irows, &icols, &irows, &b[*ilo + *ilo * b_dim1], ldb, &	    work[itau], &a[*ilo + *ilo * a_dim1], lda, &work[iwrk], &i__1, &	    ierr);/*     Initialize VL and/or VR *//*     (Workspace: need N, prefer N*NB) */    if (ilvl) {	zlaset_("Full", n, n, &c_b1, &c_b2, &vl[vl_offset], ldvl);	if (irows > 1) {	    i__1 = irows - 1;	    i__2 = irows - 1;	    zlacpy_("L", &i__1, &i__2, &b[*ilo + 1 + *ilo * b_dim1], ldb, &vl[		    *ilo + 1 + *ilo * vl_dim1], ldvl);	}	i__1 = *lwork + 1 - iwrk;	zungqr_(&irows, &irows, &irows, &vl[*ilo + *ilo * vl_dim1], ldvl, &		work[itau], &work[iwrk], &i__1, &ierr);    }    if (ilvr) {	zlaset_("Full", n, n, &c_b1, &c_b2, &vr[vr_offset], ldvr);    }/*     Reduce to generalized Hessenberg form *//*     (Workspace: none needed) */    if (ilv || ! wantsn) {/*        Eigenvectors requested -- work on whole matrix. */	zgghrd_(jobvl, jobvr, n, ilo, ihi, &a[a_offset], lda, &b[b_offset], 		ldb, &vl[vl_offset], ldvl, &vr[vr_offset], ldvr, &ierr);    } else {	zgghrd_("N", "N", &irows, &c__1, &irows, &a[*ilo + *ilo * a_dim1], 		lda, &b[*ilo + *ilo * b_dim1], ldb, &vl[vl_offset], ldvl, &vr[		vr_offset], ldvr, &ierr);    }/*     Perform QZ algorithm (Compute eigenvalues, and optionally, the *//*     Schur forms and Schur vectors) */

/* Subroutine */ int zspt01_(char *uplo, integer *n, doublecomplex *a, 	doublecomplex *afac, integer *ipiv, doublecomplex *c__, integer *ldc, 	doublereal *rwork, doublereal *resid){    /* System generated locals */    integer c_dim1, c_offset, i__1, i__2, i__3, i__4, i__5;    doublecomplex z__1;    /* Local variables */    static integer info, i__, j;    extern logical lsame_(char *, char *);    static doublereal anorm;    static integer jc;    extern doublereal dlamch_(char *);    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);    extern doublereal zlansp_(char *, char *, integer *, doublecomplex *, 	    doublereal *);    extern /* Subroutine */ int zlavsp_(char *, char *, char *, integer *, 	    integer *, doublecomplex *, integer *, doublecomplex *, integer *,	     integer *);    extern doublereal zlansy_(char *, char *, integer *, doublecomplex *, 	    integer *, doublereal *);    static doublereal eps;#define c___subscr(a_1,a_2) (a_2)*c_dim1 + a_1#define c___ref(a_1,a_2) c__[c___subscr(a_1,a_2)]/*  -- LAPACK test routine (version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          September 30, 1994       Purpose       =======       ZSPT01 reconstructs a symmetric indefinite packed matrix A from its       diagonal pivoting factorization A = U*D*U' or A = L*D*L' and computes       the residual          norm( C - A ) / ( N * norm(A) * EPS ),       where C is the reconstructed matrix and EPS is the machine epsilon.       Arguments       ==========       UPLO    (input) CHARACTER*1               Specifies whether the upper or lower triangular part of the               Hermitian matrix A is stored:               = 'U':  Upper triangular               = 'L':  Lower triangular       N       (input) INTEGER               The order of the matrix A.  N >= 0.       A       (input) COMPLEX*16 array, dimension (N*(N+1)/2)               The original symmetric matrix A, stored as a packed               triangular matrix.       AFAC    (input) COMPLEX*16 array, dimension (N*(N+1)/2)               The factored form of the matrix A, stored as a packed               triangular matrix.  AFAC contains the block diagonal matrix D               and the multipliers used to obtain the factor L or U from the               L*D*L' or U*D*U' factorization as computed by ZSPTRF.       IPIV    (input) INTEGER array, dimension (N)               The pivot indices from ZSPTRF.       C       (workspace) COMPLEX*16 array, dimension (LDC,N)       LDC     (integer) INTEGER               The leading dimension of the array C.  LDC >= max(1,N).       RWORK   (workspace) DOUBLE PRECISION array, dimension (N)       RESID   (output) DOUBLE PRECISION               If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )               If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )       =====================================================================          Quick exit if N = 0.          Parameter adjustments */    --a;    --afac;    --ipiv;    c_dim1 = *ldc;    c_offset = 1 + c_dim1 * 1;    c__ -= c_offset;    --rwork;    /* Function Body */    if (*n <= 0) {	*resid = 0.;	return 0;    }//.........这里部分代码省略.........

/* Subroutine */ int zneigh_(doublereal *rnorm, integer *n, doublecomplex *	h__, integer *ldh, doublecomplex *ritz, doublecomplex *bounds, 	doublecomplex *q, integer *ldq, doublecomplex *workl, doublereal *	rwork, integer *ierr){    /* System generated locals */    integer h_dim1, h_offset, q_dim1, q_offset, i__1;    doublereal d__1;    /* Local variables */    static integer j;    static real t0, t1;    static doublecomplex vl[1];    static doublereal temp;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zmout_(integer *, integer *, integer 	    *, doublecomplex *, integer *, integer *, char *, ftnlen), zvout_(	    integer *, integer *, doublecomplex *, integer *, char *, ftnlen);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);    extern /* Subroutine */ int second_(real *);    static logical select[1];    static integer msglvl;    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, ftnlen), 	    zlahqr_(logical *, logical *, integer *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, integer *,	     doublecomplex *, integer *, integer *), ztrevc_(char *, char *, 	    logical *, integer *, doublecomplex *, integer *, doublecomplex *,	     integer *, doublecomplex *, integer *, integer *, integer *, 	    doublecomplex *, doublereal *, integer *, ftnlen, ftnlen), 	    zdscal_(integer *, doublereal *, doublecomplex *, integer *), 	    zlaset_(char *, integer *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *, integer *, ftnlen);/*     %----------------------------------------------------% *//*     | Include files for debugging and timing information | *//*     %----------------------------------------------------% *//* /SCCS Information: @(#) *//* FILE: debug.h   SID: 2.3   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %---------------------------------% *//*     | See debug.doc for documentation | *//*     %---------------------------------% *//*     %------------------% *//*     | Scalar Arguments | *//*     %------------------% *//*     %--------------------------------% *//*     | See stat.doc for documentation | *//*     %--------------------------------% *//* /SCCS Information: @(#) *//* FILE: stat.h   SID: 2.2   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %-----------------% *//*     | Array Arguments | *//*     %-----------------% *//*     %------------% *//*     | Parameters | *//*     %------------% *//*     %------------------------% *//*     | Local Scalars & Arrays | *//*     %------------------------% *//*     %----------------------% *//*     | External Subroutines | *//*     %----------------------% *//*     %--------------------% *//*     | External Functions | *//*     %--------------------% *//*     %-----------------------% *//*     | Executable Statements | *//*     %-----------------------% *//*     %-------------------------------% *//*     | Initialize timing statistics  | *//*     | & message level for debugging | *//*     %-------------------------------% */    /* Parameter adjustments */    --rwork;    --workl;    --bounds;    --ritz;    h_dim1 = *ldh;//.........这里部分代码省略.........

/* Subroutine */ int zpteqr_(char *compz, integer *n, doublereal *d__, 	doublereal *e, doublecomplex *z__, integer *ldz, doublereal *work, 	integer *info, ftnlen compz_len){    /* System generated locals */    integer z_dim1, z_offset, i__1;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    static doublecomplex c__[1]	/* was [1][1] */;    static integer i__;    static doublecomplex vt[1]	/* was [1][1] */;    static integer nru;    extern logical lsame_(char *, char *, ftnlen, ftnlen);    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);    static integer icompz;    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *, 	    ftnlen), dpttrf_(integer *, doublereal *, doublereal *, integer *)	    , zbdsqr_(char *, integer *, integer *, integer *, integer *, 	    doublereal *, doublereal *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublereal *, integer *, ftnlen);/*  -- LAPACK routine (version 3.0) -- *//*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., *//*     Courant Institute, Argonne National Lab, and Rice University *//*     October 31, 1999 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZPTEQR computes all eigenvalues and, optionally, eigenvectors of a *//*  symmetric positive definite tridiagonal matrix by first factoring the *//*  matrix using DPTTRF and then calling ZBDSQR to compute the singular *//*  values of the bidiagonal factor. *//*  This routine computes the eigenvalues of the positive definite *//*  tridiagonal matrix to high relative accuracy.  This means that if the *//*  eigenvalues range over many orders of magnitude in size, then the *//*  small eigenvalues and corresponding eigenvectors will be computed *//*  more accurately than, for example, with the standard QR method. *//*  The eigenvectors of a full or band positive definite Hermitian matrix *//*  can also be found if ZHETRD, ZHPTRD, or ZHBTRD has been used to *//*  reduce this matrix to tridiagonal form.  (The reduction to *//*  tridiagonal form, however, may preclude the possibility of obtaining *//*  high relative accuracy in the small eigenvalues of the original *//*  matrix, if these eigenvalues range over many orders of magnitude.) *//*  Arguments *//*  ========= *//*  COMPZ   (input) CHARACTER*1 *//*          = 'N':  Compute eigenvalues only. *//*          = 'V':  Compute eigenvectors of original Hermitian *//*                  matrix also.  Array Z contains the unitary matrix *//*                  used to reduce the original matrix to tridiagonal *//*                  form. *//*          = 'I':  Compute eigenvectors of tridiagonal matrix also. *//*  N       (input) INTEGER *//*          The order of the matrix.  N >= 0. *//*  D       (input/output) DOUBLE PRECISION array, dimension (N) *//*          On entry, the n diagonal elements of the tridiagonal matrix. *//*          On normal exit, D contains the eigenvalues, in descending *//*          order. *//*  E       (input/output) DOUBLE PRECISION array, dimension (N-1) *//*          On entry, the (n-1) subdiagonal elements of the tridiagonal *//*          matrix. *//*          On exit, E has been destroyed. *//*  Z       (input/output) COMPLEX*16 array, dimension (LDZ, N) *//*          On entry, if COMPZ = 'V', the unitary matrix used in the *//*          reduction to tridiagonal form. *//*          On exit, if COMPZ = 'V', the orthonormal eigenvectors of the *//*          original Hermitian matrix; *//*          if COMPZ = 'I', the orthonormal eigenvectors of the *//*          tridiagonal matrix. *//*          If INFO > 0 on exit, Z contains the eigenvectors associated *//*          with only the stored eigenvalues. *//*          If  COMPZ = 'N', then Z is not referenced. *//*  LDZ     (input) INTEGER *//*          The leading dimension of the array Z.  LDZ >= 1, and if *//*          COMPZ = 'V' or 'I', LDZ >= max(1,N). *//*  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N) *//*  INFO    (output) INTEGER *///.........这里部分代码省略.........

/* Subroutine */ int zsteqr_(char *compz, integer *n, doublereal *d__, 	doublereal *e, doublecomplex *z__, integer *ldz, doublereal *work, 	integer *info){    /* System generated locals */    integer z_dim1, z_offset, i__1, i__2;    doublereal d__1, d__2;    /* Builtin functions */    double sqrt(doublereal), d_sign(doublereal *, doublereal *);    /* Local variables */    doublereal b, c__, f, g;    integer i__, j, k, l, m;    doublereal p, r__, s;    integer l1, ii, mm, lm1, mm1, nm1;    doublereal rt1, rt2, eps;    integer lsv;    doublereal tst, eps2;    integer lend, jtot;    extern /* Subroutine */ int dlae2_(doublereal *, doublereal *, doublereal 	    *, doublereal *, doublereal *);    extern logical lsame_(char *, char *);    doublereal anorm;    extern /* Subroutine */ int zlasr_(char *, char *, char *, integer *, 	    integer *, doublereal *, doublereal *, doublecomplex *, integer *), zswap_(integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *), dlaev2_(doublereal *, 	    doublereal *, doublereal *, doublereal *, doublereal *, 	    doublereal *, doublereal *);    integer lendm1, lendp1;    extern doublereal dlapy2_(doublereal *, doublereal *), dlamch_(char *);    integer iscale;    extern /* Subroutine */ int dlascl_(char *, integer *, integer *, 	    doublereal *, doublereal *, integer *, integer *, doublereal *, 	    integer *, integer *);    doublereal safmin;    extern /* Subroutine */ int dlartg_(doublereal *, doublereal *, 	    doublereal *, doublereal *, doublereal *);    doublereal safmax;    extern /* Subroutine */ int xerbla_(char *, integer *);    extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *);    extern /* Subroutine */ int dlasrt_(char *, integer *, doublereal *, 	    integer *);    integer lendsv;    doublereal ssfmin;    integer nmaxit, icompz;    doublereal ssfmax;    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZSTEQR computes all eigenvalues and, optionally, eigenvectors of a *//*  symmetric tridiagonal matrix using the implicit QL or QR method. *//*  The eigenvectors of a full or band complex Hermitian matrix can also *//*  be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this *//*  matrix to tridiagonal form. *//*  Arguments *//*  ========= *//*  COMPZ   (input) CHARACTER*1 *//*          = 'N':  Compute eigenvalues only. *//*          = 'V':  Compute eigenvalues and eigenvectors of the original *//*                  Hermitian matrix.  On entry, Z must contain the *//*                  unitary matrix used to reduce the original matrix *//*                  to tridiagonal form. *//*          = 'I':  Compute eigenvalues and eigenvectors of the *//*                  tridiagonal matrix.  Z is initialized to the identity *//*                  matrix. *//*  N       (input) INTEGER *//*          The order of the matrix.  N >= 0. *//*  D       (input/output) DOUBLE PRECISION array, dimension (N) *//*          On entry, the diagonal elements of the tridiagonal matrix. *//*          On exit, if INFO = 0, the eigenvalues in ascending order. *//*  E       (input/output) DOUBLE PRECISION array, dimension (N-1) *//*          On entry, the (n-1) subdiagonal elements of the tridiagonal *//*          matrix. *//*          On exit, E has been destroyed. *//*  Z       (input/output) COMPLEX*16 array, dimension (LDZ, N) *//*          On entry, if  COMPZ = 'V', then Z contains the unitary *//*          matrix used in the reduction to tridiagonal form. *//*          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the *//*          orthonormal eigenvectors of the original Hermitian matrix, *//*          and if COMPZ = 'I', Z contains the orthonormal eigenvectors *///.........这里部分代码省略.........

/* Subroutine */ int zlatm4_(integer *itype, integer *n, integer *nz1, 	integer *nz2, logical *rsign, doublereal *amagn, doublereal *rcond, 	doublereal *triang, integer *idist, integer *iseed, doublecomplex *a, 	integer *lda){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3, i__4;    doublereal d__1, d__2;    doublecomplex z__1, z__2;    /* Local variables */    integer i__, k, jc, jd, jr, kbeg, isdb, kend, isde, klen;    doublereal alpha;    doublecomplex ctemp;    extern doublereal dlaran_(integer *);    extern /* Double Complex */ void zlarnd_(doublecomplex *, integer *, 	    integer *);    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);/*  -- LAPACK auxiliary test routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLATM4 generates basic square matrices, which may later be *//*  multiplied by others in order to produce test matrices.  It is *//*  intended mainly to be used to test the generalized eigenvalue *//*  routines. *//*  It first generates the diagonal and (possibly) subdiagonal, *//*  according to the value of ITYPE, NZ1, NZ2, RSIGN, AMAGN, and RCOND. *//*  It then fills in the upper triangle with random numbers, if TRIANG is *//*  non-zero. *//*  Arguments *//*  ========= *//*  ITYPE   (input) INTEGER *//*          The "type" of matrix on the diagonal and sub-diagonal. *//*          If ITYPE < 0, then type abs(ITYPE) is generated and then *//*             swapped end for end (A(I,J) := A'(N-J,N-I).)  See also *//*             the description of AMAGN and RSIGN. *//*          Special types: *//*          = 0:  the zero matrix. *//*          = 1:  the identity. *//*          = 2:  a transposed Jordan block. *//*          = 3:  If N is odd, then a k+1 x k+1 transposed Jordan block *//*                followed by a k x k identity block, where k=(N-1)/2. *//*                If N is even, then k=(N-2)/2, and a zero diagonal entry *//*                is tacked onto the end. *//*          Diagonal types.  The diagonal consists of NZ1 zeros, then *//*             k=N-NZ1-NZ2 nonzeros.  The subdiagonal is zero.  ITYPE *//*             specifies the nonzero diagonal entries as follows: *//*          = 4:  1, ..., k *//*          = 5:  1, RCOND, ..., RCOND *//*          = 6:  1, ..., 1, RCOND *//*          = 7:  1, a, a^2, ..., a^(k-1)=RCOND *//*          = 8:  1, 1-d, 1-2*d, ..., 1-(k-1)*d=RCOND *//*          = 9:  random numbers chosen from (RCOND,1) *//*          = 10: random numbers with distribution IDIST (see ZLARND.) *//*  N       (input) INTEGER *//*          The order of the matrix. *//*  NZ1     (input) INTEGER *//*          If abs(ITYPE) > 3, then the first NZ1 diagonal entries will *//*          be zero. *//*  NZ2     (input) INTEGER *//*          If abs(ITYPE) > 3, then the last NZ2 diagonal entries will *//*          be zero. *//*  RSIGN   (input) LOGICAL *//*          = .TRUE.:  The diagonal and subdiagonal entries will be *//*                     multiplied by random numbers of magnitude 1. *//*          = .FALSE.: The diagonal and subdiagonal entries will be *//*                     left as they are (usually non-negative real.) *//*  AMAGN   (input) DOUBLE PRECISION *//*          The diagonal and subdiagonal entries will be multiplied by *//*          AMAGN. *//*  RCOND   (input) DOUBLE PRECISION *//*          If abs(ITYPE) > 4, then the smallest diagonal entry will be *//*          RCOND.  RCOND must be between 0 and 1. *//*  TRIANG  (input) DOUBLE PRECISION *//*          The entries above the diagonal will be random numbers with *//*          magnitude bounded by TRIANG (i.e., random numbers multiplied *///.........这里部分代码省略.........

/* Subroutine */ int zunt01_(char *rowcol, integer *m, integer *n, 	doublecomplex *u, integer *ldu, doublecomplex *work, integer *lwork, 	doublereal *rwork, doublereal *resid){    /* System generated locals */    integer u_dim1, u_offset, i__1, i__2;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1, z__2;    /* Builtin functions */    double d_imag(doublecomplex *);    /* Local variables */    integer i__, j, k;    doublereal eps;    doublecomplex tmp;    extern logical lsame_(char *, char *);    integer mnmin;    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern /* Subroutine */ int zherk_(char *, char *, integer *, integer *, 	    doublereal *, doublecomplex *, integer *, doublereal *, 	    doublecomplex *, integer *);    extern doublereal dlamch_(char *);    integer ldwork;    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);    char transu[1];    extern doublereal zlansy_(char *, char *, integer *, doublecomplex *, 	    integer *, doublereal *);/*  -- LAPACK test routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZUNT01 checks that the matrix U is unitary by computing the ratio *//*     RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R', *//*  or *//*     RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'. *//*  Alternatively, if there isn't sufficient workspace to form *//*  I - U*U' or I - U'*U, the ratio is computed as *//*     RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R', *//*  or *//*     RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'. *//*  where EPS is the machine precision.  ROWCOL is used only if m = n; *//*  if m > n, ROWCOL is assumed to be 'C', and if m < n, ROWCOL is *//*  assumed to be 'R'. *//*  Arguments *//*  ========= *//*  ROWCOL  (input) CHARACTER *//*          Specifies whether the rows or columns of U should be checked *//*          for orthogonality.  Used only if M = N. *//*          = 'R':  Check for orthogonal rows of U *//*          = 'C':  Check for orthogonal columns of U *//*  M       (input) INTEGER *//*          The number of rows of the matrix U. *//*  N       (input) INTEGER *//*          The number of columns of the matrix U. *//*  U       (input) COMPLEX*16 array, dimension (LDU,N) *//*          The unitary matrix U.  U is checked for orthogonal columns *//*          if m > n or if m = n and ROWCOL = 'C'.  U is checked for *//*          orthogonal rows if m < n or if m = n and ROWCOL = 'R'. *//*  LDU     (input) INTEGER *//*          The leading dimension of the array U.  LDU >= max(1,M). *//*  WORK    (workspace) COMPLEX*16 array, dimension (LWORK) *//*  LWORK   (input) INTEGER *//*          The length of the array WORK.  For best performance, LWORK *//*          should be at least N*N if ROWCOL = 'C' or M*M if *//*          ROWCOL = 'R', but the test will be done even if LWORK is 0. *//*  RWORK   (workspace) DOUBLE PRECISION array, dimension (min(M,N)) *//*          Used only if LWORK is large enough to use the Level 3 BLAS *//*          code. *//*  RESID   (output) DOUBLE PRECISION *//*          RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or *//*          RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'. *//*  ===================================================================== *///.........这里部分代码省略.........

/* Subroutine */ int zhgeqz_(char *job, char *compq, char *compz, integer *n, 	integer *ilo, integer *ihi, doublecomplex *a, integer *lda, 	doublecomplex *b, integer *ldb, doublecomplex *alpha, doublecomplex *	beta, doublecomplex *q, integer *ldq, doublecomplex *z__, integer *	ldz, doublecomplex *work, integer *lwork, doublereal *rwork, integer *	info){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, 	    z_offset, i__1, i__2, i__3, i__4, i__5, i__6;    doublereal d__1, d__2, d__3, d__4, d__5, d__6;    doublecomplex z__1, z__2, z__3, z__4, z__5, z__6;    /* Builtin functions */    double z_abs(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    double d_imag(doublecomplex *);    void z_div(doublecomplex *, doublecomplex *, doublecomplex *), pow_zi(	    doublecomplex *, doublecomplex *, integer *), z_sqrt(	    doublecomplex *, doublecomplex *);    /* Local variables */    static doublereal absb, atol, btol, temp, opst;    extern /* Subroutine */ int zrot_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublecomplex *);    static doublereal temp2, c__;    static integer j;    static doublecomplex s, t;    extern logical lsame_(char *, char *);    static doublecomplex ctemp;    static integer iiter, ilast, jiter;    static doublereal anorm;    static integer maxit;    static doublereal bnorm;    static doublecomplex shift;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *);    static doublereal tempr;    static doublecomplex ctemp2, ctemp3;    static logical ilazr2;    static integer jc, in;    static doublereal ascale, bscale;    static doublecomplex u12;    extern doublereal dlamch_(char *);    static integer jr, nq;    static doublecomplex signbc;    static integer nz;    static doublereal safmin;    extern /* Subroutine */ int xerbla_(char *, integer *);    static doublecomplex eshift;    static logical ilschr;    static integer icompq, ilastm;    static doublecomplex rtdisc;    static integer ischur;    extern doublereal zlanhs_(char *, integer *, doublecomplex *, integer *, 	    doublereal *);    static logical ilazro;    static integer icompz, ifirst;    extern /* Subroutine */ int zlartg_(doublecomplex *, doublecomplex *, 	    doublereal *, doublecomplex *, doublecomplex *);    static integer ifrstm;    extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);    static integer istart;    static logical lquery;    static doublecomplex ad11, ad12, ad21, ad22;    static integer jch;    static logical ilq, ilz;    static doublereal ulp;    static doublecomplex abi22;#define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1#define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]#define b_subscr(a_1,a_2) (a_2)*b_dim1 + a_1#define b_ref(a_1,a_2) b[b_subscr(a_1,a_2)]#define q_subscr(a_1,a_2) (a_2)*q_dim1 + a_1#define q_ref(a_1,a_2) q[q_subscr(a_1,a_2)]#define z___subscr(a_1,a_2) (a_2)*z_dim1 + a_1#define z___ref(a_1,a_2) z__[z___subscr(a_1,a_2)]/*  -- LAPACK routine (instrumented to count operations, version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          June 30, 1999          ----------------------- Begin Timing Code ------------------------          Common block to return operation count and iteration count          ITCNT is initialized to 0, OPS is only incremented          OPST is used to accumulate small contributions to OPS          to avoid roundoff error          ------------------------ End Timing Code -------------------------       Purpose       =======       ZHGEQZ implements a single-shift version of the QZ   //.........这里部分代码省略.........

/* Subroutine */ int zget36_(doublereal *rmax, integer *lmax, integer *ninfo, 	integer *knt, integer *nin){    /* System generated locals */    integer i__1, i__2, i__3, i__4;    /* Builtin functions */    integer s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 	    e_rsle(void);    /* Local variables */    static doublecomplex diag[10];    static integer ifst, ilst;    static doublecomplex work[200];    static integer info1, info2, i__, j, n;    static doublecomplex q[100]	/* was [10][10] */, ctemp;    extern /* Subroutine */ int zhst01_(integer *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublereal *, doublereal *);    static doublereal rwork[10];    static doublecomplex t1[100]	/* was [10][10] */, t2[100]	/* 	    was [10][10] */;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    extern doublereal dlamch_(char *);    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *), 	    zlaset_(char *, integer *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *, integer *);    static doublereal result[2];    extern /* Subroutine */ int ztrexc_(char *, integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *, integer *, integer *, 	    integer *);    static doublereal eps, res;    static doublecomplex tmp[100]	/* was [10][10] */;    /* Fortran I/O blocks */    static cilist io___2 = { 0, 0, 0, 0, 0 };    static cilist io___7 = { 0, 0, 0, 0, 0 };#define q_subscr(a_1,a_2) (a_2)*10 + a_1 - 11#define q_ref(a_1,a_2) q[q_subscr(a_1,a_2)]#define t1_subscr(a_1,a_2) (a_2)*10 + a_1 - 11#define t1_ref(a_1,a_2) t1[t1_subscr(a_1,a_2)]#define t2_subscr(a_1,a_2) (a_2)*10 + a_1 - 11#define t2_ref(a_1,a_2) t2[t2_subscr(a_1,a_2)]#define tmp_subscr(a_1,a_2) (a_2)*10 + a_1 - 11#define tmp_ref(a_1,a_2) tmp[tmp_subscr(a_1,a_2)]/*  -- LAPACK test routine (version 3.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          September 30, 1994       Purpose       =======       ZGET36 tests ZTREXC, a routine for reordering diagonal entries of a       matrix in complex Schur form. Thus, ZLAEXC computes a unitary matrix       Q such that          Q' * T1 * Q  = T2       and where one of the diagonal blocks of T1 (the one at row IFST) has       been moved to position ILST.       The test code verifies that the residual Q'*T1*Q-T2 is small, that T2       is in Schur form, and that the final position of the IFST block is       ILST.       The test matrices are read from a file with logical unit number NIN.       Arguments       ==========       RMAX    (output) DOUBLE PRECISION               Value of the largest test ratio.       LMAX    (output) INTEGER               Example number where largest test ratio achieved.       NINFO   (output) INTEGER               Number of examples where INFO is nonzero.       KNT     (output) INTEGER               Total number of examples tested.       NIN     (input) INTEGER               Input logical unit number.       ===================================================================== */    eps = dlamch_("P");    *rmax = 0.;//.........这里部分代码省略.........

/* Subroutine */ int zlqt02_(integer *m, integer *n, integer *k, 	doublecomplex *a, doublecomplex *af, doublecomplex *q, doublecomplex *	l, integer *lda, doublecomplex *tau, doublecomplex *work, integer *	lwork, doublereal *rwork, doublereal *result){    /* System generated locals */    integer a_dim1, a_offset, af_dim1, af_offset, l_dim1, l_offset, q_dim1, 	    q_offset, i__1;    /* Builtin functions */    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);    /* Local variables */    doublereal eps;    integer info;    doublereal resid, anorm;    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *), zherk_(char *, char *, integer *, 	    integer *, doublereal *, doublecomplex *, integer *, doublereal *, 	     doublecomplex *, integer *);    extern doublereal dlamch_(char *), zlange_(char *, integer *, 	    integer *, doublecomplex *, integer *, doublereal *);    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *), 	    zlaset_(char *, integer *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *, integer *);    extern doublereal zlansy_(char *, char *, integer *, doublecomplex *, 	    integer *, doublereal *);    extern /* Subroutine */ int zunglq_(integer *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *, integer *);/*  -- LAPACK test routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLQT02 tests ZUNGLQ, which generates an m-by-n matrix Q with *//*  orthonornmal rows that is defined as the product of k elementary *//*  reflectors. *//*  Given the LQ factorization of an m-by-n matrix A, ZLQT02 generates *//*  the orthogonal matrix Q defined by the factorization of the first k *//*  rows of A; it compares L(1:k,1:m) with A(1:k,1:n)*Q(1:m,1:n)', and *//*  checks that the rows of Q are orthonormal. *//*  Arguments *//*  ========= *//*  M       (input) INTEGER *//*          The number of rows of the matrix Q to be generated.  M >= 0. *//*  N       (input) INTEGER *//*          The number of columns of the matrix Q to be generated. *//*          N >= M >= 0. *//*  K       (input) INTEGER *//*          The number of elementary reflectors whose product defines the *//*          matrix Q. M >= K >= 0. *//*  A       (input) COMPLEX*16 array, dimension (LDA,N) *//*          The m-by-n matrix A which was factorized by ZLQT01. *//*  AF      (input) COMPLEX*16 array, dimension (LDA,N) *//*          Details of the LQ factorization of A, as returned by ZGELQF. *//*          See ZGELQF for further details. *//*  Q       (workspace) COMPLEX*16 array, dimension (LDA,N) *//*  L       (workspace) COMPLEX*16 array, dimension (LDA,M) *//*  LDA     (input) INTEGER *//*          The leading dimension of the arrays A, AF, Q and L. LDA >= N. *//*  TAU     (input) COMPLEX*16 array, dimension (M) *//*          The scalar factors of the elementary reflectors corresponding *//*          to the LQ factorization in AF. *//*  WORK    (workspace) COMPLEX*16 array, dimension (LWORK) *//*  LWORK   (input) INTEGER *//*          The dimension of the array WORK. *//*  RWORK   (workspace) DOUBLE PRECISION array, dimension (M) *//*  RESULT  (output) DOUBLE PRECISION array, dimension (2) *//*          The test ratios: *//*          RESULT(1) = norm( L - A*Q' ) / ( N * norm(A) * EPS ) *//*          RESULT(2) = norm( I - Q*Q' ) / ( N * EPS ) *///.........这里部分代码省略.........

/* Subroutine */ int ztgsja_(char *jobu, char *jobv, char *jobq, integer *m, 	integer *p, integer *n, integer *k, integer *l, doublecomplex *a, 	integer *lda, doublecomplex *b, integer *ldb, doublereal *tola, 	doublereal *tolb, doublereal *alpha, doublereal *beta, doublecomplex *	u, integer *ldu, doublecomplex *v, integer *ldv, doublecomplex *q, 	integer *ldq, doublecomplex *work, integer *ncycle, integer *info){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, u_dim1, 	    u_offset, v_dim1, v_offset, i__1, i__2, i__3, i__4;    doublereal d__1;    doublecomplex z__1;    /* Builtin functions */    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j;    doublereal a1, b1, a3, b3;    doublecomplex a2, b2;    doublereal csq, csu, csv;    doublecomplex snq;    doublereal rwk;    doublecomplex snu, snv;    extern /* Subroutine */ int zrot_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublecomplex *);    doublereal gamma;    extern logical lsame_(char *, char *);    logical initq, initu, initv, wantq, upper;    doublereal error, ssmin;    logical wantu, wantv;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zlags2_(logical *, doublereal *, 	    doublecomplex *, doublereal *, doublereal *, doublecomplex *, 	    doublereal *, doublereal *, doublecomplex *, doublereal *, 	    doublecomplex *, doublereal *, doublecomplex *);    integer kcycle;    extern /* Subroutine */ int dlartg_(doublereal *, doublereal *, 	    doublereal *, doublereal *, doublereal *), xerbla_(char *, 	    integer *), zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *), zlapll_(integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *, doublereal *), zlaset_(	    char *, integer *, integer *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *);/*  -- LAPACK routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZTGSJA computes the generalized singular value decomposition (GSVD) *//*  of two complex upper triangular (or trapezoidal) matrices A and B. *//*  On entry, it is assumed that matrices A and B have the following *//*  forms, which may be obtained by the preprocessing subroutine ZGGSVP *//*  from a general M-by-N matrix A and P-by-N matrix B: *//*               N-K-L  K    L *//*     A =    K ( 0    A12  A13 ) if M-K-L >= 0; *//*            L ( 0     0   A23 ) *//*        M-K-L ( 0     0    0  ) *//*             N-K-L  K    L *//*     A =  K ( 0    A12  A13 ) if M-K-L < 0; *//*        M-K ( 0     0   A23 ) *//*             N-K-L  K    L *//*     B =  L ( 0     0   B13 ) *//*        P-L ( 0     0    0  ) *//*  where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular *//*  upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, *//*  otherwise A23 is (M-K)-by-L upper trapezoidal. *//*  On exit, *//*         U'*A*Q = D1*( 0 R ),    V'*B*Q = D2*( 0 R ), *//*  where U, V and Q are unitary matrices, Z' denotes the conjugate *//*  transpose of Z, R is a nonsingular upper triangular matrix, and D1 *//*  and D2 are ``diagonal'' matrices, which are of the following *//*  structures: *//*  If M-K-L >= 0, *//*                      K  L *//*         D1 =     K ( I  0 ) *//*                  L ( 0  C ) *//*              M-K-L ( 0  0 ) *//*                     K  L *//*         D2 = L   ( 0  S ) *///.........这里部分代码省略.........

/* Subroutine */ int zchkqp_(logical *dotype, integer *nm, integer *mval, 	integer *nn, integer *nval, doublereal *thresh, logical *tsterr, 	doublecomplex *a, doublecomplex *copya, doublereal *s, doublereal *	copys, doublecomplex *tau, doublecomplex *work, doublereal *rwork, 	integer *iwork, integer *nout){    /* Initialized data */    static integer iseedy[4] = { 1988,1989,1990,1991 };    /* Format strings */    static char fmt_9999[] = "(/002 M =/002,i5,/002, N =/002,i5,/002, type"	    " /002,i2,/002, test /002,i2,/002, ratio =/002,g12.5)";    /* System generated locals */    integer i__1, i__2, i__3, i__4;    doublereal d__1;    /* Builtin functions */    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);    /* Local variables */    integer i__, k, m, n, im, in, lda;    doublereal eps;    integer mode, info;    char path[3];    integer ilow, nrun;    extern /* Subroutine */ int alahd_(integer *, char *);    integer ihigh, nfail, iseed[4], imode, mnmin, istep, nerrs, lwork;    extern doublereal zqpt01_(integer *, integer *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zqrt11_(integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *), zqrt12_(integer *, integer *, doublecomplex *, 	    integer *, doublereal *, doublecomplex *, integer *, doublereal *)	    , dlamch_(char *);    extern /* Subroutine */ int dlaord_(char *, integer *, doublereal *, 	    integer *), alasum_(char *, integer *, integer *, integer 	    *, integer *), zgeqpf_(integer *, integer *, 	    doublecomplex *, integer *, integer *, doublecomplex *, 	    doublecomplex *, doublereal *, integer *), zlacpy_(char *, 	    integer *, integer *, doublecomplex *, integer *, doublecomplex *, 	     integer *), zlaset_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *), zlatms_(integer *, integer *, char *, integer *, char *, 	    doublereal *, integer *, doublereal *, doublereal *, integer *, 	    integer *, char *, doublecomplex *, integer *, doublecomplex *, 	    integer *);    doublereal result[3];    extern /* Subroutine */ int zerrqp_(char *, integer *);    /* Fortran I/O blocks */    static cilist io___24 = { 0, 0, 0, fmt_9999, 0 };/*  -- LAPACK test routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZCHKQP tests ZGEQPF. *//*  Arguments *//*  ========= *//*  DOTYPE  (input) LOGICAL array, dimension (NTYPES) *//*          The matrix types to be used for testing.  Matrices of type j *//*          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = *//*          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. *//*  NM      (input) INTEGER *//*          The number of values of M contained in the vector MVAL. *//*  MVAL    (input) INTEGER array, dimension (NM) *//*          The values of the matrix row dimension M. *//*  NN      (input) INTEGER *//*          The number of values of N contained in the vector NVAL. *//*  NVAL    (input) INTEGER array, dimension (NN) *//*          The values of the matrix column dimension N. *//*  THRESH  (input) DOUBLE PRECISION *//*          The threshold value for the test ratios.  A result is *//*          included in the output file if RESULT >= THRESH.  To have *//*          every test ratio printed, use THRESH = 0. *//*  TSTERR  (input) LOGICAL *//*          Flag that indicates whether error exits are to be tested. *//*  A       (workspace) COMPLEX*16 array, dimension (MMAX*NMAX) *//*          where MMAX is the maximum value of M in MVAL and NMAX is the *///.........这里部分代码省略.........

/* Subroutine */ int zgrqts_(integer *m, integer *p, integer *n,                             doublecomplex *a, doublecomplex *af, doublecomplex *q, doublecomplex *                             r__, integer *lda, doublecomplex *taua, doublecomplex *b,                             doublecomplex *bf, doublecomplex *z__, doublecomplex *t,                             doublecomplex *bwk, integer *ldb, doublecomplex *taub, doublecomplex *                             work, integer *lwork, doublereal *rwork, doublereal *result){    /* System generated locals */    integer a_dim1, a_offset, af_dim1, af_offset, b_dim1, b_offset, bf_dim1,            bf_offset, bwk_dim1, bwk_offset, q_dim1, q_offset, r_dim1,            r_offset, t_dim1, t_offset, z_dim1, z_offset, i__1, i__2;    doublereal d__1;    doublecomplex z__1;    /* Local variables */    static integer info;    static doublereal unfl, resid, anorm, bnorm;    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *,                                       integer *, doublecomplex *, doublecomplex *, integer *,                                       doublecomplex *, integer *, doublecomplex *, doublecomplex *,                                       integer *), zherk_(char *, char *, integer *,                                               integer *, doublereal *, doublecomplex *, integer *, doublereal *,                                               doublecomplex *, integer *);    extern doublereal dlamch_(char *), zlange_(char *, integer *,            integer *, doublecomplex *, integer *, doublereal *),                    zlanhe_(char *, char *, integer *, doublecomplex *, integer *,                            doublereal *);    extern /* Subroutine */ int zggrqf_(integer *, integer *, integer *,                                        doublecomplex *, integer *, doublecomplex *, doublecomplex *,                                        integer *, doublecomplex *, doublecomplex *, integer *, integer *)    , zlacpy_(char *, integer *, integer *, doublecomplex *, integer *              , doublecomplex *, integer *), zlaset_(char *, integer *,                      integer *, doublecomplex *, doublecomplex *, doublecomplex *,                      integer *), zungqr_(integer *, integer *, integer *,                                          doublecomplex *, integer *, doublecomplex *, doublecomplex *,                                          integer *, integer *), zungrq_(integer *, integer *, integer *,                                                  doublecomplex *, integer *, doublecomplex *, doublecomplex *,                                                  integer *, integer *);    static doublereal ulp;#define q_subscr(a_1,a_2) (a_2)*q_dim1 + a_1#define q_ref(a_1,a_2) q[q_subscr(a_1,a_2)]#define r___subscr(a_1,a_2) (a_2)*r_dim1 + a_1#define r___ref(a_1,a_2) r__[r___subscr(a_1,a_2)]#define z___subscr(a_1,a_2) (a_2)*z_dim1 + a_1#define z___ref(a_1,a_2) z__[z___subscr(a_1,a_2)]#define af_subscr(a_1,a_2) (a_2)*af_dim1 + a_1#define af_ref(a_1,a_2) af[af_subscr(a_1,a_2)]#define bf_subscr(a_1,a_2) (a_2)*bf_dim1 + a_1#define bf_ref(a_1,a_2) bf[bf_subscr(a_1,a_2)]    /*  -- LAPACK test routine (version 3.0) --           Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,           Courant Institute, Argonne National Lab, and Rice University           September 30, 1994        Purpose        =======        ZGRQTS tests ZGGRQF, which computes the GRQ factorization of an        M-by-N matrix A and a P-by-N matrix B: A = R*Q and B = Z*T*Q.        Arguments        =========        M       (input) INTEGER                The number of rows of the matrix A.  M >= 0.        P       (input) INTEGER                The number of rows of the matrix B.  P >= 0.        N       (input) INTEGER                The number of columns of the matrices A and B.  N >= 0.        A       (input) COMPLEX*16 array, dimension (LDA,N)                The M-by-N matrix A.        AF      (output) COMPLEX*16 array, dimension (LDA,N)                Details of the GRQ factorization of A and B, as returned                by ZGGRQF, see CGGRQF for further details.        Q       (output) COMPLEX*16 array, dimension (LDA,N)                The N-by-N unitary matrix Q.        R       (workspace) COMPLEX*16 array, dimension (LDA,MAX(M,N))        LDA     (input) INTEGER                The leading dimension of the arrays A, AF, R and Q.                LDA >= max(M,N).        TAUA    (output) COMPLEX*16 array, dimension (min(M,N))                The scalar factors of the elementary reflectors, as returned                by DGGQRC.        B       (input) COMPLEX*16 array, dimension (LDB,N)                On entry, the P-by-N matrix A.//.........这里部分代码省略.........