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

static void zget_wv(double complex *w, double complex *v,                    double complex *cache, double complex *fvohalf,                    double complex *vooo, double complex *vv_op,                    double complex *t1Thalf, double complex *t2T,                    int nocc, int nvir, int a, int b, int c, int *idx){        const double complex D0 = 0;        const double complex D1 = 1;        const double complex DN1 =-1;        const char TRANS_N = 'N';        const int nmo = nocc + nvir;        const int noo = nocc * nocc;        const size_t nooo = nocc * noo;        const size_t nvoo = nvir * noo;        int i, j, k, n;        double complex *pt2T;        zgemm_(&TRANS_N, &TRANS_N, &noo, &nocc, &nvir,               &D1, t2T+c*nvoo, &noo, vv_op+nocc, &nmo,               &D0, cache, &noo);        zgemm_(&TRANS_N, &TRANS_N, &nocc, &noo, &nocc,               &DN1, t2T+c*nvoo+b*noo, &nocc, vooo+a*nooo, &nocc,               &D1, cache, &nocc);        pt2T = t2T + b * nvoo + a * noo;        for (n = 0, i = 0; i < nocc; i++) {        for (j = 0; j < nocc; j++) {        for (k = 0; k < nocc; k++, n++) {                w[idx[n]] += cache[n];                v[idx[n]] +=(vv_op[i*nmo+j] * t1Thalf[c*nocc+k]                           + pt2T[i*nocc+j] * fvohalf[c*nocc+k]);        } } }}

