From 6a69c33b8b0a8452d8c04a7d28564e11abea8d41 Mon Sep 17 00:00:00 2001
From: Andrew Jewett <jewett.aij@gmail.com>
Date: Wed, 29 Jan 2020 16:40:02 -0800
Subject: [PATCH] fixed run time errors (mostly in test.cpp) introduced by the
 recent changes in the Jacobi API (adding a ConstMatrix template argument)

---
 include/jacobi.hpp       |  1 +
 include/matrix_alloc.hpp | 23 +++++++++++++++++++++--
 tests/test.cpp           | 10 +++++-----
 3 files changed, 27 insertions(+), 7 deletions(-)

diff --git a/include/jacobi.hpp b/include/jacobi.hpp
index 7a1d54a..707c972 100644
--- a/include/jacobi.hpp
+++ b/include/jacobi.hpp
@@ -459,6 +459,7 @@ template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
 void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
 Init() {
   n = 0;
+  M = nullptr;
   max_idx_row = nullptr;
 }
 
diff --git a/include/matrix_alloc.hpp b/include/matrix_alloc.hpp
index d94bff7..5f9f52b 100644
--- a/include/matrix_alloc.hpp
+++ b/include/matrix_alloc.hpp
@@ -33,6 +33,24 @@ void Dealloc2D(Entry ***paaX);      //!< pointer to 2-D multidimensional array
 
 // --- IMPLEMENTATION ---
 
+/// @brief
+/// Allocate a 2-dimensional table row-major order.
+/// In this version, the the size of the array specified by 2 integer arguments.
+template<typename Entry, typename Integer>
+void Alloc2D(Integer const size[2], //!< size of the array in x,y directions
+             Entry ***paaX)         //!< pointer to 2-D multidimensional array
+{
+  assert(paaX);
+  // Allocate a conventional 2-dimensional
+  // pointer-to-a-pointer data structure.
+  *paaX = new Entry* [size[0]];
+  (*paaX)[0] = new Entry [size[0] * size[1]];
+  for(Integer iy=0; iy<size[0]; iy++)
+    (*paaX)[iy] = (*paaX)[0] + iy*size[1];
+  // The caller can access the contents of *paX using (*paaX)[i][j] notation.
+}
+
+
 template<typename Entry>
 void Alloc2D(int nrows,             //!< size of the array (outer)
              int ncolumns,          //!< size of the array (inner)
@@ -47,10 +65,11 @@ void Alloc2D(int nrows,             //!< size of the array (outer)
 template<typename Entry>
 void Dealloc2D(Entry ***paaX)        //!< pointer to 2-D multidimensional array
 {
+  //*paaX = new Entry* [size[0]];
+  //(*paaX)[0] = new Entry [size[0] * size[1]];
   if (paaX && *paaX) {
-    Entry *aX = &((*paaX)[0][0]);
+    delete [] (*paaX)[0];
     delete [] (*paaX);
-    delete [] aX;
     *paaX = nullptr;
   }
 }
diff --git a/tests/test.cpp b/tests/test.cpp
index fdc6b31..f501c41 100644
--- a/tests/test.cpp
+++ b/tests/test.cpp
@@ -209,10 +209,10 @@ void GenRandSymm(Matrix M, //<! store the matrix here
   GenRandOrth<Scalar, Matrix>(evects, n, generator);
 
   // Construct the test matrix, M, where M = Rt * D * R
-  mmult(evects, D, tmp, n);
-  for (int i = 0; i < n; i++)
-    for (int j = 0; j < n; j++)
-      evects[i][j] = evects[j][i]; //transpose "evects"
+  mmult(evects, D, tmp, n);  //tmp = Rt * D
+  for (int i = 0; i < n-1; i++)
+    for (int j = i+1; j < n; j++)
+      std::swap(evects[i][j], evects[j][i]); //transpose "evects"
   mmult(tmp, evects, M, n);  //at this point M = Rt * D * R (where "R"="evects")
   Dealloc2D(&D);
   Dealloc2D(&tmp);
@@ -304,7 +304,7 @@ void TestJacobi(int n, //<! matrix size
       //Check each eigenvector
       for (int i = 0; i < n; i++) {
         for (int a = 0; a < n; a++) {
-          test_evec[i] = 0.0;
+          test_evec[a] = 0.0;
           for (int b = 0; b < n; b++)
             test_evec[a] += M[a][b] * evects[i][b];
           assert(Similar(test_evec[a], evals[i] * evects[i][a], 1.0e-06));