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;
}

int
f2c_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;
}

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

void
zgemm(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);
}