/*! zgematrix*=_zgematrix operator */inline zgematrix& zgematrix::operator*=(const _zgematrix& mat){#ifdef  CPPL_VERBOSE  std::cerr << "# [MARK] zgematrix::operator*=(const _zgematrix&)"            << std::endl;#endif//CPPL_VERBOSE  #ifdef  CPPL_DEBUG  if(N!=mat.M){    std::cerr << "[ERROR] zgematrix::operator*=(_zgematrix&)" << std::endl              << "These two matrises can not make a product." << std::endl              << "Your input was (" << M << "x" << N << ") *= ("              << mat.M << "x" << mat.N << ")." << std::endl;    exit(1);  }#endif//CPPL_DEBUG    zgematrix newmat( M, mat.N );  zgemm_( 'N', 'N', M, mat.N, N, std::complex<double>(1.0,0.0), Array, M,          mat.Array, mat.M, std::complex<double>(0.0,0.0), newmat.array, M );    swap(*this,newmat);  mat.destroy();  return *this;}

intf2c_zgemm(char* transA, char* transB, integer* M, integer* N, integer* K,          doublecomplex* alpha,          doublecomplex* A, integer* lda,          doublecomplex* B, integer* ldb,          doublecomplex* beta,          doublecomplex* C, integer* ldc){    zgemm_(transA, transB, M, N, K,           alpha, A, lda, B, ldb, beta, C, ldc);    return 0;}

static void dot_ao_dm(double complex *vm, double complex *ao, double complex *dm,                      int nao, int nocc, int ngrids, int bgrids,                      unsigned char *non0table, int *shls_slice, int *ao_loc){        int nbox = (nao+BOXSIZE-1) / BOXSIZE;        char empty[nbox];        int has0 = VXCao_empty_blocks(empty, non0table, shls_slice, ao_loc);        const char TRANS_T = 'T';        const char TRANS_N = 'N';        const double complex Z1 = 1;        double complex beta = 0;        if (has0) {                int box_id, blen, i, j;                size_t b0;                for (box_id = 0; box_id < nbox; box_id++) {                        if (!empty[box_id]) {                                b0 = box_id * BOXSIZE;                                blen = MIN(nao-b0, BOXSIZE);                                zgemm_(&TRANS_N, &TRANS_T, &bgrids, &nocc, &blen,                                       &Z1, ao+b0*ngrids, &ngrids, dm+b0*nocc, &nocc,                                       &beta, vm, &ngrids);                                beta = 1.0;                        }                }                if (beta == 0) { // all empty                        for (i = 0; i < nocc; i++) {                                for (j = 0; j < bgrids; j++) {                                        vm[i*ngrids+j] = 0;                                }                        }                }        } else {                zgemm_(&TRANS_N, &TRANS_T, &bgrids, &nocc, &nao,                       &Z1, ao, &ngrids, dm, &nocc, &beta, vm, &ngrids);        }}

/* conj(vv[n,m]) = ao1[n,ngrids] * conj(ao2[m,ngrids]) */static void dot_ao_ao(double complex *vv, double complex *ao1, double complex *ao2,                      int nao, int ngrids, int bgrids, int hermi,                      unsigned char *non0table, int *shls_slice, int *ao_loc){        int nbox = (nao+BOXSIZE-1) / BOXSIZE;        char empty[nbox];        int has0 = VXCao_empty_blocks(empty, non0table, shls_slice, ao_loc);        const char TRANS_C = 'C';        const char TRANS_N = 'N';        const double complex Z1 = 1;        if (has0) {                int ib, jb, leni, lenj;                int j1 = nbox;                size_t b0i, b0j;                for (ib = 0; ib < nbox; ib++) {                if (!empty[ib]) {                        b0i = ib * BOXSIZE;                        leni = MIN(nao-b0i, BOXSIZE);                        if (hermi) {                                j1 = ib + 1;                        }                        for (jb = 0; jb < j1; jb++) {                        if (!empty[jb]) {                                b0j = jb * BOXSIZE;                                lenj = MIN(nao-b0j, BOXSIZE);                                zgemm_(&TRANS_C, &TRANS_N, &lenj, &leni, &bgrids, &Z1,                                       ao2+b0j*ngrids, &ngrids, ao1+b0i*ngrids, &ngrids,                                       &Z1, vv+b0i*nao+b0j, &nao);                        } }                } }        } else {                zgemm_(&TRANS_C, &TRANS_N, &nao, &nao, &bgrids,                       &Z1, ao2, &ngrids, ao1, &ngrids, &Z1, vv, &nao);        }}

PyObject* gemm(PyObject *self, PyObject *args){  Py_complex alpha;  PyArrayObject* a;  PyArrayObject* b;  Py_complex beta;  PyArrayObject* c;  char transa = 'n';  if (!PyArg_ParseTuple(args, "DOODO|c", &alpha, &a, &b, &beta, &c, &transa))    return NULL;  int m, k, lda, ldb, ldc;  if (transa == 'n')    {      m = PyArray_DIMS(a)[1];      for (int i = 2; i < PyArray_NDIM(a); i++)        m *= PyArray_DIMS(a)[i];      k = PyArray_DIMS(a)[0];      lda = MAX(1, PyArray_STRIDES(a)[0] / PyArray_STRIDES(a)[PyArray_NDIM(a) - 1]);      ldb = MAX(1, PyArray_STRIDES(b)[0] / PyArray_STRIDES(b)[1]);      ldc = MAX(1, PyArray_STRIDES(c)[0] / PyArray_STRIDES(c)[PyArray_NDIM(c) - 1]);    }  else    {      k = PyArray_DIMS(a)[1];      for (int i = 2; i < PyArray_NDIM(a); i++)        k *= PyArray_DIMS(a)[i];      m = PyArray_DIMS(a)[0];      lda = MAX(1, k);      ldb = MAX(1, PyArray_STRIDES(b)[0] / PyArray_STRIDES(b)[PyArray_NDIM(b) - 1]);      ldc = MAX(1, PyArray_STRIDES(c)[0] / PyArray_STRIDES(c)[1]);    }  int n = PyArray_DIMS(b)[0];  if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)    dgemm_(&transa, "n", &m, &n, &k,           &(alpha.real),           DOUBLEP(a), &lda,           DOUBLEP(b), &ldb,           &(beta.real),           DOUBLEP(c), &ldc);  else    zgemm_(&transa, "n", &m, &n, &k,           &alpha,           (void*)COMPLEXP(a), &lda,           (void*)COMPLEXP(b), &ldb,           &beta,           (void*)COMPLEXP(c), &ldc);  Py_RETURN_NONE;}

/*! zgematrix*_zgematrix operator */inline _zgematrix operator*(const zgematrix& matA, const _zgematrix& matB){VERBOSE_REPORT;#ifdef  CPPL_DEBUG  if(matA.n!=matB.m){    ERROR_REPORT;    std::cerr << "These two matrises can not make a product." << std::endl              << "Your input was (" << matA.m << "x" << matA.n << ") * (" << matB.m << "x" << matB.n << ")." << std::endl;    exit(1);  }#endif//CPPL_DEBUG    zgematrix newmat( matA.m, matB.n );  zgemm_( 'n', 'n', matA.m, matB.n, matA.n, comple(1.0,0.0),          matA.array, matA.m, matB.array, matB.m,          comple(0.0,0.0), newmat.array, matA.m );    matB.destroy();  return _(newmat);}

/* * transform bra (without doing conj(mo)), v_{iq} = C_{pi} v_{pq} * s1 to label AO symmetry */int RIhalfmmm_r_s1_bra_noconj(double complex *vout, double complex *vin,                              struct _AO2MOEnvs *envs, int seekdim){        switch (seekdim) {                case 1: return envs->bra_count * envs->nao;                case 2: return envs->nao * envs->nao;        }        const double complex Z0 = 0;        const double complex Z1 = 1;        const char TRANS_N = 'N';        int n2c = envs->nao;        int i_start = envs->bra_start;        int i_count = envs->bra_count;        double complex *mo_coeff = envs->mo_coeff;        zgemm_(&TRANS_N, &TRANS_N, &n2c, &i_count, &n2c,               &Z1, vin, &n2c, mo_coeff+i_start*n2c, &n2c,               &Z0, vout, &n2c);        return 0;}

voidcublasDgemm( char transA, char transB, int m, int n, int k,    double alpha, double *A, int ldA,                           double *B, int ldB,    double beta,  double *C, int ldC ){  double zalpha_[REAL_PART+IMAG_PART+1];  double zbeta_[REAL_PART+IMAG_PART+1];  double *zalpha = &(zalpha_[0]);  double *zbeta = &(zbeta_[0]);  zalpha[REAL_PART] = creal(alpha);  zalpha[IMAG_PART] = cimag(alpha);  zbeta[REAL_PART] = creal(beta);  zbeta[IMAG_PART] = cimag(beta);  zgemm_( &transA, &transB, &m, &n, &k,        zalpha, (double *) A, &ldA,  (double *) B, &ldB,        zbeta,  (double *) C, &ldC );}

/* * transform ket, s1 to label AO symmetry */int RIhalfmmm_r_s1_ket(double complex *vout, double complex *vin,                       struct _AO2MOEnvs *envs, int seekdim){        switch (seekdim) {                case 1: return envs->nao * envs->ket_count;                case 2: return envs->nao * envs->nao;        }        const double complex Z0 = 0;        const double complex Z1 = 1;        const char TRANS_T = 'T';        const char TRANS_N = 'N';        int n2c = envs->nao;        int j_start = envs->ket_start;        int j_count = envs->ket_count;        double complex *mo_coeff = envs->mo_coeff;        zgemm_(&TRANS_T, &TRANS_N, &j_count, &n2c, &n2c,               &Z1, mo_coeff+j_start*n2c, &n2c, vin, &n2c,               &Z0, vout, &j_count);        return 0;}

// Interface to lapack routine zgemm// mm: nrow of A// nn: ncol of B// kk: ncol and nrow of Cvoid lapack_zgemm(int mm, int nn, int kk,                  char transa, char transb,                  dcmplx alpha, dcmplx *AA, dcmplx *BB,                  dcmplx beta, dcmplx *CC){  int lda, ldb, ldc;    if(transa == 'N' || transa == 'n') {    lda = (1 > mm) ? 1 : mm;  } else {    lda = (1 > kk) ? 1 : kk;  }    if(transb == 'N' || transb == 'n') {    ldb = (1 > kk) ? 1 : kk;  } else {    ldb = (1 > nn) ? 1 : nn;  }    ldc = (1 > mm) ? 1 : mm;    zgemm_(&transa, &transb, &mm, &nn, &kk, &alpha, AA, &lda, BB, &ldb, &beta, CC, &ldc); }

/* Subroutine */ int zlahr2_(integer *n, integer *k, integer *nb, 	doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *t, 	integer *ldt, doublecomplex *y, integer *ldy){    /* System generated locals */    integer a_dim1, a_offset, t_dim1, t_offset, y_dim1, y_offset, i__1, i__2, 	    i__3;    doublecomplex z__1;    /* Local variables */    integer i__;    doublecomplex ei;    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 	    doublecomplex *, integer *), zgemm_(char *, char *, integer *, 	    integer *, integer *, doublecomplex *, doublecomplex *, integer *, 	     doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *), zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *), 	    zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, 	    integer *), ztrmm_(char *, char *, char *, char *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *), 	    zaxpy_(integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *), ztrmv_(char *, char *, char *, 	    integer *, doublecomplex *, integer *, doublecomplex *, integer *), zlarfg_(integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *), zlacgv_(integer *, 	    doublecomplex *, integer *), zlacpy_(char *, integer *, integer *, 	     doublecomplex *, integer *, doublecomplex *, integer *);/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLAHR2 reduces the first NB columns of A complex general n-BY-(n-k+1) *//*  matrix A so that elements below the k-th subdiagonal are zero. The *//*  reduction is performed by an unitary similarity transformation *//*  Q' * A * Q. The routine returns the matrices V and T which determine *//*  Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T. *//*  This is an auxiliary routine called by ZGEHRD. *//*  Arguments *//*  ========= *//*  N       (input) INTEGER *//*          The order of the matrix A. *//*  K       (input) INTEGER *//*          The offset for the reduction. Elements below the k-th *//*          subdiagonal in the first NB columns are reduced to zero. *//*          K < N. *//*  NB      (input) INTEGER *//*          The number of columns to be reduced. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,N-K+1) *//*          On entry, the n-by-(n-k+1) general matrix A. *//*          On exit, the elements on and above the k-th subdiagonal in *//*          the first NB columns are overwritten with the corresponding *//*          elements of the reduced matrix; the elements below the k-th *//*          subdiagonal, with the array TAU, represent the matrix Q as a *//*          product of elementary reflectors. The other columns of A are *//*          unchanged. See Further Details. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A.  LDA >= max(1,N). *//*  TAU     (output) COMPLEX*16 array, dimension (NB) *//*          The scalar factors of the elementary reflectors. See Further *//*          Details. *//*  T       (output) COMPLEX*16 array, dimension (LDT,NB) *//*          The upper triangular matrix T. *//*  LDT     (input) INTEGER *//*          The leading dimension of the array T.  LDT >= NB. *//*  Y       (output) COMPLEX*16 array, dimension (LDY,NB) *//*          The n-by-nb matrix Y. *//*  LDY     (input) INTEGER *//*          The leading dimension of the array Y. LDY >= N. *//*  Further Details *//*  =============== *//*  The matrix Q is represented as a product of nb elementary reflectors *//*     Q = H(1) H(2) . . . H(nb). *///.........这里部分代码省略.........

//.........这里部分代码省略.........    m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):         (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 )         (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 )         (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 )         (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 )         (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )         (  v1  v2  v3  v4  v5 )       where d and e denote diagonal and off-diagonal elements of B, vi       denotes an element of the vector defining H(i), and ui an element of       the vector defining G(i).       =====================================================================          Test the input parameters          Parameter adjustments */    /* Table of constant values */    static doublecomplex c_b1 = {1.,0.};    static integer c__1 = 1;    static integer c_n1 = -1;    static integer c__3 = 3;    static integer c__2 = 2;        /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;    doublereal d__1;    doublecomplex z__1;    /* Local variables */    static integer i__, j, nbmin, iinfo, minmn;    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *), zgebd2_(integer *, integer *, 	    doublecomplex *, integer *, doublereal *, doublereal *, 	    doublecomplex *, doublecomplex *, doublecomplex *, integer *);    static integer nb, nx;    static doublereal ws;    extern /* Subroutine */ int xerbla_(char *, integer *), zlabrd_(	    integer *, integer *, integer *, doublecomplex *, integer *, 	    doublereal *, doublereal *, doublecomplex *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 	    integer *, integer *, ftnlen, ftnlen);    static integer ldwrkx, ldwrky, lwkopt;    static logical lquery;#define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1#define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]    a_dim1 = *lda;    a_offset = 1 + a_dim1 * 1;    a -= a_offset;    --d__;    --e;    --tauq;    --taup;    --work;    /* Function Body */    *info = 0;/* Computing MAX */    i__1 = 1, i__2 = ilaenv_(&c__1, "ZGEBRD", " ", m, n, &c_n1, &c_n1, (

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

voidzgemm(char transa, char transb, int m, int n, int k, doublecomplex *alpha, doublecomplex *a, int lda, doublecomplex *b, int ldb, doublecomplex *beta, doublecomplex *c, int ldc){    zgemm_(&transa, &transb, &m, &n, &k, alpha, a, &lda, b, &ldb, beta, c, &ldc);}

/* Subroutine */ int zhfrk_(char *transr, char *uplo, char *trans, integer *n, 	 integer *k, doublereal *alpha, doublecomplex *a, integer *lda, 	doublereal *beta, doublecomplex *c__){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2;    doublecomplex z__1;    /* Local variables */    integer j, n1, n2, nk, info;    doublecomplex cbeta;    logical normaltransr;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *), zherk_(char *, char *, integer *, 	    integer *, doublereal *, doublecomplex *, integer *, doublereal *, 	     doublecomplex *, integer *);    integer nrowa;    logical lower;    doublecomplex calpha;    extern /* Subroutine */ int xerbla_(char *, integer *);    logical nisodd, notrans;/*  -- LAPACK routine (version 3.2)                                    -- *//*  -- Contributed by Julien Langou of the Univ. of Colorado Denver    -- *//*  -- November 2008                                                   -- *//*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- *//*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- *//*     .. *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  Level 3 BLAS like routine for C in RFP Format. *//*  ZHFRK performs one of the Hermitian rank--k operations *//*     C := alpha*A*conjg( A' ) + beta*C, *//*  or *//*     C := alpha*conjg( A' )*A + beta*C, *//*  where alpha and beta are real scalars, C is an n--by--n Hermitian *//*  matrix and A is an n--by--k matrix in the first case and a k--by--n *//*  matrix in the second case. *//*  Arguments *//*  ========== *//*  TRANSR    (input) CHARACTER. *//*          = 'N':  The Normal Form of RFP A is stored; *//*          = 'C':  The Conjugate-transpose Form of RFP A is stored. *//*  UPLO   - (input) CHARACTER. *//*           On  entry,   UPLO  specifies  whether  the  upper  or  lower *//*           triangular  part  of the  array  C  is to be  referenced  as *//*           follows: *//*              UPLO = 'U' or 'u'   Only the  upper triangular part of  C *//*                                  is to be referenced. *//*              UPLO = 'L' or 'l'   Only the  lower triangular part of  C *//*                                  is to be referenced. *//*           Unchanged on exit. *//*  TRANS  - (input) CHARACTER. *//*           On entry,  TRANS  specifies the operation to be performed as *//*           follows: *//*              TRANS = 'N' or 'n'   C := alpha*A*conjg( A' ) + beta*C. *//*              TRANS = 'C' or 'c'   C := alpha*conjg( A' )*A + beta*C. *//*           Unchanged on exit. *//*  N      - (input) INTEGER. *//*           On entry,  N specifies the order of the matrix C.  N must be *//*           at least zero. *//*           Unchanged on exit. *//*  K      - (input) INTEGER. *//*           On entry with  TRANS = 'N' or 'n',  K  specifies  the number *//*           of  columns   of  the   matrix   A,   and  on   entry   with *//*           TRANS = 'C' or 'c',  K  specifies  the number of rows of the *//*           matrix A.  K must be at least zero. *//*           Unchanged on exit. *//*  ALPHA  - (input) DOUBLE PRECISION. *///.........这里部分代码省略.........

/* Subroutine */ int zget51_(integer *itype, integer *n, doublecomplex *a, 	integer *lda, doublecomplex *b, integer *ldb, doublecomplex *u, 	integer *ldu, doublecomplex *v, integer *ldv, doublecomplex *work, 	doublereal *rwork, doublereal *result){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, u_dim1, u_offset, v_dim1, 	    v_offset, i__1, i__2, i__3, i__4, i__5;    doublereal d__1, d__2;    doublecomplex z__1;    /* Local variables */    doublereal ulp;    integer jcol;    doublereal unfl;    integer jrow, jdiag;    doublereal anorm;    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *);    doublereal wnorm;    extern doublereal dlamch_(char *), zlange_(char *, integer *, 	    integer *, doublecomplex *, integer *, doublereal *);    extern /* Subroutine */ int 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 *//*  ======= *//*       ZGET51  generally checks a decomposition of the form *//*               A = U B V* *//*       where * means conjugate transpose and U and V are unitary. *//*       Specifically, if ITYPE=1 *//*               RESULT = | A - U B V* | / ( |A| n ulp ) *//*       If ITYPE=2, then: *//*               RESULT = | A - B | / ( |A| n ulp ) *//*       If ITYPE=3, then: *//*               RESULT = | I - UU* | / ( n ulp ) *//*  Arguments *//*  ========= *//*  ITYPE   (input) INTEGER *//*          Specifies the type of tests to be performed. *//*          =1: RESULT = | A - U B V* | / ( |A| n ulp ) *//*          =2: RESULT = | A - B | / ( |A| n ulp ) *//*          =3: RESULT = | I - UU* | / ( n ulp ) *//*  N       (input) INTEGER *//*          The size of the matrix.  If it is zero, ZGET51 does nothing. *//*          It must be at least zero. *//*  A       (input) COMPLEX*16 array, dimension (LDA, N) *//*          The original (unfactored) matrix. *//*  LDA     (input) INTEGER *//*          The leading dimension of A.  It must be at least 1 *//*          and at least N. *//*  B       (input) COMPLEX*16 array, dimension (LDB, N) *//*          The factored matrix. *//*  LDB     (input) INTEGER *//*          The leading dimension of B.  It must be at least 1 *//*          and at least N. *//*  U       (input) COMPLEX*16 array, dimension (LDU, N) *//*          The unitary matrix on the left-hand side in the *//*          decomposition. *//*          Not referenced if ITYPE=2 *//*  LDU     (input) INTEGER *//*          The leading dimension of U.  LDU must be at least N and *//*          at least 1. *//*  V       (input) COMPLEX*16 array, dimension (LDV, N) *//*          The unitary matrix on the left-hand side in the *//*          decomposition. *//*          Not referenced if ITYPE=2 *//*  LDV     (input) INTEGER *///.........这里部分代码省略.........

/* 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 zlaqps_(integer *m, integer *n, integer *offset, integer 	*nb, integer *kb, doublecomplex *a, integer *lda, integer *jpvt, 	doublecomplex *tau, doublereal *vn1, doublereal *vn2, doublecomplex *	auxv, doublecomplex *f, integer *ldf){    /* System generated locals */    integer a_dim1, a_offset, f_dim1, f_offset, i__1, i__2, i__3;    doublereal d__1, d__2;    doublecomplex z__1;    /* Builtin functions */    double sqrt(doublereal);    void d_cnjg(doublecomplex *, doublecomplex *);    double z_abs(doublecomplex *);    integer i_dnnt(doublereal *);    /* Local variables */    integer j, k, rk;    doublecomplex akk;    integer pvt;    doublereal temp, temp2, tol3z;    integer itemp;    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *), zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *), 	    zswap_(integer *, doublecomplex *, integer *, doublecomplex *, 	    integer *);    extern doublereal dznrm2_(integer *, doublecomplex *, integer *), dlamch_(	    char *);    extern integer idamax_(integer *, doublereal *, integer *);    integer lsticc;    extern /* Subroutine */ int zlarfp_(integer *, doublecomplex *, 	    doublecomplex *, integer *, doublecomplex *);    integer lastrk;/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZLAQPS computes a step of QR factorization with column pivoting *//*  of a complex M-by-N matrix A by using Blas-3.  It tries to factorize *//*  NB columns from A starting from the row OFFSET+1, and updates all *//*  of the matrix with Blas-3 xGEMM. *//*  In some cases, due to catastrophic cancellations, it cannot *//*  factorize NB columns.  Hence, the actual number of factorized *//*  columns is returned in KB. *//*  Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. *//*  Arguments *//*  ========= *//*  M       (input) INTEGER *//*          The number of rows of the matrix A. M >= 0. *//*  N       (input) INTEGER *//*          The number of columns of the matrix A. N >= 0 *//*  OFFSET  (input) INTEGER *//*          The number of rows of A that have been factorized in *//*          previous steps. *//*  NB      (input) INTEGER *//*          The number of columns to factorize. *//*  KB      (output) INTEGER *//*          The number of columns actually factorized. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,N) *//*          On entry, the M-by-N matrix A. *//*          On exit, block A(OFFSET+1:M,1:KB) is the triangular *//*          factor obtained and block A(1:OFFSET,1:N) has been *//*          accordingly pivoted, but no factorized. *//*          The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has *//*          been updated. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A. LDA >= max(1,M). *//*  JPVT    (input/output) INTEGER array, dimension (N) *//*          JPVT(I) = K <==> Column K of the full matrix A has been *//*          permuted into position I in AP. *//*  TAU     (output) COMPLEX*16 array, dimension (KB) *//*          The scalar factors of the elementary reflectors. *//*  VN1     (input/output) DOUBLE PRECISION array, dimension (N) *///.........这里部分代码省略.........

