void calcstep(int m, int n, double *A, double *B, double *s, double *y,
              double *r1, double *r2, double r3, double *r4, double *dx,
              double *ds, double *dt, double *dy) {
  char Transpose = 'T';
  char Normal = 'N';
  int n1 = n + 1;
  int oneI = 1;
  double none = -1.0;
  double one = 1.0;
  int info;
  int i;

  int *myworkI;
  double *dxdt;
  double *tmp;
  double *tmpB;

  tmp = pswarm_malloc(m * sizeof(double));
  dxdt = pswarm_malloc(n1 * sizeof(double));

  memset(dxdt, 0, n1 * sizeof(double));

  dxdt[n] = 0.0;
  for (i = 0; i < m; i++) {
    tmp[i] = (r1[i] * y[i] - r4[i]) / s[i];
    dxdt[n] += tmp[i];
  }

  memcpy(dxdt, r2, n * sizeof(double));
  dgemv_(&Transpose, &m, &n, &one, A, &m, tmp, &oneI, &one, dxdt, &oneI);

  /*  dpotrs_(&Upper, &n1, &oneI, B, &n1, dxdt, &n1, &info); */

  free(tmp);

  tmpB = pswarm_malloc(n1 * n1 * sizeof(double));
  myworkI = pswarm_malloc(n1 * sizeof(int));

  memcpy(tmpB, B, n1 * n1 * sizeof(double));

  dgesv_(&n1, &oneI, tmpB, &n1, myworkI, dxdt, &n1, &info);

  memcpy(dx, dxdt, n * sizeof(double));
  *dt = dxdt[n];

  memcpy(ds, r1, m * sizeof(double));
  dgemv_(&Normal, &m, &n, &none, A, &m, dx, &oneI, &one, ds, &oneI);

  for (i = 0; i < m; i++) {
    ds[i] -= (*dt);
    dy[i] = (r4[i] - y[i] * ds[i]) / s[i];
  }

  free(myworkI);
  free(dxdt);
  free(tmpB);
}

void ProtoMol::Lapack::dgemv(char *transA, int *m, int *n, double *alpha,
                             double *A, int *lda, double *x, int *incx,
                             double *beta, double *Y, int *incY) {
  FAHCheckIn();
#if defined(HAVE_LAPACK)
  dgemv_(transA, m, n, alpha, A, lda, x, incx, beta, Y, incY);
#elif defined(HAVE_SIMTK_LAPACK)
  dgemv_(*transA, *m, *n, *alpha, A, *lda, x, *incx, *beta, Y, *incY, 1);
#elif defined(HAVE_MKL_LAPACK)
  DGEMV(transA, m, n, alpha, A, lda, x, incx, beta, Y, incY);
#else
  THROW(std::string(__func__) + " not supported");
#endif
}

/* Subroutine */ int bicgkernel_(integer *lda, integer *n, doublereal *a, 
	doublereal *p, doublereal *r__, real *s, real *q)
{
    /* System generated locals */
    integer a_dim1, a_offset;

    /* Local variables */
    static doublereal one;
    static integer incx, incy;
    static doublereal zero;
    extern /* Subroutine */ int dgemv_(char *, integer *, integer *, 
	    doublereal *, doublereal *, integer *, doublereal *, integer *, 
	    doublereal *, real *, integer *, ftnlen);


/*   BICG */
/*   in */
/*     A : column matrix, p : vector, r : vector */
/*   out */
/*     q : vector, s : vector */
/*   { */
/*     q = A * p */
/*     s = A' * r */
/*   } */

    /* Parameter adjustments */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    --p;
    --r__;

    /* Function Body */
    incx = 1;
    incy = 1;
    one = 1.;
    zero = 0.;

/*    Put A*p in q */

    dgemv_("n", n, n, &one, &a[a_offset], lda, &p[1], &incx, &zero, q, &incy, 
	    (ftnlen)1);

/*    Put A'*r in s */

    dgemv_("t", n, n, &one, &a[a_offset], lda, &r__[1], &incx, &zero, s, &
	    incy, (ftnlen)1);

    return 0;
} /* bicgkernel_ */

