/* Subroutine */ int cdrvpb_(logical *dotype, integer *nn, integer *nval, 
	integer *nrhs, real *thresh, logical *tsterr, integer *nmax, complex *
	a, complex *afac, complex *asav, complex *b, complex *bsav, complex *
	x, complex *xact, real *s, complex *work, real *rwork, integer *nout)
{
    /* Initialized data */

    static integer iseedy[4] = { 1988,1989,1990,1991 };
    static char facts[1*3] = "F" "N" "E";
    static char equeds[1*2] = "N" "Y";

    /* Format strings */
    static char fmt_9999[] = "(1x,a6,\002, UPLO='\002,a1,\002', N =\002,i5"
	    ",\002, KD =\002,i5,\002, type \002,i1,\002, test(\002,i1,\002)"
	    "=\002,g12.5)";
    static char fmt_9997[] = "(1x,a6,\002( '\002,a1,\002', '\002,a1,\002',"
	    " \002,i5,\002, \002,i5,\002, ... ), EQUED='\002,a1,\002', type"
	    " \002,i1,\002, test(\002,i1,\002)=\002,g12.5)";
    static char fmt_9998[] = "(1x,a6,\002( '\002,a1,\002', '\002,a1,\002',"
	    " \002,i5,\002, \002,i5,\002, ... ), type \002,i1,\002, test(\002"
	    ",i1,\002)=\002,g12.5)";

    /* System generated locals */
    address a__1[2];
    integer i__1, i__2, i__3, i__4, i__5, i__6, i__7[2];
    char ch__1[2];

    /* Builtin functions */
    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
    /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);

    /* Local variables */
    integer i__, k, n, i1, i2, k1, kd, nb, in, kl, iw, ku, nt, lda, ikd, nkd, 
	    ldab;
    char fact[1];
    integer ioff, mode, koff;
    real amax;
    char path[3];
    integer imat, info;
    char dist[1], uplo[1], type__[1];
    integer nrun, ifact;
    extern /* Subroutine */ int cget04_(integer *, integer *, complex *, 
	    integer *, complex *, integer *, real *, real *);
    integer nfail, iseed[4], nfact;
    extern /* Subroutine */ int cpbt01_(char *, integer *, integer *, complex 
	    *, integer *, complex *, integer *, real *, real *), 
	    cpbt02_(char *, integer *, integer *, integer *, complex *, 
	    integer *, complex *, integer *, complex *, integer *, real *, 
	    real *), cpbt05_(char *, integer *, integer *, integer *, 
	    complex *, integer *, complex *, integer *, complex *, integer *, 
	    complex *, integer *, real *, real *, real *);
    integer kdval[4];
    extern logical lsame_(char *, char *);
    char equed[1];
    integer nbmin;
    real rcond, roldc, scond;
    integer nimat;
    extern doublereal sget06_(real *, real *);
    real anorm;
    extern /* Subroutine */ int ccopy_(integer *, complex *, integer *, 
	    complex *, integer *), cpbsv_(char *, integer *, integer *, 
	    integer *, complex *, integer *, complex *, integer *, integer *);
    logical equil;
    extern /* Subroutine */ int cswap_(integer *, complex *, integer *, 
	    complex *, integer *);
    integer iuplo, izero, nerrs;
    logical zerot;
    char xtype[1];
    extern /* Subroutine */ int clatb4_(char *, integer *, integer *, integer 
	    *, char *, integer *, integer *, real *, integer *, real *, char *
), aladhd_(integer *, char *);
    extern doublereal clanhb_(char *, char *, integer *, integer *, complex *, 
	     integer *, real *), clange_(char *, integer *, 
	    integer *, complex *, integer *, real *);
    extern /* Subroutine */ int claqhb_(char *, integer *, integer *, complex 
	    *, integer *, real *, real *, real *, char *), 
	    alaerh_(char *, char *, integer *, integer *, char *, integer *, 
	    integer *, integer *, integer *, integer *, integer *, integer *, 
	    integer *, integer *), claipd_(integer *, 
	    complex *, integer *, integer *);
    logical prefac;
    real rcondc;
    logical nofact;
    char packit[1];
    integer iequed;
    extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex 
	    *, integer *, complex *, integer *), clarhs_(char *, char 
	    *, char *, char *, integer *, integer *, integer *, integer *, 
	    integer *, complex *, integer *, complex *, integer *, complex *, 
	    integer *, integer *, integer *), 
	    claset_(char *, integer *, integer *, complex *, complex *, 
	    complex *, integer *), cpbequ_(char *, integer *, integer 
	    *, complex *, integer *, real *, real *, real *, integer *), alasvm_(char *, integer *, integer *, integer *, integer 
	    *);
    real cndnum;
    extern /* Subroutine */ int clatms_(integer *, integer *, char *, integer 
	    *, char *, real *, integer *, real *, real *, integer *, integer *
, char *, complex *, integer *, complex *, integer *), cpbtrf_(char *, integer *, integer *, complex *, 
	    integer *, integer *);
