PyObject* scalapack_diagonalize_ex(PyObject *self, PyObject *args)
{
  // Expert Driver for QR algorithm
  // Computes *all* eigenvalues and eigenvectors

  PyArrayObject* a_obj; // Hamiltonian matrix
  PyArrayObject* b_obj; // overlap matrix
  PyArrayObject* adesc; // Hamintonian matrix descriptor
  PyArrayObject* z_obj; // eigenvector matrix
  PyArrayObject* w_obj; // eigenvalue array

  //PyArrayObject* H_MM_obj;
  //PyArrayObject* eps_M_obj;
  
  int ibtype  =  1; // Solve H*psi = lambda*S*psi
  int z_mycol = -1;
  int z_myrow = -1;
  int z_nprow, z_npcol;
  int z_type, z_ConTxt, z_m, z_n, z_mb, z_nb, z_rsrc, z_csrc, z_lld;
  int zdesc[9];
  int il, iu;  // not used when range = 'A' or 'V'
  int eigvalm, nz;
  static int one = 1;

  double vl, vu; // not used when range = 'A' or 'I'

  char jobz = 'V'; // eigenvectors also
  char range = 'A'; // all eigenvalues
  char uplo;
  char cmach = 'U'; // most orthogonal eigenvectors
  // char cmach = 'S'; // most acccurate eigenvalues

  int isgeneral; // flag for general diagonalize
  if (!PyArg_ParseTuple(args, "OOcOOO", &a_obj, &adesc, &uplo, &b_obj,
                        //&H_MM_obj, &eps_M_obj
                        &z_obj, &w_obj))
    return NULL;

  // adesc,
  // bdesc = adesc
  // This is generally not required, as long as the alignment properties
  // are satisfied, see pdsygvx.f. In the context of GPAW, don't see why
  // bdesc would not be equal to adesc so I am just hard-coding it in.
  int a_type   = INTP(adesc)[0];
  int a_ConTxt = INTP(adesc)[1];
  int a_m      = INTP(adesc)[2];
  int a_n      = INTP(adesc)[3];
  int a_mb     = INTP(adesc)[4];
  int a_nb     = INTP(adesc)[5];
  int a_rsrc   = INTP(adesc)[6];
  int a_csrc   = INTP(adesc)[7];
  int a_lld    = INTP(adesc)[8];

  // Note that A is symmetric, so n = a_m = a_n;
  // We do not test for that here.
  int n = a_n;
  //printf("C: n=%d", n);

  // zdesc = adesc
  // This is generally not required, as long as the alignment properties
  // are satisfied, see pdsygvx.f. In the context of GPAW, don't see why
  // zdesc would not be equal to adesc so I am just hard-coding it in.
  z_type   = a_type;
  z_ConTxt = a_ConTxt;
  z_m      = a_m;
  z_n      = a_n;
  z_mb     = a_mb;
  z_nb     = a_nb;
  z_rsrc   = a_rsrc;
  z_csrc   = a_csrc;
  z_lld    = a_lld;
  zdesc[0] = z_type;
  zdesc[1] = z_ConTxt;
  zdesc[2] = z_m;
  zdesc[3] = z_n;
  zdesc[4] = z_mb;
  zdesc[5] = z_nb;
  zdesc[6] = z_rsrc;
  zdesc[7] = z_csrc;
  zdesc[8] = z_lld;

  Cblacs_gridinfo_(z_ConTxt, &z_nprow, &z_npcol, &z_myrow, &z_mycol);

  assert(z_ConTxt != -1);

  if (z_ConTxt != -1)
    {
      if (PyArray_Check(b_obj))
        isgeneral = 1;
      else
        isgeneral = 0;

      // Convergence tolerance
      // most orthogonal eigenvectors
      double abstol = pdlamch_(&z_ConTxt, &cmach);
      double orfac = -1.0;

      // Query part, need to find the optimal size of a number of work arrays
      int info;
      int *ifail;
      ifail = GPAW_MALLOC(int, n);
      int *iclustr;
      iclustr = GPAW_MALLOC(int, 2*z_nprow*z_npcol);
      double  *gap;
      gap = GPAW_MALLOC(double, z_nprow*z_npcol);
      int querywork = -1;
      int* iwork;
      int liwork;
      int lwork;
      int lrwork;
      int i_work;
      double d_work;
      double D_work[1];
      //double* D_work = &d_work;
      double_complex c_work;

      if (a_obj->descr->type_num == PyArray_DOUBLE)
        {
	  if(!isgeneral){
	    pdsyevx_(&jobz, &range, &uplo, &n,
                     DOUBLEP(a_obj), &one, &one, INTP(adesc),
		     &vl, &vu, &il, &iu, &abstol, &eigvalm,
		     &nz, DOUBLEP(w_obj), &orfac,
                     DOUBLEP(z_obj), &one, &one, zdesc,
                     &d_work, &querywork,  &i_work, &querywork,
                     ifail, iclustr, gap, &info);
          }
	  else {
	    pdsygvx_(&ibtype, &jobz, &range, &uplo, &n,
                     DOUBLEP(a_obj), &one, &one, INTP(adesc),
                     DOUBLEP(b_obj), &one, &one, INTP(adesc),
                     &vl, &vu, &il, &iu, &abstol, &eigvalm,
		     &nz, DOUBLEP(w_obj), &orfac,
                     DOUBLEP(z_obj),  &one, &one, zdesc,
                     D_work, 
                     &querywork, &i_work, &querywork,
                     ifail, iclustr, gap, &info);
          }
          lwork = (int)(D_work[0]);
          //lwork = (int)(d_work);
          printf("abstol after=%f\n", abstol);

        }
      else
        {
	  if(!isgeneral)
	    pzheevx_(&jobz, &range, &uplo, &n,
                     (void*)COMPLEXP(a_obj),&one, &one, INTP(adesc),
		     &vl, &vu, &il, &iu, &abstol, &eigvalm,
		     &nz, DOUBLEP(w_obj), &orfac,
                     (void*)COMPLEXP(z_obj), &one, &one, zdesc,
                     &c_work, &querywork, &d_work, &querywork,
                     &i_work, &querywork,
                     ifail, iclustr, gap, &info);
	  else
	    pzhegvx_(&ibtype, &jobz, &range, &uplo, &n,
                     (void*)COMPLEXP(a_obj), &one, &one, INTP(adesc),
                     (void*)COMPLEXP(b_obj), &one, &one, INTP(adesc),
                     &vl, &vu, &il, &iu, &abstol, &eigvalm,
		     &nz, DOUBLEP(w_obj), &orfac,
                     (void*)COMPLEXP(z_obj), &one, &one, zdesc,
                     &c_work, &querywork, &d_work, &querywork,
		     &i_work, &querywork,
                     ifail, iclustr, gap, &info);
          lwork = (int)(c_work);
          lrwork = (int)(d_work);
        }
      if (info != 0)
        {
          printf("ERROR %d\n", info);
          PyErr_SetString(PyExc_RuntimeError,
                          "scalapack_diagonalize_ex error in query.");
          return NULL;
        }

      // Computation part
      lwork = lwork + (n-1)*n;
      liwork = i_work;
      iwork = GPAW_MALLOC(int, liwork);
      if (a_obj->descr->type_num == PyArray_DOUBLE)
        {
          double* work = GPAW_MALLOC(double, lwork);
	  if (!isgeneral)
	    pdsyevx_(&jobz, &range, &uplo, &n,
                     DOUBLEP(a_obj), &one, &one, INTP(adesc),
		     &vl, &vu, &il, &iu, &abstol, &eigvalm,
		     &nz, DOUBLEP(w_obj), &orfac,
                     DOUBLEP(z_obj), &one, &one, zdesc,
                     work, &lwork, iwork, &liwork,
                     ifail, iclustr, gap, &info);
	  else
	    pdsygvx_(&ibtype, &jobz, &range, &uplo, &n,
                     DOUBLEP(a_obj), &one, &one, INTP(adesc),
                     DOUBLEP(b_obj), &one, &one, INTP(adesc),
                     &vl, &vu, &il, &iu, &abstol, &eigvalm,
		     &nz, DOUBLEP(w_obj), &orfac,
                     DOUBLEP(z_obj), &one, &one,  zdesc,
                     work, &lwork,  iwork, &liwork,
                     ifail, iclustr, gap, &info);
          free(work);
        }
      else
        {
          double_complex* work = GPAW_MALLOC(double_complex, lwork);
          double* rwork = GPAW_MALLOC(double, lrwork);
	  if (!isgeneral)
	    pzheevx_(&jobz, &range, &uplo, &n,
                     (void*)COMPLEXP(a_obj), &one, &one, INTP(adesc),
		     &vl, &vu, &il, &iu, &abstol, &eigvalm,
		     &nz, DOUBLEP(w_obj), &orfac,
                     (void*)COMPLEXP(z_obj), &one, &one, zdesc, work,
                     &lwork, rwork, &lrwork,
		     iwork, &liwork,
                     ifail, iclustr, gap, &info);
	  else
	    pzhegvx_(&ibtype, &jobz, &range, &uplo, &n,
                     (void*)COMPLEXP(a_obj), &one, &one, INTP(adesc),
                     (void*)COMPLEXP(b_obj), &one, &one, INTP(adesc),
                     &vl, &vu, &il, &iu, &abstol, &eigvalm,
		     &nz, DOUBLEP(w_obj), &orfac,
                     (void*)COMPLEXP(z_obj), &one, &one, zdesc,
                     work, &lwork, rwork, &lrwork,
		     iwork, &liwork,
                     ifail, iclustr, gap, &info);
          free(rwork);
          free(work);
        }
      if (info != 0)
        {
          PyErr_SetString(PyExc_RuntimeError,
                          "scalapack_diagonalize_ex error in compute.");
          return NULL;
        }

      free(iwork);
      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;
      Py_RETURN_NONE;
    }
  else
    {
      Py_RETURN_NONE;
      //return Py_BuildValue("(OO)", Py_None, Py_None);
    }
}