void FCItdm12kern_b(double *tdm1, double *tdm2, double *bra, double *ket,
                    int bcount, int stra_id, int strb_id,
                    int norb, int na, int nb, int nlinka, int nlinkb,
                    _LinkT *clink_indexa, _LinkT *clink_indexb, int symm)
{
        const int INC1 = 1;
        const char TRANS_N = 'N';
        const char TRANS_T = 'T';
        const double D1 = 1;
        const int nnorb = norb * norb;
        double csum;
        double *buf0 = calloc(nnorb*bcount, sizeof(double));
        double *buf1 = calloc(nnorb*bcount, sizeof(double));

        csum = FCIrdm2_b_t1ci(bra, buf1, bcount, stra_id, strb_id,
                              norb, nb, nlinkb, clink_indexb);
        if (csum < CSUMTHR) { goto _normal_end; }
        csum = FCIrdm2_b_t1ci(ket, buf0, bcount, stra_id, strb_id,
                              norb, nb, nlinkb, clink_indexb);
        if (csum < CSUMTHR) { goto _normal_end; }
        dgemv_(&TRANS_N, &nnorb, &bcount, &D1, buf0, &nnorb,
               bra+stra_id*nb+strb_id, &INC1, &D1, tdm1, &INC1);
        switch (symm) {
        case PARTICLESYM:
                tril_particle_symm(tdm2, buf1, buf0, bcount, norb, D1, D1);
                break;
        default:
                dgemm_(&TRANS_N, &TRANS_T, &nnorb, &nnorb, &bcount,
                       &D1, buf0, &nnorb, buf1, &nnorb, &D1, tdm2, &nnorb);
        }
_normal_end:
        free(buf0);
        free(buf1);
}

static void make_rdm12_sf(double *rdm1, double *rdm2,
                          double *bra, double *ket, double *t1bra, double *t1ket,
                          int bcount, int stra_id, int strb_id,
                          int norb, int na, int nb)
{
        const char TRANS_N = 'N';
        const char TRANS_T = 'T';
        const int INC1 = 1;
        const double D1 = 1;
        const int nnorb = norb * norb;
        int k, l;
        size_t n;
        double *tbra = malloc(sizeof(double) * nnorb * bcount);
        double *pbra, *pt1;

        for (n = 0; n < bcount; n++) {
                pbra = tbra + n * nnorb;
                pt1 = t1bra + n * nnorb;
                for (k = 0; k < norb; k++) {
                        for (l = 0; l < norb; l++) {
                                pbra[k*norb+l] = pt1[l*norb+k];
                        }
                }
        }
        dgemm_(&TRANS_N, &TRANS_T, &nnorb, &nnorb, &bcount,
               &D1, t1ket, &nnorb, tbra, &nnorb,
               &D1, rdm2, &nnorb);

        dgemv_(&TRANS_N, &nnorb, &bcount, &D1, t1ket, &nnorb,
               bra+stra_id*nb+strb_id, &INC1, &D1, rdm1, &INC1);

        free(tbra);
}

