Index: c/_gpaw.c =================================================================== --- c/_gpaw.c (revision 5444) +++ c/_gpaw.c (working copy) @@ -76,6 +76,7 @@ PyObject* scalapack_redist(PyObject *self, PyObject *args); PyObject* scalapack_diagonalize_dc(PyObject *self, PyObject *args); PyObject* scalapack_diagonalize_ex(PyObject *self, PyObject *args); +PyObject* scalapack_diagonalize2_ex(PyObject *self, PyObject *args); PyObject* scalapack_inverse_cholesky(PyObject *self, PyObject *args); #endif @@ -146,6 +147,7 @@ {"scalapack_redist", scalapack_redist, METH_VARARGS, 0}, {"scalapack_diagonalize_dc", scalapack_diagonalize_dc, METH_VARARGS, 0}, {"scalapack_diagonalize_ex", scalapack_diagonalize_ex, METH_VARARGS, 0}, + {"scalapack_diagonalize2_ex", scalapack_diagonalize2_ex, METH_VARARGS, 0}, {"scalapack_inverse_cholesky", scalapack_inverse_cholesky, METH_VARARGS, 0}, #endif #ifdef GPAW_HPM Index: c/blacs.c =================================================================== --- c/blacs.c (revision 5444) +++ c/blacs.c (working copy) @@ -615,6 +615,21 @@ } } +PyObject* scalapack_diagonalize2_ex(PyObject *self, PyObject *args) +{ + PyArrayObject* h_obj, s_obj; + PyArrayObject* c_obj, eps_obj; + PyArrayObject* descriptor; + char uplo; + + if (!PyArg_ParseTuple(args, "OOOOOOc", &h_obj, &s_obj, &c_obj, &eps_obj, + &descriptor, &uplo)) { + return NULL; + } + + Py_RETURN_NONE; +} + PyObject* scalapack_diagonalize_ex(PyObject *self, PyObject *args) { // Expert Driver for QR algorithm @@ -635,6 +650,9 @@ int eigvalm, nz; static int one = 1; + //PyArrayObject* c_obj; + //PyArrayObject* eps_obj; + double vl, vu; // not used when range = 'A' or 'I' char jobz = 'V'; // eigenvectors also @@ -644,7 +662,8 @@ // char cmach = 'S'; // most acccurate eigenvalues int isgeneral; // flag for general diagonalize - if (!PyArg_ParseTuple(args, "OOc|O", &a_obj, &adesc, &uplo, &b_obj)) + if (!PyArg_ParseTuple(args, "OOcOOO", &a_obj, &adesc, &uplo, &b_obj, + &z_obj, &w_obj)) return NULL; // adesc, @@ -703,6 +722,8 @@ // Convergence tolerance // most orthogonal eigenvectors double abstol = pdlamch_(&z_ConTxt, &cmach); + printf("first %f\n", abstol); + fflush(stdout); // most accurate eigenvalues // double abstol = 2.0*pdlamch_(&z_ConTxt, &cmach); @@ -714,18 +735,24 @@ int z_locN = numroc_(&z_n, &z_nb, &z_mycol, &z_csrc, &z_npcol); // Eigenvectors - npy_intp z_dims[2] = {z_locM, z_locN}; + /*npy_intp z_dims[2] = {z_locM, z_locN}; if (a_obj->descr->type_num == PyArray_DOUBLE) z_obj = (PyArrayObject*)PyArray_EMPTY(2, z_dims, NPY_DOUBLE, NPY_F_CONTIGUOUS); else z_obj = (PyArrayObject*)PyArray_EMPTY(2, z_dims, NPY_CDOUBLE, NPY_F_CONTIGUOUS); + */ + //z_obj = c_obj; + // Eigenvalues, since w_obj is 1D-array NPY_F_CONTIGUOUS not needed here npy_intp w_dims[1] = {n}; - w_obj = (PyArrayObject*)PyArray_SimpleNew(1, w_dims, NPY_DOUBLE); + PyArrayObject* w1_obj; + w1_obj = (PyArrayObject*)PyArray_SimpleNew(1, w_dims, NPY_DOUBLE); + //w_obj = (PyArrayObject*)eps_obj; + // Query part, need to find the optimal size of a number of work arrays int info; int *ifail; @@ -752,7 +779,10 @@ DOUBLEP(z_obj), &one, &one, zdesc, &d_work, &querywork, &i_work, &querywork, ifail, iclustr, gap, &info); - else + else { + printf("middle %f\n", abstol); + fflush(stdout); + pdsygvx_(&ibtype, &jobz, &range, &uplo, &n, DOUBLEP(a_obj), &one, &one, INTP(adesc), DOUBLEP(b_obj), &one, &one, INTP(adesc), @@ -761,6 +791,7 @@ DOUBLEP(z_obj), &one, &one, zdesc, &d_work, &querywork, &i_work, &querywork, ifail, iclustr, gap, &info); + } lwork = (int)(d_work); } else @@ -789,6 +820,10 @@ } if (info != 0) { + printf("last %f", abstol); + fflush(stdout); + + PyErr_SetString(PyExc_RuntimeError, "scalapack_diagonalize_ex error in query."); return NULL; @@ -857,15 +892,16 @@ free(gap); free(iclustr); free(ifail); - - PyObject* values = Py_BuildValue("(OO)", w_obj, z_obj); - Py_DECREF(w_obj); - Py_DECREF(z_obj); - return values; + //PyObject* values = Py_BuildValue("(OO)", w_obj, z_obj); + //Py_DECREF(w_obj); + //Py_DECREF(z_obj); + Py_RETURN_NONE; + //return values; } else { - return Py_BuildValue("(OO)", Py_None, Py_None); + Py_RETURN_NONE; + //return Py_BuildValue("(OO)", Py_None, Py_None); } } Index: gpaw/blacs.py =================================================================== --- gpaw/blacs.py (revision 5446) +++ gpaw/blacs.py (working copy) @@ -17,6 +17,7 @@ self.bd = bd self.kpt_comm = kpt_comm self.gd = gd + self.nao = nao bcommsize = bd.comm.size gcommsize = gd.comm.size @@ -68,8 +69,19 @@ redistributor.redistribute(colS, sqrS) redistributor.redistribute(colH, sqrH) - eps_n[:], H1_MM = scalapack_diagonalize_ex(sqrH.A_mn, descriptor2, - 'U', sqrS.A_mn) + H0_MM = np.zeros_like(sqrH.A_mn) + eps_M = np.zeros(self.nao, order='F') + + x=scalapack_diagonalize_ex(sqrH.A_mn, descriptor2, + 'U', sqrS.A_mn, + H0_MM, eps_M) + #eps_n[:], H1_MM = x + + #print self.supercomm.rank, len(eps_n), len(eps_M) + + eps_n[:] = eps_M + H1_MM = H0_MM + if H1_MM is not None: sqrH.A_mn[:] = H1_MM Index: gpaw/utilities/blacs.py =================================================================== --- gpaw/utilities/blacs.py (revision 5446) +++ gpaw/utilities/blacs.py (working copy) @@ -67,7 +67,10 @@ assert uplo in ['U','L'] return _gpaw.scalapack_diagonalize_dc(a_obj, adesc, uplo) -def scalapack_diagonalize_ex(a_obj, adescriptor, uplo, b_obj=None): +def scalapack_diagonalize2_ex(H, S, C, eps, uplo): + _gpaw.scalapack_diagonalize2_ex(H, S, C, eps, uplo) + +def scalapack_diagonalize_ex(a_obj, adescriptor, uplo, b_obj, C_MM, eps_M): adesc = adescriptor.asarray() if a_obj is not None: assert a_obj.ndim == 2 @@ -81,7 +84,8 @@ assert b_obj is None assert len(adesc) == 9 assert uplo in ['U','L'] - return _gpaw.scalapack_diagonalize_ex(a_obj, adesc, uplo, b_obj) + return _gpaw.scalapack_diagonalize_ex(a_obj, adesc, uplo, b_obj, + C_MM, eps_M) def scalapack_inverse_cholesky(a_obj, adesc, uplo): if a_obj is not None: Index: customize.py =================================================================== --- customize.py (revision 5444) +++ customize.py (working copy) @@ -1,57 +1,26 @@ -#User provided customizations for the gpaw setup - -#Here, one can override the default arguments, or append own -#arguments to default ones -#To override use the form -# libraries = ['somelib','otherlib'] -#To append use the form -# libraries += ['somelib','otherlib'] - -# Valid values for scalapack are False, or True: -# False (the default) - no ScaLapack compiled in -# True - ScaLapack compiled in -#scalapack = True - -#compiler = 'mpcc' -#libraries = [] -#libraries += [] - -#library_dirs = [] -#library_dirs += [] - -#include_dirs = [] -#include_dirs += [] - -#extra_link_args = [] -#extra_link_args += [] - -#extra_compile_args = [] -#extra_compile_args += [] - -#runtime_library_dirs = [] -#runtime_library_dirs += [] - -#extra_objects = [] -#extra_objects += [] - -#define_macros = [] -#define_macros += [] - -#mpicompiler = None -#mpilinker = None -#mpi_libraries = [] -#mpi_libraries += [] - -#mpi_library_dirs = [] -#mpi_library_dirs += [] - -#mpi_include_dirs = [] -#mpi_include_dirs += [] - -#mpi_runtime_library_dirs = [] -#mpi_runtime_library_dirs += [] - -#mpi_define_macros = [] -#mpi_define_macros += [] - -#platform_id = '' +scalapack = True +compiler = 'gcc43' +libraries = [ + 'gfortran', 'scalapack', 'mpiblacsF77init', 'mpiblacs', 'scalapack', + 'acml', + 'mpi', 'mpi_f77' + ] +library_dirs = [ + '/opt/openmpi/1.3.3-1.el5.fys.gfortran43.4.3.2/lib64', + '/opt/acml/4.3.0/gfortran4364/lib', + '/opt/blacs/1.1/24.el5.fys.gfortran43.4.3.2.openmpi.1.3.3/lib64', + '/opt/scalapack/1.8.0/1.el5.fys.gfortran43.4.3.2.openmpi.1.3.3.acml.4.3.0.acml.4.3.0/lib64' + ] +include_dirs += ['/opt/openmpi/1.3.3-1.el5.fys.gfortran43.4.3.2/include'] +extra_link_args = [ + '-Wl,-rpath=/opt/openmpi/1.3.3-1.el5.fys.gfortran43.4.3.2/lib64,' + '-rpath=/opt/acml/4.3.0/gfortran4364/lib,' + '-rpath=/opt/blacs/1.1/24.el5.fys.gfortran43.4.3.2.openmpi.1.3.3/lib64,' + '-rpath=/opt/scalapack/1.8.0/1.el5.fys.gfortran43.4.3.2.openmpi.1.3.3.acml.4.3.0.acml.4.3.0/lib64' + ] +extra_compile_args = ['-O3', '-std=c99', '-funroll-all-loops', '-fPIC'] +define_macros += [('GPAW_NO_UNDERSCORE_CBLACS', '1')] +define_macros += [('GPAW_NO_UNDERSCORE_CSCALAPACK', '1')] +mpicompiler = '/opt/openmpi/1.3.3-1.el5.fys.gfortran43.4.3.2/bin/mpicc' +mpilinker = mpicompiler +platform_id = 'xeon'