doublereal dtzt02_(integer *m, integer *n, doublereal *af, integer *lda, 
	doublereal *tau, doublereal *work, integer *lwork)
{
    /* System generated locals */
    integer af_dim1, af_offset, i__1, i__2;
    doublereal ret_val;

    /* Local variables */
    static integer i__;
    static doublereal rwork[1];
    extern doublereal dlamch_(char *), dlange_(char *, integer *, 
	    integer *, doublereal *, integer *, doublereal *);
    extern /* Subroutine */ int dlaset_(char *, integer *, integer *, 
	    doublereal *, doublereal *, doublereal *, integer *), 
	    xerbla_(char *, integer *), dlatzm_(char *, integer *, 
	    integer *, doublereal *, integer *, doublereal *, doublereal *, 
	    doublereal *, integer *, doublereal *);


#define af_ref(a_1,a_2) af[(a_2)*af_dim1 + a_1]


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

    DTZT02 returns   
         || I - Q'*Q || / ( M * eps)   
    where the matrix Q is defined by the Householder transformations   
    generated by DTZRQF.   

    Arguments   
    =========   

    M       (input) INTEGER   
            The number of rows of the matrix AF.   

    N       (input) INTEGER   
            The number of columns of the matrix AF.   

    AF      (input) DOUBLE PRECISION array, dimension (LDA,N)   
            The output of DTZRQF.   

    LDA     (input) INTEGER   
            The leading dimension of the array AF.   

    TAU     (input) DOUBLE PRECISION array, dimension (M)   
            Details of the Householder transformations as returned by   
            DTZRQF.   

    WORK    (workspace) DOUBLE PRECISION array, dimension (LWORK)   

    LWORK   (input) INTEGER   
            length of WORK array. Must be >= N*N+N   

    =====================================================================   


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

    /* Function Body */
    ret_val = 0.;

    if (*lwork < *n * *n + *n) {
	xerbla_("DTZT02", &c__7);
	return ret_val;
    }

/*     Quick return if possible */

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

/*     Q := I */

    dlaset_("Full", n, n, &c_b5, &c_b6, &work[1], n);

/*     Q := P(1) * ... * P(m) * Q */

    for (i__ = *m; i__ >= 1; --i__) {
	i__1 = *n - *m + 1;
	dlatzm_("Left", &i__1, n, &af_ref(i__, *m + 1), lda, &tau[i__], &work[
		i__], &work[*m + 1], n, &work[*n * *n + 1]);
/* L10: */
    }

/*     Q := P(m) * ... * P(1) * Q */

    i__1 = *m;
//.........这里部分代码省略.........

 int zlaqgb_(int *m, int *n, int *kl, int *ku, 
	 doublecomplex *ab, int *ldab, double *r__, double *c__, 
	double *rowcnd, double *colcnd, double *amax, char *equed)
{
    /* System generated locals */
    int ab_dim1, ab_offset, i__1, i__2, i__3, i__4, i__5, i__6;
    double d__1;
    doublecomplex z__1;

    /* Local variables */
    int i__, j;
    double cj, large, small;
    extern double dlamch_(char *);


/*  -- LAPACK auxiliary routine (version 3.2) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

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

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

/*  ZLAQGB equilibrates a general M by N band matrix A with KL */
/*  subdiagonals and KU superdiagonals using the row and scaling factors */
/*  in the vectors R and C. */

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

/*  M       (input) INTEGER */
/*          The number of rows of the matrix A.  M >= 0. */

/*  N       (input) INTEGER */
/*          The number of columns of the matrix A.  N >= 0. */

/*  KL      (input) INTEGER */
/*          The number of subdiagonals within the band of A.  KL >= 0. */

/*  KU      (input) INTEGER */
/*          The number of superdiagonals within the band of A.  KU >= 0. */

/*  AB      (input/output) COMPLEX*16 array, dimension (LDAB,N) */
/*          On entry, the matrix A in band storage, in rows 1 to KL+KU+1. */
/*          The j-th column of A is stored in the j-th column of the */
/*          array AB as follows: */
/*          AB(ku+1+i-j,j) = A(i,j) for MAX(1,j-ku)<=i<=MIN(m,j+kl) */

/*          On exit, the equilibrated matrix, in the same storage format */
/*          as A.  See EQUED for the form of the equilibrated matrix. */

/*  LDAB    (input) INTEGER */
/*          The leading dimension of the array AB.  LDA >= KL+KU+1. */

/*  R       (input) DOUBLE PRECISION array, dimension (M) */
/*          The row scale factors for A. */

/*  C       (input) DOUBLE PRECISION array, dimension (N) */
/*          The column scale factors for A. */

/*  ROWCND  (input) DOUBLE PRECISION */
/*          Ratio of the smallest R(i) to the largest R(i). */

/*  COLCND  (input) DOUBLE PRECISION */
/*          Ratio of the smallest C(i) to the largest C(i). */

/*  AMAX    (input) DOUBLE PRECISION */
/*          Absolute value of largest matrix entry. */

/*  EQUED   (output) CHARACTER*1 */
/*          Specifies the form of equilibration that was done. */
/*          = 'N':  No equilibration */
/*          = 'R':  Row equilibration, i.e., A has been premultiplied by */
/*                  diag(R). */
/*          = 'C':  Column equilibration, i.e., A has been postmultiplied */
/*                  by diag(C). */
/*          = 'B':  Both row and column equilibration, i.e., A has been */
/*                  replaced by diag(R) * A * diag(C). */

/*  Internal Parameters */
/*  =================== */

/*  THRESH is a threshold value used to decide if row or column scaling */
/*  should be done based on the ratio of the row or column scaling */
/*  factors.  If ROWCND < THRESH, row scaling is done, and if */
/*  COLCND < THRESH, column scaling is done. */

/*  LARGE and SMALL are threshold values used to decide if row scaling */
/*  should be done based on the absolute size of the largest matrix */
/*  element.  If AMAX > LARGE or AMAX < SMALL, row scaling is done. */

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

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            /* Permute the global matrix GA for symbfact() */
            for (i = 0; i < colptr[n]; ++i) {
                irow = rowind[i];
                rowind[i] = perm_r[irow];
            }
        } else if ( !factored && Fact != SamePattern_SameRowPerm ) {
            /* Get a new perm_r[] */
            if ( job == 5 ) {
                /* Allocate storage for scaling factors. */
                if ( !(R1 = (double *) SUPERLU_MALLOC(m * sizeof(double))) )
                    ABORT("SUPERLU_MALLOC fails for R1[]");
                if ( !(C1 = (double *) SUPERLU_MALLOC(n * sizeof(double))) )
                    ABORT("SUPERLU_MALLOC fails for C1[]");
            }

            if ( !iam ) {
                /* Process 0 finds a row permutation for large diagonal. */
                dldperm(job, m, nnz, colptr, rowind, a_GA, perm_r, R1, C1);

                MPI_Bcast( perm_r, m, mpi_int_t, 0, grid->comm );
                if ( job == 5 && Equil ) {
                    MPI_Bcast( R1, m, MPI_DOUBLE, 0, grid->comm );
                    MPI_Bcast( C1, n, MPI_DOUBLE, 0, grid->comm );
                }
            } else {
                MPI_Bcast( perm_r, m, mpi_int_t, 0, grid->comm );
                if ( job == 5 && Equil ) {
                    MPI_Bcast( R1, m, MPI_DOUBLE, 0, grid->comm );
                    MPI_Bcast( C1, n, MPI_DOUBLE, 0, grid->comm );
                }
            }

#if ( PRNTlevel>=2 )
            dmin = dlamch_("Overflow");
            dsum = 0.0;
            dprod = 1.0;
#endif
            if ( job == 5 ) {
                if ( Equil ) {
                    for (i = 0; i < n; ++i) {
                        R1[i] = exp(R1[i]);
                        C1[i] = exp(C1[i]);
                    }

                    /* Permute the global matrix GA for symbfact(). */
                    for (j = 0; j < n; ++j) {
                        for (i = colptr[j]; i < colptr[j+1]; ++i) {
                            irow = rowind[i];
                            rowind[i] = perm_r[irow];
#if ( PRNTlevel>=2 )
                            if ( rowind[i] == j ) /* New diagonal */
                                dprod *= fabs(a[i]);
#endif
                        }
                    }

                    /* Scale the distributed matrix */
                    irow = fst_row;
                    for (j = 0; j < m_loc; ++j) {
                        for (i = rowptr[j]; i < rowptr[j+1]; ++i) {
                            icol = colind[i];
                            a[i] *= R1[irow] * C1[icol];
                        }
                        ++irow;
                    }

/*<        >*/
/* Subroutine */ int dsapps_(integer *n, integer *kev, integer *np,
        doublereal *shift, doublereal *v, integer *ldv, doublereal *h__,
        integer *ldh, doublereal *resid, doublereal *q, integer *ldq,
        doublereal *workd)
{
    /* Initialized data */

    static logical first = TRUE_;

    /* System generated locals */
    integer h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
            i__3, i__4;
    doublereal d__1, d__2;

    /* Local variables */
    doublereal c__, f, g;
    integer i__, j;
    doublereal r__, s, a1, a2, a3, a4;
/*  static real t0, t1; */
    integer jj;
    doublereal big;
    integer iend, itop;
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
            integer *), dgemv_(char *, integer *, integer *, doublereal *,
            doublereal *, integer *, doublereal *, integer *, doublereal *,
            doublereal *, integer *, ftnlen), dcopy_(integer *, doublereal *,
            integer *, doublereal *, integer *), daxpy_(integer *, doublereal
            *, doublereal *, integer *, doublereal *, integer *);
    extern doublereal dlamch_(char *, ftnlen);
    extern /* Subroutine */ int second_(real *);
    static doublereal epsmch;
    integer istart, kplusp /*, msglvl */;
    extern /* Subroutine */ int dlacpy_(char *, integer *, integer *,
            doublereal *, integer *, doublereal *, integer *, ftnlen),
            dlartg_(doublereal *, doublereal *, doublereal *, doublereal *,
            doublereal *), dlaset_(char *, integer *, integer *, doublereal *,
             doublereal *, doublereal *, integer *, ftnlen);


/*     %----------------------------------------------------% */
/*     | Include files for debugging and timing information | */
/*     %----------------------------------------------------% */

/*<       include   'debug.h' >*/
/*<       include   'stat.h' >*/

/* \SCCS Information: @(#) */
/* FILE: debug.h   SID: 2.3   DATE OF SID: 11/16/95   RELEASE: 2 */

/*     %---------------------------------% */
/*     | See debug.doc for documentation | */
/*     %---------------------------------% */
/*<        >*/
/*<       integer    kev, ldh, ldq, ldv, n, np >*/

/*     %------------------% */
/*     | Scalar Arguments | */
/*     %------------------% */

/*     %--------------------------------% */
/*     | See stat.doc for documentation | */
/*     %--------------------------------% */

/* \SCCS Information: @(#) */
/* FILE: stat.h   SID: 2.2   DATE OF SID: 11/16/95   RELEASE: 2 */

/*<       save       t0, t1, t2, t3, t4, t5 >*/

/*<       integer    nopx, nbx, nrorth, nitref, nrstrt >*/
/*<        >*/
/*<        >*/

/*     %-----------------% */
/*     | Array Arguments | */
/*     %-----------------% */

/*<        >*/

/*     %------------% */
/*     | Parameters | */
/*     %------------% */

/*<        >*/
/*<       parameter (one = 1.0D+0, zero = 0.0D+0) >*/

/*     %---------------% */
/*     | Local Scalars | */
/*     %---------------% */

/*<       integer    i, iend, istart, itop, j, jj, kplusp, msglvl >*/
/*<       logical    first >*/
/*<        >*/
/*<       save       epsmch, first >*/


/*     %----------------------% */
/*     | External Subroutines | */
/*     %----------------------% */

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

/* Subroutine */ int dlanv2_(doublereal *a, doublereal *b, doublereal *c__, 
	doublereal *d__, doublereal *rt1r, doublereal *rt1i, doublereal *rt2r,
	 doublereal *rt2i, doublereal *cs, doublereal *sn)
{
/*  -- LAPACK driver routine (version 3.1) --   
       Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..   
       November 2006   


    Purpose   
    =======   

    DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric   
    matrix in standard form:   

         [ A  B ] = [ CS -SN ] [ AA  BB ] [ CS  SN ]   
         [ C  D ]   [ SN  CS ] [ CC  DD ] [-SN  CS ]   

    where either   
    1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or   
    2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex   
    conjugate eigenvalues.   

    Arguments   
    =========   

    A       (input/output) DOUBLE PRECISION   
    B       (input/output) DOUBLE PRECISION   
    C       (input/output) DOUBLE PRECISION   
    D       (input/output) DOUBLE PRECISION   
            On entry, the elements of the input matrix.   
            On exit, they are overwritten by the elements of the   
            standardised Schur form.   

    RT1R    (output) DOUBLE PRECISION   
    RT1I    (output) DOUBLE PRECISION   
    RT2R    (output) DOUBLE PRECISION   
    RT2I    (output) DOUBLE PRECISION   
            The real and imaginary parts of the eigenvalues. If the   
            eigenvalues are a complex conjugate pair, RT1I > 0.   

    CS      (output) DOUBLE PRECISION   
    SN      (output) DOUBLE PRECISION   
            Parameters of the rotation matrix.   

    Further Details   
    ===============   

    Modified by V. Sima, Research Institute for Informatics, Bucharest,   
    Romania, to reduce the risk of cancellation errors,   
    when computing real eigenvalues, and to ensure, if possible, that   
    abs(RT1R) >= abs(RT2R).   

    ===================================================================== */
    /* Table of constant values */
    static doublereal c_b4 = 1.;
    
    /* System generated locals */
    doublereal d__1, d__2;
    /* Builtin functions */
    double d_sign(doublereal *, doublereal *), sqrt(doublereal);
    /* Local variables */
    static doublereal p, z__, aa, bb, cc, dd, cs1, sn1, sab, sac, eps, tau, 
	    temp, scale, bcmax, bcmis, sigma;
    extern doublereal dlapy2_(doublereal *, doublereal *), dlamch_(char *);




    eps = dlamch_("P");
    if (*c__ == 0.) {
	*cs = 1.;
	*sn = 0.;
	goto L10;

    } else if (*b == 0.) {

/*        Swap rows and columns */

	*cs = 0.;
	*sn = 1.;
	temp = *d__;
	*d__ = *a;
	*a = temp;
	*b = -(*c__);
	*c__ = 0.;
	goto L10;
    } else if (*a - *d__ == 0. && d_sign(&c_b4, b) != d_sign(&c_b4, c__)) {
	*cs = 1.;
	*sn = 0.;
	goto L10;
    } else {

	temp = *a - *d__;
	p = temp * .5;
/* Computing MAX */
	d__1 = abs(*b), d__2 = abs(*c__);
	bcmax = max(d__1,d__2);
/* Computing MIN */
	d__1 = abs(*b), d__2 = abs(*c__);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    /* System generated locals */
    integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, qe_dim1, 
	    qe_offset, s_dim1, s_offset, t_dim1, t_offset, z_dim1, z_offset, 
	    i__1, i__2, i__3, i__4, i__5, i__6, i__7;
    doublereal d__1, d__2;
    doublecomplex z__1, z__2, z__3;

    /* Builtin functions */
    double d_sign(doublereal *, doublereal *), z_abs(doublecomplex *);
    void d_cnjg(doublecomplex *, doublecomplex *);
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

    /* Local variables */
    static integer iadd, ierr, nmax, i__, j, n;
    static logical badnn;
    static doublereal rmagn[4];
    static doublecomplex ctemp;
    extern /* Subroutine */ int zget52_(logical *, integer *, doublecomplex *,
	     integer *, doublecomplex *, integer *, doublecomplex *, integer *
	    , doublecomplex *, doublecomplex *, doublecomplex *, doublereal *,
	     doublereal *);
    static integer nmats, jsize;
    extern /* Subroutine */ int zggev_(char *, char *, integer *, 
	    doublecomplex *, integer *, doublecomplex *, integer *, 
	    doublecomplex *, doublecomplex *, doublecomplex *, integer *, 
	    doublecomplex *, integer *, doublecomplex *, integer *, 
	    doublereal *, integer *);
    static integer nerrs, jtype, n1;
    extern /* Subroutine */ int dlabad_(doublereal *, doublereal *), zlatm4_(
	    integer *, integer *, integer *, integer *, logical *, doublereal 
	    *, doublereal *, doublereal *, integer *, integer *, 
	    doublecomplex *, integer *);
    static integer jc, nb, in;
    extern doublereal dlamch_(char *);
    static integer jr;
    extern /* Subroutine */ int zunm2r_(char *, char *, integer *, integer *, 
	    integer *, doublecomplex *, integer *, doublecomplex *, 
	    doublecomplex *, integer *, doublecomplex *, integer *);
    static doublereal safmin, safmax;
    static integer ioldsd[4];
    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
	    integer *, integer *, ftnlen, ftnlen);
    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 *);
    static integer minwrk, maxwrk;
    static doublereal ulpinv;
    static integer mtypes, ntestt;
    static doublereal ulp;

    /* Fortran I/O blocks */
    static cilist io___40 = { 0, 0, 0, fmt_9999, 0 };
    static cilist io___42 = { 0, 0, 0, fmt_9999, 0 };
    static cilist io___43 = { 0, 0, 0, fmt_9998, 0 };
    static cilist io___44 = { 0, 0, 0, fmt_9998, 0 };
    static cilist io___45 = { 0, 0, 0, fmt_9999, 0 };
    static cilist io___46 = { 0, 0, 0, fmt_9999, 0 };
    static cilist io___47 = { 0, 0, 0, fmt_9999, 0 };
    static cilist io___48 = { 0, 0, 0, fmt_9997, 0 };

//.........这里部分代码省略.........
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

/*     Quick exit if N = 0 or NRHS = 0. */

    /* Parameter adjustments */
    ab_dim1 = *ldab;
    ab_offset = 1 + ab_dim1;
    ab -= ab_offset;
    b_dim1 = *ldb;
    b_offset = 1 + b_dim1;
    b -= b_offset;
    x_dim1 = *ldx;
    x_offset = 1 + x_dim1;
    x -= x_offset;
    xact_dim1 = *ldxact;
    xact_offset = 1 + xact_dim1;
    xact -= xact_offset;
    --ferr;
    --berr;
    --reslts;

    /* Function Body */
    if (*n <= 0 || *nrhs <= 0) {
	reslts[1] = 0.;
	reslts[2] = 0.;
	return 0;
    }

    eps = dlamch_("Epsilon");
    unfl = dlamch_("Safe minimum");
    ovfl = 1. / unfl;
    notran = lsame_(trans, "N");
/* Computing MIN */
    i__1 = *kl + *ku + 2, i__2 = *n + 1;
    nz = min(i__1,i__2);

/*     Test 1:  Compute the maximum of */
/*        norm(X - XACT) / ( norm(X) * FERR ) */
/*     over all the vectors X and XACT using the infinity-norm. */

    errbnd = 0.;
    i__1 = *nrhs;
    for (j = 1; j <= i__1; ++j) {
	imax = idamax_(n, &x[j * x_dim1 + 1], &c__1);
/* Computing MAX */
	d__2 = (d__1 = x[imax + j * x_dim1], abs(d__1));
	xnorm = max(d__2,unfl);
	diff = 0.;
	i__2 = *n;
	for (i__ = 1; i__ <= i__2; ++i__) {
/* Computing MAX */
	    d__2 = diff, d__3 = (d__1 = x[i__ + j * x_dim1] - xact[i__ + j * 
		    xact_dim1], abs(d__1));
	    diff = max(d__2,d__3);
/* L10: */
	}

	if (xnorm > 1.) {
	    goto L20;
	} else if (diff <= ovfl * xnorm) {
	    goto L20;

/* Subroutine */ int zlaqhp_(char *uplo, integer *n, doublecomplex *ap, 
	doublereal *s, doublereal *scond, doublereal *amax, char *equed)
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4;
    doublereal d__1;
    doublecomplex z__1;

    /* Local variables */
    integer i__, j, jc;
    doublereal cj, large;
    extern logical lsame_(char *, char *);
    doublereal small;
    extern doublereal dlamch_(char *);


/*  -- LAPACK auxiliary routine (version 3.1) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

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

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

/*  ZLAQHP equilibrates a Hermitian matrix A using the scaling factors */
/*  in the vector S. */

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

/*  AP      (input/output) COMPLEX*16 array, dimension (N*(N+1)/2) */
/*          On entry, the upper or lower triangle of the Hermitian matrix */
/*          A, packed columnwise in a linear array.  The j-th column of A */
/*          is stored in the array AP as follows: */
/*          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
/*          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */

/*          On exit, the equilibrated matrix:  diag(S) * A * diag(S), in */
/*          the same storage format as A. */

/*  S       (input) DOUBLE PRECISION array, dimension (N) */
/*          The scale factors for A. */

/*  SCOND   (input) DOUBLE PRECISION */
/*          Ratio of the smallest S(i) to the largest S(i). */

/*  AMAX    (input) DOUBLE PRECISION */
/*          Absolute value of largest matrix entry. */

/*  EQUED   (output) CHARACTER*1 */
/*          Specifies whether or not equilibration was done. */
/*          = 'N':  No equilibration. */
/*          = 'Y':  Equilibration was done, i.e., A has been replaced by */
/*                  diag(S) * A * diag(S). */

/*  Internal Parameters */
/*  =================== */

/*  THRESH is a threshold value used to decide if scaling should be done */
/*  based on the ratio of the scaling factors.  If SCOND < THRESH, */
/*  scaling is done. */

/*  LARGE and SMALL are threshold values used to decide if scaling should */
/*  be done based on the absolute size of the largest matrix element. */
/*  If AMAX > LARGE or AMAX < SMALL, scaling is done. */

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

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

/*     Quick return if possible */

    /* Parameter adjustments */
    --s;
    --ap;

    /* Function Body */
    if (*n <= 0) {
	*(unsigned char *)equed = 'N';
//.........这里部分代码省略.........

 int dgerfs_(char *trans, int *n, int *nrhs, 
	double *a, int *lda, double *af, int *ldaf, int *
	ipiv, double *b, int *ldb, double *x, int *ldx, 
	double *ferr, double *berr, double *work, int *iwork, 
	int *info)
{
    /* System generated locals */
    int a_dim1, a_offset, af_dim1, af_offset, b_dim1, b_offset, x_dim1, 
	    x_offset, i__1, i__2, i__3;
    double d__1, d__2, d__3;

    /* Local variables */
    int i__, j, k;
    double s, xk;
    int nz;
    double eps;
    int kase;
    double safe1, safe2;
    extern int lsame_(char *, char *);
    extern  int dgemv_(char *, int *, int *, 
	    double *, double *, int *, double *, int *, 
	    double *, double *, int *);
    int isave[3];
    extern  int dcopy_(int *, double *, int *, 
	    double *, int *), daxpy_(int *, double *, 
	    double *, int *, double *, int *);
    int count;
    extern  int dlacn2_(int *, double *, double *, 
	     int *, double *, int *, int *);
    extern double dlamch_(char *);
    double safmin;
    extern  int xerbla_(char *, int *), dgetrs_(
	    char *, int *, int *, double *, int *, int *, 
	    double *, int *, int *);
    int notran;
    char transt[1];
    double lstres;


/*  -- LAPACK routine (version 3.2) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH. */

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

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

/*  DGERFS improves the computed solution to a system of linear */
/*  equations and provides error bounds and backward error estimates for */
/*  the solution. */

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

/*  TRANS   (input) CHARACTER*1 */
/*          Specifies the form of the system of equations: */
/*          = 'N':  A * X = B     (No transpose) */
/*          = 'T':  A**T * X = B  (Transpose) */
/*          = 'C':  A**H * X = B  (Conjugate transpose = Transpose) */

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

/*  NRHS    (input) INTEGER */
/*          The number of right hand sides, i.e., the number of columns */
/*          of the matrices B and X.  NRHS >= 0. */

/*  A       (input) DOUBLE PRECISION array, dimension (LDA,N) */
/*          The original N-by-N matrix A. */

/*  LDA     (input) INTEGER */
/*          The leading dimension of the array A.  LDA >= MAX(1,N). */

/*  AF      (input) DOUBLE PRECISION array, dimension (LDAF,N) */
/*          The factors L and U from the factorization A = P*L*U */
/*          as computed by DGETRF. */

/*  LDAF    (input) INTEGER */
/*          The leading dimension of the array AF.  LDAF >= MAX(1,N). */

/*  IPIV    (input) INTEGER array, dimension (N) */
/*          The pivot indices from DGETRF; for 1<=i<=N, row i of the */
/*          matrix was interchanged with row IPIV(i). */

/*  B       (input) DOUBLE PRECISION array, dimension (LDB,NRHS) */
/*          The right hand side matrix B. */

/*  LDB     (input) INTEGER */
/*          The leading dimension of the array B.  LDB >= MAX(1,N). */

/*  X       (input/output) DOUBLE PRECISION array, dimension (LDX,NRHS) */
/*          On entry, the solution matrix X, as computed by DGETRS. */
/*          On exit, the improved solution matrix X. */

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

/*<    >*/
/* Subroutine */ int dtgsen_(integer *ijob, logical *wantq, logical *wantz,
        logical *select, integer *n, doublereal *a, integer *lda, doublereal *
        b, integer *ldb, doublereal *alphar, doublereal *alphai, doublereal *
        beta, doublereal *q, integer *ldq, doublereal *z__, integer *ldz,
        integer *m, doublereal *pl, doublereal *pr, doublereal *dif,
        doublereal *work, integer *lwork, integer *iwork, integer *liwork,
        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;
    doublereal d__1;

    /* Builtin functions */
    double sqrt(doublereal), d_sign(doublereal *, doublereal *);

    /* Local variables */
    integer i__, k, n1, n2, kk, ks, mn2, ijb;
    doublereal eps;
    integer kase;
    logical pair;
    integer ierr;
    doublereal dsum;
    logical swap;
    extern /* Subroutine */ int dlag2_(doublereal *, integer *, doublereal *,
            integer *, doublereal *, doublereal *, doublereal *, doublereal *,
             doublereal *, doublereal *);
    logical wantd;
    integer lwmin;
    logical wantp, wantd1, wantd2;
    extern doublereal dlamch_(char *, ftnlen);
    doublereal dscale;
    extern /* Subroutine */ int dlacon_(integer *, doublereal *, doublereal *,
             integer *, doublereal *, integer *);
    doublereal rdscal;
    extern /* Subroutine */ int dlacpy_(char *, integer *, integer *,
            doublereal *, integer *, doublereal *, integer *, ftnlen),
            xerbla_(char *, integer *, ftnlen), dtgexc_(logical *, logical *,
            integer *, doublereal *, integer *, doublereal *, integer *,
            doublereal *, integer *, doublereal *, integer *, integer *,
            integer *, doublereal *, integer *, integer *), dlassq_(integer *,
             doublereal *, integer *, doublereal *, doublereal *);
    integer liwmin;
    extern /* Subroutine */ int dtgsyl_(char *, integer *, integer *, integer
            *, doublereal *, integer *, doublereal *, integer *, doublereal *,
             integer *, doublereal *, integer *, doublereal *, integer *,
            doublereal *, integer *, doublereal *, doublereal *, doublereal *,
             integer *, integer *, integer *, ftnlen);
    doublereal smlnum;
    logical lquery;


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

/*     .. Scalar Arguments .. */
/*<       LOGICAL            WANTQ, WANTZ >*/
/*<    >*/
/*<       DOUBLE PRECISION   PL, PR >*/
/*     .. */
/*     .. Array Arguments .. */
/*<       LOGICAL            SELECT( * ) >*/
/*<       INTEGER            IWORK( * ) >*/
/*<    >*/
/*     .. */

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

/*  DTGSEN reorders the generalized real Schur decomposition of a real */
/*  matrix pair (A, B) (in terms of an orthonormal equivalence trans- */
/*  formation Q' * (A, B) * Z), so that a selected cluster of eigenvalues */
/*  appears in the leading diagonal blocks of the upper quasi-triangular */
/*  matrix A and the upper triangular B. The leading columns of Q and */
/*  Z form orthonormal bases of the corresponding left and right eigen- */
/*  spaces (deflating subspaces). (A, B) must be in generalized real */
/*  Schur canonical form (as returned by DGGES), i.e. A is block upper */
/*  triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper */
/*  triangular. */

/*  DTGSEN also computes the generalized eigenvalues */

/*              w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j) */

/*  of the reordered matrix pair (A, B). */

/*  Optionally, DTGSEN computes the estimates of reciprocal condition */
/*  numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11), */
/*  (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s) */
/*  between the matrix pairs (A11, B11) and (A22,B22) that correspond to */
/*  the selected cluster and the eigenvalues outside the cluster, resp., */
/*  and norms of "projections" onto left and right eigenspaces w.r.t. */
/*  the selected cluster in the (1,1)-block. */

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

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

//.........这里部分代码省略.........
 *                    been completed, but the factor U is exactly   
 *                    singular, so the solution and error bounds   
 *                    could not be computed.   
 *              = A->ncol+1: U is nonsingular, but RCOND is less than machine
 *                    precision, meaning that the matrix is singular to
 *                    working precision. Nevertheless, the solution and
 *                    error bounds are computed because there are a number
 *                    of situations where the computed solution can be more
 *                    accurate than the value of RCOND would suggest.   
 *              > A->ncol+1: number of bytes allocated when memory allocation
 *                    failure occurred, plus A->ncol.
 *
 */

    NCformat  *Astore;
    DNformat  *Bstore, *Xstore;
    double    *Bmat, *Xmat;
    int       ldb, ldx, nrhs;
    SuperMatrix *AA; /* A in NC format used by the factorization routine.*/
    SuperMatrix AC; /* Matrix postmultiplied by Pc */
    int       colequ, equil, dofact, notran, rowequ;
    char      norm[1];
    trans_t   trant;
    int       i, j, info1;
    double    amax, anorm, bignum, smlnum, colcnd, rowcnd, rcmax, rcmin;
    int       n, relax, panel_size;
    Gstat_t   Gstat;
    double    t0;      /* temporary time */
    double    *utime;
    flops_t   *ops, flopcnt;
   
    /* External functions */
    extern double dlangs(char *, SuperMatrix *);
    extern double dlamch_(char *);

    Astore = A->Store;
    Bstore = B->Store;
    Xstore = X->Store;
    Bmat   = Bstore->nzval;
    Xmat   = Xstore->nzval;
    n      = A->ncol;
    ldb    = Bstore->lda;
    ldx    = Xstore->lda;
    nrhs   = B->ncol;
    pdgstrf_options->perm_c = perm_c;
    pdgstrf_options->perm_r = perm_r;

    *info = 0;
    dofact = (pdgstrf_options->fact == DOFACT);
    equil = (pdgstrf_options->fact == EQUILIBRATE);
    notran = (pdgstrf_options->trans == NOTRANS);
    if (dofact || equil) {
	*equed = NOEQUIL;
	rowequ = FALSE;
	colequ = FALSE;
    } else {
	rowequ = (*equed == ROW) || (*equed == BOTH);
	colequ = (*equed == COL) || (*equed == BOTH);
	smlnum = dlamch_("Safe minimum");
	bignum = 1. / smlnum;
    }

    /* ------------------------------------------------------------
       Test the input parameters.
       ------------------------------------------------------------*/
    if ( nprocs <= 0 ) *info = -1;

/* Subroutine */ int zchkhb_(integer *nsizes, integer *nn, integer *nwdths, 
	integer *kk, integer *ntypes, logical *dotype, integer *iseed, 
	doublereal *thresh, integer *nounit, doublecomplex *a, integer *lda, 
	doublereal *sd, doublereal *se, doublecomplex *u, integer *ldu, 
	doublecomplex *work, integer *lwork, doublereal *rwork, doublereal *
	result, integer *info)
{
    /* Initialized data */

    static integer ktype[15] = { 1,2,4,4,4,4,4,5,5,5,5,5,8,8,8 };
    static integer kmagn[15] = { 1,1,1,1,1,2,3,1,1,1,2,3,1,2,3 };
    static integer kmode[15] = { 0,0,4,3,1,4,4,4,3,1,4,4,0,0,0 };

    /* Format strings */
    static char fmt_9999[] = "(\002 ZCHKHB: \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[] = "(/1x,a3,\002 -- Complex Hermitian Banded Tridi"
	    "agonal Reduction Routines\002)";
    static char fmt_9997[] = "(\002 Matrix types (see DCHK23 for details):"
	    " \002)";
    static char fmt_9996[] = "(/\002 Special Matrices:\002,/\002  1=Zero mat"
	    "rix.                        \002,\002  5=Diagonal: clustered ent"
	    "ries.\002,/\002  2=Identity matrix.                    \002,\002"
	    "  6=Diagonal: large, evenly spaced.\002,/\002  3=Diagonal: evenl"
	    "y spaced entries.    \002,\002  7=Diagonal: small, evenly spaced."
	    "\002,/\002  4=Diagonal: geometr. spaced entries.\002)";
    static char fmt_9995[] = "(\002 Dense \002,a,\002 Banded Matrices:\002,"
	    "/\002  8=Evenly spaced eigenvals.            \002,\002 12=Small,"
	    " evenly spaced eigenvals.\002,/\002  9=Geometrically spaced eige"
	    "nvals.     \002,\002 13=Matrix with random O(1) entries.\002,"
	    "/\002 10=Clustered eigenvalues.              \002,\002 14=Matrix"
	    " with large random entries.\002,/\002 11=Large, evenly spaced ei"
	    "genvals.     \002,\002 15=Matrix with small random entries.\002)";
    static char fmt_9994[] = "(/\002 Tests performed:   (S is Tridiag,  U "
	    "is \002,a,\002,\002,/20x,a,\002 means \002,a,\002.\002,/\002 UPL"
	    "O='U':\002,/\002  1= | A - U S U\002,a1,\002 | / ( |A| n ulp )  "
	    "   \002,\002  2= | I - U U\002,a1,\002 | / ( n ulp )\002,/\002 U"
	    "PLO='L':\002,/\002  3= | A - U S U\002,a1,\002 | / ( |A| n ulp )"
	    "     \002,\002  4= | I - U U\002,a1,\002 | / ( n ulp )\002)";
    static char fmt_9993[] = "(\002 N=\002,i5,\002, K=\002,i4,\002, seed="
	    "\002,4(i4,\002,\002),\002 type \002,i2,\002, test(\002,i2,\002)"
	    "=\002,g10.3)";

    /* System generated locals */
    integer a_dim1, a_offset, u_dim1, u_offset, i__1, i__2, i__3, i__4, i__5, 
	    i__6, i__7;
    doublereal d__1;
    doublecomplex z__1;

    /* Local variables */
    integer i__, j, k, n, jc, jr;
    doublereal ulp, cond;
    integer jcol, kmax, nmax;
    doublereal unfl, ovfl, temp1;
    logical badnn;
    integer imode, iinfo;
    extern /* Subroutine */ int zhbt21_(char *, integer *, integer *, integer 
	    *, doublecomplex *, integer *, doublereal *, doublereal *, 
	    doublecomplex *, integer *, doublecomplex *, doublereal *, 
	    doublereal *);
    doublereal aninv, anorm;
    integer nmats, jsize, nerrs, itype, jtype, ntest;
    logical badnnb;
    extern doublereal dlamch_(char *);
    integer idumma[1];
    integer ioldsd[4];
    extern /* Subroutine */ int dlasum_(char *, integer *, integer *, integer 
	    *);
    integer jwidth;
    extern /* Subroutine */ int zhbtrd_(char *, char *, integer *, integer *, 
	    doublecomplex *, integer *, doublereal *, doublereal *, 
	    doublecomplex *, integer *, doublecomplex *, 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, ulpinv;
    extern /* Subroutine */ int zlatms_(integer *, integer *, char *, integer 
	    *, char *, doublereal *, integer *, doublereal *, doublereal *, 
	    integer *, integer *, char *, doublecomplex *, integer *, 
	    doublecomplex *, integer *);
    integer mtypes, ntestt;

    /* Fortran I/O blocks */
    static cilist io___36 = { 0, 0, 0, fmt_9999, 0 };
    static cilist io___37 = { 0, 0, 0, fmt_9999, 0 };
    static cilist io___40 = { 0, 0, 0, fmt_9999, 0 };
    static cilist io___41 = { 0, 0, 0, fmt_9998, 0 };
    static cilist io___42 = { 0, 0, 0, fmt_9997, 0 };
    static cilist io___43 = { 0, 0, 0, fmt_9996, 0 };
    static cilist io___44 = { 0, 0, 0, fmt_9995, 0 };
    static cilist io___45 = { 0, 0, 0, fmt_9994, 0 };
    static cilist io___46 = { 0, 0, 0, fmt_9993, 0 };

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

/* Subroutine */ int dlalsd_(char *uplo, integer *smlsiz, integer *n, integer 
	*nrhs, doublereal *d__, doublereal *e, doublereal *b, integer *ldb, 
	doublereal *rcond, integer *rank, doublereal *work, integer *iwork, 
	integer *info)
{
    /* System generated locals */
    integer b_dim1, b_offset, i__1, i__2;
    doublereal d__1;

    /* Builtin functions */
    double log(doublereal), d_sign(doublereal *, doublereal *);

    /* Local variables */
    static integer difl, difr, perm, nsub;
    extern /* Subroutine */ int drot_(integer *, doublereal *, integer *, 
	    doublereal *, integer *, doublereal *, doublereal *);
    static integer nlvl, sqre, bxst, c__, i__, j, k;
    static doublereal r__;
    static integer s, u;
    extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *, 
	    integer *, doublereal *, doublereal *, integer *, doublereal *, 
	    integer *, doublereal *, doublereal *, integer *);
    static integer z__;
    extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 
	    doublereal *, integer *);
    static integer poles, sizei, nsize, nwork, 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;
    extern /* Subroutine */ int dlalsa_(integer *, integer *, integer *, 
	    integer *, doublereal *, integer *, doublereal *, integer *, 
	    doublereal *, integer *, doublereal *, integer *, doublereal *, 
	    doublereal *, doublereal *, doublereal *, integer *, integer *, 
	    integer *, integer *, doublereal *, doublereal *, doublereal *, 
	    doublereal *, integer *, integer *);
    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 dlacpy_(char *, integer *, integer *, 
	    doublereal *, integer *, doublereal *, integer *), 
	    dlartg_(doublereal *, doublereal *, doublereal *, doublereal *, 
	    doublereal *), dlaset_(char *, integer *, integer *, doublereal *,
	     doublereal *, doublereal *, integer *), xerbla_(char *, 
	    integer *);
    static integer givcol;
    extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *);
    extern /* Subroutine */ int dlasrt_(char *, integer *, doublereal *, 
	    integer *);
    static doublereal orgnrm;
    static integer givnum, givptr, nm1, smlszp, st1;
    static doublereal eps;
    static integer iwk;
    static doublereal tol;


#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]


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

    DLALSD 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.   
    It could conceivably fail on hexadecimal or decimal machines   
    without guard digits, but we know of none.   

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

/* Subroutine */ int zchkbl_(integer *nin, integer *nout)
{
    /* Format strings */
    static char fmt_9999[] = "(1x,\002.. test output of ZGEBAL .. \002)";
    static char fmt_9998[] = "(1x,\002value of largest test error           "
	    " = \002,d12.3)";
    static char fmt_9997[] = "(1x,\002example number where info is not zero "
	    " = \002,i4)";
    static char fmt_9996[] = "(1x,\002example number where ILO or IHI wrong "
	    " = \002,i4)";
    static char fmt_9995[] = "(1x,\002example number having largest error   "
	    " = \002,i4)";
    static char fmt_9994[] = "(1x,\002number of examples where info is not 0"
	    " = \002,i4)";
    static char fmt_9993[] = "(1x,\002total number of examples tested       &qu