/*! dgematrix*_dcovector operator */
inline _dcovector operator*(const dgematrix& mat, const _dcovector& vec)
{
#ifdef  CPPL_VERBOSE
  std::cerr << "# [MARK] operator*(const dgematrix&, const _dcovector&)"
            << std::endl;
#endif//CPPL_VERBOSE

#ifdef  CPPL_DEBUG
  if(mat.N!=vec.L){
    std::cerr << "[ERROR] operator*(const dgematrix&, const _dcovector&)"
              << std::endl
              << "These matrix and vector can not make a product." << std::endl
              << "Your input was (" << mat.M << "x" << mat.N << ") * ("
              << vec.L << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dcovector newvec(mat.M);
  dgemv_( 'N', mat.M, mat.N, 1.0, mat.Array, mat.M,
          vec.Array, 1, 0.0, newvec.array, 1 );
  
  vec.destroy();
  return _(newvec);
}

/*
 * 2pdm kernel for  beta^i beta_j | ci0 >
 */
void FCIrdm12kern_b(double *rdm1, double *rdm2, double *bra, double *ket,
                    int bcount, int stra_id, int strb_id,
                    int norb, int na, int nb, int nlinka, int nlinkb,
                    _LinkT *clink_indexa, _LinkT *clink_indexb, int symm)
{
        const int INC1 = 1;
        const char UP = 'U';
        const char TRANS_N = 'N';
        const char TRANS_T = 'T';
        const double D1 = 1;
        const int nnorb = norb * norb;
        double csum;
        double *buf = calloc(nnorb*bcount, sizeof(double));

        csum = FCIrdm2_b_t1ci(ket, buf, bcount, stra_id, strb_id,
                              norb, nb, nlinkb, clink_indexb);
        if (csum > CSUMTHR) {
                dgemv_(&TRANS_N, &nnorb, &bcount, &D1, buf, &nnorb,
                       ket+stra_id*nb+strb_id, &INC1, &D1, rdm1, &INC1);
                switch (symm) {
                case BRAKETSYM:
                        dsyrk_(&UP, &TRANS_N, &nnorb, &bcount,
                               &D1, buf, &nnorb, &D1, rdm2, &nnorb);
                        break;
                case PARTICLESYM:
                        tril_particle_symm(rdm2, buf, buf, bcount, norb, 1, 1);
                        break;
                default:
                        dgemm_(&TRANS_N, &TRANS_T, &nnorb, &nnorb, &bcount,
                               &D1, buf, &nnorb, buf, &nnorb,
                               &D1, rdm2, &nnorb);
                }
        }
        free(buf);
}

/** Calculates w = G * s.
 */
void calc_w(int m, int p, double** G, double* s, double* w)
{
    double a = 1.0;
    int inc = 1;
    double b = 0.0;

    dgemv_(&noT, &m, &p, &a, G[0], &m, s, &inc, &b, w, &inc);
}

void THBlas_(gemv)(char trans, int64_t m, int64_t n, real alpha, real *a, int64_t lda, real *x, int64_t incx, real beta, real *y, int64_t incy)
{
  if(n == 1)
    lda = m;

#if defined(USE_BLAS) && (defined(TH_REAL_IS_DOUBLE) || defined(TH_REAL_IS_FLOAT))
  if( (m <= INT_MAX) && (n <= INT_MAX) && (lda <= INT_MAX) &&
      (incx > 0) && (incx <= INT_MAX) &&
      (incy > 0) && (incy <= INT_MAX) )
  {
    THArgCheck(lda >= THMax(1, m), 6,
      "lda should be at least max(1, m=%d), but have %d", m, lda);
    int i_m = (int)m;
    int i_n = (int)n;
    int i_lda = (int)lda;
    int i_incx = (int)incx;
    int i_incy = (int)incy;

#if defined(TH_REAL_IS_DOUBLE)
    dgemv_(&trans, &i_m, &i_n, &alpha, a, &i_lda, x, &i_incx, &beta, y, &i_incy);
#else
    sgemv_(&trans, &i_m, &i_n, &alpha, a, &i_lda, x, &i_incx, &beta, y, &i_incy);
#endif
    return;
  }
#endif
  {
    int64_t i, j;

    if( (trans == 'T') || (trans == 't') )
    {
      for(i = 0; i < n; i++)
      {
        real sum = 0;
        real *row_ = a+lda*i;
        for(j = 0; j < m; j++)
          sum += x[j*incx]*row_[j];
	if (beta == 0)
	  y[i*incy] = alpha*sum;
	else
	  y[i*incy] = beta*y[i*incy] + alpha*sum;
      }
    }
    else
    {
      if(beta != 1)
        THBlas_(scal)(m, beta, y, incy);

      for(j = 0; j < n; j++)
      {
        real *column_ = a+lda*j;
        real z = alpha*x[j*incx];
        for(i = 0; i < m; i++)
          y[i*incy] += z*column_[i];
      }
    }
  }
}

 /* ---- PBC functions ---- */
void DoCorrectionTable(double *ChebyshevWeightSource, double *ChebyshevWeightField, int n, int2 dof, double Len, double alpha, int lpbc, kernel_t kernel, double *Tkz) {

  kfun_t kfun = kernel.kfun;
  double homogen = kernel.homogen;
  
  char kpbcFilename[50];
  sprintf(kpbcFilename, "Kn%dpbc%da%.1f.out", n, lpbc, alpha);
  FILE *kfile = fopen(kpbcFilename, "r");

  // Create correction table
  if (kfile == NULL) {
    printf("kpbc file: %s does not exist. Creating now ...\n", kpbcFilename);
    CreateTableCorrection(kfun, n, dof, alpha, lpbc, Tkz, kpbcFilename);
    kfile = fopen(kpbcFilename, "r");
    assert(kfile != NULL);
  }
  else
    printf("kpbc file exits. Reading now ...\n");
  
  // Read and scale for elastic constants
  int i, j=0;
  double c3Read[3];
  for (i=0; i<3; i++)
    j += fscanf(kfile, "%lf", c3Read+i);
  assert (j == 3);

  int n3 = n*n*n;
  int n3f = n3 * dof.f, n3s = n3 * dof.s;
  int dof2n6 = n3f * n3s;
  double *KPBC = (double *) malloc( dof2n6 * sizeof(double) );
  for (i=0; i<dof2n6; i++) 
    j += fscanf(kfile, "%lf", KPBC+i);
  fclose(kfile);
  assert( j-3 == dof2n6 );
  
  // Compute stress from PBC
  int incr = 1;
  double beta = 0;
  char trans = 'n';
  double scale = pow(1/Len, homogen); // Scale for the box size

  // check if parameters are exactly scaled
  if ( fabs(AnisoParameters->c3[0]/c3Read[0] -
	    AnisoParameters->c3[1]/c3Read[1]) <1e-6 &&   
       fabs(AnisoParameters->c3[2]/c3Read[2] -
	    AnisoParameters->c3[1]/c3Read[1]) <1e-6 )
      
    scale *= AnisoParameters->c3[0]/c3Read[0];
  else
    printf("Error: elastic constants do not match (scale).\n");

  dgemv_(&trans, &n3f, &n3s, &scale, KPBC,
	 &n3f, ChebyshevWeightSource, &incr, &beta,
	 ChebyshevWeightField, &incr);
    
  free(KPBC), KPBC=NULL;

}

void DenseGenMatrix::mult ( double beta,  double y[], int incy,
				double alpha, double x[], int incx )
{
  char fortranTrans = 'T';
  int n = mStorage->n, m = mStorage->m;
  
  dgemv_( &fortranTrans, &n, &m, &alpha, &mStorage->M[0][0], &n,
	  x, &incx, &beta, y, &incy );
}

/* y := alpha*A*x + beta*y.
 * INPUT
 *  m : the number of rows   of the matrix A, that is, y[m]
 *  n : the number of colums of the matrix A, that is, x[n]
 *      A is 'm by n' matrix.
 *  dgemv_() is implemented in FORTRAN, so that,
 *  for 'N' case,
 *     y[i] = sum_j A[i,j] x[j]
 *          = sum_J a[I+m*J] * x[J], where I:=i-1, etc.
 *          = sum_J a[J*m+I] * x[J]
 *          = A_C[J,I] * x[J]
 *  for 'T' case,
 *     y[i] = sum_j A[j,i] x[j]
 *          = sum_J a[J+n*I] * x[J], where I:=i-1, etc.
 *          = sum_J a[I*n+J] * x[J]
 *          = A_C[I,J] * x[J]
 *  NOTE, in this case, m and n also should be exchanged...
 */
void dgemv_wrap (int m, int n, double alpha, double *a,
		 double *x, double beta,
		 double *y)
{
  char trans = 'T'; /* fortran's memory allocation is transposed */
  int one = 1;

  dgemv_ (&trans, &n, &m, &alpha, a, &n, x, &one, &beta, y, &one);
}

/*
 * _spin0 assumes the strict symmetry on alpha and beta electrons
 */
void FCIrdm12kern_spin0(double *rdm1, double *rdm2, double *bra, double *ket,
                        int bcount, int stra_id, int strb_id,
                        int norb, int na, int nb, int nlinka, int nlinkb,
                        _LinkT *clink_indexa, _LinkT *clink_indexb, int symm)
{
        if (stra_id < strb_id) {
                return;
        }
        const int INC1 = 1;
        const char UP = 'U';
        const char TRANS_N = 'N';
        const char TRANS_T = 'T';
        const double D1 = 1;
        const double D2 = 2;
        const int nnorb = norb * norb;
        int fill0, fill1, i;
        double csum = 0;
        double *buf = calloc(nnorb * na, sizeof(double));

        if (strb_id+bcount <= stra_id) {
                fill0 = bcount;
                fill1 = bcount;
                csum = FCIrdm2_b_t1ci(ket, buf, fill0, stra_id, strb_id,
                                      norb, na, nlinka, clink_indexa)
                     + FCIrdm2_a_t1ci(ket, buf, fill1, stra_id, strb_id,
                                      norb, na, nlinka, clink_indexa);
        } else if (stra_id >= strb_id) {
                fill0 = stra_id - strb_id;
                fill1 = stra_id - strb_id + 1;
                csum = FCIrdm2_b_t1ci(ket, buf, fill0, stra_id, strb_id,
                                      norb, na, nlinka, clink_indexa)
                     + FCIrdm2_a_t1ci(ket, buf, fill1, stra_id, strb_id,
                                      norb, na, nlinka, clink_indexa);
        }
        if (csum > CSUMTHR) {
                dgemv_(&TRANS_N, &nnorb, &fill1, &D2, buf, &nnorb,
                       ket+stra_id*na+strb_id, &INC1, &D1, rdm1, &INC1);

                for (i = fill0*nnorb; i < fill1*nnorb; i++) {
                        buf[i] *= SQRT2;
                }
                switch (symm) {
                case BRAKETSYM:
                        dsyrk_(&UP, &TRANS_N, &nnorb, &fill1,
                               &D2, buf, &nnorb, &D1, rdm2, &nnorb);
                        break;
                case PARTICLESYM:
                        tril_particle_symm(rdm2, buf, buf, fill1, norb, D2, D1);
                        break;
                default:
                        dgemm_(&TRANS_N, &TRANS_T, &nnorb, &nnorb, &fill1,
                               &D2, buf, &nnorb, buf, &nnorb,
                               &D1, rdm2, &nnorb);
                }
        }
        free(buf);
}

INLINE void do_dgemv(double *a, double *x, double *y, int iterations, int *limit, int *lda)
{
  REGISTER int i = 0;
  extern int dgemv_();

  for (;i<iterations;i++)
    {
      dgemv_(foo,limit,limit,&dalpha,a,lda,x,&stride,&dbeta,y,&stride);
    }
}

void EucQuadratic::HessianEta(Variable *x, Vector *etax, Vector *xix) const
{
	const double *v = etax->ObtainReadData();
	double *xixTV = xix->ObtainWriteEntireData();

	char *transn = const_cast<char *> ("n");
	integer N = Dim, inc = 1;
	double two = 2, zero = 0;
	dgemv_(transn, &N, &N, &two, A, &N, const_cast<double *> (v), &inc, &zero, xixTV, &inc);
};

PyObject* gemv(PyObject *self, PyObject *args)
{
  Py_complex alpha;
  PyArrayObject* a;
  PyArrayObject* x;
  Py_complex beta;
  PyArrayObject* y;
  char trans = 't';
  if (!PyArg_ParseTuple(args, "DOODO|c", &alpha, &a, &x, &beta, &y, &trans))
    return NULL;

  int m, n, lda, itemsize, incx, incy;

  if (trans == 'n')
    {
      m = PyArray_DIMS(a)[1];
      for (int i = 2; i < PyArray_NDIM(a); i++)
        m *= PyArray_DIMS(a)[i];
      n = PyArray_DIMS(a)[0];
      lda = MAX(1, m);
    }
  else
    {
      n = PyArray_DIMS(a)[0];
      for (int i = 1; i < PyArray_NDIM(a)-1; i++)
        n *= PyArray_DIMS(a)[i];
      m = PyArray_DIMS(a)[PyArray_NDIM(a)-1];
      lda = MAX(1, m);
    }

  if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)
    itemsize = sizeof(double);
  else
    itemsize = sizeof(double_complex);

  incx = PyArray_STRIDES(x)[0]/itemsize;
  incy = 1;

  if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)
    dgemv_(&trans, &m, &n,
           &(alpha.real),
           DOUBLEP(a), &lda,
           DOUBLEP(x), &incx,
           &(beta.real),
           DOUBLEP(y), &incy);
  else
    zgemv_(&trans, &m, &n,
           &alpha,
           (void*)COMPLEXP(a), &lda,
           (void*)COMPLEXP(x), &incx,
           &beta,
           (void*)COMPLEXP(y), &incy);
  Py_RETURN_NONE;
}

	void EucQuadratic::HessianEta(Variable *x, Vector *etax, Vector *xix) const
	{
		const double *v = etax->ObtainReadData();
		double *xixTV = xix->ObtainWriteEntireData();

		char *transn = const_cast<char *> ("n");
		integer N = Dim, inc = 1;
		double two = 2, zero = 0;
		// xixTV <- 2 * A * v, details: http://www.netlib.org/lapack/explore-html/dc/da8/dgemv_8f.html
		dgemv_(transn, &N, &N, &two, A, &N, const_cast<double *> (v), &inc, &zero, xixTV, &inc);
	};

