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

/*! return a transposed column vector */inline _zcovector t(const zrovector& rovec){VERBOSE_REPORT;  zcovector covec(rovec.l);  zcopy_(rovec.l, rovec.array, 1, covec.array, 1);    return _(covec);}

intf2c_zcopy(integer* N,          doublecomplex* X, integer* incX,          doublecomplex* Y, integer* incY){    zcopy_(N, X, incX, Y, incY);    return 0;}

static PyObject* diag(PyObject *self, PyObject *args){    PyObject *F;    matrix *d=NULL;    cholmod_factor *L;#if PY_MAJOR_VERSION >= 3    const char *descr;#else    char *descr;#endif    int k, strt, incx=1, incy, nrows, ncols;    if (!set_options()) return NULL;    if (!PyArg_ParseTuple(args, "O", &F)) return NULL;#if PY_MAJOR_VERSION >= 3    if (!PyCapsule_CheckExact(F) || !(descr = PyCapsule_GetName(F)))        err_CO("F");    if (strncmp(descr, "CHOLMOD FACTOR", 14))        PY_ERR_TYPE("F is not a CHOLMOD factor");    L = (cholmod_factor *) PyCapsule_GetPointer(F, descr);#else    if (!PyCObject_Check(F)) err_CO("F");    descr = PyCObject_GetDesc(F);    if (!descr || strncmp(descr, "CHOLMOD FACTOR", 14))        PY_ERR_TYPE("F is not a CHOLMOD factor");    L = (cholmod_factor *) PyCObject_AsVoidPtr(F);#endif    /* Check factorization */    if (L->xtype == CHOLMOD_PATTERN  || L->minor<L->n || !L->is_ll        || !L->is_super)        PY_ERR(PyExc_ValueError, "F must be a nonsingular supernodal "            "Cholesky factor");    if (!(d = Matrix_New(L->n,1,L->xtype == CHOLMOD_REAL ? DOUBLE :        COMPLEX))) return PyErr_NoMemory();    strt = 0;    for (k=0; k<L->nsuper; k++){	/* x[L->px[k], .... ,L->px[k+1]-1] is a dense lower-triangular	 * nrowx times ncols matrix.  We copy its diagonal to	 * d[strt, ..., strt+ncols-1] */        ncols = (int)((int_t *) L->super)[k+1] -            ((int_t *) L->super)[k];        nrows = (int)((int_t *) L->pi)[k+1] - ((int_t *) L->pi)[k];        incy = nrows+1;        if (MAT_ID(d) == DOUBLE)	    dcopy_(&ncols, ((double *) L->x) + ((int_t *) L->px)[k],                &incy, MAT_BUFD(d)+strt, &incx);        else	    zcopy_(&ncols, ((double complex *) L->x) + ((int_t *) L->px)[k],                &incy, MAT_BUFZ(d)+strt, &incx);        strt += ncols;    }    return (PyObject *)d;}

/*! zhematrix copy constructor */inline zhematrix::zhematrix(const zhematrix& mat)  : m(n){VERBOSE_REPORT;  //////// initialize ////////  n =mat.n;  array =new comple[n*n];  darray =new comple*[n];  for(int i=0; i<n; i++){ darray[i] =&array[i*n]; }    //////// copy ////////  zcopy_(n*n, mat.array, 1, array, 1);}

/*! zgematrix copy constructor */inline zgematrix::zgematrix(const zgematrix& mat){VERBOSE_REPORT;  //////// initialize ////////  m =mat.m;  n =mat.n;  array =new comple[m*n];  darray =new comple*[n];  for(int i=0; i<n; i++){ darray[i] =&array[i*m]; }  //////// copy ////////  zcopy_(m*n, mat.array, 1, array, 1);}

/*! return a transposed column vector */inline _zcovector t(const zrovector& rovec){#ifdef  CPPL_VERBOSE  std::cerr << "# [MARK] t(const zrovector&)"            << std::endl;#endif//CPPL_VERBOSE    zcovector covec(rovec.L);  zcopy_(rovec.L, rovec.Array, 1, covec.array, 1);    return _(covec);}

void parallelStorageDTCCplx::getField(){    // If process has a part of the field and is stores the outGrid copy relevant field info directly to the out_grid    if(toOutGrid_)    {        for(int kk = 0; kk < toOutGrid_->opSz_[2]; ++kk )        {            for(int jj = 0; jj < toOutGrid_->opSz_[1]; ++jj)            {                zcopy_(toOutGrid_->opSz_[0], &grid_->point(toOutGrid_->loc_[0]+jj*toOutGrid_->addVec1_[0]+kk*toOutGrid_->addVec2_[0], toOutGrid_->loc_[1]+jj*toOutGrid_->addVec1_[1]+kk*toOutGrid_->addVec2_[1], toOutGrid_->loc_[2]+jj*toOutGrid_->addVec1_[2]+kk*toOutGrid_->addVec2_[2]), toOutGrid_->stride_, &outGrid_->point(toOutGrid_->locOutGrid_[0]+jj*toOutGrid_->addVec1_[0]+kk*toOutGrid_->addVec2_[0], toOutGrid_->locOutGrid_[1]+jj*toOutGrid_->addVec1_[1]+kk*toOutGrid_->addVec2_[1], toOutGrid_->locOutGrid_[2]+jj*toOutGrid_->addVec1_[2]+kk*toOutGrid_->addVec2_[2]), toOutGrid_->strideOutGrid_);            }        }    }    // If the process is a slave process not holding the outGrid then copy the field information to a vector and send it to master    if(slave_)    {        for(int kk = 0; kk < slave_->opSz_[2]; ++kk )            for(int jj = 0; jj < slave_->opSz_[1]; ++jj)                zcopy_(slave_->opSz_[0], &grid_->point(slave_->loc_[0]+jj*slave_->addVec1_[0]+kk*slave_->addVec2_[0], slave_->loc_[1]+jj*slave_->addVec1_[1]+kk*slave_->addVec2_[1], slave_->loc_[2]+jj*slave_->addVec1_[2]+kk*slave_->addVec2_[2]), slave_->stride_, &scratch_[ slave_->opSz_[0]*(jj + kk*slave_->opSz_[1]) ], 1);        gridComm_->send(slave_->masterProc_, gridComm_->cantorTagGen(gridComm_->rank(), slave_->masterProc_, 1, 0), scratch_);    }    // If master then for each slave recv the information and copy it to outGrid    if(masterBool_)    {        for(auto & slave : master_)        {            gridComm_->recv(slave->slaveProc_, gridComm_->cantorTagGen(slave->slaveProc_, gridComm_->rank(), 1, 0), scratch_);            for(int kk = 0; kk < slave->sz_[2]; ++kk)            {                for(int jj = 0; jj < slave->sz_[1]; ++jj)                {                    zcopy_(slave->sz_[0], &scratch_[(jj + slave->sz_[1] * kk) * slave->sz_[0] ], 1, &outGrid_->point(slave->addVec1_[0]*jj+slave->addVec2_[0]*kk+slave->loc_[0], slave->addVec1_[1]*jj+slave->addVec2_[1]*kk+slave->loc_[1], slave->addVec1_[2]*jj+slave->addVec2_[2]*kk+slave->loc_[2]), slave->stride_);                }            }        }    }    return;}

/*! zgbmatrix copy constructor */inline zgbmatrix::zgbmatrix(const zgbmatrix& mat){VERBOSE_REPORT;  //////// initialize ////////  m =mat.m;  n =mat.n;  kl =mat.kl;  ku =mat.ku;  array =new comple[(kl+ku+1)*n];  darray =new comple*[n];  for(int i=0; i<n; i++){ darray[i] =&array[i*(kl+ku+1)]; }  //////// copy ////////  zcopy_((kl+ku+1)*n, mat.array, 1, array, 1);}

/*! make a deep copy of the matrix */inline void zgbmatrix::copy(const zgbmatrix& mat){VERBOSE_REPORT;  m =mat.m;  n =mat.n;  kl =mat.kl;  ku =mat.ku;  delete [] array;  array =new comple[(mat.kl+mat.ku+1)*mat.n];  delete [] darray;  darray =new comple*[n];  for(int i=0; i<n; i++){ darray[i] =&array[i*(kl+ku+1)]; }    zcopy_((mat.kl+mat.ku+1)*mat.n, mat.array, 1, array, 1);}

void parallelStorageFreqDTCCplx::fieldIn(cplx* fftFact){    if(!fieldInFreq_)        return;    // Copy the field information into a vector    for(int jj = 0; jj < fieldInFreq_->sz_[2]; ++jj)    {        for(int ii = 0; ii < fieldInFreq_->sz_[1]; ++ii)        {            zcopy_(fieldInFreq_->sz_[0], &grid_->point(fieldInFreq_->loc_[0]+ii*fieldInFreq_->addVec1_[0]+jj*fieldInFreq_->addVec2_[0], fieldInFreq_->loc_[1]+ii*fieldInFreq_->addVec1_[1]+jj*fieldInFreq_->addVec2_[1], fieldInFreq_->loc_[2]+ii*fieldInFreq_->addVec1_[2]+jj*fieldInFreq_->addVec2_[2]), fieldInFreq_->stride_, &fIn_[ (ii*fieldInFreq_->sz_[2] + jj)*fieldInFreq_->sz_[0]]  , 1 );        }    }    // Take an outer product of the prefactor vector and the field vectors to get the discrete Fourier Transform at all points    zgerc_(nfreq_, fieldInFreq_->sz_[0]*fieldInFreq_->sz_[1]*fieldInFreq_->sz_[2], ONE_, fftFact, 1, fIn_.data(), 1, outGrid_->data(), nfreq_);}

/*! zhematrix copy constructor */inline zhematrix::zhematrix(const zhematrix& mat)  : n(N), array(Array), darray(Darray){#ifdef  CPPL_VERBOSE  std::cerr << "# [MARK] zhematrix::zhematrix(const zhematrix&)"            << std::endl;#endif//CPPL_VERBOSE    //////// initialize ////////  N =mat.N;  Array =new std::complex<double>[N*N];  Darray =new std::complex<double>*[N];  for(int i=0; i<N; i++){ Darray[i] =&Array[i*N]; }    //////// copy ////////  zcopy_(N*N, mat.Array, 1, Array, 1);  #ifdef  CPPL_DEBUG  std::cerr << "# [NOTE] zhematrix::zhematrix(const zhematrix&) "            << "A new matrix at " << Array << " has been made." << std::endl;#endif//CPPL_DEBUG}

//.........这里部分代码省略.........    (x,y) is taken to be |x| + |y|.       =====================================================================          Decode and test the input parameters          Parameter adjustments */    /* Table of constant values */    static doublecomplex c_b2 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer t_dim1, t_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, 	    i__2, i__3, i__4, i__5;    doublereal d__1, d__2, d__3;    doublecomplex z__1, z__2;    /* Builtin functions */    double d_imag(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    static logical allv;    static doublereal unfl, ovfl, smin;    static logical over;    static integer i__, j, k;    static doublereal scale;    extern logical lsame_(char *, char *);    static doublereal remax;    static logical leftv, bothv;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *);    static logical somev;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), dlabad_(doublereal *, doublereal *);    static integer ii, ki;    extern doublereal dlamch_(char *);    static integer is;    extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_(	    integer *, doublereal *, doublecomplex *, integer *);    extern integer izamax_(integer *, doublecomplex *, integer *);    static logical rightv;    extern doublereal dzasum_(integer *, doublecomplex *, integer *);    static doublereal smlnum;    extern /* Subroutine */ int zlatrs_(char *, char *, char *, char *, 	    integer *, doublecomplex *, integer *, doublecomplex *, 	    doublereal *, doublereal *, integer *);    static doublereal ulp;#define t_subscr(a_1,a_2) (a_2)*t_dim1 + a_1#define t_ref(a_1,a_2) t[t_subscr(a_1,a_2)]#define vl_subscr(a_1,a_2) (a_2)*vl_dim1 + a_1#define vl_ref(a_1,a_2) vl[vl_subscr(a_1,a_2)]#define vr_subscr(a_1,a_2) (a_2)*vr_dim1 + a_1#define vr_ref(a_1,a_2) vr[vr_subscr(a_1,a_2)]    --select;    t_dim1 = *ldt;    t_offset = 1 + t_dim1 * 1;    t -= t_offset;    vl_dim1 = *ldvl;    vl_offset = 1 + vl_dim1 * 1;    vl -= vl_offset;    vr_dim1 = *ldvr;    vr_offset = 1 + vr_dim1 * 1;    vr -= vr_offset;

 int zggglm_(int *n, int *m, int *p, 	doublecomplex *a, int *lda, doublecomplex *b, int *ldb, 	doublecomplex *d__, doublecomplex *x, doublecomplex *y, doublecomplex 	*work, int *lwork, int *info){    /* System generated locals */    int a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4;    doublecomplex z__1;    /* Local variables */    int i__, nb, np, nb1, nb2, nb3, nb4, lopt;    extern  int zgemv_(char *, int *, int *, 	    doublecomplex *, doublecomplex *, int *, doublecomplex *, 	    int *, doublecomplex *, doublecomplex *, int *), 	    zcopy_(int *, doublecomplex *, int *, doublecomplex *, 	    int *), xerbla_(char *, int *);    extern int ilaenv_(int *, char *, char *, int *, int *, 	    int *, int *);    extern  int zggqrf_(int *, int *, int *, 	    doublecomplex *, int *, doublecomplex *, doublecomplex *, 	    int *, doublecomplex *, doublecomplex *, int *, int *)	    ;    int lwkmin, lwkopt;    int lquery;    extern  int zunmqr_(char *, char *, int *, int *, 	    int *, doublecomplex *, int *, doublecomplex *, 	    doublecomplex *, int *, doublecomplex *, int *, int *), zunmrq_(char *, char *, int *, int *, 	    int *, doublecomplex *, int *, doublecomplex *, 	    doublecomplex *, int *, doublecomplex *, int *, int *), ztrtrs_(char *, char *, char *, int *, 	    int *, doublecomplex *, int *, doublecomplex *, int *, 	     int *);/*  -- LAPACK driver routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZGGGLM solves a general Gauss-Markov linear model (GLM) problem: *//*          minimize || y ||_2   subject to   d = A*x + B*y *//*              x *//*  where A is an N-by-M matrix, B is an N-by-P matrix, and d is a *//*  given N-vector. It is assumed that M <= N <= M+P, and *//*             rank(A) = M    and    rank( A B ) = N. *//*  Under these assumptions, the constrained equation is always *//*  consistent, and there is a unique solution x and a minimal 2-norm *//*  solution y, which is obtained using a generalized QR factorization *//*  of the matrices (A, B) given by *//*     A = Q*(R),   B = Q*T*Z. *//*           (0) *//*  In particular, if matrix B is square nonsingular, then the problem *//*  GLM is equivalent to the following weighted linear least squares *//*  problem *//*               minimize || inv(B)*(d-A*x) ||_2 *//*                   x *//*  where inv(B) denotes the inverse of B. *//*  Arguments *//*  ========= *//*  N       (input) INTEGER *//*          The number of rows of the matrices A and B.  N >= 0. *//*  M       (input) INTEGER *//*          The number of columns of the matrix A.  0 <= M <= N. *//*  P       (input) INTEGER *//*          The number of columns of the matrix B.  P >= N-M. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,M) *//*          On entry, the N-by-M matrix A. *//*          On exit, the upper triangular part of the array A contains *//*          the M-by-M upper triangular matrix R. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A. LDA >= MAX(1,N). *//*  B       (input/output) COMPLEX*16 array, dimension (LDB,P) *//*          On entry, the N-by-P matrix B. *//*          On exit, if N <= P, the upper triangle of the subarray *//*          B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T; *//*          if N > P, the elements on and above the (N-P)th subdiagonal *//*          contain the N-by-P upper trapezoidal matrix T. *//*  LDB     (input) INTEGER *///.........这里部分代码省略.........

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

/* Subroutine */ int znaitr_(integer *ido, char *bmat, integer *n, integer *k,	 integer *np, integer *nb, doublecomplex *resid, doublereal *rnorm, 	doublecomplex *v, integer *ldv, doublecomplex *h__, integer *ldh, 	integer *ipntr, doublecomplex *workd, integer *info, ftnlen bmat_len){    /* Initialized data */    static logical first = TRUE_;    /* System generated locals */    integer h_dim1, h_offset, v_dim1, v_offset, i__1, i__2, i__3;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1;    /* Builtin functions */    double d_imag(doublecomplex *), sqrt(doublereal);    /* Local variables */    static integer i__, j;    static real t0, t1, t2, t3, t4, t5;    static integer jj, ipj, irj, ivj;    static doublereal ulp, tst1;    static integer ierr, iter;    static doublereal unfl, ovfl;    static integer itry;    static doublereal temp1;    static logical orth1, orth2, step3, step4;    static doublereal betaj;    static integer infol;    static doublecomplex cnorm;    extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    static doublereal rtemp[2];    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *, ftnlen);    static doublereal wnorm;    extern /* Subroutine */ int dvout_(integer *, integer *, doublereal *, 	    integer *, char *, ftnlen), zcopy_(integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *), ivout_(integer *, integer 	    *, integer *, integer *, char *, ftnlen), zaxpy_(integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *), zmout_(integer *, integer *, integer *, doublecomplex 	    *, integer *, integer *, char *, ftnlen), zvout_(integer *, 	    integer *, doublecomplex *, integer *, char *, ftnlen);    extern doublereal dlapy2_(doublereal *, doublereal *);    extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *);    static doublereal rnorm1;    extern /* Subroutine */ int zgetv0_(integer *, char *, integer *, logical 	    *, integer *, integer *, doublecomplex *, integer *, 	    doublecomplex *, doublereal *, integer *, doublecomplex *, 	    integer *, ftnlen);    extern doublereal dlamch_(char *, ftnlen);    extern /* Subroutine */ int second_(real *), zdscal_(integer *, 	    doublereal *, doublecomplex *, integer *);    static logical rstart;    static integer msglvl;    static doublereal smlnum;    extern doublereal zlanhs_(char *, integer *, doublecomplex *, integer *, 	    doublecomplex *, ftnlen);    extern /* Subroutine */ int zlascl_(char *, integer *, integer *, 	    doublereal *, doublereal *, integer *, integer *, doublecomplex *,	     integer *, integer *, ftnlen);/*     %----------------------------------------------------% *//*     | Include files for debugging and timing information | *//*     %----------------------------------------------------% *//* /SCCS Information: @(#) *//* FILE: debug.h   SID: 2.3   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %---------------------------------% *//*     | See debug.doc for documentation | *//*     %---------------------------------% *//*     %------------------% *//*     | Scalar Arguments | *//*     %------------------% *//*     %--------------------------------% *//*     | See stat.doc for documentation | *//*     %--------------------------------% *//* /SCCS Information: @(#) *//* FILE: stat.h   SID: 2.2   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %-----------------% *//*     | Array Arguments | *//*     %-----------------% *//*     %------------% *//*     | Parameters | *//*     %------------% *///.........这里部分代码省略.........

//.........这里部分代码省略.........            On entry, the (m - 1) x n matrix C2 if SIDE = 'L', or the               m x (n - 1) matrix C2 if SIDE = 'R'.               On exit, rows 2:m of P*C if SIDE = 'L', or columns 2:m of C*P               if SIDE = 'R'.       LDC     (input) INTEGER               The leading dimension of the arrays C1 and C2.               LDC >= max(1,M).       WORK    (workspace) COMPLEX*16 array, dimension                           (N) if SIDE = 'L'                           (M) if SIDE = 'R'       =====================================================================          Parameter adjustments          Function Body */    /* Table of constant values */    static doublecomplex c_b1 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer c1_dim1, c1_offset, c2_dim1, c2_offset, i__1;    doublecomplex z__1;    /* Local variables */    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zgerc_(integer *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *), 	    zgeru_(integer *, integer *, doublecomplex *, doublecomplex *, 	    integer *, doublecomplex *, integer *, doublecomplex *, integer *)	    , zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, 	    integer *), zaxpy_(integer *, doublecomplex *, doublecomplex *, 	    integer *, doublecomplex *, integer *), zlacgv_(integer *, 	    doublecomplex *, integer *);#define V(I) v[(I)-1]#define WORK(I) work[(I)-1]#define C2(I,J) c2[(I)-1 + ((J)-1)* ( *ldc)]#define C1(I,J) c1[(I)-1 + ((J)-1)* ( *ldc)]    if (min(*m,*n) == 0 || tau->r == 0. && tau->i == 0.) {	return 0;    }    if (lsame_(side, "L")) {/*        w :=  conjg( C1 + v' * C2 ) */	zcopy_(n, &C1(1,1), ldc, &WORK(1), &c__1);	zlacgv_(n, &WORK(1), &c__1);	i__1 = *m - 1;	zgemv_("Conjugate transpose", &i__1, n, &c_b1, &C2(1,1), ldc, &		V(1), incv, &c_b1, &WORK(1), &c__1);/*        [ C1 ] := [ C1 ] - tau* [ 1 ] * w'             [ C2 ]    [ C2 ]        [ v ] */	zlacgv_(n, &WORK(1), &c__1);	z__1.r = -tau->r, z__1.i = -tau->i;	zaxpy_(n, &z__1, &WORK(1), &c__1, &C1(1,1), ldc);	i__1 = *m - 1;	z__1.r = -tau->r, z__1.i = -tau->i;	zgeru_(&i__1, n, &z__1, &V(1), incv, &WORK(1), &c__1, &C2(1,1), 		ldc);    } else if (lsame_(side, "R")) {/*        w := C1 + C2 * v */	zcopy_(m, &C1(1,1), &c__1, &WORK(1), &c__1);	i__1 = *n - 1;	zgemv_("No transpose", m, &i__1, &c_b1, &C2(1,1), ldc, &V(1), 		incv, &c_b1, &WORK(1), &c__1);/*        [ C1, C2 ] := [ C1, C2 ] - tau* w * [ 1 , v'] */	z__1.r = -tau->r, z__1.i = -tau->i;	zaxpy_(m, &z__1, &WORK(1), &c__1, &C1(1,1), &c__1);	i__1 = *n - 1;	z__1.r = -tau->r, z__1.i = -tau->i;	zgerc_(m, &i__1, &z__1, &WORK(1), &c__1, &V(1), incv, &C2(1,1), 		ldc);    }    return 0;/*     End of ZLATZM */} /* zlatzm_ */

/* Subroutine */ int zdrvpb_(logical *dotype, integer *nn, integer *nval, 	integer *nrhs, doublereal *thresh, logical *tsterr, integer *nmax, 	doublecomplex *a, doublecomplex *afac, doublecomplex *asav, 	doublecomplex *b, doublecomplex *bsav, doublecomplex *x, 	doublecomplex *xact, doublereal *s, doublecomplex *work, doublereal *	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 */    static integer ldab;    static char fact[1];    static integer ioff, mode, koff;    static doublereal amax;    static char path[3];    static integer imat, info;    static char dist[1], uplo[1], type__[1];    static integer nrun, i__, k, n, ifact, nfail, iseed[4], nfact;    extern doublereal dget06_(doublereal *, doublereal *);    static integer kdval[4];    extern logical lsame_(char *, char *);    static char equed[1];    static integer nbmin;    static doublereal rcond, roldc, scond;    static integer nimat;    static doublereal anorm;    extern /* Subroutine */ int zget04_(integer *, integer *, doublecomplex *,	     integer *, doublecomplex *, integer *, doublereal *, doublereal *	    );    static logical equil;    extern /* Subroutine */ int zpbt01_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublereal *, doublereal *), zpbt02_(char *, integer *, 	    integer *, integer *, doublecomplex *, integer *, doublecomplex *,	     integer *, doublecomplex *, integer *, doublereal *, doublereal *	    ), zpbt05_(char *, integer *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublereal *, doublereal *, doublereal *);    static integer iuplo, izero, i1, i2, k1, nerrs;    static logical zerot;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zpbsv_(char *, integer *, integer *, 	    integer *, doublecomplex *, integer *, doublecomplex *, integer *,	     integer *), zswap_(integer *, doublecomplex *, integer *,	     doublecomplex *, integer *);    static char xtype[1];    extern /* Subroutine */ int zlatb4_(char *, integer *, integer *, integer 	    *, char *, integer *, integer *, doublereal *, integer *, 	    doublereal *, char *), aladhd_(integer *, 	    char *);    static integer kd, nb, in, kl;    extern /* Subroutine */ int alaerh_(char *, char *, integer *, integer *, 	    char *, integer *, integer *, integer *, integer *, integer *, 	    integer *, integer *, integer *, integer *);    static logical prefac;    static integer iw, ku, nt;    static doublereal rcondc;    static logical nofact;    static char packit[1];    static integer iequed;    extern doublereal zlanhb_(char *, char *, integer *, integer *, 	    doublecomplex *, integer *, doublereal *), 	    zlange_(char *, integer *, integer *, doublecomplex *, integer *, 	    doublereal *);    extern /* Subroutine */ int zlaqhb_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublereal *, doublereal *, 	    doublereal *, char *), alasvm_(char *, integer *, 	    integer *, integer *, integer *);    static doublereal cndnum;    extern /* Subroutine */ int zlaipd_(integer *, doublecomplex *, integer *,	     integer *);    static doublereal ainvnm;    extern /* Subroutine */ int xlaenv_(integer *, integer *), zlacpy_(char *,	     integer *, integer *, doublecomplex *, integer *, doublecomplex *//.........这里部分代码省略.........

/* Subroutine */ int zchkgt_(logical *dotype, integer *nn, integer *nval, 	integer *nns, integer *nsval, doublereal *thresh, logical *tsterr, 	doublecomplex *a, doublecomplex *af, doublecomplex *b, doublecomplex *	x, doublecomplex *xact, doublecomplex *work, doublereal *rwork, 	integer *iwork, integer *nout){    /* Initialized data */    static integer iseedy[4] = { 0,0,0,1 };    static char transs[1*3] = "N" "T" "C";    /* Format strings */    static char fmt_9999[] = "(12x,/002N =/002,i5,/002,/002,10x,/002 type"	    " /002,i2,/002, test(/002,i2,/002) = /002,g12.5)";    static char fmt_9997[] = "(/002 NORM ='/002,a1,/002', N =/002,i5,/002"	    ",/002,10x,/002 type /002,i2,/002, test(/002,i2,/002) = /002,g12."	    "5)";    static char fmt_9998[] = "(/002 TRANS='/002,a1,/002', N =/002,i5,/002, N"	    "RHS=/002,i3,/002, type /002,i2,/002, test(/002,i2,/002) = /002,g"	    "12.5)";    /* System generated locals */    integer i__1, i__2, i__3, i__4, i__5;    doublereal d__1, d__2;    /* Builtin functions */    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);    /* Local variables */    integer i__, j, k, m, n;    doublecomplex z__[3];    integer in, kl, ku, ix, lda;    doublereal cond;    integer mode, koff, imat, info;    char path[3], dist[1];    integer irhs, nrhs;    char norm[1], type__[1];    integer nrun;    extern /* Subroutine */ int alahd_(integer *, char *);    integer nfail, iseed[4];    extern doublereal dget06_(doublereal *, doublereal *);    doublereal rcond;    integer nimat;    doublereal anorm;    integer itran;    extern /* Subroutine */ int zget04_(integer *, integer *, doublecomplex *, 	     integer *, doublecomplex *, integer *, doublereal *, doublereal *);    char trans[1];    integer izero, nerrs;    extern /* Subroutine */ int zgtt01_(integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublereal *, doublereal *), zgtt02_(char *, integer *, 	     integer *, doublecomplex *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublereal *, doublereal *), zgtt05_(char *, integer *, 	    integer *, doublecomplex *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublereal *, 	    doublereal *);    logical zerot;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zlatb4_(char *, integer *, integer *, 	     integer *, char *, integer *, integer *, doublereal *, integer *, 	     doublereal *, char *), alaerh_(char *, 	    char *, integer *, integer *, char *, integer *, integer *, 	    integer *, integer *, integer *, integer *, integer *, integer *, 	    integer *);    doublereal rcondc, rcondi;    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *), alasum_(char *, integer *, integer *, 	     integer *, integer *);    doublereal rcondo, ainvnm;    logical trfcon;    extern /* Subroutine */ int zerrge_(char *, integer *);    extern doublereal zlangt_(char *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *);    extern /* Subroutine */ int zlagtm_(char *, integer *, integer *, 	    doublereal *, doublecomplex *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *, doublereal *, doublecomplex *, 	    integer *), zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern doublereal dzasum_(integer *, doublecomplex *, integer *);    extern /* Subroutine */ int zgtcon_(char *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *, 	    doublereal *, doublereal *, doublecomplex *, integer *), 	    zlatms_(integer *, integer *, char *, integer *, char *, 	    doublereal *, integer *, doublereal *, doublereal *, integer *, 	    integer *, char *, doublecomplex *, integer *, doublecomplex *, 	    integer *), zlarnv_(integer *, integer *, 	    integer *, doublecomplex *);    doublereal result[7];    extern /* Subroutine */ int zgtrfs_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublereal *, doublereal *, doublecomplex *, doublereal *, 	    integer *), zgttrf_(integer *, doublecomplex *, //.........这里部分代码省略.........

/* Subroutine */ int zgetv0_(integer *ido, char *bmat, integer *itry, logical 	*initv, integer *n, integer *j, doublecomplex *v, integer *ldv, 	doublecomplex *resid, doublereal *rnorm, integer *ipntr, 	doublecomplex *workd, integer *ierr, ftnlen bmat_len){    /* Initialized data */    static logical inits = TRUE_;    /* System generated locals */    integer v_dim1, v_offset, i__1, i__2;    doublereal d__1, d__2;    doublecomplex z__1;    /* Local variables */    static real t0, t1, t2, t3;    static integer jj, iter;    static logical orth;    static integer iseed[4], idist;    static doublecomplex cnorm;    extern /* Double Complex */ void zdotc_(doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    static logical first;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *, ftnlen), 	    dvout_(integer *, integer *, doublereal *, integer *, char *, 	    ftnlen), zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zvout_(integer *, integer *, 	    doublecomplex *, integer *, char *, ftnlen);    extern doublereal dlapy2_(doublereal *, doublereal *), dznrm2_(integer *, 	    doublecomplex *, integer *);    static doublereal rnorm0;    extern /* Subroutine */ int arscnd_(real *);    static integer msglvl;    extern /* Subroutine */ int zlarnv_(integer *, integer *, integer *, 	    doublecomplex *);/*     %----------------------------------------------------% *//*     | Include files for debugging and timing information | *//*     %----------------------------------------------------% *//* /SCCS Information: @(#) *//* FILE: debug.h   SID: 2.3   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %---------------------------------% *//*     | See debug.doc for documentation | *//*     %---------------------------------% *//*     %------------------% *//*     | Scalar Arguments | *//*     %------------------% *//*     %--------------------------------% *//*     | See stat.doc for documentation | *//*     %--------------------------------% *//* /SCCS Information: @(#) *//* FILE: stat.h   SID: 2.2   DATE OF SID: 11/16/95   RELEASE: 2 *//*     %-----------------% *//*     | Array Arguments | *//*     %-----------------% *//*     %------------% *//*     | Parameters | *//*     %------------% *//*     %------------------------% *//*     | Local Scalars & Arrays | *//*     %------------------------% *//*     %----------------------% *//*     | External Subroutines | *//*     %----------------------% *//*     %--------------------% *//*     | External Functions | *//*     %--------------------% *//*     %-----------------% *//*     | Data Statements | *//*     %-----------------% */    /* Parameter adjustments */    --workd;    --resid;    v_dim1 = *ldv;    v_offset = 1 + v_dim1;    v -= v_offset;    --ipntr;//.........这里部分代码省略.........

//.........这里部分代码省略.........    RWORK   (workspace) DOUBLE PRECISION array, dimension (N)       INFO    (output) INTEGER               = 0:  successful exit               < 0:  if INFO = -i, the i-th argument had an illegal value               > 0:  if INFO = i, and i is                     <= N:  the leading minor of order i of A is                            not positive definite, so the factorization                            could not be completed, and the solution has not                            been computed. RCOND = 0 is returned.                     = N+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.       =====================================================================          Test the input parameters.          Parameter adjustments */    /* Table of constant values */    static integer c__1 = 1;        /* System generated locals */    integer b_dim1, b_offset, x_dim1, x_offset, i__1;    /* Local variables */    extern logical lsame_(char *, char *);    static doublereal anorm;    extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 	    doublereal *, integer *), zcopy_(integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *);    extern doublereal dlamch_(char *);    static logical nofact;    extern /* Subroutine */ int xerbla_(char *, integer *);    extern doublereal zlanht_(char *, integer *, doublereal *, doublecomplex *	    );    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *), 	    zptcon_(integer *, doublereal *, doublecomplex *, doublereal *, 	    doublereal *, doublereal *, integer *), zptrfs_(char *, integer *,	     integer *, doublereal *, doublecomplex *, doublereal *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublereal *, doublereal *, doublecomplex *, 	    doublereal *, integer *), zpttrf_(integer *, doublereal *,	     doublecomplex *, integer *), zpttrs_(char *, integer *, integer *	    , doublereal *, doublecomplex *, doublecomplex *, integer *, 	    integer *);    --d__;    --e;    --df;    --ef;    b_dim1 = *ldb;    b_offset = 1 + b_dim1 * 1;    b -= b_offset;    x_dim1 = *ldx;    x_offset = 1 + x_dim1 * 1;    x -= x_offset;    --ferr;    --berr;    --work;

/* Subroutine */ int zlarz_(char *side, integer *m, integer *n, integer *l, 	doublecomplex *v, integer *incv, doublecomplex *tau, doublecomplex *	c__, integer *ldc, doublecomplex *work){    /* System generated locals */    integer c_dim1, c_offset;    doublecomplex z__1;    /* Local variables *//*  -- LAPACK routine (version 3.2) -- *//*     November 2006 *//*  Purpose *//*  ======= *//*  ZLARZ applies a complex elementary reflector H to a complex *//*  M-by-N matrix C, from either the left or the right. H is represented *//*  in the form *//*        H = I - tau * v * v' *//*  where tau is a complex scalar and v is a complex vector. *//*  If tau = 0, then H is taken to be the unit matrix. *//*  To apply H' (the conjugate transpose of H), supply conjg(tau) instead *//*  tau. *//*  H is a product of k elementary reflectors as returned by ZTZRZF. *//*  Arguments *//*  ========= *//*  SIDE    (input) CHARACTER*1 *//*          = 'L': form  H * C *//*          = 'R': form  C * H *//*  M       (input) INTEGER *//*          The number of rows of the matrix C. *//*  N       (input) INTEGER *//*          The number of columns of the matrix C. *//*  L       (input) INTEGER *//*          The number of entries of the vector V containing *//*          the meaningful part of the Householder vectors. *//*          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0. *//*  V       (input) COMPLEX*16 array, dimension (1+(L-1)*abs(INCV)) *//*          The vector v in the representation of H as returned by *//*          ZTZRZF. V is not used if TAU = 0. *//*  INCV    (input) INTEGER *//*          The increment between elements of v. INCV <> 0. *//*  TAU     (input) COMPLEX*16 *//*          The value tau in the representation of H. *//*  C       (input/output) COMPLEX*16 array, dimension (LDC,N) *//*          On entry, the M-by-N matrix C. *//*          On exit, C is overwritten by the matrix H * C if SIDE = 'L', *//*          or C * H if SIDE = 'R'. *//*  LDC     (input) INTEGER *//*          The leading dimension of the array C. LDC >= max(1,M). *//*  WORK    (workspace) COMPLEX*16 array, dimension *//*                         (N) if SIDE = 'L' *//*                      or (M) if SIDE = 'R' *//*  Further Details *//*  =============== *//*  Based on contributions by *//*    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA *//*  ===================================================================== */    /* Parameter adjustments */    --v;    c_dim1 = *ldc;    c_offset = 1 + c_dim1;    c__ -= c_offset;    --work;    /* Function Body */    if (lsame_(side, "L")) {/*        Form  H * C */	if (tau->r != 0. || tau->i != 0.) {/*           w( 1:n ) = conjg( C( 1, 1:n ) ) */	    zcopy_(n, &c__[c_offset], ldc, &work[1], &c__1);	    zlacgv_(n, &work[1], &c__1);/*           w( 1:n ) = conjg( w( 1:n ) + C( m-l+1:m, 1:n )' * v( 1:l ) ) *///.........这里部分代码省略.........

/* Subroutine */ int zgtsvx_(char *fact, char *trans, integer *n, integer *	nrhs, doublecomplex *dl, doublecomplex *d__, doublecomplex *du, 	doublecomplex *dlf, doublecomplex *df, doublecomplex *duf, 	doublecomplex *du2, integer *ipiv, doublecomplex *b, integer *ldb, 	doublecomplex *x, integer *ldx, doublereal *rcond, doublereal *ferr, 	doublereal *berr, doublecomplex *work, doublereal *rwork, integer *	info){    /* System generated locals */    integer b_dim1, b_offset, x_dim1, x_offset, i__1;    /* Local variables */    char norm[1];    extern logical lsame_(char *, char *);    doublereal anorm;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    extern doublereal dlamch_(char *);    logical nofact;    extern /* Subroutine */ int xerbla_(char *, integer *);    extern doublereal zlangt_(char *, integer *, doublecomplex *, 	    doublecomplex *, doublecomplex *);    logical notran;    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *), 	    zgtcon_(char *, integer *, doublecomplex *, doublecomplex *, 	    doublecomplex *, doublecomplex *, integer *, doublereal *, 	    doublereal *, doublecomplex *, integer *), zgtrfs_(char *, 	     integer *, integer *, doublecomplex *, doublecomplex *, 	    doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublereal *, 	    doublecomplex *, doublereal *, integer *), zgttrf_(	    integer *, doublecomplex *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *, integer *), zgttrs_(char *, integer *, 	     integer *, doublecomplex *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, integer *);/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZGTSVX uses the LU factorization to compute the solution to a complex *//*  system of linear equations A * X = B, A**T * X = B, or A**H * X = B, *//*  where A is a tridiagonal matrix of order N and X and B are N-by-NRHS *//*  matrices. *//*  Error bounds on the solution and a condition estimate are also *//*  provided. *//*  Description *//*  =========== *//*  The following steps are performed: *//*  1. If FACT = 'N', the LU decomposition is used to factor the matrix A *//*     as A = L * U, where L is a product of permutation and unit lower *//*     bidiagonal matrices and U is upper triangular with nonzeros in *//*     only the main diagonal and first two superdiagonals. *//*  2. If some U(i,i)=0, so that U is exactly singular, then the routine *//*     returns with INFO = i. Otherwise, the factored form of A is used *//*     to estimate the condition number of the matrix A.  If the *//*     reciprocal of the condition number is less than machine precision, *//*     INFO = N+1 is returned as a warning, but the routine still goes on *//*     to solve for X and compute error bounds as described below. *//*  3. The system of equations is solved for X using the factored form *//*     of A. *//*  4. Iterative refinement is applied to improve the computed solution *//*     matrix and calculate error bounds and backward error estimates *//*     for it. *//*  Arguments *//*  ========= *//*  FACT    (input) CHARACTER*1 *//*          Specifies whether or not the factored form of A has been *//*          supplied on entry. *//*          = 'F':  DLF, DF, DUF, DU2, and IPIV contain the factored form *//*                  of A; DL, D, DU, DLF, DF, DUF, DU2 and IPIV will not *//*                  be modified. *//*          = 'N':  The matrix will be copied to DLF, DF, and DUF *//*                  and factored. *//*  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) *///.........这里部分代码省略.........

/* Subroutine */ int zlarfb_(char *side, char *trans, char *direct, char *	storev, integer *m, integer *n, integer *k, doublecomplex *v, integer 	*ldv, doublecomplex *t, integer *ldt, doublecomplex *c__, integer *	ldc, doublecomplex *work, integer *ldwork){    /* System generated locals */    integer c_dim1, c_offset, t_dim1, t_offset, v_dim1, v_offset, work_dim1, 	    work_offset, i__1, i__2, i__3, i__4, i__5;    doublecomplex z__1, z__2;    /* Builtin functions */    void d_cnjg(doublecomplex *, doublecomplex *);    /* Local variables */    integer i__, j;    extern logical lsame_(char *, char *);    integer lastc;    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *);    integer lastv;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), ztrmm_(char *, char *, char *, char *, integer *, integer *, doublecomplex *, doublecomplex *, integer 	    *, doublecomplex *, integer *);    extern integer ilazlc_(integer *, integer *, doublecomplex *, integer *);    extern /* Subroutine */ int zlacgv_(integer *, doublecomplex *, integer *)	    ;    extern integer ilazlr_(integer *, integer *, doublecomplex *, integer *);    char transt[1];/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLARFB applies a complex block reflector H or its transpose H' to a *//*  complex M-by-N matrix C, from either the left or the right. *//*  Arguments *//*  ========= *//*  SIDE    (input) CHARACTER*1 *//*          = 'L': apply H or H' from the Left *//*          = 'R': apply H or H' from the Right *//*  TRANS   (input) CHARACTER*1 *//*          = 'N': apply H (No transpose) *//*          = 'C': apply H' (Conjugate transpose) *//*  DIRECT  (input) CHARACTER*1 *//*          Indicates how H is formed from a product of elementary *//*          reflectors *//*          = 'F': H = H(1) H(2) . . . H(k) (Forward) *//*          = 'B': H = H(k) . . . H(2) H(1) (Backward) *//*  STOREV  (input) CHARACTER*1 *//*          Indicates how the vectors which define the elementary *//*          reflectors are stored: *//*          = 'C': Columnwise *//*          = 'R': Rowwise *//*  M       (input) INTEGER *//*          The number of rows of the matrix C. *//*  N       (input) INTEGER *//*          The number of columns of the matrix C. *//*  K       (input) INTEGER *//*          The order of the matrix T (= the number of elementary *//*          reflectors whose product defines the block reflector). *//*  V       (input) COMPLEX*16 array, dimension *//*                                (LDV,K) if STOREV = 'C' *//*                                (LDV,M) if STOREV = 'R' and SIDE = 'L' *//*                                (LDV,N) if STOREV = 'R' and SIDE = 'R' *//*          The matrix V. See further details. *//*  LDV     (input) INTEGER *//*          The leading dimension of the array V. *//*          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); *//*          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); *//*          if STOREV = 'R', LDV >= K. *//*  T       (input) COMPLEX*16 array, dimension (LDT,K) *//*          The triangular K-by-K matrix T in the representation of the *//*          block reflector. *//*  LDT     (input) INTEGER *//*          The leading dimension of the array T. LDT >= K. *///.........这里部分代码省略.........

/* Subroutine */int zpbsvx_(char *fact, char *uplo, integer *n, integer *kd, integer *nrhs, doublecomplex *ab, integer *ldab, doublecomplex *afb, integer *ldafb, char *equed, doublereal *s, doublecomplex *b, integer *ldb, doublecomplex *x, integer *ldx, doublereal *rcond, doublereal * ferr, doublereal *berr, doublecomplex *work, doublereal *rwork, integer *info){    /* System generated locals */    integer ab_dim1, ab_offset, afb_dim1, afb_offset, b_dim1, b_offset, x_dim1, x_offset, i__1, i__2, i__3, i__4, i__5;    doublereal d__1, d__2;    doublecomplex z__1;    /* Local variables */    integer i__, j, j1, j2;    doublereal amax, smin, smax;    extern logical lsame_(char *, char *);    doublereal scond, anorm;    logical equil, rcequ, upper;    extern /* Subroutine */    int zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, integer *);    extern doublereal dlamch_(char *);    logical nofact;    extern /* Subroutine */    int xerbla_(char *, integer *);    extern doublereal zlanhb_(char *, char *, integer *, integer *, doublecomplex *, integer *, doublereal *);    doublereal bignum;    extern /* Subroutine */    int zlaqhb_(char *, integer *, integer *, doublecomplex *, integer *, doublereal *, doublereal *, doublereal *, char *);    integer infequ;    extern /* Subroutine */    int zpbcon_(char *, integer *, integer *, doublecomplex *, integer *, doublereal *, doublereal *, doublecomplex *, doublereal *, integer *), zlacpy_(char *, integer *, integer *, doublecomplex *, integer *, doublecomplex * , integer *), zpbequ_(char *, integer *, integer *, doublecomplex *, integer *, doublereal *, doublereal *, doublereal *, integer *), zpbrfs_(char *, integer *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, integer *, doublecomplex *, integer *, doublecomplex *, integer * , doublereal *, doublereal *, doublecomplex *, doublereal *, integer *), zpbtrf_(char *, integer *, integer *, doublecomplex *, integer *, integer *);    doublereal smlnum;    extern /* Subroutine */    int zpbtrs_(char *, integer *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, integer *, integer *);    /* -- LAPACK driver routine (version 3.4.1) -- */    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */    /* April 2012 */    /* .. Scalar Arguments .. */    /* .. */    /* .. Array Arguments .. */    /* .. */    /* ===================================================================== */    /* .. Parameters .. */    /* .. */    /* .. Local Scalars .. */    /* .. */    /* .. External Functions .. */    /* .. */    /* .. External Subroutines .. */    /* .. */    /* .. Intrinsic Functions .. */    /* .. */    /* .. Executable Statements .. */    /* Parameter adjustments */    ab_dim1 = *ldab;    ab_offset = 1 + ab_dim1;    ab -= ab_offset;    afb_dim1 = *ldafb;    afb_offset = 1 + afb_dim1;    afb -= afb_offset;    --s;    b_dim1 = *ldb;    b_offset = 1 + b_dim1;    b -= b_offset;    x_dim1 = *ldx;    x_offset = 1 + x_dim1;    x -= x_offset;    --ferr;    --berr;    --work;    --rwork;    /* Function Body */    *info = 0;    nofact = lsame_(fact, "N");    equil = lsame_(fact, "E");    upper = lsame_(uplo, "U");    if (nofact || equil)    {        *(unsigned char *)equed = 'N';        rcequ = FALSE_;    }    else    {        rcequ = lsame_(equed, "Y");        smlnum = dlamch_("Safe minimum");        bignum = 1. / smlnum;    }    /* Test the input parameters. */    if (! nofact && ! equil && ! lsame_(fact, "F"))    {        *info = -1;    }    else if (! upper && ! lsame_(uplo, "L"))    {        *info = -2;    }    else if (*n < 0)    {        *info = -3;    }    else if (*kd < 0)    {        *info = -4;    }//.........这里部分代码省略.........

//.........这里部分代码省略.........            On exit, if INFO = 0, WORK(1) returns the optimal LWORK.       LWORK   (input) INTEGER               The dimension of the array WORK. LWORK >= max(1,N+M+P).               For optimum performance, LWORK >= M+min(N,P)+max(N,P)*NB,               where NB is an upper bound for the optimal blocksizes for               ZGEQRF, CGERQF, ZUNMQR and CUNMRQ.       INFO    (output) INTEGER               = 0:  successful exit.               < 0:  if INFO = -i, the i-th argument had an illegal value.       ===================================================================          Test the input parameters          Parameter adjustments          Function Body */    /* Table of constant values */    static doublecomplex c_b2 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4;    doublereal d__1;    doublecomplex z__1;    /* Local variables */    static integer lopt, i;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *), 	    zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, 	    integer *), ztrsv_(char *, char *, char *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    static integer np;    extern /* Subroutine */ int xerbla_(char *, integer *), zggqrf_(	    integer *, integer *, integer *, doublecomplex *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    doublecomplex *, integer *, integer *), zunmqr_(char *, char *, 	    integer *, integer *, integer *, doublecomplex *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, integer *), zunmrq_(char *, char *, 	    integer *, integer *, integer *, doublecomplex *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, integer *);#define D(I) d[(I)-1]#define X(I) x[(I)-1]#define Y(I) y[(I)-1]#define WORK(I) work[(I)-1]#define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]#define B(I,J) b[(I)-1 + ((J)-1)* ( *ldb)]    *info = 0;    np = min(*n,*p);    if (*n < 0) {	*info = -1;    } else if (*m < 0 || *m > *n) {	*info = -2;    } else if (*p < 0 || *p < *n - *m) {	*info = -3;

/* Subroutine */ int zlacn2_(integer *n, doublecomplex *v, doublecomplex *x, 	doublereal *est, integer *kase, integer *isave){    /* System generated locals */    integer i__1, i__2, i__3;    doublereal d__1, d__2;    doublecomplex z__1;    /* Builtin functions */    double z_abs(doublecomplex *), d_imag(doublecomplex *);    /* Local variables */    integer i__;    doublereal temp, absxi;    integer jlast;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *);    extern integer izmax1_(integer *, doublecomplex *, integer *);    extern doublereal dzsum1_(integer *, doublecomplex *, integer *), dlamch_(	    char *);    doublereal safmin, altsgn, estold;/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLACN2 estimates the 1-norm of a square, complex matrix A. *//*  Reverse communication is used for evaluating matrix-vector products. *//*  Arguments *//*  ========= *//*  N      (input) INTEGER *//*         The order of the matrix.  N >= 1. *//*  V      (workspace) COMPLEX*16 array, dimension (N) *//*         On the final return, V = A*W,  where  EST = norm(V)/norm(W) *//*         (W is not returned). *//*  X      (input/output) COMPLEX*16 array, dimension (N) *//*         On an intermediate return, X should be overwritten by *//*               A * X,   if KASE=1, *//*               A' * X,  if KASE=2, *//*         where A' is the conjugate transpose of A, and ZLACN2 must be *//*         re-called with all the other parameters unchanged. *//*  EST    (input/output) DOUBLE PRECISION *//*         On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be *//*         unchanged from the previous call to ZLACN2. *//*         On exit, EST is an estimate (a lower bound) for norm(A). *//*  KASE   (input/output) INTEGER *//*         On the initial call to ZLACN2, KASE should be 0. *//*         On an intermediate return, KASE will be 1 or 2, indicating *//*         whether X should be overwritten by A * X  or A' * X. *//*         On the final return from ZLACN2, KASE will again be 0. *//*  ISAVE  (input/output) INTEGER array, dimension (3) *//*         ISAVE is used to save variables between calls to ZLACN2 *//*  Further Details *//*  ======= ======= *//*  Contributed by Nick Higham, University of Manchester. *//*  Originally named CONEST, dated March 16, 1988. *//*  Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of *//*  a real or complex matrix, with applications to condition estimation", *//*  ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. *//*  Last modified:  April, 1999 *//*  This is a thread safe version of ZLACON, which uses the array ISAVE *//*  in place of a SAVE statement, as follows: *//*     ZLACON     ZLACN2 *//*      JUMP     ISAVE(1) *//*      J        ISAVE(2) *//*      ITER     ISAVE(3) *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Intrinsic Functions .. *///.........这里部分代码省略.........

//.........这里部分代码省略........./*     Quick return if possible */    if (*n == 0 || *nrhs == 0) {	i__1 = *nrhs;	for (j = 1; j <= i__1; ++j) {	    ferr[j] = 0.;	    berr[j] = 0.;	}	return 0;    }/*     NZ = maximum number of nonzero elements in each row of A, plus 1 */    nz = *n + 1;    eps = dlamch_("Epsilon");    safmin = dlamch_("Safe minimum");    safe1 = nz * safmin;    safe2 = safe1 / eps;/*     Do for each right hand side */    i__1 = *nrhs;    for (j = 1; j <= i__1; ++j) {	count = 1;	lstres = 3.;L20:/*        Loop until stopping criterion is satisfied. *//*        Compute residual R = B - A * X */	zcopy_(n, &b[j * b_dim1 + 1], &c__1, &work[1], &c__1);	z__1.r = -1., z__1.i = -0.;	zhemv_(uplo, n, &z__1, &a[a_offset], lda, &x[j * x_dim1 + 1], &c__1, &		c_b1, &work[1], &c__1);/*        Compute componentwise relative backward error from formula *//*        max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) ) *//*        where abs(Z) is the componentwise absolute value of the matrix *//*        or vector Z.  If the i-th component of the denominator is less *//*        than SAFE2, then SAFE1 is added to the i-th components of the *//*        numerator and denominator before dividing. */	i__2 = *n;	for (i__ = 1; i__ <= i__2; ++i__) {	    i__3 = i__ + j * b_dim1;	    rwork[i__] = (d__1 = b[i__3].r, abs(d__1)) + (d__2 = d_imag(&b[		    i__ + j * b_dim1]), abs(d__2));	}/*        Compute abs(A)*abs(X) + abs(B). */	if (upper) {	    i__2 = *n;	    for (k = 1; k <= i__2; ++k) {		s = 0.;		i__3 = k + j * x_dim1;		xk = (d__1 = x[i__3].r, abs(d__1)) + (d__2 = d_imag(&x[k + j *			 x_dim1]), abs(d__2));		i__3 = k - 1;		for (i__ = 1; i__ <= i__3; ++i__) {		    i__4 = i__ + k * a_dim1;

/* Subroutine */ int zgbrfs_(char *trans, integer *n, integer *kl, integer *	ku, integer *nrhs, doublecomplex *ab, integer *ldab, doublecomplex *	afb, integer *ldafb, integer *ipiv, doublecomplex *b, integer *ldb, 	doublecomplex *x, integer *ldx, doublereal *ferr, doublereal *berr, 	doublecomplex *work, doublereal *rwork, integer *info){    /* System generated locals */    integer ab_dim1, ab_offset, afb_dim1, afb_offset, b_dim1, b_offset, 	    x_dim1, x_offset, i__1, i__2, i__3, i__4, i__5, i__6, i__7;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1;    /* Builtin functions */    double d_imag(doublecomplex *);    /* Local variables */    integer i__, j, k;    doublereal s;    integer kk;    doublereal xk;    integer nz;    doublereal eps;    integer kase;    doublereal safe1, safe2;    extern logical lsame_(char *, char *);    integer isave[3];    extern /* Subroutine */ int zgbmv_(char *, integer *, integer *, integer *, integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *);    integer count;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zaxpy_(integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *), zlacn2_(	    integer *, doublecomplex *, doublecomplex *, doublereal *, 	    integer *, integer *);    extern doublereal dlamch_(char *);    doublereal safmin;    extern /* Subroutine */ int xerbla_(char *, integer *);    logical notran;    char transn[1], transt[1];    doublereal lstres;    extern /* Subroutine */ int zgbtrs_(char *, integer *, integer *, integer 	    *, integer *, doublecomplex *, integer *, integer *, 	    doublecomplex *, integer *, integer *);/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH. *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZGBRFS improves the computed solution to a system of linear *//*  equations when the coefficient matrix is banded, 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) *//*  N       (input) INTEGER *//*          The order 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. *//*  NRHS    (input) INTEGER *//*          The number of right hand sides, i.e., the number of columns *//*          of the matrices B and X.  NRHS >= 0. *//*  AB      (input) COMPLEX*16 array, dimension (LDAB,N) *//*          The original band matrix A, stored 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(n,j+kl). *//*  LDAB    (input) INTEGER *//*          The leading dimension of the array AB.  LDAB >= KL+KU+1. *//*  AFB     (input) COMPLEX*16 array, dimension (LDAFB,N) *//*          Details of the LU factorization of the band matrix A, as *//*          computed by ZGBTRF.  U is stored as an upper triangular band *//*          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and *///.........这里部分代码省略.........

/* 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 zchkhp_(logical *dotype, integer *nn, integer *nval, 	integer *nns, integer *nsval, doublereal *thresh, logical *tsterr, 	integer *nmax, doublecomplex *a, doublecomplex *afac, doublecomplex *	ainv, doublecomplex *b, doublecomplex *x, doublecomplex *xact, 	doublecomplex *work, doublereal *rwork, integer *iwork, integer *nout){    /* Initialized data */    static integer iseedy[4] = { 1988,1989,1990,1991 };    static char uplos[1*2] = "U" "L";    /* Format strings */    static char fmt_9999[] = "(/002 UPLO = '/002,a1,/002', N =/002,i5,/002, "	    "type /002,i2,/002, test /002,i2,/002, ratio =/002,g12.5)";    static char fmt_9998[] = "(/002 UPLO = '/002,a1,/002', N =/002,i5,/002, "	    "NRHS=/002,i3,/002, type /002,i2,/002, test(/002,i2,/002) =/002,g"	    "12.5)";    /* System generated locals */    integer i__1, i__2, i__3, i__4, i__5;    /* Builtin functions */    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);    /* Local variables */    integer i__, j, k, n, i1, i2, in, kl, ku, nt, lda, npp, ioff, mode, imat, 	    info;    char path[3], dist[1];    integer irhs, nrhs;    char uplo[1], type__[1];    integer nrun;    extern /* Subroutine */ int alahd_(integer *, char *);    integer nfail, iseed[4];    extern doublereal dget06_(doublereal *, doublereal *);    extern logical lsame_(char *, char *);    doublereal rcond;    integer nimat;    doublereal anorm;    extern /* Subroutine */ int zget04_(integer *, integer *, doublecomplex *, 	     integer *, doublecomplex *, integer *, doublereal *, doublereal *), zhpt01_(char *, integer *, doublecomplex *, doublecomplex *, 	    integer *, doublecomplex *, integer *, doublereal *, doublereal *);    integer iuplo, izero, nerrs;    extern /* Subroutine */ int zppt02_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublereal *, doublereal *), zppt03_(char *, 	    integer *, doublecomplex *, doublecomplex *, doublecomplex *, 	    integer *, doublereal *, doublereal *, doublereal *);    logical zerot;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zppt05_(char *, integer *, integer *, 	     doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *, doublereal *, doublereal *, 	     doublereal *);    char xtype[1];    extern /* Subroutine */ int zlatb4_(char *, integer *, integer *, integer 	    *, char *, integer *, integer *, doublereal *, integer *, 	    doublereal *, char *), alaerh_(char *, 	    char *, integer *, integer *, char *, integer *, integer *, 	    integer *, integer *, integer *, integer *, integer *, integer *, 	    integer *);    doublereal rcondc;    char packit[1];    extern /* Subroutine */ int alasum_(char *, integer *, integer *, integer 	    *, integer *);    doublereal cndnum;    extern /* Subroutine */ int zlaipd_(integer *, doublecomplex *, integer *, 	     integer *);    logical trfcon;    extern doublereal zlanhp_(char *, char *, integer *, doublecomplex *, 	    doublereal *);    extern /* Subroutine */ int zhpcon_(char *, integer *, doublecomplex *, 	    integer *, doublereal *, doublereal *, doublecomplex *, integer *), zlacpy_(char *, integer *, integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *), zlarhs_(char *, 	    char *, char *, char *, integer *, integer *, integer *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, integer 	    *, doublecomplex *, integer *, integer *, integer *), zlatms_(integer *, integer *, char *, 	    integer *, char *, doublereal *, integer *, doublereal *, 	    doublereal *, integer *, integer *, char *, doublecomplex *, 	    integer *, doublecomplex *, integer *), 	    zhprfs_(char *, integer *, integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublereal *, doublereal *, 	    doublecomplex *, doublereal *, integer *), zhptrf_(char *, 	     integer *, doublecomplex *, integer *, integer *);    doublereal result[8];    extern /* Subroutine */ int zhptri_(char *, integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *), zhptrs_(char *, 	    integer *, integer *, doublecomplex *, integer *, doublecomplex *, 	     integer *, integer *), zerrsy_(char *, integer *)	    ;    /* Fortran I/O blocks */    static cilist io___38 = { 0, 0, 0, fmt_9999, 0 };    static cilist io___41 = { 0, 0, 0, fmt_9998, 0 };    static cilist io___43 = { 0, 0, 0, fmt_9999, 0 };//.........这里部分代码省略.........