//.........这里部分代码省略.........

 int chgeqz_(char *job, char *compq, char *compz, int *n, 
	int *ilo, int *ihi, complex *h__, int *ldh, complex *t, 
	int *ldt, complex *alpha, complex *beta, complex *q, int *ldq, 
	 complex *z__, int *ldz, complex *work, int *lwork, float *
	rwork, int *info)
{
    /* System generated locals */
    int h_dim1, h_offset, q_dim1, q_offset, t_dim1, t_offset, z_dim1, 
	    z_offset, i__1, i__2, i__3, i__4, i__5, i__6;
    float r__1, r__2, r__3, r__4, r__5, r__6;
    complex q__1, q__2, q__3, q__4, q__5, q__6;

    /* Builtin functions */
    double c_abs(complex *);
    void r_cnjg(complex *, complex *);
    double r_imag(complex *);
    void c_div(complex *, complex *, complex *), pow_ci(complex *, complex *, 
	    int *), c_sqrt(complex *, complex *);

    /* Local variables */
    float c__;
    int j;
    complex s, t1;
    int jc, in;
    complex u12;
    int jr;
    complex ad11, ad12, ad21, ad22;
    int jch;
    int ilq, ilz;
    float ulp;
    complex abi22;
    float absb, atol, btol, temp;
    extern  int crot_(int *, complex *, int *, 
	    complex *, int *, float *, complex *);
    float temp2;
    extern  int cscal_(int *, complex *, complex *, 
	    int *);
    extern int lsame_(char *, char *);
    complex ctemp;
    int iiter, ilast, jiter;
    float anorm, bnorm;
    int maxit;
    complex shift;
    float tempr;
    complex ctemp2, ctemp3;
    int ilazr2;
    float ascale, bscale;
    complex signbc;
    extern double slamch_(char *), clanhs_(char *, int *, 
	    complex *, int *, float *);
    extern  int claset_(char *, int *, int *, complex 
	    *, complex *, complex *, int *), clartg_(complex *, 
	    complex *, float *, complex *, complex *);
    float safmin;
    extern  int xerbla_(char *, int *);
    complex eshift;
    int ilschr;
    int icompq, ilastm;
    complex rtdisc;
    int ischur;
    int ilazro;
    int icompz, ifirst, ifrstm, istart;
    int lquery;


/*  -- 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 .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

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

/*  CHGEQZ computes the eigenvalues of a complex matrix pair (H,T), */
/*  where H is an upper Hessenberg matrix and T is upper triangular, */
/*  using the single-shift QZ method. */
/*  Matrix pairs of this type are produced by the reduction to */
/*  generalized upper Hessenberg form of a complex matrix pair (A,B): */

/*     A = Q1*H*Z1**H,  B = Q1*T*Z1**H, */

/*  as computed by CGGHRD. */

/*  If JOB='S', then the Hessenberg-triangular pair (H,T) is */
/*  also reduced to generalized Schur form, */

/*     H = Q*S*Z**H,  T = Q*P*Z**H, */

/*  where Q and Z are unitary matrices and S and P are upper triangular. */

/*  Optionally, the unitary matrix Q from the generalized Schur */
/*  factorization may be postmultiplied into an input matrix Q1, and the */
/*  unitary matrix Z may be postmultiplied into an input matrix Z1. */
/*  If Q1 and Z1 are the unitary matrices from CGGHRD that reduced */
/*  the matrix pair (A,B) to generalized Hessenberg form, then the output */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
/*          The leading dimension of the arrays A and AF. */

/*  TAU     (input) COMPLEX array, dimension (M) */
/*          Details of the  Householder transformations as returned by */
/*          CTZRQF. */

/*  WORK    (workspace) COMPLEX array, dimension (LWORK) */

/*  LWORK   (input) INTEGER */
/*          The length of the array WORK.  LWORK >= m*n + m. */

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

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Local Arrays .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

    /* Parameter adjustments */
    af_dim1 = *lda;
    af_offset = 1 + af_dim1;
    af -= af_offset;
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    --tau;
    --work;

    /* Function Body */
    ret_val = 0.f;

    if (*lwork < *m * *n + *m) {
	this_xerbla_("CTZT01", &c__8);
	return ret_val;
    }

/*     Quick return if possible */

    if (*m <= 0 || *n <= 0) {
	return ret_val;
    }

    norma = clange_("One-norm", m, n, &a[a_offset], lda, rwork);

/*     Copy upper triangle R */

    claset_("Full", m, n, &c_b6, &c_b6, &work[1], m);
    i__1 = *m;
    for (j = 1; j <= i__1; ++j) {
	i__2 = j;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    i__3 = (j - 1) * *m + i__;
	    i__4 = i__ + j * af_dim1;
	    work[i__3].r = af[i__4].r, work[i__3].i = af[i__4].i;
/* L10: */
	}
/* L20: */
    }

/*     R = R * P(1) * ... *P(m) */

    i__1 = *m;
    for (i__ = 1; i__ <= i__1; ++i__) {
	i__2 = *n - *m + 1;
	clatzm_("Right", &i__, &i__2, &af[i__ + (*m + 1) * af_dim1], lda, &
		tau[i__], &work[(i__ - 1) * *m + 1], &work[*m * *m + 1], m, &
		work[*m * *n + 1]);
/* L30: */
    }

/*     R = R - A */

    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	caxpy_(m, &c_b15, &a[i__ * a_dim1 + 1], &c__1, &work[(i__ - 1) * *m + 
		1], &c__1);
/* L40: */
    }

    ret_val = clange_("One-norm", m, n, &work[1], m, rwork);

    ret_val /= slamch_("Epsilon") * (real) max(*m,*n);
    if (norma != 0.f) {
	ret_val /= norma;
    }

    return ret_val;

/*     End of CTZT01 */

} /* ctzt01_ */

//.........这里部分代码省略.........
				    ;
			    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(real));
			    e_wsfe();
			    ++nfail;
			}
/* L70: */
		    }
		    nrun += nt;
		}