GURLS_EXPORT void gemv(const CBLAS_TRANSPOSE TransA,
          const int M, const int N, const double alpha, const double *A,
          const int lda, const double *X, const int incX,
          const double beta, double *Y, const int incY)
{
    char transA = BlasUtils::charValue(TransA);

    dgemv_(&transA, const_cast<int*>(&M), const_cast<int*>(&N),
          const_cast<double*>(&alpha), const_cast<double*>(A), const_cast<int*>(&lda),
          const_cast<double*>(X), const_cast<int*>(&incX), const_cast<double*>(&beta),
          const_cast<double*>(Y), const_cast<int*>(&incY));
}

int
f2c_dgemv(char* trans, integer* M, integer* N,
          doublereal* alpha,
          doublereal* A, integer* lda,
          doublereal* X, integer* incX,
          doublereal* beta,
          doublereal* Y, integer* incY)
{
    dgemv_(trans, M, N,
           alpha, A, lda, X, incX, beta, Y, incY);
    return 0;
}

/*
 * INPUT
 *  *user_data = (double *) mat
 */
static void
atimes_by_matrix (int n, const double *x, double *b, void *user_data)
{
  char trans = 'T'; /* fortran's memory allocation is transposed */
  int i_1 = 1;
  double d_1 = 1.0;
  double d_0 = 0.0;

  double *mat = (double *)user_data;
  dgemv_ (&trans, &n, &n, &d_1, mat, &n,
	  x, &i_1,
	  &d_0, b, &i_1);
}