Fixed a bug in a LAPACK wrapper

Author: Anton Leykin

M2/Macaulay2/e/lapack.cpp

@@ -2166,13 +2166,17 @@ bool Lapack::eigenvectors(const DMatRR *A,
     }
   else
     {
-
-      ///// YYY REVIEW THIS...!
       // Make the complex arrays of eigvals and eigvecs
       eigvals->resize(size, 1);
       eigvecs->resize(size, size);
       //      DMatCC::ElementType *elems = eigvals->rowMajorArray();
       double* eigenLoc = eigen;
+      for (int i = 0; i < size; ++i) {
+        for (int j = 0; j < size; ++j) {
+            std::cout << eigen[i + j * size] << " ";
+        }
+        std::cout << std::endl;
+      }
       for (int j = 0; j < size; j++, eigenLoc += size)
         {
           eigvals->ring().set_from_doubles(eigvals->entry(j,0), real[j], imag[j]);
@@ -2191,9 +2195,9 @@ bool Lapack::eigenvectors(const DMatRR *A,
               for (int i = 0; i < size; ++i)
                 {
                   eigvecs->ring().set_from_doubles(eigvecs->entry(i,j),
-                                                   eigen[i], eigen[size + i]);
+                                                   eigenLoc[i], eigenLoc[size + i]);
                   eigvecs->ring().set_from_doubles(eigvecs->entry(i,j+1),
-                                                   eigen[i], - eigen[size + i]);
+                                                   eigenLoc[i], - eigenLoc[size + i]);
                 }
             }
         }

M2/Macaulay2/tests/normal/000-core.m2

@@ -1394,8 +1394,14 @@ M = matrix{{1.0,1.0},{0.0,1.0}}
 eigenvalues M
 eigenvectors M
 
+M0 = matrix{{1_RR,-1},{1,1}};
+M = M0++M0;
+(D, P) = eigenvectors M
+
+
 M = matrix{{1.0, 2.0}, {2.0, 1.0}}
 eigenvectors(M, Hermitian=>true)
+assert( 1e-10 > norm(M*P - P * diagonalMatrix(D)))
 
 M = matrix{{1.0, 2.0}, {5.0, 7.0}}
 (eigvals, eigvecs) = eigenvectors M