//.........这里部分代码省略.........	    nrow = nsupr - nsupc;	    solve_ops += 4 * nsupc * (nsupc - 1) * nrhs;	    solve_ops += 8 * nrow * nsupc * nrhs;	    	    if ( nsupc == 1 ) {		for (j = 0; j < nrhs; j++) {		    rhs_work = &Bmat[j*ldb];	    	    luptr = L_NZ_START(fsupc);		    for (iptr=istart+1; iptr < L_SUB_START(fsupc+1); iptr++){			irow = L_SUB(iptr);			++luptr;			zz_mult(&temp_comp, &rhs_work[fsupc], &Lval[luptr]);			z_sub(&rhs_work[irow], &rhs_work[irow], &temp_comp);		    }		}	    } else {	    	luptr = L_NZ_START(fsupc);#ifdef USE_VENDOR_BLAS#ifdef _CRAY		ftcs1 = _cptofcd("L", strlen("L"));		ftcs2 = _cptofcd("N", strlen("N"));		ftcs3 = _cptofcd("U", strlen("U"));		CTRSM( ftcs1, ftcs1, ftcs2, ftcs3, &nsupc, &nrhs, &alpha,		       &Lval[luptr], &nsupr, &Bmat[fsupc], &ldb);				CGEMM( ftcs2, ftcs2, &nrow, &nrhs, &nsupc, &alpha, 			&Lval[luptr+nsupc], &nsupr, &Bmat[fsupc], &ldb, 			&beta, &work[0], &n );#else		ztrsm_("L", "L", "N", "U", &nsupc, &nrhs, &alpha,		       &Lval[luptr], &nsupr, &Bmat[fsupc], &ldb);				zgemm_( "N", "N", &nrow, &nrhs, &nsupc, &alpha, 			&Lval[luptr+nsupc], &nsupr, &Bmat[fsupc], &ldb, 			&beta, &work[0], &n );#endif		for (j = 0; j < nrhs; j++) {		    rhs_work = &Bmat[j*ldb];		    work_col = &work[j*n];		    iptr = istart + nsupc;		    for (i = 0; i < nrow; i++) {			irow = L_SUB(iptr);			z_sub(&rhs_work[irow], &rhs_work[irow], &work_col[i]);			work_col[i].r = 0.0;	                work_col[i].i = 0.0;			iptr++;		    }		}#else				for (j = 0; j < nrhs; j++) {		    rhs_work = &Bmat[j*ldb];		    zlsolve (nsupr, nsupc, &Lval[luptr], &rhs_work[fsupc]);		    zmatvec (nsupr, nrow, nsupc, &Lval[luptr+nsupc],			    &rhs_work[fsupc], &work[0] );		    iptr = istart + nsupc;		    for (i = 0; i < nrow; i++) {			irow = L_SUB(iptr);			z_sub(&rhs_work[irow], &rhs_work[irow], &work[i]);			work[i].r = 0.;	                work[i].i = 0.;			iptr++;		    }		}#endif		    