/*              --- Test CPTSVX --- */

		if (ifact > 1) {

/*                 Initialize D( N+1:2*N ) and E( N+1:2*N ) to zero. */

		    i__3 = n - 1;
		    for (i__ = 1; i__ <= i__3; ++i__) {
			d__[n + i__] = 0.f;
			i__4 = n + i__;
			e[i__4].r = 0.f, e[i__4].i = 0.f;
/* L80: */
		    }
		    if (n > 0) {
			d__[n + n] = 0.f;
		    }
		}

		claset_("Full", &n, nrhs, &c_b62, &c_b62, &x[1], &lda);

/*              Solve the system and compute the condition number and */
/*              error bounds using CPTSVX. */

		s_copy(srnamc_1.srnamt, "CPTSVX", (ftnlen)32, (ftnlen)6);
		cptsvx_(fact, &n, nrhs, &d__[1], &e[1], &d__[n + 1], &e[n + 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 CPTSVX. */

		if (info != izero) {
		    alaerh_(path, "CPTSVX", &info, &izero, fact, &n, &n, &
			    c__1, &c__1, nrhs, &imat, &nfail, &nerrs, nout);
		}
		if (izero == 0) {
		    if (ifact == 2) {

/*                    Check the factorization by computing the ratio */
/*                       norm(L*D*L' - A) / (n * norm(A) * EPS ) */

			k1 = 1;
			cptt01_(&n, &d__[1], &e[1], &d__[n + 1], &e[n + 1], &
				work[1], result);
		    } else {
			k1 = 2;
		    }

/*                 Compute the residual in the solution. */

		    clacpy_("Full", &n, nrhs, &b[1], &lda, &work[1], &lda);
		    cptt02_("Lower", &n, nrhs, &d__[1], &e[1], &x[1], &lda, &

//.........这里部分代码省略.........
    integer izero, nerrs;
    logical zerot;
    char xtype[1];
    extern /* Subroutine */ int clatb4_(char *, integer *, integer *, integer 
	    *, char *, integer *, integer *, real *, integer *, real *, char *
), aladhd_(integer *, char *);
    extern doublereal clangb_(char *, integer *, integer *, integer *, 
	    complex *, integer *, real *), clange_(char *, integer *, 
	    integer *, complex *, integer *, real *);
    extern /* Subroutine */ int claqgb_(integer *, integer *, integer *, 
	    integer *, complex *, integer *, real *, real *, real *, real *, 
	    real *, char *), alaerh_(char *, char *, integer *, 
	    integer *, char *, integer *, integer *, integer *, integer *, 
	    integer *, integer *, integer *, integer *, integer *);
    logical prefac;
    real colcnd;
    extern doublereal clantb_(char *, char *, char *, integer *, integer *, 
	    complex *, integer *, real *);
    extern /* Subroutine */ int cgbequ_(integer *, integer *, integer *, 
	    integer *, complex *, integer *, real *, real *, real *, real *, 
	    real *, integer *);
    real rcondc;
    extern doublereal slamch_(char *);
    logical nofact;
    extern /* Subroutine */ int cgbtrf_(integer *, integer *, integer *, 
	    integer *, complex *, integer *, integer *, integer *);
    integer iequed;
    extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex 
	    *, integer *, complex *, integer *);
    real rcondi;
    extern /* Subroutine */ int clarhs_(char *, char *, char *, char *, 
	    integer *, integer *, integer *, integer *, integer *, complex *, 
	    integer *, complex *, integer *, complex *, integer *, integer *, 
	    integer *), claset_(char *, 
	    integer *, integer *, complex *, complex *, complex *, integer *), alasvm_(char *, integer *, integer *, integer *, integer 
	    *);
    real cndnum, anormi, rcondo, ainvnm;
    extern /* Subroutine */ int cgbtrs_(char *, integer *, integer *, integer 
	    *, integer *, complex *, integer *, integer *, complex *, integer 
	    *, integer *), clatms_(integer *, integer *, char *, 
	    integer *, char *, real *, integer *, real *, real *, integer *, 
	    integer *, char *, complex *, integer *, complex *, integer *);
    logical trfcon;
    real anormo, rowcnd;
    extern /* Subroutine */ int cgbsvx_(char *, char *, integer *, integer *, 
	    integer *, integer *, complex *, integer *, complex *, integer *, 
	    integer *, char *, real *, real *, complex *, integer *, complex *
, integer *, real *, real *, real *, complex *, real *, integer *), xlaenv_(integer *, integer *);
    real anrmpv;
    extern /* Subroutine */ int cerrvx_(char *, integer *);
    real result[7], rpvgrw;

    /* Fortran I/O blocks */
    static cilist io___26 = { 0, 0, 0, fmt_9999, 0 };
    static cilist io___27 = { 0, 0, 0, fmt_9998, 0 };
    static cilist io___65 = { 0, 0, 0, fmt_9997, 0 };
    static cilist io___73 = { 0, 0, 0, fmt_9995, 0 };
    static cilist io___74 = { 0, 0, 0, fmt_9996, 0 };
    static cilist io___75 = { 0, 0, 0, fmt_9995, 0 };
    static cilist io___76 = { 0, 0, 0, fmt_9996, 0 };
    static cilist io___77 = { 0, 0, 0, fmt_9995, 0 };
    static cilist io___78 = { 0, 0, 0, fmt_9996, 0 };
    static cilist io___79 = { 0, 0, 0, fmt_9995, 0 };
    static cilist io___80 = { 0, 0, 0, fmt_9996, 0 };

/* Subroutine */
int cgelsd_(integer *m, integer *n, integer *nrhs, complex * a, integer *lda, complex *b, integer *ldb, real *s, real *rcond, integer *rank, complex *work, integer *lwork, real *rwork, integer * iwork, integer *info)
{
    /* System generated locals */
    integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4;
    /* Builtin functions */
    double log(doublereal);
    /* Local variables */
    integer ie, il, mm;
    real eps, anrm, bnrm;
    integer itau, nlvl, iascl, ibscl;
    real sfmin;
    integer minmn, maxmn, itaup, itauq, mnthr, nwork;
    extern /* Subroutine */
    int cgebrd_(integer *, integer *, complex *, integer *, real *, real *, complex *, complex *, complex *, integer *, integer *), slabad_(real *, real *);
    extern real clange_(char *, integer *, integer *, complex *, integer *, real *);
    extern /* Subroutine */
    int cgelqf_(integer *, integer *, complex *, integer *, complex *, complex *, integer *, integer *), clalsd_( char *, integer *, integer *, integer *, real *, real *, complex * , integer *, real *, integer *, complex *, real *, integer *, integer *), clascl_(char *, integer *, integer *, real *, real *, integer *, integer *, complex *, integer *, integer *), cgeqrf_(integer *, integer *, complex *, integer *, complex *, complex *, integer *, integer *);
    extern real slamch_(char *);
    extern /* Subroutine */
    int clacpy_(char *, integer *, integer *, complex *, integer *, complex *, integer *), claset_(char *, integer *, integer *, complex *, complex *, complex *, integer *), xerbla_(char *, integer *);
    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *);
    real bignum;
    extern /* Subroutine */
    int slascl_(char *, integer *, integer *, real *, real *, integer *, integer *, real *, integer *, integer *), cunmbr_(char *, char *, char *, integer *, integer *, integer *, complex *, integer *, complex *, complex *, integer *, complex *, integer *, integer *), slaset_( char *, integer *, integer *, real *, real *, real *, integer *), cunmlq_(char *, char *, integer *, integer *, integer *, complex *, integer *, complex *, complex *, integer *, complex *, integer *, integer *);
    integer ldwork;
    extern /* Subroutine */
    int cunmqr_(char *, char *, integer *, integer *, integer *, complex *, integer *, complex *, complex *, integer *, complex *, integer *, integer *);
    integer liwork, minwrk, maxwrk;
    real smlnum;
    integer lrwork;
    logical lquery;
    integer nrwork, smlsiz;
    /* -- LAPACK driver routine (version 3.4.0) -- */
    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
    /* November 2011 */
    /* .. Scalar Arguments .. */
    /* .. */
    /* .. Array Arguments .. */
    /* .. */
    /* ===================================================================== */
    /* .. Parameters .. */
    /* .. */
    /* .. Local Scalars .. */
    /* .. */
    /* .. External Subroutines .. */
    /* .. */
    /* .. External Functions .. */
    /* .. */
    /* .. Intrinsic Functions .. */
    /* .. */
    /* .. Executable Statements .. */
    /* Test the input arguments. */
    /* Parameter adjustments */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    b_dim1 = *ldb;
    b_offset = 1 + b_dim1;
    b -= b_offset;
    --s;
    --work;
    --rwork;
    --iwork;
    /* Function Body */
    *info = 0;
    minmn = min(*m,*n);
    maxmn = max(*m,*n);
    lquery = *lwork == -1;
    if (*m < 0)
    {
        *info = -1;
    }
    else if (*n < 0)
    {
        *info = -2;
    }
    else if (*nrhs < 0)
    {
        *info = -3;
    }
    else if (*lda < max(1,*m))
    {
        *info = -5;
    }
    else if (*ldb < max(1,maxmn))
    {
        *info = -7;
    }
    /* Compute workspace. */
    /* (Note: Comments in the code beginning "Workspace:" describe the */
    /* minimal amount of workspace needed at that point in the code, */
    /* as well as the preferred amount for good performance. */
    /* NB refers to the optimal block size for the immediately */
    /* following subroutine, as returned by ILAENV.) */
    if (*info == 0)
    {
        minwrk = 1;
        maxwrk = 1;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
         =8                      random symmetric   
         =9                      random nonsymmetric   
         =10                     random bidiagonal (log. distrib.) */

	    if (mtypes > 16) {
		goto L100;
	    }

	    itype = ktype[jtype - 1];
	    imode = kmode[jtype - 1];

/*           Compute norm */

	    switch (kmagn[jtype - 1]) {
		case 1:  goto L40;
		case 2:  goto L50;
		case 3:  goto L60;
	    }

L40:
	    anorm = 1.f;
	    goto L70;

L50:
	    anorm = rtovfl * ulp * amninv;
	    goto L70;

L60:
	    anorm = rtunfl * max(m,n) * ulpinv;
	    goto L70;

L70:

	    claset_("Full", lda, &n, &c_b1, &c_b1, &a[a_offset], lda);
	    iinfo = 0;
	    cond = ulpinv;

	    bidiag = FALSE_;
	    if (itype == 1) {

/*              Zero matrix */

		iinfo = 0;

	    } else if (itype == 2) {

/*              Identity */

		i__3 = mnmin;
		for (jcol = 1; jcol <= i__3; ++jcol) {
		    i__4 = a_subscr(jcol, jcol);
		    a[i__4].r = anorm, a[i__4].i = 0.f;
/* L80: */
		}

	    } else if (itype == 4) {

/*              Diagonal Matrix, [Eigen]values Specified */

		clatms_(&mnmin, &mnmin, "S", &iseed[1], "N", &rwork[1], &
			imode, &cond, &anorm, &c__0, &c__0, "N", &a[a_offset],
			 lda, &work[1], &iinfo);

	    } else if (itype == 5) {

/*              Symmetric, eigenvalues specified */

/* Subroutine */ int claqr2_(logical *wantt, logical *wantz, integer *n, 
	integer *ktop, integer *kbot, integer *nw, complex *h__, integer *ldh, 
	 integer *iloz, integer *ihiz, complex *z__, integer *ldz, integer *
	ns, integer *nd, complex *sh, complex *v, integer *ldv, integer *nh, 
	complex *t, integer *ldt, integer *nv, complex *wv, integer *ldwv, 
	complex *work, integer *lwork)
{
    /* System generated locals */
    integer h_dim1, h_offset, t_dim1, t_offset, v_dim1, v_offset, wv_dim1, 
	    wv_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4;
    real r__1, r__2, r__3, r__4, r__5, r__6;
    complex q__1, q__2;

    /* Builtin functions */
    double r_imag(complex *);
    void r_cnjg(complex *, complex *);

    /* Local variables */
    integer i__, j;
    complex s;
    integer jw;
    real foo;
    integer kln;
    complex tau;
    integer knt;
    real ulp;
    integer lwk1, lwk2;
    complex beta;
    integer kcol, info, ifst, ilst, ltop, krow;
    extern /* Subroutine */ int clarf_(char *, integer *, integer *, complex *
, integer *, complex *, complex *, integer *, complex *), 
	    cgemm_(char *, char *, integer *, integer *, integer *, complex *, 
	     complex *, integer *, complex *, integer *, complex *, complex *, 
	     integer *), ccopy_(integer *, complex *, integer 
	    *, complex *, integer *);
    integer infqr, kwtop;
    extern /* Subroutine */ int slabad_(real *, real *), cgehrd_(integer *, 
	    integer *, integer *, complex *, integer *, complex *, complex *, 
	    integer *, integer *), clarfg_(integer *, complex *, complex *, 
	    integer *, complex *);
    extern doublereal slamch_(char *);
    extern /* Subroutine */ int clahqr_(logical *, logical *, integer *, 
	    integer *, integer *, complex *, integer *, complex *, integer *, 
	    integer *, complex *, integer *, integer *), clacpy_(char *, 
	    integer *, integer *, complex *, integer *, complex *, integer *), claset_(char *, integer *, integer *, complex *, complex 
	    *, complex *, integer *);
    real safmin, safmax;
    extern /* Subroutine */ int ctrexc_(char *, integer *, complex *, integer 
	    *, complex *, integer *, integer *, integer *, integer *),
	     cunmhr_(char *, char *, integer *, integer *, integer *, integer 
	    *, complex *, integer *, complex *, complex *, integer *, complex 
	    *, integer *, integer *);
    real smlnum;
    integer lwkopt;


/*  -- LAPACK auxiliary routine (version 3.2.1)                        -- */
/*     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd.. */
/*  -- April 2009                                                      -- */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*     This subroutine is identical to CLAQR3 except that it avoids */
/*     recursion by calling CLAHQR instead of CLAQR4. */


/*     ****************************************************************** */
/*     Aggressive early deflation: */

/*     This subroutine accepts as input an upper Hessenberg matrix */
/*     H and performs an unitary similarity transformation */
/*     designed to detect and deflate fully converged eigenvalues from */
/*     a trailing principal submatrix.  On output H has been over- */
/*     written by a new Hessenberg matrix that is a perturbation of */
/*     an unitary similarity transformation of H.  It is to be */
/*     hoped that the final version of H has many zero subdiagonal */
/*     entries. */

/*     ****************************************************************** */
/*     WANTT   (input) LOGICAL */
/*          If .TRUE., then the Hessenberg matrix H is fully updated */
/*          so that the triangular Schur factor may be */
/*          computed (in cooperation with the calling subroutine). */
/*          If .FALSE., then only enough of H is updated to preserve */
/*          the eigenvalues. */

/*     WANTZ   (input) LOGICAL */
/*          If .TRUE., then the unitary matrix Z is updated so */
/*          so that the unitary Schur factor may be computed */
/*          (in cooperation with the calling subroutine). */
/*          If .FALSE., then Z is not referenced. */

/*     N       (input) INTEGER */
/*          The order of the matrix H and (if WANTZ is .TRUE.) the */
/*          order of the unitary matrix Z. */

/*     KTOP    (input) INTEGER */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    static integer c__2 = 2;
    
    /* System generated locals */
    integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4;
    real r__1, r__2;
    complex q__1;
    /* Builtin functions */
    double c_abs(complex *);
    /* Local variables */
    static real anrm, bnrm, smin, smax;
    static integer i__, j, iascl, ibscl;
    extern /* Subroutine */ int ccopy_(integer *, complex *, integer *, 
	    complex *, integer *);
    static integer ismin, ismax;
    static complex c1, c2;
    extern /* Subroutine */ int ctrsm_(char *, char *, char *, char *, 
	    integer *, integer *, complex *, complex *, integer *, complex *, 
	    integer *), claic1_(integer *, 
	    integer *, complex *, real *, complex *, complex *, real *, 
	    complex *, complex *);
    static real wsize;
    static complex s1, s2;
    extern /* Subroutine */ int cgeqp3_(integer *, integer *, complex *, 
	    integer *, integer *, complex *, complex *, integer *, real *, 
	    integer *);
    static integer nb;
    extern /* Subroutine */ int slabad_(real *, real *);
    extern doublereal clange_(char *, integer *, integer *, complex *, 
	    integer *, real *);
    static integer mn;
    extern /* Subroutine */ int clascl_(char *, integer *, integer *, real *, 
	    real *, integer *, integer *, complex *, integer *, integer *);
    extern doublereal slamch_(char *);
    extern /* Subroutine */ int claset_(char *, integer *, integer *, complex 
	    *, complex *, complex *, integer *), xerbla_(char *, 
	    integer *);
    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
	    integer *, integer *, ftnlen, ftnlen);
    static real bignum;
    static integer nb1, nb2, nb3, nb4;
    extern /* Subroutine */ int cunmqr_(char *, char *, integer *, integer *, 
	    integer *, complex *, integer *, complex *, complex *, integer *, 
	    complex *, integer *, integer *);
    static real sminpr, smaxpr, smlnum;
    extern /* Subroutine */ int cunmrz_(char *, char *, integer *, integer *, 
	    integer *, integer *, complex *, integer *, complex *, complex *, 
	    integer *, complex *, integer *, integer *);
    static integer lwkopt;
    static logical lquery;
    extern /* Subroutine */ int ctzrzf_(integer *, integer *, complex *, 
	    integer *, complex *, complex *, integer *, integer *);
#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)]


    a_dim1 = *lda;
    a_offset = 1 + a_dim1 * 1;
    a -= a_offset;
    b_dim1 = *ldb;
    b_offset = 1 + b_dim1 * 1;
    b -= b_offset;
    --jpvt;
    --work;
    --rwork;

//.........这里部分代码省略.........
    icols = *n + 1 - ilo;
    itau = iwork;
    iwork = itau + irows;
    i__1 = *lwork + 1 - iwork;
    cgeqrf_(&irows, &icols, &b[ilo + ilo * b_dim1], ldb, &work[itau], &work[
                iwork], &i__1, &iinfo);
    if (iinfo >= 0) {
        /* Computing MAX */
        i__3 = iwork;
        i__1 = lwkopt, i__2 = (integer) work[i__3].r + iwork - 1;
        lwkopt = max(i__1,i__2);
    }
    if (iinfo != 0) {
        *info = *n + 2;
        goto L10;
    }

    i__1 = *lwork + 1 - iwork;
    cunmqr_("L", "C", &irows, &icols, &irows, &b[ilo + ilo * b_dim1], ldb, &
            work[itau], &a[ilo + ilo * a_dim1], lda, &work[iwork], &i__1, &
            iinfo);
    if (iinfo >= 0) {
        /* Computing MAX */
        i__3 = iwork;
        i__1 = lwkopt, i__2 = (integer) work[i__3].r + iwork - 1;
        lwkopt = max(i__1,i__2);
    }
    if (iinfo != 0) {
        *info = *n + 3;
        goto L10;
    }

    if (ilvsl) {
        claset_("Full", n, n, &c_b1, &c_b2, &vsl[vsl_offset], ldvsl);
        i__1 = irows - 1;
        i__2 = irows - 1;
        clacpy_("L", &i__1, &i__2, &b[ilo + 1 + ilo * b_dim1], ldb, &vsl[ilo
                + 1 + ilo * vsl_dim1], ldvsl);
        i__1 = *lwork + 1 - iwork;
        cungqr_(&irows, &irows, &irows, &vsl[ilo + ilo * vsl_dim1], ldvsl, &
                work[itau], &work[iwork], &i__1, &iinfo);
        if (iinfo >= 0) {
            /* Computing MAX */
            i__3 = iwork;
            i__1 = lwkopt, i__2 = (integer) work[i__3].r + iwork - 1;
            lwkopt = max(i__1,i__2);
        }
        if (iinfo != 0) {
            *info = *n + 4;
            goto L10;
        }
    }

    if (ilvsr) {
        claset_("Full", n, n, &c_b1, &c_b2, &vsr[vsr_offset], ldvsr);
    }

    /*     Reduce to generalized Hessenberg form */

    cgghrd_(jobvsl, jobvsr, n, &ilo, &ihi, &a[a_offset], lda, &b[b_offset],
            ldb, &vsl[vsl_offset], ldvsl, &vsr[vsr_offset], ldvsr, &iinfo);
    if (iinfo != 0) {
        *info = *n + 5;
        goto L10;
    }

//.........这里部分代码省略.........
/*     Determine EPS and the norm of A. */

    eps = slamch_("Epsilon");
    anorm = clanhp_("1", uplo, n, &a[1], &rwork[1]);

/*     Check the imaginary parts of the diagonal elements and return with */
/*     an error code if any are nonzero. */

    jc = 1;
    if (lsame_(uplo, "U")) {
	i__1 = *n;
	for (j = 1; j <= i__1; ++j) {
	    if (r_imag(&afac[jc]) != 0.f) {
		*resid = 1.f / eps;
		return 0;
	    }
	    jc = jc + j + 1;
/* L10: */
	}
    } else {
	i__1 = *n;
	for (j = 1; j <= i__1; ++j) {
	    if (r_imag(&afac[jc]) != 0.f) {
		*resid = 1.f / eps;
		return 0;
	    }
	    jc = jc + *n - j + 1;
/* L20: */
	}
    }

/*     Initialize C to the identity matrix. */

    claset_("Full", n, n, &c_b1, &c_b2, &c__[c_offset], ldc);

/*     Call CLAVHP to form the product D * U' (or D * L' ). */

    clavhp_(uplo, "Conjugate", "Non-unit", n, n, &afac[1], &ipiv[1], &c__[
	    c_offset], ldc, &info);

/*     Call CLAVHP again to multiply by U ( or L ). */

    clavhp_(uplo, "No transpose", "Unit", n, n, &afac[1], &ipiv[1], &c__[
	    c_offset], ldc, &info);

/*     Compute the difference  C - A . */

    if (lsame_(uplo, "U")) {
	jc = 0;
	i__1 = *n;
	for (j = 1; j <= i__1; ++j) {
	    i__2 = j - 1;
	    for (i__ = 1; i__ <= i__2; ++i__) {
		i__3 = i__ + j * c_dim1;
		i__4 = i__ + j * c_dim1;
		i__5 = jc + i__;
		q__1.r = c__[i__4].r - a[i__5].r, q__1.i = c__[i__4].i - a[
			i__5].i;
		c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
/* L30: */
	    }
	    i__2 = j + j * c_dim1;
	    i__3 = j + j * c_dim1;
	    i__4 = jc + j;
	    r__1 = a[i__4].r;
	    q__1.r = c__[i__3].r - r__1, q__1.i = c__[i__3].i;

//.........这里部分代码省略.........
       Function Body */
    /* Table of constant values */
    static complex c_b1 = {0.f,0.f};
    static complex c_b2 = {1.f,0.f};
    static integer c__1 = 1;
    
    /* System generated locals */
    integer ab_dim1, ab_offset, c_dim1, c_offset, pt_dim1, pt_offset, q_dim1, 
	    q_offset, i__1, i__2, i__3, i__4, i__5, i__6, i__7;
    complex q__1, q__2, q__3;
    /* Builtin functions */
    void r_cnjg(complex *, complex *);
    double c_abs(complex *);
    /* Local variables */
    static integer inca;
    static real abst;
    extern /* Subroutine */ int crot_(integer *, complex *, integer *, 
	    complex *, integer *, real *, complex *);
    static integer i, j, l;
    static complex t;
    extern /* Subroutine */ int cscal_(integer *, complex *, complex *, 
	    integer *);
    extern logical lsame_(char *, char *);
    static logical wantb, wantc;
    static integer minmn;
    static logical wantq;
    static integer j1, j2, kb;
    static complex ra;
    static real rc;
    static integer kk;
    static complex rb;
    static integer ml, nr, mu;
    static complex rs;
    extern /* Subroutine */ int claset_(char *, integer *, integer *, complex 
	    *, complex *, complex *, integer *), clartg_(complex *, 
	    complex *, real *, complex *, complex *), xerbla_(char *, integer 
	    *), clargv_(integer *, complex *, integer *, complex *, 
	    integer *, real *, integer *), clartv_(integer *, complex *, 
	    integer *, complex *, integer *, real *, complex *, integer *);
    static integer kb1, ml0;
    static logical wantpt;
    static integer mu0, klm, kun, nrt, klu1;



#define D(I) d[(I)-1]
#define E(I) e[(I)-1]
#define WORK(I) work[(I)-1]
#define RWORK(I) rwork[(I)-1]

#define AB(I,J) ab[(I)-1 + ((J)-1)* ( *ldab)]
#define Q(I,J) q[(I)-1 + ((J)-1)* ( *ldq)]
#define PT(I,J) pt[(I)-1 + ((J)-1)* ( *ldpt)]
#define C(I,J) c[(I)-1 + ((J)-1)* ( *ldc)]

    wantb = lsame_(vect, "B");
    wantq = lsame_(vect, "Q") || wantb;
    wantpt = lsame_(vect, "P") || wantb;
    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) {

//.........这里部分代码省略.........
       Parameter adjustments */
    l_dim1 = *lda;
    l_offset = 1 + l_dim1 * 1;
    l -= l_offset;
    q_dim1 = *lda;
    q_offset = 1 + q_dim1 * 1;
    q -= q_offset;
    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;
    --rwork;
    --result;

    /* Function Body */
    minmn = min(*m,*n);
    eps = slamch_("Epsilon");

/*     Copy the matrix A to the array AF. */

    clacpy_("Full", m, n, &a[a_offset], lda, &af[af_offset], lda);

/*     Factorize the matrix A in the array AF. */

    s_copy(srnamc_1.srnamt, "CGEQLF", (ftnlen)6, (ftnlen)6);
    cgeqlf_(m, n, &af[af_offset], lda, &tau[1], &work[1], lwork, &info);

/*     Copy details of Q */

    claset_("Full", m, m, &c_b1, &c_b1, &q[q_offset], lda);
    if (*m >= *n) {
	if (*n < *m && *n > 0) {
	    i__1 = *m - *n;
	    clacpy_("Full", &i__1, n, &af[af_offset], lda, &q_ref(1, *m - *n 
		    + 1), lda);
	}
	if (*n > 1) {
	    i__1 = *n - 1;
	    i__2 = *n - 1;
	    clacpy_("Upper", &i__1, &i__2, &af_ref(*m - *n + 1, 2), lda, &
		    q_ref(*m - *n + 1, *m - *n + 2), lda);
	}
    } else {
	if (*m > 1) {
	    i__1 = *m - 1;
	    i__2 = *m - 1;
	    clacpy_("Upper", &i__1, &i__2, &af_ref(1, *n - *m + 2), lda, &
		    q_ref(1, 2), lda);
	}
    }

/*     Generate the m-by-m matrix Q */

    s_copy(srnamc_1.srnamt, "CUNGQL", (ftnlen)6, (ftnlen)6);
    cungql_(m, m, &minmn, &q[q_offset], lda, &tau[1], &work[1], lwork, &info);

/*     Copy L */

    claset_("Full", m, n, &c_b12, &c_b12, &l[l_offset], lda);
    if (*m >= *n) {
	if (*n > 0) {
	    clacpy_("Lower", n, n, &af_ref(*m - *n + 1, 1), lda, &l_ref(*m - *

//.........这里部分代码省略.........
    *ninfo = 0;

/*     Read input data until N=0 */

L10:
    io___2.ciunit = *nin;
    s_rsle(&io___2);
    do_lio(&c__3, &c__1, (char *)&n, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&ifst, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&ilst, (ftnlen)sizeof(integer));
    e_rsle();
    if (n == 0) {
	return 0;
    }
    ++(*knt);
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	io___7.ciunit = *nin;
	s_rsle(&io___7);
	i__2 = n;
	for (j = 1; j <= i__2; ++j) {
	    do_lio(&c__6, &c__1, (char *)&tmp_ref(i__, j), (ftnlen)sizeof(
		    complex));
	}
	e_rsle();
/* L20: */
    }
    clacpy_("F", &n, &n, tmp, &c__10, t1, &c__10);
    clacpy_("F", &n, &n, tmp, &c__10, t2, &c__10);
    res = 0.f;

/*     Test without accumulating Q */

    claset_("Full", &n, &n, &c_b1, &c_b2, q, &c__10);
    ctrexc_("N", &n, t1, &c__10, q, &c__10, &ifst, &ilst, &info1);
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	i__2 = n;
	for (j = 1; j <= i__2; ++j) {
	    i__3 = q_subscr(i__, j);
	    if (i__ == j && (q[i__3].r != 1.f || q[i__3].i != 0.f)) {
		res += 1.f / eps;
	    }
	    i__3 = q_subscr(i__, j);
	    if (i__ != j && (q[i__3].r != 0.f || q[i__3].i != 0.f)) {
		res += 1.f / eps;
	    }
/* L30: */
	}
/* L40: */
    }

/*     Test with accumulating Q */

    claset_("Full", &n, &n, &c_b1, &c_b2, q, &c__10);
    ctrexc_("V", &n, t2, &c__10, q, &c__10, &ifst, &ilst, &info2);

/*     Compare T1 with T2 */

    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	i__2 = n;
	for (j = 1; j <= i__2; ++j) {
	    i__3 = t1_subscr(i__, j);
	    i__4 = t2_subscr(i__, j);
	    if (t1[i__3].r != t2[i__4].r || t1[i__3].i != t2[i__4].i) {

/* Subroutine */ int cgrqts_(integer *m, integer *p, integer *n, complex *a, 
	complex *af, complex *q, complex *r__, integer *lda, complex *taua, 
	complex *b, complex *bf, complex *z__, complex *t, complex *bwk, 
	integer *ldb, complex *taub, complex *work, integer *lwork, real *
	rwork, real *result)
{
    /* System generated locals */
    integer a_dim1, a_offset, af_dim1, af_offset, r_dim1, r_offset, q_dim1, 
	    q_offset, b_dim1, b_offset, bf_dim1, bf_offset, t_dim1, t_offset, 
	    z_dim1, z_offset, bwk_dim1, bwk_offset, i__1, i__2;
    real r__1;
    complex q__1;

    /* Local variables */
    static integer info;
    static real unfl;
    extern /* Subroutine */ int cgemm_(char *, char *, integer *, integer *, 
	    integer *, complex *, complex *, integer *, complex *, integer *, 
	    complex *, complex *, integer *), cherk_(char *, 
	    char *, integer *, integer *, real *, complex *, integer *, real *
	    , complex *, integer *);
    static real resid, anorm, bnorm;
    extern doublereal clange_(char *, integer *, integer *, complex *, 
	    integer *, real *), clanhe_(char *, char *, integer *, 
	    complex *, integer *, real *), slamch_(char *);
    extern /* Subroutine */ int cggrqf_(integer *, integer *, integer *, 
	    complex *, integer *, complex *, complex *, integer *, complex *, 
	    complex *, integer *, integer *), clacpy_(char *, integer *, 
	    integer *, complex *, integer *, complex *, integer *), 
	    claset_(char *, integer *, integer *, complex *, complex *, 
	    complex *, integer *), cungqr_(integer *, integer *, 
	    integer *, complex *, integer *, complex *, complex *, integer *, 
	    integer *), cungrq_(integer *, integer *, integer *, complex *, 
	    integer *, complex *, complex *, integer *, integer *);
    static real 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   
    =======   

    CGRQTS tests CGGRQF, 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 array, dimension (LDA,N)   
            The M-by-N matrix A.   

    AF      (output) COMPLEX array, dimension (LDA,N)   
            Details of the GRQ factorization of A and B, as returned   
            by CGGRQF, see CGGRQF for further details.   

    Q       (output) COMPLEX array, dimension (LDA,N)   
            The N-by-N unitary matrix Q.   

    R       (workspace) COMPLEX 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 array, dimension (min(M,N))   
            The scalar factors of the elementary reflectors, as returned   
            by SGGQRC.   

    B       (input) COMPLEX array, dimension (LDB,N)   
            On entry, the P-by-N matrix A.   

    BF      (output) COMPLEX array, dimension (LDB,N)   
            Details of the GQR factorization of A and B, as returned   
            by CGGRQF, see CGGRQF for further details.   

//.........这里部分代码省略.........