/* Subroutine */
int zlaev2_(doublecomplex *a, doublecomplex *b, doublecomplex *c__, doublereal *rt1, doublereal *rt2, doublereal *cs1, doublecomplex *sn1)
{
    /* System generated locals */
    doublereal d__1, d__2, d__3;
    doublecomplex z__1, z__2;
    /* Builtin functions */
    double z_abs(doublecomplex *);
    void d_cnjg(doublecomplex *, doublecomplex *);
    /* Local variables */
    doublereal t;
    doublecomplex w;
    extern /* Subroutine */
    int dlaev2_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *);
    /* -- LAPACK auxiliary routine (version 3.4.2) -- */
    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
    /* September 2012 */
    /* .. Scalar Arguments .. */
    /* .. */
    /* ===================================================================== */
    /* .. Parameters .. */
    /* .. */
    /* .. Local Scalars .. */
    /* .. */
    /* .. External Subroutines .. */
    /* .. */
    /* .. Intrinsic Functions .. */
    /* .. */
    /* .. Executable Statements .. */
    if (z_abs(b) == 0.)
    {
        w.r = 1.;
        w.i = 0.; // , expr subst
    }
    else
    {
        d_cnjg(&z__2, b);
        d__1 = z_abs(b);
        z__1.r = z__2.r / d__1;
        z__1.i = z__2.i / d__1; // , expr subst
        w.r = z__1.r;
        w.i = z__1.i; // , expr subst
    }
    d__1 = a->r;
    d__2 = z_abs(b);
    d__3 = c__->r;
    dlaev2_(&d__1, &d__2, &d__3, rt1, rt2, cs1, &t);
    z__1.r = t * w.r;
    z__1.i = t * w.i; // , expr subst
    sn1->r = z__1.r, sn1->i = z__1.i;
    return 0;
    /* End of ZLAEV2 */
}

/* Subroutine */ int zrotg_(doublecomplex *ca, doublecomplex *cb, doublereal *
	c, doublecomplex *s)
{
    /* System generated locals */
    doublereal d__1, d__2;
    doublecomplex z__1, z__2, z__3, z__4;

    /* Builtin functions */
    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
    double sqrt(doublereal);
    void d_cnjg(doublecomplex *, doublecomplex *);

    /* Local variables */
    static doublereal norm;
    extern doublereal cdabs_(doublecomplex *);
    static doublecomplex alpha;
    static doublereal scale;

    if (cdabs_(ca) != 0.) {
	goto L10;
    }
    *c = 0.;
    s->r = 1., s->i = 0.;
    ca->r = cb->r, ca->i = cb->i;
    goto L20;
L10:
    scale = cdabs_(ca) + cdabs_(cb);
    z__2.r = scale, z__2.i = 0.;
    z_div(&z__1, ca, &z__2);
/* Computing 2nd power */
    d__1 = cdabs_(&z__1);
    z__4.r = scale, z__4.i = 0.;
    z_div(&z__3, cb, &z__4);
/* Computing 2nd power */
    d__2 = cdabs_(&z__3);
    norm = scale * sqrt(d__1 * d__1 + d__2 * d__2);
    d__1 = cdabs_(ca);
    z__1.r = ca->r / d__1, z__1.i = ca->i / d__1;
    alpha.r = z__1.r, alpha.i = z__1.i;
    *c = cdabs_(ca) / norm;
    d_cnjg(&z__3, cb);
    z__2.r = alpha.r * z__3.r - alpha.i * z__3.i, z__2.i = alpha.r * z__3.i + 
	    alpha.i * z__3.r;
    z__1.r = z__2.r / norm, z__1.i = z__2.i / norm;
    s->r = z__1.r, s->i = z__1.i;
    z__1.r = norm * alpha.r, z__1.i = norm * alpha.i;
    ca->r = z__1.r, ca->i = z__1.i;
L20:
    return 0;
} /* zrotg_ */

/* 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 */
//.........这里部分代码省略.........