/* Subroutine */ int zchkgk_(integer *nin, integer *nout){    /* Format strings */    static char fmt_9999[] = "(1x,/002.. test output of ZGGBAK .. /002)";    static char fmt_9998[] = "(/002 value of largest test error             "	    "     =/002,d12.3)";    static char fmt_9997[] = "(/002 example number where ZGGBAL info is not "	    "0    =/002,i4)";    static char fmt_9996[] = "(/002 example number where ZGGBAK(L) info is n"	    "ot 0 =/002,i4)";    static char fmt_9995[] = "(/002 example number where ZGGBAK(R) info is n"	    "ot 0 =/002,i4)";    static char fmt_9994[] = "(/002 example number having largest error     "	    "     =/002,i4)";    static char fmt_9992[] = "(/002 number of examples where info is not 0  "	    "     =/002,i4)";    static char fmt_9991[] = "(/002 total number of examples tested         "	    "     =/002,i4)";    /* System generated locals */    integer i__1, i__2, i__3, i__4;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1, z__2;    /* Builtin functions */    integer s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 	    e_rsle(void);    double d_imag(doublecomplex *);    integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen);    /* Local variables */    doublecomplex a[2500]	/* was [50][50] */, b[2500]	/* was [50][	    50] */, e[2500]	/* was [50][50] */, f[2500]	/* was [50][	    50] */;    integer i__, j, m, n;    doublecomplex af[2500]	/* was [50][50] */, bf[2500]	/* was [50][	    50] */, vl[2500]	/* was [50][50] */, vr[2500]	/* was [50][	    50] */;    integer ihi, ilo;    doublereal eps;    doublecomplex vlf[2500]	/* was [50][50] */;    integer knt;    doublecomplex vrf[2500]	/* was [50][50] */;    integer info, lmax[4];    doublereal rmax, vmax;    doublecomplex work[2500]	/* was [50][50] */;    integer ninfo;    doublereal anorm, bnorm;    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *);    doublereal rwork[300];    extern doublereal dlamch_(char *);    doublereal lscale[50];    extern /* Subroutine */ int zggbak_(char *, char *, integer *, integer *, 	    integer *, doublereal *, doublereal *, integer *, doublecomplex *, 	     integer *, integer *), zggbal_(char *, integer *, 	     doublecomplex *, integer *, doublecomplex *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, integer *);    doublereal rscale[50];    extern doublereal zlange_(char *, integer *, integer *, doublecomplex *, 	    integer *, doublereal *);    extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 	    doublecomplex *, integer *, doublecomplex *, integer *);    /* Fortran I/O blocks */    static cilist io___6 = { 0, 0, 0, 0, 0 };    static cilist io___10 = { 0, 0, 0, 0, 0 };    static cilist io___13 = { 0, 0, 0, 0, 0 };    static cilist io___15 = { 0, 0, 0, 0, 0 };    static cilist io___17 = { 0, 0, 0, 0, 0 };    static cilist io___35 = { 0, 0, 0, fmt_9999, 0 };    static cilist io___36 = { 0, 0, 0, fmt_9998, 0 };    static cilist io___37 = { 0, 0, 0, fmt_9997, 0 };    static cilist io___38 = { 0, 0, 0, fmt_9996, 0 };    static cilist io___39 = { 0, 0, 0, fmt_9995, 0 };    static cilist io___40 = { 0, 0, 0, fmt_9994, 0 };    static cilist io___41 = { 0, 0, 0, fmt_9992, 0 };    static cilist io___42 = { 0, 0, 0, fmt_9991, 0 };/*  -- LAPACK test routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZCHKGK tests ZGGBAK, a routine for backward balancing  of *//*  a matrix pair (A, B). *//*  Arguments *//*  ========= *//*  NIN     (input) INTEGER *///.........这里部分代码省略.........

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

