Fix arpack solver (#468)

* arpackmm: improve command line.

* [BUG FIX] arpackmm: run arpack with real or imag shift.

* [BUG FIX] arpackmm: make sure the restart file is written.

* [BUG FIX] arpackmm: rewrite the test.

  - Do NOT use eval in bash script which does NOT return command exit code.

  - Use a dedicated maxResNorm to test the residual (no point to deduce
    it from tol).

* Update .gitignore.

* Update CHANGELOG.

.github/workflows/emulated.yml

@@ -1,4 +1,4 @@
-name: arch-emu
+name: arpack-ng-emu
 on:
   workflow_dispatch:
   push:
@@ -110,4 +110,4 @@ jobs:
       - name: test
         run: |
           cd ${GITHUB_WORKSPACE}/build
-          ctest . || ctest . --rerun-failed --output-on-failure
+          CTEST_OUTPUT_ON_FAILURE=1 make all test

.gitignore

@@ -20,6 +20,8 @@ parpack*.pc
 arpackdef.h
 arpackicb.h
 tstAutotoolsInstall.sh
+cmake/arpackng-config-version.cmake
+cmake/arpackng-config.cmake
 
 # Generated by `make`
 .dirstamp

CHANGES

@@ -1,3 +1,8 @@
+[ Franck Houssen ]
+ * [BUG FIX] Fix arpack solver to handle real and imag shifts.
+ * [BUG FIX] Make sure the restart file of the arpack solver i written.
+ * [BUG FIX] Change arpack solver API.
+
 [ Henri Menke ]
  * [BUG FIX] Add missing stdexcept header
 

CMakeLists.txt

@@ -704,13 +704,13 @@ function(build_tests)
       configure_file(EXAMPLES/MATRIX_MARKET/B.mtx  ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/B.mtx)
       configure_file(EXAMPLES/MATRIX_MARKET/Bz.mtx ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/Bz.mtx)
       configure_file(EXAMPLES/MATRIX_MARKET/arpackmm.sh ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/arpackmm.sh)
-      add_test(NAME arpackmm_tst WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} COMMAND ${BASH_PROGRAM} -eu arpackmm.sh)
+      add_test(NAME arpackmm_tst WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} COMMAND ${BASH_PROGRAM} arpackmm.sh)
       configure_file(EXAMPLES/MATRIX_MARKET/issue401.mtx ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/issue401.mtx)
       configure_file(EXAMPLES/MATRIX_MARKET/issue401.sh ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/issue401.sh)
-      add_test(NAME issue401_tst WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} COMMAND ${BASH_PROGRAM} -eu issue401.sh)
+      add_test(NAME issue401_tst WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} COMMAND ${BASH_PROGRAM} issue401.sh)
       configure_file(EXAMPLES/MATRIX_MARKET/issue215.mtx ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/issue215.mtx)
       configure_file(EXAMPLES/MATRIX_MARKET/issue215.sh ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/issue215.sh)
-      add_test(NAME issue215_tst WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} COMMAND ${BASH_PROGRAM} -eu issue215.sh)
+      add_test(NAME issue215_tst WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} COMMAND ${BASH_PROGRAM} issue215.sh)
     endif()
 
     if (PYTHON3)

EXAMPLES/MATRIX_MARKET/arpackSolver.hpp

@@ -176,6 +176,7 @@ class arpackSolver {
     int createMatrix(string const & fileName, Eigen::SparseMatrix<RC> & M) {
       // Read matrix from file.
 
+      if (fileName.empty()) { cerr << "Error: matrix file missing" << endl; return 1; }
       a_uint n = 0, m = 0;
       vector<a_uint> i, j;
       vector<RC> Mij;
@@ -255,7 +256,8 @@ class arpackSolver {
       mode = 0;
       if (stdPb) {
         mode = 1;
-        if (shiftReal && !shiftImag) {
+        // CAUTION: back transform must be done only if mode = 1.
+        if (shiftReal || shiftImag) {
           EM I(A.rows(), A.cols());
           I.setIdentity();
           RC sigma; makeSigma(sigma);
@@ -264,6 +266,7 @@ class arpackSolver {
         }
       }
       else {
+        // CAUTION: back transform must NOT be done if mode > 1.
         mode = 2;
         if (shiftReal || shiftImag) mode = 3;
       }
@@ -294,10 +297,7 @@ class arpackSolver {
       return 0;
     };
 
-    int checkEigVec(EM const & A, EM const * B = nullptr, double const * diffTol = nullptr) {
-      stdPb = !B ? true : false;
-      double dTol = !diffTol ? sqrt(tol) : *diffTol;
-
+    int checkEigVec(EM const & A, EM const * B = nullptr, double const maxResNorm = 1.e-3) const {
       // Check eigen vectors.
 
       string rs = schur ? "Schur" : "Ritz";
@@ -322,23 +322,25 @@ class arpackSolver {
         EigVecZ left = A.template cast<complex<double>>() * V;
         EigVecZ right = stdPb ? V : B->template cast<complex<double>>() * V;
         right *= lambda;
-        EigVecZ diff = left - right;
-        if (diff.norm() > dTol) {
-          cerr << endl << "Error: bad vector " << setw(3) << i << " (norm " << V.norm() << "):" << endl;
+        EigVecZ residual = left - right;
+        if (residual.norm() > maxResNorm) {
+          cerr << endl << "Error: bad eigen value " << i << " (norm " << std::norm(lambda) << "):" << endl;
+          cerr << endl << lambda << endl;
+          cerr << endl << "Error: bad eigen vector " << i << " (norm " << V.norm() << "):" << endl;
           cerr << endl << V << endl;
-          cerr << endl << "Error: left side (A*V - norm " << left.norm() << "):" << endl;
+          cerr << endl << "Error: left side (A*V has norm " << left.norm() << "):" << endl;
           cerr << endl << left << endl;
-          cerr << endl << "Error: right side (lambda*" << (stdPb ? "" : "B*") << "V - norm " << right.norm() << "):" << endl;
+          cerr << endl << "Error: right side (lambda*" << (stdPb ? "" : "B*") << "V has norm " << right.norm() << "):" << endl;
           cerr << endl << right << endl;
-          cerr << endl << "Error: diff (norm " << diff.norm() << ", tol " << dTol << "):" << endl;
-          cerr << endl << diff << endl;
+          cerr << endl << "Error: residual (norm " << residual.norm() << ", maxResNorm " << maxResNorm << "):" << endl;
+          cerr << endl << residual << endl;
           return 1;
         }
         else {
           if (verbose >= 1) {
             cout << endl << "arpackSolver:" << endl;
             cout << endl << rs << " value/vector " << setw(3) << i << ": check OK";
-            cout << ", diff (norm " << diff.norm() << ", tol " << dTol << ")" << endl;
+            cout << ", residual (norm " << residual.norm() << ", maxResNorm " << maxResNorm << ")" << endl;
           }
         }
       }
@@ -667,6 +669,7 @@ class arpackSolver {
       if (ofs.is_open()) {
         ofs << nbDim << endl;
         for (a_int n = 0; rv && n < nbDim; n++) ofs << rv[n] << endl;
+        ofs.close(); // Make sure the file is written.
       }
       return 0;
     };