Use syevd/heevd in JD stepper, workaround lwork size calculation

src/JDWavefunctionStepper.C

@@ -257,7 +257,7 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
         // q is (2n,2n)
         DoubleMatrix q(h.context(),h.n(),h.n(),h.nb(),h.nb());
         tmap_["jd_syev"].start();
-        h.syev('l',w,q);
+        h.syevd('l',w,q);
         tmap_["jd_syev"].stop();
 
         // compute the first n eigenvectors and store in wf
@@ -332,7 +332,7 @@ void JDWavefunctionStepper::update(Wavefunction& dwf)
         // q is (2n,2n)
         ComplexMatrix q(h.context(),h.n(),h.n(),h.nb(),h.nb());
         tmap_["jd_heev"].start();
-        h.heev('l',w,q);
+        h.heevd('l',w,q);
         tmap_["jd_heev"].stop();
 
         // compute the first n eigenvectors and store in wf

src/Matrix.C

@@ -2819,6 +2819,8 @@ void DoubleMatrix::syev(char uplo, valarray<double>& w, DoubleMatrix& z)
             &info);
 
     lwork = (int) (tmpwork + 1);
+    // set lwork to max value among all tasks
+    ctxt_.imax(1,1,&lwork,1);
     double* work=new double[lwork];
     pdsyev(&jobz, &uplo, &m_, val, &ione, &ione, desc_, &w[0],
             z.val, &ione, &ione, z.desc_, work, &lwork,
@@ -2876,7 +2878,10 @@ void DoubleMatrix::syevd(char uplo, valarray<double>& w, DoubleMatrix& z)
             &tmpiwork, &liwork, &info);
 
     lwork = (int) (tmpwork + 1);
+    // set lwork to max value among all tasks
+    ctxt_.imax(1,1,&lwork,1);
     double* work=new double[lwork];
+
     liwork = tmpiwork;
     int* iwork = new int[liwork];
     pdsyevd(&jobz, &uplo, &m_, val, &ione, &ione, desc_, &w[0],