//.........这里部分代码省略.........    On entry:                       On exit:          a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66          a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *          a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *       Array elements marked * are not used by the routine.       Contributed by       Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989       =====================================================================          Test the input parameters.          Parameter adjustments */    /* Table of constant values */    static doublecomplex c_b1 = {1.,0.};    static integer c__1 = 1;    static integer c_n1 = -1;    static doublereal c_b21 = -1.;    static doublereal c_b22 = 1.;    static integer c__33 = 33;        /* System generated locals */    integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4, i__5, i__6;    doublecomplex z__1;    /* Local variables */    static doublecomplex work[1056]	/* was [33][32] */;    static integer i__, j;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *), zherk_(char *, char *, integer *, 	    integer *, doublereal *, doublecomplex *, integer *, doublereal *,	     doublecomplex *, integer *);    static integer i2, i3;    extern /* Subroutine */ int ztrsm_(char *, char *, char *, char *, 	    integer *, integer *, doublecomplex *, doublecomplex *, integer *,	     doublecomplex *, integer *), 	    zpbtf2_(char *, integer *, integer *, doublecomplex *, integer *, 	    integer *);    static integer ib, nb, ii, jj;    extern /* Subroutine */ int zpotf2_(char *, integer *, doublecomplex *, 	    integer *, integer *), xerbla_(char *, integer *);    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 	    integer *, integer *, ftnlen, ftnlen);#define work_subscr(a_1,a_2) (a_2)*33 + a_1 - 34#define work_ref(a_1,a_2) work[work_subscr(a_1,a_2)]#define ab_subscr(a_1,a_2) (a_2)*ab_dim1 + a_1#define ab_ref(a_1,a_2) ab[ab_subscr(a_1,a_2)]    ab_dim1 = *ldab;    ab_offset = 1 + ab_dim1 * 1;    ab -= ab_offset;    /* Function Body */    *info = 0;    if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {	*info = -1;    } else if (*n < 0) {	*info = -2;

 int zgetrf_(int *m, int *n, doublecomplex *a, 	int *lda, int *ipiv, int *info){    /* System generated locals */    int a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;    doublecomplex z__1;    /* Local variables */    int i__, j, jb, nb, iinfo;    extern  int zgemm_(char *, char *, int *, int *, 	    int *, doublecomplex *, doublecomplex *, int *, 	    doublecomplex *, int *, doublecomplex *, doublecomplex *, 	    int *), ztrsm_(char *, char *, char *, char *, 	     int *, int *, doublecomplex *, doublecomplex *, int *, doublecomplex *, int *), 	    zgetf2_(int *, int *, doublecomplex *, int *, int 	    *, int *), xerbla_(char *, int *);    extern int ilaenv_(int *, char *, char *, int *, int *, 	    int *, int *);    extern  int zlaswp_(int *, doublecomplex *, int *, 	     int *, int *, int *, int *);/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZGETRF computes an LU factorization of a general M-by-N matrix A *//*  using partial pivoting with row interchanges. *//*  The factorization has the form *//*     A = P * L * U *//*  where P is a permutation matrix, L is lower triangular with unit *//*  diagonal elements (lower trapezoidal if m > n), and U is upper *//*  triangular (upper trapezoidal if m < n). *//*  This is the right-looking Level 3 BLAS version of the algorithm. *//*  Arguments *//*  ========= *//*  M       (input) INTEGER *//*          The number of rows of the matrix A.  M >= 0. *//*  N       (input) INTEGER *//*          The number of columns of the matrix A.  N >= 0. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,N) *//*          On entry, the M-by-N matrix to be factored. *//*          On exit, the factors L and U from the factorization *//*          A = P*L*U; the unit diagonal elements of L are not stored. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A.  LDA >= MAX(1,M). *//*  IPIV    (output) INTEGER array, dimension (MIN(M,N)) *//*          The pivot indices; for 1 <= i <= MIN(M,N), row i of the *//*          matrix was interchanged with row IPIV(i). *//*  INFO    (output) INTEGER *//*          = 0:  successful exit *//*          < 0:  if INFO = -i, the i-th argument had an illegal value *//*          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization *//*                has been completed, but the factor U is exactly *//*                singular, and division by zero will occur if it is used *//*                to solve a system of equations. *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. Executable Statements .. *//*     Test the input parameters. */    /* Parameter adjustments */    a_dim1 = *lda;    a_offset = 1 + a_dim1;    a -= a_offset;    --ipiv;    /* Function Body */    *info = 0;    if (*m < 0) {//.........这里部分代码省略.........

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

//.........这里部分代码省略.........            interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.       W       (workspace) COMPLEX*16 array, dimension (LDW,NB)       LDW     (input) INTEGER               The leading dimension of the array W.  LDW >= max(1,N).       INFO    (output) INTEGER               = 0: successful exit               > 0: if INFO = k, D(k,k) is exactly zero.  The factorization                    has been completed, but the block diagonal matrix D is                    exactly singular.       =====================================================================          Parameter adjustments */    /* Table of constant values */    static doublecomplex c_b1 = {1.,0.};    static integer c__1 = 1;        /* System generated locals */    integer a_dim1, a_offset, w_dim1, w_offset, i__1, i__2, i__3, i__4, i__5;    doublereal d__1, d__2, d__3, d__4;    doublecomplex z__1, z__2, z__3, z__4;    /* Builtin functions */    double sqrt(doublereal), d_imag(doublecomplex *);    void d_cnjg(doublecomplex *, doublecomplex *), z_div(doublecomplex *, 	    doublecomplex *, doublecomplex *);    /* Local variables */    static integer imax, jmax, j, k;    static doublereal t, alpha;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 	    integer *, doublecomplex *, doublecomplex *, integer *, 	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 	    integer *);    static integer kstep;    extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 	    doublecomplex *, doublecomplex *, integer *, doublecomplex *, 	    integer *, doublecomplex *, doublecomplex *, integer *);    static doublereal r1;    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 	    doublecomplex *, integer *), zswap_(integer *, doublecomplex *, 	    integer *, doublecomplex *, integer *);    static doublecomplex d11, d21, d22;    static integer jb, jj, kk, jp, kp;    static doublereal absakk;    static integer kw;    extern /* Subroutine */ int zdscal_(integer *, doublereal *, 	    doublecomplex *, integer *);    static doublereal colmax;    extern /* Subroutine */ int zlacgv_(integer *, doublecomplex *, integer *)	    ;    extern integer izamax_(integer *, doublecomplex *, integer *);    static doublereal rowmax;    static integer kkw;#define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1#define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]#define w_subscr(a_1,a_2) (a_2)*w_dim1 + a_1#define w_ref(a_1,a_2) w[w_subscr(a_1,a_2)]    a_dim1 = *lda;    a_offset = 1 + a_dim1 * 1;    a -= a_offset;

 int zpotrf_(char *uplo, int *n, doublecomplex *a, 	int *lda, int *info){    /* System generated locals */    int a_dim1, a_offset, i__1, i__2, i__3, i__4;    doublecomplex z__1;    /* Local variables */    int j, jb, nb;    extern int lsame_(char *, char *);    extern  int zgemm_(char *, char *, int *, int *, 	    int *, doublecomplex *, doublecomplex *, int *, 	    doublecomplex *, int *, doublecomplex *, doublecomplex *, 	    int *), zherk_(char *, char *, int *, 	    int *, double *, doublecomplex *, int *, double *, 	     doublecomplex *, int *);    int upper;    extern  int ztrsm_(char *, char *, char *, char *, 	    int *, int *, doublecomplex *, doublecomplex *, int *, 	     doublecomplex *, int *), 	    zpotf2_(char *, int *, doublecomplex *, int *, int *), xerbla_(char *, int *);    extern int ilaenv_(int *, char *, char *, int *, int *, 	    int *, int *);/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  ZPOTRF computes the Cholesky factorization of a complex Hermitian *//*  positive definite matrix A. *//*  The factorization has the form *//*     A = U**H * U,  if UPLO = 'U', or *//*     A = L  * L**H,  if UPLO = 'L', *//*  where U is an upper triangular matrix and L is lower triangular. *//*  This is the block version of the algorithm, calling Level 3 BLAS. *//*  Arguments *//*  ========= *//*  UPLO    (input) CHARACTER*1 *//*          = 'U':  Upper triangle of A is stored; *//*          = 'L':  Lower triangle of A is stored. *//*  N       (input) INTEGER *//*          The order of the matrix A.  N >= 0. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,N) *//*          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading *//*          N-by-N upper triangular part of A contains the upper *//*          triangular part of the matrix A, and the strictly lower *//*          triangular part of A is not referenced.  If UPLO = 'L', the *//*          leading N-by-N lower triangular part of A contains the lower *//*          triangular part of the matrix A, and the strictly upper *//*          triangular part of A is not referenced. *//*          On exit, if INFO = 0, the factor U or L from the Cholesky *//*          factorization A = U**H*U or A = L*L**H. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A.  LDA >= MAX(1,N). *//*  INFO    (output) INTEGER *//*          = 0:  successful exit *//*          < 0:  if INFO = -i, the i-th argument had an illegal value *//*          > 0:  if INFO = i, the leading minor of order i is not *//*                positive definite, and the factorization could not be *//*                completed. *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. Executable Statements .. *//*     Test the input parameters. */    /* Parameter adjustments */    a_dim1 = *lda;    a_offset = 1 + a_dim1;    a -= a_offset;//.........这里部分代码省略.........

/*<       SUBROUTINE ZGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO ) >*//* Subroutine */ int zgehrd_(integer *n, integer *ilo, integer *ihi,        doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *        work, integer *lwork, integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3, i__4;    doublecomplex z__1;    /* Local variables */    integer i__;    doublecomplex t[4160]       /* was [65][64] */;    integer ib;    doublecomplex ei;    integer nb, nh, nx=0, iws, nbmin, iinfo;    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *,            integer *, doublecomplex *, doublecomplex *, integer *,            doublecomplex *, integer *, doublecomplex *, doublecomplex *,            integer *, ftnlen, ftnlen), zgehd2_(integer *, integer *, integer            *, doublecomplex *, integer *, doublecomplex *, doublecomplex *,            integer *), xerbla_(char *, integer *, ftnlen);    extern integer ilaenv_(integer *, char *, char *, integer *, integer *,            integer *, integer *, ftnlen, ftnlen);    extern /* Subroutine */ int zlarfb_(char *, char *, char *, char *,            integer *, integer *, integer *, doublecomplex *, integer *,            doublecomplex *, integer *, doublecomplex *, integer *,            doublecomplex *, integer *, ftnlen, ftnlen, ftnlen, ftnlen),            zlahrd_(integer *, integer *, integer *, doublecomplex *, integer            *, doublecomplex *, doublecomplex *, integer *, doublecomplex *,            integer *);    integer ldwork, lwkopt;    logical lquery;/*  -- LAPACK routine (version 3.0) -- *//*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., *//*     Courant Institute, Argonne National Lab, and Rice University *//*     June 30, 1999 *//*     .. Scalar Arguments .. *//*<       INTEGER            IHI, ILO, INFO, LDA, LWORK, N >*//*     .. *//*     .. Array Arguments .. *//*<       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * ) >*//*     .. *//*  Purpose *//*  ======= *//*  ZGEHRD 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 <= N, if N > 0; ILO=1 and IHI=0, if N=0. *//*  A       (input/output) COMPLEX*16 array, dimension (LDA,N) *//*          On entry, the N-by-N general matrix to be reduced. *//*          On exit, the upper triangle and the first subdiagonal of A *//*          are overwritten with the upper Hessenberg matrix H, and the *//*          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). Elements 1:ILO-1 and IHI:N-1 of TAU are set to *//*          zero. *//*  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK) *//*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. *//*  LWORK   (input) INTEGER *//*          The length of the array WORK.  LWORK >= max(1,N). *//*          For optimum performance LWORK >= N*NB, where NB is the *//*          optimal blocksize. *//*          If LWORK = -1, then a workspace query is assumed; the routine *//*          only calculates the optimal size of the WORK array, returns *//*          this value as the first entry of the WORK array, and no error *//*          message related to LWORK is issued by XERBLA. *//*  INFO    (output) INTEGER *//*          = 0:  successful exit *//*          < 0:  if INFO = -i, the i-th argument had an illegal value. *//*  Further Details *///.........这里部分代码省略.........