/* Double Complex */ VOID zlatm2_(doublecomplex * ret_val, integer *m, 
	integer *n, integer *i, integer *j, integer *kl, integer *ku, integer 
	*idist, integer *iseed, doublecomplex *d, integer *igrade, 
	doublecomplex *dl, doublecomplex *dr, integer *ipvtng, integer *iwork,
	 doublereal *sparse)
{
    /* System generated locals */
    integer i__1, i__2;
    doublecomplex z__1, z__2, z__3;

    /* Builtin functions */
    void z_div(doublecomplex *, doublecomplex *, doublecomplex *), d_cnjg(
	    doublecomplex *, doublecomplex *);

    /* Local variables */
    static integer isub, jsub;
    static doublecomplex ctemp;
    extern doublereal dlaran_(integer *);
    extern /* Double Complex */ VOID zlarnd_(doublecomplex *, integer *, 
	    integer *);


/*  -- LAPACK auxiliary test routine (version 2.0) --   
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
       Courant Institute, Argonne National Lab, and Rice University   
       February 29, 1992   





    Purpose   
    =======   

       ZLATM2 returns the (I,J) entry of a random matrix of dimension   
       (M, N) described by the other paramters. It is called by the   
       ZLATMR routine in order to build random test matrices. No error   
       checking on parameters is done, because this routine is called in 
  
       a tight loop by ZLATMR which has already checked the parameters.   

       Use of ZLATM2 differs from CLATM3 in the order in which the random 
  
       number generator is called to fill in random matrix entries.   
       With ZLATM2, the generator is called to fill in the pivoted matrix 
  
       columnwise. With ZLATM3, the generator is called to fill in the   
       matrix columnwise, after which it is pivoted. Thus, ZLATM3 can   
       be used to construct random matrices which differ only in their   
       order of rows and/or columns. ZLATM2 is used to construct band   
       matrices while avoiding calling the random number generator for   
       entries outside the band (and therefore generating random numbers 
  

       The matrix whose (I,J) entry is returned is constructed as   
       follows (this routine only computes one entry):   

         If I is outside (1..M) or J is outside (1..N), return zero   
            (this is convenient for generating matrices in band format). 
  

         Generate a matrix A with random entries of distribution IDIST.   

         Set the diagonal to D.   

         Grade the matrix, if desired, from the left (by DL) and/or   
            from the right (by DR or DL) as specified by IGRADE.   

         Permute, if desired, the rows and/or columns as specified by   
            IPVTNG and IWORK.   

         Band the matrix to have lower bandwidth KL and upper   
            bandwidth KU.   

         Set random entries to zero as specified by SPARSE.   

    Arguments   
    =========   

    M      - INTEGER   
             Number of rows of matrix. Not modified.   

    N      - INTEGER   
             Number of columns of matrix. Not modified.   

    I      - INTEGER   
             Row of entry to be returned. Not modified.   

    J      - INTEGER   
             Column of entry to be returned. Not modified.   

    KL     - INTEGER   
             Lower bandwidth. Not modified.   

    KU     - INTEGER   
             Upper bandwidth. Not modified.   

    IDIST  - INTEGER   
             On entry, IDIST specifies the type of distribution to be   
             used to generate a random matrix .   
//.........这里部分代码省略.........

/* Subroutine */ int zlatm6_(integer *type__, integer *n, doublecomplex *a, 
	integer *lda, doublecomplex *b, doublecomplex *x, integer *ldx, 
	doublecomplex *y, integer *ldy, doublecomplex *alpha, doublecomplex *
	beta, doublecomplex *wx, doublecomplex *wy, doublereal *s, doublereal 
	*dif)
{
    /* System generated locals */
    integer a_dim1, a_offset, b_dim1, b_offset, x_dim1, x_offset, y_dim1, 
	    y_offset, i__1, i__2, i__3;
    doublereal d__1, d__2;
    doublecomplex z__1, z__2, z__3, z__4;

    /* Builtin functions */
    void d_cnjg(doublecomplex *, doublecomplex *);
    double z_abs(doublecomplex *), sqrt(doublereal);

    /* Local variables */
    integer i__, j;
    doublecomplex z__[64]	/* was [8][8] */;
    integer info;
    doublecomplex work[26];
    doublereal rwork[50];
    extern /* Subroutine */ int zlakf2_(integer *, integer *, doublecomplex *, 
	     integer *, doublecomplex *, doublecomplex *, doublecomplex *, 
	    doublecomplex *, integer *), zgesvd_(char *, char *, integer *, 
	    integer *, doublecomplex *, integer *, doublereal *, 
	    doublecomplex *, integer *, doublecomplex *, integer *, 
	    doublecomplex *, integer *, doublereal *, integer *), zlacpy_(char *, integer *, integer *, doublecomplex *, 
	    integer *, doublecomplex *, integer *);


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

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

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

/*  ZLATM6 generates test matrices for the generalized eigenvalue */
/*  problem, their corresponding right and left eigenvector matrices, */
/*  and also reciprocal condition numbers for all eigenvalues and */
/*  the reciprocal condition numbers of eigenvectors corresponding to */
/*  the 1th and 5th eigenvalues. */

/*  Test Matrices */
/*  ============= */

/*  Two kinds of test matrix pairs */
/*           (A, B) = inverse(YH) * (Da, Db) * inverse(X) */
/*  are used in the tests: */

/*  Type 1: */
/*     Da = 1+a   0    0    0    0    Db = 1   0   0   0   0 */
/*           0   2+a   0    0    0         0   1   0   0   0 */
/*           0    0   3+a   0    0         0   0   1   0   0 */
/*           0    0    0   4+a   0         0   0   0   1   0 */
/*           0    0    0    0   5+a ,      0   0   0   0   1 */
/*  and Type 2: */
/*     Da = 1+i   0    0       0       0    Db = 1   0   0   0   0 */
/*           0   1-i   0       0       0         0   1   0   0   0 */
/*           0    0    1       0       0         0   0   1   0   0 */
/*           0    0    0 (1+a)+(1+b)i  0         0   0   0   1   0 */
/*           0    0    0       0 (1+a)-(1+b)i,   0   0   0   0   1 . */

/*  In both cases the same inverse(YH) and inverse(X) are used to compute */
/*  (A, B), giving the exact eigenvectors to (A,B) as (YH, X): */

/*  YH:  =  1    0   -y    y   -y    X =  1   0  -x  -x   x */
/*          0    1   -y    y   -y         0   1   x  -x  -x */
/*          0    0    1    0    0         0   0   1   0   0 */
/*          0    0    0    1    0         0   0   0   1   0 */
/*          0    0    0    0    1,        0   0   0   0   1 , where */

/*  a, b, x and y will have all values independently of each other. */

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

/*  TYPE    (input) INTEGER */
/*          Specifies the problem type (see futher details). */

/*  N       (input) INTEGER */
/*          Size of the matrices A and B. */

/*  A       (output) COMPLEX*16 array, dimension (LDA, N). */
/*          On exit A N-by-N is initialized according to TYPE. */

/*  LDA     (input) INTEGER */
/*          The leading dimension of A and of B. */

/*  B       (output) COMPLEX*16 array, dimension (LDA, N). */
/*          On exit B N-by-N is initialized according to TYPE. */

/*  X       (output) COMPLEX*16 array, dimension (LDX, N). */
/*          On exit X is the N-by-N matrix of right eigenvectors. */
//.........这里部分代码省略.........

/* Subroutine */
int zhptri_(char *uplo, integer *n, doublecomplex *ap, integer *ipiv, doublecomplex *work, integer *info)
{
    /* System generated locals */
    integer i__1, i__2, i__3;
    doublereal d__1;
    doublecomplex z__1, z__2;
    /* Builtin functions */
    double z_abs(doublecomplex *);
    void d_cnjg(doublecomplex *, doublecomplex *);
    /* Local variables */
    doublereal d__;
    integer j, k;
    doublereal t, ak;
    integer kc, kp, kx, kpc, npp;
    doublereal akp1;
    doublecomplex temp, akkp1;
    extern logical lsame_(char *, char *);
    extern /* Double Complex */
    VOID zdotc_f2c_(doublecomplex *, integer *, doublecomplex *, integer *, doublecomplex *, integer *);
    integer kstep;
    logical upper;
    extern /* Subroutine */
    int zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, integer *), zhpmv_(char *, integer *, doublecomplex *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, doublecomplex *, integer *), zswap_( integer *, doublecomplex *, integer *, doublecomplex *, integer *) , xerbla_(char *, integer *);
    integer kcnext;
    /* -- LAPACK computational 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 Functions .. */
    /* .. */
    /* .. External Subroutines .. */
    /* .. */
    /* .. Intrinsic Functions .. */
    /* .. */
    /* .. Executable Statements .. */
    /* Test the input parameters. */
    /* Parameter adjustments */
    --work;
    --ipiv;
    --ap;
    /* Function Body */
    *info = 0;
    upper = lsame_(uplo, "U");
    if (! upper && ! lsame_(uplo, "L"))
    {
        *info = -1;
    }
    else if (*n < 0)
    {
        *info = -2;
    }
    if (*info != 0)
    {
        i__1 = -(*info);
        xerbla_("ZHPTRI", &i__1);
        return 0;
    }
    /* Quick return if possible */
    if (*n == 0)
    {
        return 0;
    }
    /* Check that the diagonal matrix D is nonsingular. */
    if (upper)
    {
        /* Upper triangular storage: examine D from bottom to top */
        kp = *n * (*n + 1) / 2;
        for (*info = *n;
                *info >= 1;
                --(*info))
        {
            i__1 = kp;
            if (ipiv[*info] > 0 && (ap[i__1].r == 0. && ap[i__1].i == 0.))
            {
                return 0;
            }
            kp -= *info;
            /* L10: */
        }
    }
    else
    {
        /* Lower triangular storage: examine D from top to bottom. */
        kp = 1;
        i__1 = *n;
        for (*info = 1;
                *info <= i__1;
                ++(*info))
        {
            i__2 = kp;
//.........这里部分代码省略.........

/* Subroutine */ int zlaptm_(char *uplo, integer *n, integer *nrhs, 
	doublereal *alpha, doublereal *d__, doublecomplex *e, doublecomplex *
	x, integer *ldx, doublereal *beta, doublecomplex *b, integer *ldb)
{
    /* System generated locals */
    integer b_dim1, b_offset, x_dim1, x_offset, i__1, i__2, i__3, i__4, i__5, 
	    i__6, i__7, i__8, i__9;
    doublecomplex z__1, z__2, z__3, z__4, z__5, z__6, z__7;

    /* Builtin functions */
    void d_cnjg(doublecomplex *, doublecomplex *);

    /* Local variables */
    static integer i__, j;
    extern logical lsame_(char *, char *);


#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 x_subscr(a_1,a_2) (a_2)*x_dim1 + a_1
#define x_ref(a_1,a_2) x[x_subscr(a_1,a_2)]


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


    Purpose   
    =======   

    ZLAPTM multiplies an N by NRHS matrix X by a Hermitian tridiagonal   
    matrix A and stores the result in a matrix B.  The operation has the   
    form   

       B := alpha * A * X + beta * B   

    where alpha may be either 1. or -1. and beta may be 0., 1., or -1.   

    Arguments   
    =========   

    UPLO    (input) CHARACTER   
            Specifies whether the superdiagonal or the subdiagonal of the   
            tridiagonal matrix A is stored.   
            = 'U':  Upper, E is the superdiagonal of A.   
            = 'L':  Lower, E is the subdiagonal of A.   

    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 X and B.   

    ALPHA   (input) DOUBLE PRECISION   
            The scalar alpha.  ALPHA must be 1. or -1.; otherwise,   
            it is assumed to be 0.   

    D       (input) DOUBLE PRECISION array, dimension (N)   
            The n diagonal elements of the tridiagonal matrix A.   

    E       (input) COMPLEX*16 array, dimension (N-1)   
            The (n-1) subdiagonal or superdiagonal elements of A.   

    X       (input) COMPLEX*16 array, dimension (LDX,NRHS)   
            The N by NRHS matrix X.   

    LDX     (input) INTEGER   
            The leading dimension of the array X.  LDX >= max(N,1).   

    BETA    (input) DOUBLE PRECISION   
            The scalar beta.  BETA must be 0., 1., or -1.; otherwise,   
            it is assumed to be 1.   

    B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)   
            On entry, the N by NRHS matrix B.   
            On exit, B is overwritten by the matrix expression   
            B := alpha * A * X + beta * B.   

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

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


       Parameter adjustments */
    --d__;
    --e;
    x_dim1 = *ldx;
    x_offset = 1 + x_dim1 * 1;
    x -= x_offset;
    b_dim1 = *ldb;
    b_offset = 1 + b_dim1 * 1;
    b -= b_offset;

    /* Function Body */
    if (*n == 0) {
	return 0;
//.........这里部分代码省略.........

/* Subroutine */ int zgehd2_(integer *n, integer *ilo, integer *ihi,
                             doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *
                             work, integer *info)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2, i__3;
    doublecomplex z__1;

    /* Builtin functions */
    void d_cnjg(doublecomplex *, doublecomplex *);

    /* Local variables */
    integer i__;
    doublecomplex alpha;
    extern /* Subroutine */ int zlarf_(char *, integer *, integer *,
                                       doublecomplex *, integer *, doublecomplex *, doublecomplex *,
                                       integer *, doublecomplex *), xerbla_(char *, integer *), zlarfg_(integer *, doublecomplex *, doublecomplex *,
                                               integer *, doublecomplex *);


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

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

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

    /*  ZGEHD2 reduces a complex general matrix A to upper Hessenberg form H */
    /*  by a unitary similarity transformation:  Q' * A * Q = H . */

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

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

    /*  ILO     (input) INTEGER */
    /*  IHI     (input) INTEGER */
    /*          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 ZGEBAL; otherwise they should be */
    /*          set to 1 and N respectively. See Further Details. */
    /*          1 <= ILO <= IHI <= max(1,N). */

    /*  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 */
    /*          elements below the first subdiagonal, with the array TAU, */
    /*          represent the unitary matrix Q as a product of elementary */
    /*          reflectors. See Further Details. */

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

    /*  TAU     (output) COMPLEX*16 array, dimension (N-1) */
    /*          The scalar factors of the elementary reflectors (see Further */
    /*          Details). */

    /*  WORK    (workspace) COMPLEX*16 array, dimension (N) */

    /*  INFO    (output) INTEGER */
    /*          = 0:  successful exit */
    /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */

    /*  Further Details */
    /*  =============== */

    /*  The matrix Q is represented as a product of (ihi-ilo) elementary */
    /*  reflectors */

    /*     Q = H(ilo) H(ilo+1) . . . H(ihi-1). */

    /*  Each H(i) has the form */

    /*     H(i) = I - tau * v * v' */

    /*  where tau is a complex scalar, and v is a complex vector with */
    /*  v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on */
    /*  exit in A(i+2:ihi,i), and tau in TAU(i). */

    /*  The contents of A are illustrated by the following example, with */
    /*  n = 7, ilo = 2 and ihi = 6: */

    /*  on entry,                        on exit, */

    /*  ( a   a   a   a   a   a   a )    (  a   a   h   h   h   h   a ) */
    /*  (     a   a   a   a   a   a )    (      a   h   h   h   h   a ) */
    /*  (     a   a   a   a   a   a )    (      h   h   h   h   h   h ) */
    /*  (     a   a   a   a   a   a )    (      v2  h   h   h   h   h ) */
    /*  (     a   a   a   a   a   a )    (      v2  v3  h   h   h   h ) */
    /*  (     a   a   a   a   a   a )    (      v2  v3  v4  h   h   h ) */
    /*  (                         a )    (                          a ) */

    /*  where a denotes an element of the original matrix A, h denotes a */
//.........这里部分代码省略.........

/* Subroutine */ int zhpr_(char *uplo, integer *n, doublereal *alpha, 
	doublecomplex *x, integer *incx, doublecomplex *ap, ftnlen uplo_len)
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5;
    doublereal d__1;
    doublecomplex z__1, z__2;

    /* Builtin functions */
    void d_cnjg(doublecomplex *, doublecomplex *);

    /* Local variables */
    static integer i__, j, k, kk, ix, jx, kx, info;
    static doublecomplex temp;
    extern logical lsame_(char *, char *, ftnlen, ftnlen);
    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);

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

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

/*  ZHPR    performs the hermitian rank 1 operation */

/*     A := alpha*x*conjg( x' ) + A, */

/*  where alpha is a real scalar, x is an n element vector and A is an */
/*  n by n hermitian matrix, supplied in packed form. */

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

/*  UPLO   - CHARACTER*1. */
/*           On entry, UPLO specifies whether the upper or lower */
/*           triangular part of the matrix A is supplied in the packed */
/*           array AP as follows: */

/*              UPLO = 'U' or 'u'   The upper triangular part of A is */
/*                                  supplied in AP. */

/*              UPLO = 'L' or 'l'   The lower triangular part of A is */
/*                                  supplied in AP. */

/*           Unchanged on exit. */

/*  N      - INTEGER. */
/*           On entry, N specifies the order of the matrix A. */
/*           N must be at least zero. */
/*           Unchanged on exit. */

/*  ALPHA  - DOUBLE PRECISION. */
/*           On entry, ALPHA specifies the scalar alpha. */
/*           Unchanged on exit. */

/*  X      - COMPLEX*16       array of dimension at least */
/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
/*           Before entry, the incremented array X must contain the n */
/*           element vector x. */
/*           Unchanged on exit. */

/*  INCX   - INTEGER. */
/*           On entry, INCX specifies the increment for the elements of */
/*           X. INCX must not be zero. */
/*           Unchanged on exit. */

/*  AP     - COMPLEX*16       array of DIMENSION at least */
/*           ( ( n*( n + 1 ) )/2 ). */
/*           Before entry with  UPLO = 'U' or 'u', the array AP must */
/*           contain the upper triangular part of the hermitian matrix */
/*           packed sequentially, column by column, so that AP( 1 ) */
/*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
/*           and a( 2, 2 ) respectively, and so on. On exit, the array */
/*           AP is overwritten by the upper triangular part of the */
/*           updated matrix. */
/*           Before entry with UPLO = 'L' or 'l', the array AP must */
/*           contain the lower triangular part of the hermitian matrix */
/*           packed sequentially, column by column, so that AP( 1 ) */
/*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
/*           and a( 3, 1 ) respectively, and so on. On exit, the array */
/*           AP is overwritten by the lower triangular part of the */
/*           updated matrix. */
/*           Note that the imaginary parts of the diagonal elements need */
/*           not be set, they are assumed to be zero, and on exit they */
/*           are set to zero. */


/*  Level 2 Blas routine. */

/*  -- Written on 22-October-1986. */
/*     Jack Dongarra, Argonne National Lab. */
/*     Jeremy Du Croz, Nag Central Office. */
/*     Sven Hammarling, Nag Central Office. */
/*     Richard Hanson, Sandia National Labs. */


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

/* Subroutine */
int zhpmv_(char *uplo, integer *n, doublecomplex *alpha, doublecomplex *ap, doublecomplex *x, integer *incx, doublecomplex * beta, doublecomplex *y, integer *incy)
{
    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5;
    doublereal d__1;
    doublecomplex z__1, z__2, z__3, z__4;
    /* Builtin functions */
    void d_cnjg(doublecomplex *, doublecomplex *);
    /* Local variables */
    integer info;
    doublecomplex temp1, temp2;
    integer i__, j, k;
    extern logical lsame_(char *, char *);
    integer kk, ix, iy, jx, jy, kx, ky;
    extern /* Subroutine */
    int xerbla_(char *, integer *);
    /* .. Scalar Arguments .. */
    /* .. Array Arguments .. */
    /* .. */
    /* Purpose */
    /* ======= */
    /* ZHPMV performs the matrix-vector operation */
    /* y := alpha*A*x + beta*y, */
    /* where alpha and beta are scalars, x and y are n element vectors and */
    /* A is an n by n hermitian matrix, supplied in packed form. */
    /* Parameters */
    /* ========== */
    /* UPLO - CHARACTER*1. */
    /* On entry, UPLO specifies whether the upper or lower */
    /* triangular part of the matrix A is supplied in the packed */
    /* array AP as follows: */
    /* UPLO = 'U' or 'u' The upper triangular part of A is */
    /* supplied in AP. */
    /* UPLO = 'L' or 'l' The lower triangular part of A is */
    /* supplied in AP. */
    /* Unchanged on exit. */
    /* N - INTEGER. */
    /* On entry, N specifies the order of the matrix A. */
    /* N must be at least zero. */
    /* Unchanged on exit. */
    /* ALPHA - COMPLEX*16 . */
    /* On entry, ALPHA specifies the scalar alpha. */
    /* Unchanged on exit. */
    /* AP - COMPLEX*16 array of DIMENSION at least */
    /* ( ( n*( n + 1 ) )/2 ). */
    /* Before entry with UPLO = 'U' or 'u', the array AP must */
    /* contain the upper triangular part of the hermitian matrix */
    /* packed sequentially, column by column, so that AP( 1 ) */
    /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
    /* and a( 2, 2 ) respectively, and so on. */
    /* Before entry with UPLO = 'L' or 'l', the array AP must */
    /* contain the lower triangular part of the hermitian matrix */
    /* packed sequentially, column by column, so that AP( 1 ) */
    /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
    /* and a( 3, 1 ) respectively, and so on. */
    /* Note that the imaginary parts of the diagonal elements need */
    /* not be set and are assumed to be zero. */
    /* Unchanged on exit. */
    /* X - COMPLEX*16 array of dimension at least */
    /* ( 1 + ( n - 1 )*f2c_abs( INCX ) ). */
    /* Before entry, the incremented array X must contain the n */
    /* element vector x. */
    /* Unchanged on exit. */
    /* INCX - INTEGER. */
    /* On entry, INCX specifies the increment for the elements of */
    /* X. INCX must not be zero. */
    /* Unchanged on exit. */
    /* BETA - COMPLEX*16 . */
    /* On entry, BETA specifies the scalar beta. When BETA is */
    /* supplied as zero then Y need not be set on input. */
    /* Unchanged on exit. */
    /* Y - COMPLEX*16 array of dimension at least */
    /* ( 1 + ( n - 1 )*f2c_abs( INCY ) ). */
    /* Before entry, the incremented array Y must contain the n */
    /* element vector y. On exit, Y is overwritten by the updated */
    /* vector y. */
    /* INCY - INTEGER. */
    /* On entry, INCY specifies the increment for the elements of */
    /* Y. INCY must not be zero. */
    /* Unchanged on exit. */
    /* Level 2 Blas routine. */
    /* -- Written on 22-October-1986. */
    /* Jack Dongarra, Argonne National Lab. */
    /* Jeremy Du Croz, Nag Central Office. */
    /* Sven Hammarling, Nag Central Office. */
    /* Richard Hanson, Sandia National Labs. */
    /* .. Parameters .. */
    /* .. Local Scalars .. */
    /* .. External Functions .. */
    /* .. External Subroutines .. */
    /* .. Intrinsic Functions .. */
    /* .. */
    /* .. Executable Statements .. */
    /* Test the input parameters. */
    /* Parameter adjustments */
    --y;
    --x;
    --ap;
    /* Function Body */
//.........这里部分代码省略.........

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

/*                    generate plane rotation to annihilate a(i+ml-1,i) */
/*                    within the band, and apply rotation from the left */

			zlartg_(&ab[*ku + ml - 1 + i__ * ab_dim1], &ab[*ku + 
				ml + i__ * ab_dim1], &rwork[i__ + ml - 1], &
				work[i__ + ml - 1], &ra);
			i__3 = *ku + ml - 1 + i__ * ab_dim1;
			ab[i__3].r = ra.r, ab[i__3].i = ra.i;
			if (i__ < *n) {
/* Computing MIN */
			    i__4 = *ku + ml - 2, i__5 = *n - i__;
			    i__3 = min(i__4,i__5);
			    i__6 = *ldab - 1;
			    i__7 = *ldab - 1;
			    zrot_(&i__3, &ab[*ku + ml - 2 + (i__ + 1) * 
				    ab_dim1], &i__6, &ab[*ku + ml - 1 + (i__ 
				    + 1) * ab_dim1], &i__7, &rwork[i__ + ml - 
				    1], &work[i__ + ml - 1]);
			}
		    }
		    ++nr;
		    j1 -= kb1;
		}

		if (wantq) {

/*                 accumulate product of plane rotations in Q */

		    i__3 = j2;
		    i__4 = kb1;
		    for (j = j1; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) 
			    {
			d_cnjg(&z__1, &work[j]);
			zrot_(m, &q[(j - 1) * q_dim1 + 1], &c__1, &q[j * 
				q_dim1 + 1], &c__1, &rwork[j], &z__1);
		    }
		}

		if (wantc) {

/*                 apply plane rotations to C */

		    i__4 = j2;
		    i__3 = kb1;
		    for (j = j1; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) 
			    {
			zrot_(ncc, &c__[j - 1 + c_dim1], ldc, &c__[j + c_dim1]
, ldc, &rwork[j], &work[j]);
		    }
		}

		if (j2 + kun > *n) {

/*                 adjust J2 to keep within the bounds of the matrix */

		    --nr;
		    j2 -= kb1;
		}

		i__3 = j2;
		i__4 = kb1;
		for (j = j1; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) {

/*                 create nonzero element a(j-1,j+ku) above the band */
/*                 and store it in WORK(n+1:2*n) */

/* Subroutine */ int ztrsv_(char *uplo, char *trans, char *diag, integer *n, 
	doublecomplex *a, integer *lda, doublecomplex *x, integer *incx)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
    doublecomplex z__1, z__2, z__3;

    /* Builtin functions */
    void z_div(doublecomplex *, doublecomplex *, doublecomplex *), d_cnjg(
	    doublecomplex *, doublecomplex *);

    /* Local variables */
    integer i__, j, ix, jx, kx, info;
    doublecomplex temp;
    extern logical lsame_(char *, char *);
    extern /* Subroutine */ int xerbla_(char *, integer *);
    logical noconj, nounit;

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

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

/*  ZTRSV  solves one of the systems of equations */

/*     A*x = b,   or   A'*x = b,   or   conjg( A' )*x = b, */

/*  where b and x are n element vectors and A is an n by n unit, or */
/*  non-unit, upper or lower triangular matrix. */

/*  No test for singularity or near-singularity is included in this */
/*  routine. Such tests must be performed before calling this routine. */

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

/*  UPLO   - CHARACTER*1. */
/*           On entry, UPLO specifies whether the matrix is an upper or */
/*           lower triangular matrix as follows: */

/*              UPLO = 'U' or 'u'   A is an upper triangular matrix. */

/*              UPLO = 'L' or 'l'   A is a lower triangular matrix. */

/*           Unchanged on exit. */

/*  TRANS  - CHARACTER*1. */
/*           On entry, TRANS specifies the equations to be solved as */
/*           follows: */

/*              TRANS = 'N' or 'n'   A*x = b. */

/*              TRANS = 'T' or 't'   A'*x = b. */

/*              TRANS = 'C' or 'c'   conjg( A' )*x = b. */

/*           Unchanged on exit. */

/*  DIAG   - CHARACTER*1. */
/*           On entry, DIAG specifies whether or not A is unit */
/*           triangular as follows: */

/*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. */

/*              DIAG = 'N' or 'n'   A is not assumed to be unit */
/*                                  triangular. */

/*           Unchanged on exit. */

/*  N      - INTEGER. */
/*           On entry, N specifies the order of the matrix A. */
/*           N must be at least zero. */
/*           Unchanged on exit. */

/*  A      - COMPLEX*16       array of DIMENSION ( LDA, n ). */
/*           Before entry with  UPLO = 'U' or 'u', the leading n by n */
/*           upper triangular part of the array A must contain the upper */
/*           triangular matrix and the strictly lower triangular part of */
/*           A is not referenced. */
/*           Before entry with UPLO = 'L' or 'l', the leading n by n */
/*           lower triangular part of the array A must contain the lower */
/*           triangular matrix and the strictly upper triangular part of */
/*           A is not referenced. */
/*           Note that when  DIAG = 'U' or 'u', the diagonal elements of */
/*           A are not referenced either, but are assumed to be unity. */
/*           Unchanged on exit. */

/*  LDA    - INTEGER. */
/*           On entry, LDA specifies the first dimension of A as declared */
/*           in the calling (sub) program. LDA must be at least */
/*           max( 1, n ). */
/*           Unchanged on exit. */

/*  X      - COMPLEX*16       array of dimension at least */
/*           ( 1 + ( n - 1 )*abs( INCX ) ). */
/*           Before entry, the incremented array X must contain the n */
/*           element right-hand side vector b. On exit, X is overwritten */
//.........这里部分代码省略.........

/* Subroutine */ int zptrfs_(char *uplo, integer *n, integer *nrhs, 
	doublereal *d__, doublecomplex *e, doublereal *df, doublecomplex *ef, 
	doublecomplex *b, integer *ldb, doublecomplex *x, integer *ldx, 
	doublereal *ferr, doublereal *berr, doublecomplex *work, doublereal *
	rwork, integer *info)
{
    /* System generated locals */
    integer b_dim1, b_offset, x_dim1, x_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, d__7, d__8, d__9, d__10, 
	    d__11, d__12;
    doublecomplex z__1, z__2, z__3;

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

    /* Local variables */
    integer i__, j;
    doublereal s;
    doublecomplex bi, cx, dx, ex;
    integer ix, nz;
    doublereal eps, safe1, safe2;
    extern logical lsame_(char *, char *);
    integer count;
    logical upper;
    extern /* Subroutine */ int zaxpy_(integer *, doublecomplex *, 
	    doublecomplex *, integer *, doublecomplex *, integer *);
    extern doublereal dlamch_(char *);
    extern integer idamax_(integer *, doublereal *, integer *);
    doublereal safmin;
    extern /* Subroutine */ int xerbla_(char *, integer *);
    doublereal lstres;
    extern /* Subroutine */ int zpttrs_(char *, integer *, integer *, 
	    doublereal *, doublecomplex *, doublecomplex *, integer *, 
	    integer *);


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

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

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

/*  ZPTRFS improves the computed solution to a system of linear */
/*  equations when the coefficient matrix is Hermitian positive definite */
/*  and tridiagonal, and provides error bounds and backward error */
/*  estimates for the solution. */

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

/*  UPLO    (input) CHARACTER*1 */
/*          Specifies whether the superdiagonal or the subdiagonal of the */
/*          tridiagonal matrix A is stored and the form of the */
/*          factorization: */
/*          = 'U':  E is the superdiagonal of A, and A = U**H*D*U; */
/*          = 'L':  E is the subdiagonal of A, and A = L*D*L**H. */
/*          (The two forms are equivalent if A is real.) */

/*  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 matrix B.  NRHS >= 0. */

/*  D       (input) DOUBLE PRECISION array, dimension (N) */
/*          The n real diagonal elements of the tridiagonal matrix A. */

/*  E       (input) COMPLEX*16 array, dimension (N-1) */
/*          The (n-1) off-diagonal elements of the tridiagonal matrix A */
/*          (see UPLO). */

/*  DF      (input) DOUBLE PRECISION array, dimension (N) */
/*          The n diagonal elements of the diagonal matrix D from */
/*          the factorization computed by ZPTTRF. */

/*  EF      (input) COMPLEX*16 array, dimension (N-1) */
/*          The (n-1) off-diagonal elements of the unit bidiagonal */
/*          factor U or L from the factorization computed by ZPTTRF */
/*          (see UPLO). */

/*  B       (input) COMPLEX*16 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) COMPLEX*16 array, dimension (LDX,NRHS) */
/*          On entry, the solution matrix X, as computed by ZPTTRS. */
/*          On exit, the improved solution matrix X. */

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

/* Subroutine */ int zgemv_(char *trans, integer *m, integer *n,
                            doublecomplex *alpha, doublecomplex *a, integer *lda, doublecomplex *
                            x, integer *incx, doublecomplex *beta, doublecomplex *y, integer *
                            incy, ftnlen trans_len)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
    doublecomplex z__1, z__2, z__3;

    /* Builtin functions */
    void d_cnjg(doublecomplex *, doublecomplex *);

    /* Local variables */
    static integer i__, j, ix, iy, jx, jy, kx, ky, info;
    static doublecomplex temp;
    static integer lenx, leny;
    extern logical lsame_(char *, char *, ftnlen, ftnlen);
    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
    static logical noconj;

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

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

    /*  ZGEMV  performs one of the matrix-vector operations */

    /*     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,   or */

    /*     y := alpha*conjg( A' )*x + beta*y, */

    /*  where alpha and beta are scalars, x and y are vectors and A is an */
    /*  m by n matrix. */

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

    /*  TRANS  - CHARACTER*1. */
    /*           On entry, TRANS specifies the operation to be performed as */
    /*           follows: */

    /*              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y. */

    /*              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y. */

    /*              TRANS = 'C' or 'c'   y := alpha*conjg( A' )*x + beta*y. */

    /*           Unchanged on exit. */

    /*  M      - INTEGER. */
    /*           On entry, M specifies the number of rows of the matrix A. */
    /*           M must be at least zero. */
    /*           Unchanged on exit. */

    /*  N      - INTEGER. */
    /*           On entry, N specifies the number of columns of the matrix A. */
    /*           N must be at least zero. */
    /*           Unchanged on exit. */

    /*  ALPHA  - COMPLEX*16      . */
    /*           On entry, ALPHA specifies the scalar alpha. */
    /*           Unchanged on exit. */

    /*  A      - COMPLEX*16       array of DIMENSION ( LDA, n ). */
    /*           Before entry, the leading m by n part of the array A must */
    /*           contain the matrix of coefficients. */
    /*           Unchanged on exit. */

    /*  LDA    - INTEGER. */
    /*           On entry, LDA specifies the first dimension of A as declared */
    /*           in the calling (sub) program. LDA must be at least */