Commit 070d1372 authored by Jens Jørgen Mortensen's avatar Jens Jørgen Mortensen

Merge branch 'no-blas' into 'master'

Clean up BLAS things

See merge request !591
parents 3e5c90f3 d2e93d14
Pipeline #106623916 passed with stage
in 3 minutes and 23 seconds
......@@ -30,20 +30,13 @@ PyObject* symmetrize_wavefunction(PyObject *self, PyObject *args);
PyObject* symmetrize_return_index(PyObject *self, PyObject *args);
PyObject* symmetrize_with_index(PyObject *self, PyObject *args);
PyObject* map_k_points(PyObject *self, PyObject *args);
PyObject* scal(PyObject *self, PyObject *args);
PyObject* mmm(PyObject *self, PyObject *args);
PyObject* tetrahedron_weight(PyObject *self, PyObject *args);
#ifndef GPAW_WITHOUT_BLAS
PyObject* mmm(PyObject *self, PyObject *args);
PyObject* gemm(PyObject *self, PyObject *args);
PyObject* gemv(PyObject *self, PyObject *args);
PyObject* axpy(PyObject *self, PyObject *args);
PyObject* czher(PyObject *self, PyObject *args);
PyObject* rk(PyObject *self, PyObject *args);
PyObject* r2k(PyObject *self, PyObject *args);
PyObject* dotc(PyObject *self, PyObject *args);
PyObject* dotu(PyObject *self, PyObject *args);
PyObject* multi_dotu(PyObject *self, PyObject *args);
PyObject* multi_axpy(PyObject *self, PyObject *args);
PyObject* NewLocalizedFunctionsObject(PyObject *self, PyObject *args);
#endif
PyObject* NewOperatorObject(PyObject *self, PyObject *args);
PyObject* NewWOperatorObject(PyObject *self, PyObject *args);
PyObject* NewSplineObject(PyObject *self, PyObject *args);
......@@ -51,8 +44,6 @@ PyObject* NewTransformerObject(PyObject *self, PyObject *args);
PyObject* pc_potential(PyObject *self, PyObject *args);
PyObject* add_to_density(PyObject *self, PyObject *args);
PyObject* utilities_gaussian_wave(PyObject *self, PyObject *args);
PyObject* utilities_vdot(PyObject *self, PyObject *args);
PyObject* utilities_vdot_self(PyObject *self, PyObject *args);
PyObject* errorfunction(PyObject *self, PyObject *args);
PyObject* cerf(PyObject *self, PyObject *args);
PyObject* pack(PyObject *self, PyObject *args);
......@@ -69,7 +60,6 @@ PyObject* tci_overlap(PyObject *self, PyObject *args);
PyObject *pwlfc_expand(PyObject *self, PyObject *args);
PyObject *pw_insert(PyObject *self, PyObject *args);
PyObject *pw_precond(PyObject *self, PyObject *args);
PyObject* overlap(PyObject *self, PyObject *args);
PyObject* vdw(PyObject *self, PyObject *args);
PyObject* vdw2(PyObject *self, PyObject *args);
PyObject* spherical_harmonics(PyObject *self, PyObject *args);
......@@ -165,28 +155,19 @@ static PyMethodDef functions[] = {
{"symmetrize_return_index", symmetrize_return_index, METH_VARARGS, 0},
{"symmetrize_with_index", symmetrize_with_index, METH_VARARGS, 0},
{"map_k_points", map_k_points, METH_VARARGS, 0},
{"scal", scal, METH_VARARGS, 0},
{"mmm", mmm, METH_VARARGS, 0},
{"tetrahedron_weight", tetrahedron_weight, METH_VARARGS, 0},
#ifndef GPAW_WITHOUT_BLAS
{"mmm", mmm, METH_VARARGS, 0},
{"gemm", gemm, METH_VARARGS, 0},
{"gemv", gemv, METH_VARARGS, 0},
{"axpy", axpy, METH_VARARGS, 0},
{"czher", czher, METH_VARARGS, 0},
{"rk", rk, METH_VARARGS, 0},
{"r2k", r2k, METH_VARARGS, 0},
{"dotc", dotc, METH_VARARGS, 0},
{"dotu", dotu, METH_VARARGS, 0},
{"multi_dotu", multi_dotu, METH_VARARGS, 0},
{"multi_axpy", multi_axpy, METH_VARARGS, 0},
{"LocalizedFunctions", NewLocalizedFunctionsObject, METH_VARARGS, 0},
#endif
{"Operator", NewOperatorObject, METH_VARARGS, 0},
{"WOperator", NewWOperatorObject, METH_VARARGS, 0},
{"Spline", NewSplineObject, METH_VARARGS, 0},
{"Transformer", NewTransformerObject, METH_VARARGS, 0},
{"add_to_density", add_to_density, METH_VARARGS, 0},
{"utilities_gaussian_wave", utilities_gaussian_wave, METH_VARARGS, 0},
{"utilities_vdot", utilities_vdot, METH_VARARGS, 0},
{"utilities_vdot_self", utilities_vdot_self, METH_VARARGS, 0},
{"eed_region", exterior_electron_density_region, METH_VARARGS, 0},
{"plane_wave_grid", plane_wave_grid, METH_VARARGS, 0},
{"pwlfc_expand", pwlfc_expand, METH_VARARGS, 0},
......@@ -202,7 +183,6 @@ static PyMethodDef functions[] = {
{"XCFunctional", NewXCFunctionalObject, METH_VARARGS, 0},
{"lxcXCFunctional", NewlxcXCFunctionalObject, METH_VARARGS, 0},
{"lxcXCFuncNum", lxcXCFuncNum, METH_VARARGS, 0},
{"overlap", overlap, METH_VARARGS, 0},
{"tci_overlap", tci_overlap, METH_VARARGS, 0},
{"vdw", vdw, METH_VARARGS, 0},
{"vdw2", vdw2, METH_VARARGS, 0},
......@@ -297,7 +277,6 @@ extern PyTypeObject GPAW_MPI_Request_type;
#endif
extern PyTypeObject LFCType;
extern PyTypeObject LocalizedFunctionsType;
extern PyTypeObject OperatorType;
extern PyTypeObject WOperatorType;
extern PyTypeObject SplineType;
......@@ -362,8 +341,6 @@ static PyObject* moduleinit(void)
if (PyType_Ready(&LFCType) < 0)
return NULL;
if (PyType_Ready(&LocalizedFunctionsType) < 0)
return NULL;
if (PyType_Ready(&OperatorType) < 0)
return NULL;
if (PyType_Ready(&WOperatorType) < 0)
......@@ -395,7 +372,6 @@ static PyObject* moduleinit(void)
#endif
Py_INCREF(&LFCType);
Py_INCREF(&LocalizedFunctionsType);
Py_INCREF(&OperatorType);
Py_INCREF(&WOperatorType);
Py_INCREF(&SplineType);
......
......@@ -3,6 +3,7 @@
* Copyright (C) 2007 CSC - IT Center for Science Ltd.
* Please see the accompanying LICENSE file for further information. */
#ifndef GPAW_WITHOUT_BLAS
#include <Python.h>
#define PY_ARRAY_UNIQUE_SYMBOL GPAW_ARRAY_API
#define NO_IMPORT_ARRAY
......@@ -10,39 +11,18 @@
#include "extensions.h"
#ifdef GPAW_NO_UNDERSCORE_BLAS
# define dscal_ dscal
# define zscal_ zscal
# define daxpy_ daxpy
# define zaxpy_ zaxpy
# define dsyrk_ dsyrk
# define zher_ zher
# define zherk_ zherk
# define dsyr2k_ dsyr2k
# define zher2k_ zher2k
# define dgemm_ dgemm
# define zgemm_ zgemm
# define dgemv_ dgemv
# define zgemv_ zgemv
# define ddot_ ddot
#endif
void dscal_(int*n, double* alpha, double* x, int* incx);
void zscal_(int*n, void* alpha, void* x, int* incx);
void daxpy_(int* n, double* alpha,
double* x, int *incx,
double* y, int *incy);
void zaxpy_(int* n, void* alpha,
void* x, int *incx,
void* y, int *incy);
void dsyrk_(char *uplo, char *trans, int *n, int *k,
double *alpha, double *a, int *lda, double *beta,
double *c, int *ldc);
void zher_(char *uplo, int *n,
double *alpha, void *x, int *incx,
void *a, int *lda);
void zherk_(char *uplo, char *trans, int *n, int *k,
double *alpha, void *a, int *lda,
double *beta,
......@@ -63,34 +43,7 @@ void zgemm_(char *transa, char *transb, int *m, int * n,
int *k, void *alpha, void *a, int *lda,
void *b, int *ldb, void *beta,
void *c, int *ldc);
void dgemv_(char *trans, int *m, int * n,
double *alpha, double *a, int *lda,
double *x, int *incx, double *beta,
double *y, int *incy);
void zgemv_(char *trans, int *m, int * n,
void *alpha, void *a, int *lda,
void *x, int *incx, void *beta,
void *y, int *incy);
double ddot_(int *n, void *dx, int *incx, void *dy, int *incy);
PyObject* scal(PyObject *self, PyObject *args)
{
Py_complex alpha;
PyArrayObject* x;
if (!PyArg_ParseTuple(args, "DO", &alpha, &x))
return NULL;
int n = PyArray_DIMS(x)[0];
for (int d = 1; d < PyArray_NDIM(x); d++)
n *= PyArray_DIMS(x)[d];
int incx = 1;
if (PyArray_DESCR(x)->type_num == NPY_DOUBLE)
dscal_(&n, &(alpha.real), DOUBLEP(x), &incx);
else
zscal_(&n, &alpha, (void*)COMPLEXP(x), &incx);
Py_RETURN_NONE;
}
PyObject* gemm(PyObject *self, PyObject *args)
{
......@@ -187,106 +140,6 @@ PyObject* mmm(PyObject *self, PyObject *args)
}
PyObject* gemv(PyObject *self, PyObject *args)
{
Py_complex alpha;
PyArrayObject* a;
PyArrayObject* x;
Py_complex beta;
PyArrayObject* y;
char t = 't';
char* trans = &t;
if (!PyArg_ParseTuple(args, "DOODO|s", &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;
}
PyObject* axpy(PyObject *self, PyObject *args)
{
Py_complex alpha;
PyArrayObject* x;
PyArrayObject* y;
if (!PyArg_ParseTuple(args, "DOO", &alpha, &x, &y))
return NULL;
int n = PyArray_DIMS(x)[0];
for (int d = 1; d < PyArray_NDIM(x); d++)
n *= PyArray_DIMS(x)[d];
int incx = 1;
int incy = 1;
if (PyArray_DESCR(x)->type_num == NPY_DOUBLE)
daxpy_(&n, &(alpha.real),
DOUBLEP(x), &incx,
DOUBLEP(y), &incy);
else
zaxpy_(&n, &alpha,
(void*)COMPLEXP(x), &incx,
(void*)COMPLEXP(y), &incy);
Py_RETURN_NONE;
}
PyObject* czher(PyObject *self, PyObject *args)
{
double alpha;
PyArrayObject* x;
PyArrayObject* a;
if (!PyArg_ParseTuple(args, "dOO", &alpha, &x, &a))
return NULL;
int n = PyArray_DIMS(x)[0];
for (int d = 1; d < PyArray_NDIM(x); d++)
n *= PyArray_DIMS(x)[d];
int incx = 1;
int lda = MAX(1, n);
zher_("l", &n, &(alpha),
(void*)COMPLEXP(x), &incx,
(void*)COMPLEXP(a), &lda);
Py_RETURN_NONE;
}
PyObject* rk(PyObject *self, PyObject *args)
{
double alpha;
......@@ -352,153 +205,4 @@ PyObject* r2k(PyObject *self, PyObject *args)
(void*)COMPLEXP(c), &ldc);
Py_RETURN_NONE;
}
PyObject* dotc(PyObject *self, PyObject *args)
{
PyArrayObject* a;
PyArrayObject* b;
if (!PyArg_ParseTuple(args, "OO", &a, &b))
return NULL;
int n = PyArray_DIMS(a)[0];
for (int i = 1; i < PyArray_NDIM(a); i++)
n *= PyArray_DIMS(a)[i];
int incx = 1;
int incy = 1;
if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)
{
double result;
result = ddot_(&n, (void*)DOUBLEP(a),
&incx, (void*)DOUBLEP(b), &incy);
return PyFloat_FromDouble(result);
}
else
{
double_complex* ap = COMPLEXP(a);
double_complex* bp = COMPLEXP(b);
double_complex z = 0.0;
for (int i = 0; i < n; i++)
z += conj(ap[i]) * bp[i];
return PyComplex_FromDoubles(creal(z), cimag(z));
}
}
PyObject* dotu(PyObject *self, PyObject *args)
{
PyArrayObject* a;
PyArrayObject* b;
if (!PyArg_ParseTuple(args, "OO", &a, &b))
return NULL;
int n = PyArray_DIMS(a)[0];
for (int i = 1; i < PyArray_NDIM(a); i++)
n *= PyArray_DIMS(a)[i];
int incx = 1;
int incy = 1;
if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)
{
double result;
result = ddot_(&n, (void*)DOUBLEP(a),
&incx, (void*)DOUBLEP(b), &incy);
return PyFloat_FromDouble(result);
}
else
{
double_complex* ap = COMPLEXP(a);
double_complex* bp = COMPLEXP(b);
double_complex z = 0.0;
for (int i = 0; i < n; i++)
z += ap[i] * bp[i];
return PyComplex_FromDoubles(creal(z), cimag(z));
}
}
PyObject* multi_dotu(PyObject *self, PyObject *args)
{
PyArrayObject* a;
PyArrayObject* b;
PyArrayObject* c;
if (!PyArg_ParseTuple(args, "OOO", &a, &b, &c))
return NULL;
int n0 = PyArray_DIMS(a)[0];
int n = PyArray_DIMS(a)[1];
for (int i = 2; i < PyArray_NDIM(a); i++)
n *= PyArray_DIMS(a)[i];
int incx = 1;
int incy = 1;
if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)
{
double *ap = DOUBLEP(a);
double *bp = DOUBLEP(b);
double *cp = DOUBLEP(c);
for (int i = 0; i < n0; i++)
{
cp[i] = ddot_(&n, (void*)ap,
&incx, (void*)bp, &incy);
ap += n;
bp += n;
}
}
else
{
double_complex* ap = COMPLEXP(a);
double_complex* bp = COMPLEXP(b);
double_complex* cp = COMPLEXP(c);
for (int i = 0; i < n0; i++)
{
cp[i] = 0.0;
for (int j = 0; j < n; j++)
cp[i] += ap[j] * bp[j];
ap += n;
bp += n;
}
}
Py_RETURN_NONE;
}
PyObject* multi_axpy(PyObject *self, PyObject *args)
{
PyArrayObject* alpha;
PyArrayObject* x;
PyArrayObject* y;
if (!PyArg_ParseTuple(args, "OOO", &alpha, &x, &y))
return NULL;
int n0 = PyArray_DIMS(x)[0];
int n = PyArray_DIMS(x)[1];
for (int d = 2; d < PyArray_NDIM(x); d++)
n *= PyArray_DIMS(x)[d];
int incx = 1;
int incy = 1;
if (PyArray_DESCR(alpha)->type_num == NPY_DOUBLE)
{
if (PyArray_DESCR(x)->type_num == NPY_CDOUBLE)
n *= 2;
double *ap = DOUBLEP(alpha);
double *xp = DOUBLEP(x);
double *yp = DOUBLEP(y);
for (int i = 0; i < n0; i++)
{
daxpy_(&n, &ap[i],
(void*)xp, &incx,
(void*)yp, &incy);
xp += n;
yp += n;
}
}
else
{
double_complex *ap = COMPLEXP(alpha);
double_complex *xp = COMPLEXP(x);
double_complex *yp = COMPLEXP(y);
for (int i = 0; i < n0; i++)
{
zaxpy_(&n, (void*)(&ap[i]),
(void*)xp, &incx,
(void*)yp, &incy);
xp += n;
yp += n;
}
}
Py_RETURN_NONE;
}
#endif
......@@ -28,9 +28,7 @@
# define Py_RETURN_NONE return Py_INCREF(Py_None), Py_None
#endif
#define INLINE inline
static INLINE void* gpaw_malloc(size_t n)
static inline void* gpaw_malloc(size_t n)
{
void* p = malloc(n);
assert(p != NULL);
......
......@@ -3,25 +3,10 @@
* Please see the accompanying LICENSE file for further information. */
#include "extensions.h"
#include "localized_functions.h"
#include "bmgs/bmgs.h"
#include "spline.h"
#include <complex.h>
#ifdef GPAW_NO_UNDERSCORE_BLAS
# define dgemv_ dgemv
# define dgemm_ dgemm
#endif
int dgemv_(char *trans, int *m, int * n,
double *alpha, double *a, int *lda,
double *x, int *incx, double *beta,
double *y, int *incy);
int dgemm_(char *transa, char *transb, int *m, int * n,
int *k, const double *alpha, double *a, int *lda,
double *b, int *ldb, double *beta,
double *c, int *ldc);
// +-----------n
// +----m +----m | +----c+m |
......@@ -32,16 +17,16 @@ int dgemm_(char *transa, char *transb, int *m, int * n,
// | |
// 0-----------+
void cut(const double* a, const int n[3], const int c[3],
const double* v,
double* b, const int m[3])
const double* v,
double* b, const int m[3])
{
a += c[2] + (c[1] + c[0] * n[1]) * n[2];
for (int i0 = 0; i0 < m[0]; i0++)
{
for (int i1 = 0; i1 < m[1]; i1++)
{
for (int i2 = 0; i2 < m[2]; i2++)
b[i2] = v[i2] * a[i2];
for (int i2 = 0; i2 < m[2]; i2++)
b[i2] = v[i2] * a[i2];
a += n[2];
b += m[2];
v += m[2];
......@@ -206,168 +191,3 @@ PyObject *tci_overlap(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
PyObject * overlap(PyObject* self, PyObject *args)
{
PyObject* lfs_b_obj;
PyArrayObject* m_b_obj;
PyArrayObject* phase_bk_obj;
PyArrayObject* vt_sG_obj;
PyArrayObject* Vt_skmm_obj;
if (!PyArg_ParseTuple(args, "OOOOO", &lfs_b_obj, &m_b_obj, &phase_bk_obj,
&vt_sG_obj, &Vt_skmm_obj))
return NULL;
int nk = PyArray_DIMS(phase_bk_obj)[1];
int nm = PyArray_DIMS(Vt_skmm_obj)[2];
int nspins = PyArray_DIMS(vt_sG_obj)[0];
const long *m_b = LONGP(m_b_obj);
const double complex *phase_bk = COMPLEXP(phase_bk_obj);
const double *vt_sG = DOUBLEP(vt_sG_obj);
double *Vt_smm = 0;
double complex *Vt_skmm = 0;
if (nk == 0)
Vt_smm = DOUBLEP(Vt_skmm_obj);
else
Vt_skmm = COMPLEXP(Vt_skmm_obj);
int nb = PyList_Size(lfs_b_obj);
int nmem = 0;
double* a1 = 0;
for (int b1 = 0; b1 < nb; b1++)
{
const LocalizedFunctionsObject* lf1 =
(const LocalizedFunctionsObject*)PyList_GetItem(lfs_b_obj, b1);
int m1 = m_b[b1];
int nao1 = lf1->nf;
double* f1 = lf1->f;
double* vt1 = GPAW_MALLOC(double, lf1->ng0 * nspins);
for (int s = 0; s < nspins; s++)
bmgs_cut(vt_sG + s * lf1->ng, lf1->size, lf1->start,
vt1 + s * lf1->ng0, lf1->size0);
for (int b2 = b1; b2 < nb; b2++)
{
const LocalizedFunctionsObject* lf2 =
(const LocalizedFunctionsObject*)PyList_GetItem(lfs_b_obj, b2);
int beg[3];
int end[3];
int size[3];
int beg1[3];
int beg2[3];
bool overlap = true;
for (int c = 0; c < 3; c++)
{
beg[c] = MAX(lf1->start[c], lf2->start[c]);
end[c] = MIN(lf1->start[c] + lf1->size0[c],
lf2->start[c] + lf2->size0[c]);
size[c] = end[c] - beg[c];
if (size[c] <= 0)
{
overlap = false;
continue;
}
beg1[c] = beg[c] - lf1->start[c];
beg2[c] = beg[c] - lf2->start[c];
}
int nao2 = lf2->nf;
if (overlap)
{
int ng = size[0] * size[1] * size[2];
int n = ng * (nao1 + nao2) + nao1 * nao2;
if (n > nmem)
{
if (nmem != 0)
free(a1);
nmem = n;
a1 = GPAW_MALLOC(double, nmem);
}
double* a2 = a1 + ng * nao1;
double* H = a2 + ng * nao2;
double* f2 = lf2->f;
double* vt2 = lf2->w;
double dv = lf1->dv;
int m2 = m_b[b2];
if (b2 > b1)
for (int i = 0; i < nao2; i++)
bmgs_cut(f2 + i * lf2->ng0, lf2->size0, beg2,
a2 + i * ng, size);
else
a2 = f2;
for (int s = 0; s < nspins; s++)
{
if (b2 > b1)
{
bmgs_cut(vt1 + s * lf1->ng0, lf1->size0, beg1, vt2, size);
for (int i = 0; i < nao1; i++)
cut(f1 + i * lf1->ng0, lf1->size0, beg1, vt2,
a1 + i * ng, size);
}
else
{
for (int i1 = 0; i1 < nao1; i1++)
for (int g = 0; g < ng; g++)
a1[i1 * ng + g] = (vt1[g + s * lf1->ng0] *
f1[i1 * ng + g]);
}
double zero = 0.0;
dgemm_("t", "n", &nao2, &nao1, &ng, &dv,
a2, &ng, a1, &ng, &zero, H, &nao2);
if (nk == 0)
{
double* Vt_mm = (Vt_smm + s * nm * nm + m1 + m2 * nm);
if (b2 == b1)
for (int i1 = 0; i1 < nao1; i1++)
for (int i2 = i1; i2 < nao2; i2++)
Vt_mm[i1 + i2 * nm] += H[i2 + i1 * nao2];
else if (m1 == m2)
for (int i1 = 0; i1 < nao1; i1++)
for (int i2 = i1; i2 < nao2; i2++)
Vt_mm[i1 + i2 * nm] += (H[i2 + i1 * nao2] +
H[i1 + i2 * nao2]);
else
for (int ii = 0, i1 = 0; i1 < nao1; i1++)
for (int i2 = 0; i2 < nao2; i2++, ii++)
Vt_mm[i1 + i2 * nm] += H[ii];
}
else
for (int k = 0; k < nk; k++)
{
double complex* Vt_mm = (Vt_skmm +
(s * nk + k) * nm * nm +
m1 + m2 * nm);
if (b2 == b1)
for (int i1 = 0; i1 < nao1; i1++)
for (int i2 = i1; i2 < nao2; i2++)
Vt_mm[i1 + i2 * nm] += H[i2 + i1 * nao2];
else
{
double complex phase = \
(phase_bk[b1 * nk + k] *
conj(phase_bk[b2 * nk + k]));
if (m1 == m2)
for (int i1 = 0; i1 < nao1; i1++)
for (int i2 = i1; i2 < nao2; i2++)
Vt_mm[i1 + i2 * nm] += \
(phase * H[i2 + i1 * nao2] +
conj(phase) * H[i1 + i2 * nao2]);
else
for (int ii = 0, i1 = 0; i1 < nao1; i1++)
for (int i2 = 0; i2 < nao2; i2++, ii++)
Vt_mm[i1 + i2 * nm] += phase * H[ii];
}
}
}
}
}
free(vt1);
}
if (nmem != 0)
free(a1);
Py_RETURN_NONE;
}
This diff is collapsed.
This diff is collapsed.
/* Copyright (C) 2003-2007 CAMP
* Copyright (C) 2007-2008 CAMd
* Please see the accompanying LICENSE file for further information. */
#include <Python.h>
typedef struct
{
PyObject_HEAD
double dv; // volume per grid point
int size[3]; // dimensions of big box
int start[3]; // corner of small box
int size0[3]; // dimensions of small box
int ng; // number of grid points in big box
int ng0; // number of grid points in small box
int nf; // number of localized functions
int nfd; // number of derivatives: zero or 3*nf
// pointers to size0 arrays:
double* f; // localized functions
double* fd; // xyz-derivatives of localized functions
double* w; // work array for one double or double complex array
} LocalizedFunctionsObject;
......@@ -417,54 +417,6 @@ PyObject* utilities_gaussian_wave(PyObject *self, PyObject *args)
}
/* vdot
*
* If a and b are input vectors,
* a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + ...
* is returned.
*/
PyObject* utilities_vdot(PyObject *self, PyObject *args)
{