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'
