diff --git a/.github/workflows/emulated.yml b/.github/workflows/emulated.yml index ff299179..dc15ef4b 100644 --- a/.github/workflows/emulated.yml +++ b/.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 diff --git a/.gitignore b/.gitignore index cad9087a..ccd95057 100644 --- a/.gitignore +++ b/.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 diff --git a/CHANGES b/CHANGES index 7d16caa4..632b12fc 100644 --- a/CHANGES +++ b/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 diff --git a/CMakeLists.txt b/CMakeLists.txt index 6cc4e25a..0ac2000b 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/EXAMPLES/MATRIX_MARKET/arpackSolver.hpp b/EXAMPLES/MATRIX_MARKET/arpackSolver.hpp index 9704a1cd..53b89c4c 100644 --- a/EXAMPLES/MATRIX_MARKET/arpackSolver.hpp +++ b/EXAMPLES/MATRIX_MARKET/arpackSolver.hpp @@ -176,6 +176,7 @@ class arpackSolver { int createMatrix(string const & fileName, Eigen::SparseMatrix & M) { // Read matrix from file. + if (fileName.empty()) { cerr << "Error: matrix file missing" << endl; return 1; } a_uint n = 0, m = 0; vector i, j; vector 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>() * V; EigVecZ right = stdPb ? V : B->template cast>() * 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; }; diff --git a/EXAMPLES/MATRIX_MARKET/arpackmm.cpp b/EXAMPLES/MATRIX_MARKET/arpackmm.cpp index 05d26645..72526c46 100644 --- a/EXAMPLES/MATRIX_MARKET/arpackmm.cpp +++ b/EXAMPLES/MATRIX_MARKET/arpackmm.cpp @@ -1,3 +1,6 @@ +#include +#include +#include #include #include "arpackSolver.hpp" @@ -9,8 +12,8 @@ using namespace std; class options { public: options() { - fileA = "A.mtx"; - fileB = "N.A."; // Not available. + fileA = ""; + fileB = ""; dense = false; denseRR = true; nbEV = 1; @@ -28,6 +31,7 @@ class options { invert = false; // Eigen value invertion: look for 1./lambda instead of lambda. tol = 1.e-06; + maxResNorm = 1.e-3; maxIt = 100; schur = false; // Compute Ritz vectors. slv = "BiCG"; @@ -44,27 +48,39 @@ class options { }; int readCmdLine(int argc, char** argv) { - // Check for command line independent parameters. + // Convert command line to stream. + stringstream ssIP; // Independent parameters. for (int a = 1; argv && a < argc; a++) { - string clo = argv[a]; // Command line option. - if (clo == "--help" || clo == "-h") return usage(0); + ssIP << argv[a] << " "; + } + + // Check for command line independent parameters. + + string clo; // Command line option. + while (ssIP >> clo) { + if (clo == "--help" || clo == "-h") return usage(); + if (clo.find("--") == std::string::npos) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + if (clo == "--A") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - fileA = argv[a]; + fileA = clo; } if (clo == "--dense") { dense = true; - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - string rr(argv[a]); + string rr = clo; if (rr != "true" && rr != "false") { cerr << "Error: bad " << clo << " - bad argument" << endl; return usage(); @@ -72,12 +88,12 @@ class options { denseRR = (rr == "true") ? true : false; } if (clo == "--nbEV") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream nEV(argv[a]); + stringstream nEV(clo); nEV >> nbEV; if (!nEV) { cerr << "Error: bad " << clo << " - bad argument" << endl; @@ -85,10 +101,6 @@ class options { } nbCV = 2 * nbEV + 1; } - if (clo == "--genPb") { - stdPb = false; - fileB = "B.mtx"; - } if (clo == "--nonSymPb") symPb = false; if (clo == "--cpxPb") { symPb = false; @@ -96,12 +108,12 @@ class options { } if (clo == "--simplePrec") simplePrec = true; if (clo == "--mag") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - mag = argv[a]; // small mag (likely poor perf) <=> large mag + invert + mag = clo; // small mag (likely poor perf) <=> large mag + invert // (likely good perf). bool ok = (mag == "LM" || mag == "SM" || mag == "LR" || mag == "SR" || mag == "LA" || mag == "SA" || mag == "LI" || mag == "SI") @@ -114,12 +126,12 @@ class options { } if (clo == "--shiftReal") { shiftReal = true; - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream s(argv[a]); + stringstream s(clo); s >> sigmaReal; if (!s) { cerr << "Error: bad " << clo << " - bad argument" << endl; @@ -128,12 +140,12 @@ class options { } if (clo == "--shiftImag") { shiftImag = true; - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream s(argv[a]); + stringstream s(clo); s >> sigmaImag; if (!s) { cerr << "Error: bad " << clo << " - bad argument" << endl; @@ -142,25 +154,38 @@ class options { } if (clo == "--invert") invert = true; if (clo == "--tol") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream t(argv[a]); + stringstream t(clo); t >> tol; if (!t) { cerr << "Error: bad " << clo << " - bad argument" << endl; return usage(); } } + if (clo == "--maxResNorm") { + ssIP >> clo; + if (!ssIP) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream mrn(clo); + mrn >> maxResNorm; + if (!mrn) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + } if (clo == "--maxIt") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream mi(argv[a]); + stringstream mi(clo); mi >> maxIt; if (!mi) { cerr << "Error: bad " << clo << " - bad argument" << endl; @@ -171,21 +196,20 @@ class options { schur = true; } if (clo == "--slv") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - slv = argv[a]; + slv = clo; } if (clo == "--slvItrTol") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream t(argv[a]); - double tol = 0.; + stringstream t(clo); t >> slvItrTol; if (!t) { cerr << "Error: bad " << clo << " - bad argument" << endl; @@ -193,12 +217,12 @@ class options { } } if (clo == "--slvItrMaxIt") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream mi(argv[a]); + stringstream mi(clo); int maxIt = 0; mi >> slvItrMaxIt; if (!mi) { @@ -207,12 +231,12 @@ class options { } } if (clo == "--slvItrPC") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream pc(argv[a]); + stringstream pc(clo); pc >> slvItrPC; if (!pc) { cerr << "Error: bad " << clo << " - bad argument" << endl; @@ -220,12 +244,12 @@ class options { } } if (clo == "--slvDrtPivot") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream pv(argv[a]); + stringstream pv(clo); pv >> slvDrtPivot; if (!pv) { cerr << "Error: bad " << clo << " - bad argument" << endl; @@ -233,12 +257,12 @@ class options { } } if (clo == "--slvDrtOffset") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream of(argv[a]); + stringstream of(clo); of >> slvDrtOffset; if (!of) { cerr << "Error: bad " << clo << " - bad argument" << endl; @@ -246,12 +270,12 @@ class options { } } if (clo == "--slvDrtScale") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream sc(argv[a]); + stringstream sc(clo); sc >> slvDrtScale; if (!sc) { cerr << "Error: bad " << clo << " - bad argument" << endl; @@ -260,12 +284,12 @@ class options { } if (clo == "--noCheck") check = false; if (clo == "--verbose") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream vb(argv[a]); + stringstream vb(clo); vb >> verbose; if (!vb) { cerr << "Error: bad " << clo << " - bad argument" << endl; @@ -273,12 +297,12 @@ class options { } } if (clo == "--debug") { - a++; - if (a >= argc) { + ssIP >> clo; + if (!ssIP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream dbg(argv[a]); + stringstream dbg(clo); dbg >> debug; if (!dbg) { cerr << "Error: bad " << clo << " - bad argument" << endl; @@ -290,19 +314,29 @@ class options { debug, debug, debug, debug, debug); } if (clo == "--restart") restart = true; + + // Skip dependent parameters. + if (clo == "--nbCV") { ssIP >> clo; continue; } + if (clo == "--B") { ssIP >> clo; continue; } } - // Check for command line dependent parameters. + // Convert command line to stream. + stringstream ssDP; // Dependent parameters. for (int a = 1; argv && a < argc; a++) { - string clo = argv[a]; // Command line option. + ssDP << argv[a] << " "; + } + + // Check for command line dependent parameters. + + while (ssDP >> clo) { if (clo == "--nbCV") { - a++; - if (a >= argc) { + ssDP >> clo; + if (!ssDP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - stringstream nCV(argv[a]); + stringstream nCV(clo); nCV >> nbCV; if (!nCV) { cerr << "Error: bad " << clo << " - bad argument" << endl; @@ -310,27 +344,20 @@ class options { } } if (clo == "--B") { - a++; - if (a >= argc) { + ssDP >> clo; + if (!ssDP) { cerr << "Error: bad " << clo << " - need argument" << endl; return usage(); } - fileB = argv[a]; + fileB = clo; + stdPb = false; // Generalized problem. } } - // Sanity checks. - - if (!stdPb && fileB.empty()) { - cerr << "Error: generalized problem without B matrix" - << ", specify --B XX with XX being a matrix market file" << endl; - return usage(); - } - return 0; }; - int usage(int rc = 1) { + int usage() { cout << "Usage: running arpack with matrix market files to check for eigen " "values/vectors." << endl; @@ -375,8 +402,6 @@ class options { cout << " default: 1" << endl; cout << " --nbCV: number of columns of the matrix V." << endl; cout << " default: 2*nbEV+1" << endl; - cout << " --genPb: generalized problem." << endl; - cout << " default: standard problem" << endl; cout << " --nonSymPb: non symmetric problem (<=> use dn[ae]upd)." << endl; cout << " default: symmetric problem (<=> use ds[ae]upd)" @@ -409,6 +434,8 @@ class options { cout << " default: no invert" << endl; cout << " --tol T: tolerance T." << endl; cout << " default: 1.e-06" << endl; + cout << " --maxResNorm R: maximum residual norm R." << endl; + cout << " default: 1.e-3" << endl; cout << " --maxIt M: maximum iterations M." << endl; cout << " default: 100" << endl; cout << " --schur: compute Schur vectors." << endl; @@ -501,8 +528,7 @@ class options { "computed during a previous run." << endl; cout << " default: false" << endl; - if (rc == 0) exit(0); - return rc; + return 1; }; friend ostream& operator<<(ostream& ostr, options const& opt); @@ -523,6 +549,7 @@ class options { // lambda+sigma instead of lambda. bool invert; // Eigen value invertion: look for 1./lambda instead of lambda. double tol; + double maxResNorm; int maxIt; bool schur; string slv; @@ -539,7 +566,8 @@ class options { }; ostream& operator<<(ostream& ostr, options const& opt) { - ostr << "OPT: A " << opt.fileA << ", B " << opt.fileB; + ostr << "OPT: A " << opt.fileA; + if (!opt.fileB.empty()) ostr << ", B " << opt.fileB; if (opt.dense && opt.denseRR) ostr << ", dense yes (RR true)"; else if (opt.dense && !opt.denseRR) @@ -557,7 +585,7 @@ ostream& operator<<(ostream& ostr, options const& opt) { ostr << ", shiftImag " << (opt.shiftImag ? "yes" : "no") << ", sigmaImag " << opt.sigmaImag; ostr << ", invert " << (opt.invert ? "yes" : "no") << ", tol " << opt.tol - << ", maxIt " << opt.maxIt; + << ", maxResNorm " << opt.maxResNorm << ", maxIt " << opt.maxIt; ostr << ", " << (opt.schur ? "Schur" : "Ritz") << " vectors" << endl; ostr << "OPT: slv " << opt.slv; bool itrSlv = true; // Use iterative solvers. @@ -670,7 +698,7 @@ int itrSolve(options& opt, output& out, double const& slvItrILUDropTol, return rc; } if (opt.check) { - rc = as.checkEigVec(A, opt.stdPb ? nullptr : &B); + rc = as.checkEigVec(A, opt.stdPb ? nullptr : &B, opt.maxResNorm); if (rc != 0) { cerr << "Error: check KO" << endl; return rc; @@ -757,7 +785,7 @@ int drtSolve(options& opt, output& out) { return rc; } if (opt.check) { - rc = as.checkEigVec(A, opt.stdPb ? nullptr : &B); + rc = as.checkEigVec(A, opt.stdPb ? nullptr : &B, opt.maxResNorm); if (rc != 0) { cerr << "Error: check KO" << endl; return rc; @@ -833,7 +861,7 @@ int itrSolve(options& opt, output& out) { return rc; } -int main(int argc, char** argv) { +int run(int argc, char** argv) { // Check for options. options opt; @@ -987,6 +1015,17 @@ int main(int argc, char** argv) { return 0; } +int main(int argc, char** argv) { + int rc = 1; + try { + rc = run(argc, argv); + } + catch (exception& e) { + cout << "Error - exception:" << e.what() << endl; + } + return rc; +} + // Local Variables: // mode: c++ // c-file-style:"stroustrup" diff --git a/EXAMPLES/MATRIX_MARKET/arpackmm.sh b/EXAMPLES/MATRIX_MARKET/arpackmm.sh index 8b338126..6291209a 100755 --- a/EXAMPLES/MATRIX_MARKET/arpackmm.sh +++ b/EXAMPLES/MATRIX_MARKET/arpackmm.sh @@ -1,90 +1,94 @@ #!/bin/bash -eu -export CMD="./arpackmm --help" # For coverage -echo "$CMD" -eval "$CMD &> arpackmm.run.log" -echo "" -echo "========================================================================================" -echo "" +trap 'catch' ERR +catch() { + grep OPT arpackmm.run.log + grep Error arpackmm.run.log + exit 1 +} -for eigPb in "--A As.mtx" "--nonSymPb --A An.mtx" "--nonSymPb --cpxPb --A Az.mtx --B Bz.mtx" +# For all these eigen problems, the first eigen value is about 382 or (382, 0.). + +for stdPb in "--A As.mtx" "--nonSymPb --A An.mtx" "--nonSymPb --cpxPb --A Az.mtx" do - for genPb in "" "--genPb" + # Choose B matrix according to A. + export fileB="" + if [[ "$stdPb" == *cpxPb* ]]; then + export fileB="--B Bz.mtx" + else + export fileB="--B B.mtx" + fi + + for genPb in "" "$fileB" do - for smallMag in "" "LA" "--mag SM --noCheck" # SM is known to be difficult to converge. + # SM may not converge so we can not apply checks. + # LM + invert is equivalent to SM and does converge. + # Note: this is expected as "power-like" methods are designed to find largest eigen values. + for magOpt in "--mag LM" "--mag SM --noCheck" do - export shiftOpt="" - if [[ "$eigPb" == *nonSymPb* ]]; then - if [[ "$genPb" == *genPb* ]]; then - continue # Skip to ensure stable test: tricky to convergence. - else - export shiftOpt="--shiftReal 100.0 --shiftImag 100.0" - fi + export mrn="--maxResNorm 1.e-1" # Relax residual check to get stable tests. + + # Shift slightly to avoid the zero-vector starting problem. + export shiftZV="" + if [[ "$stdPb" == *cpxPb* ]]; then + export shiftZV="--shiftReal 1.0 --shiftImag 1.0" + else + export shiftZV="--shiftReal 1.0" + fi + + # Shift according to the estimation of the eigen value we may have. + export shiftEV="" + if [[ "$stdPb" == *cpxPb* ]]; then + export shiftEV="--shiftReal 380.0 --shiftImag 1.0" else - if [[ "$genPb" == *genPb* ]]; then - continue # Skip to ensure stable test: tricky to convergence. - else - export shiftOpt="--shiftReal 100.0" - fi + export shiftEV="--shiftReal 380.0" fi - for shiftRI in "" "$shiftOpt" + for shiftRI in "$shiftZV" "$shiftEV" do for invert in "" "--invert" do - for tol in "" "--tol 1.e-5" + # Choose solver according to the sym/non-sym type of the problem. + export cgSlv="" + if [[ "$stdPb" == *nonSymPb* ]]; then + export cgSlv="BiCG" + else + export cgSlv="CG" + fi + + for slv in "--slv $cgSlv" "--slv LU" do - for slv in "--slv BiCG --slvItrTol 1.e-06 --slvItrMaxIt 150" "--slv CG --slvItrTol 1.e-06 --slvItrMaxIt 150" \ - "--slv BiCG --slvItrPC ILU" "--slv CG --slvItrPC ILU#1.e-06#2" \ - "--slv LU" "--slv QR --slvDrtPivot 1.e-06" \ - "--slv LLT" "--slv LLT --slvDrtOffset 0." \ - "--slv LDLT" "--slv LDLT --slvDrtScale 1." + for rs in "" "--schur" do - for rs in "" "--schur" + for dsPrec in "" "--simplePrec" do - for dsPrec in "" "--simplePrec" + for dsMat in "" "--dense true" do - for dsMat in "" "--dense false" "--dense true" - do - export extraGenPb="" - if [[ "$genPb" == *genPb* ]]; then - export extraGenPb="$shiftOpt" # Force shift if genPb. - fi - - if [[ "$slv" == *" CG "* ]]; then - if [[ "$eigPb" == *nonSymPb* ]]; then - continue # Skip CG that could fail (CG is meant to deal with sym matrices). - fi + # Skip not supported cases. + if [[ "$dsMat" == *dense* ]]; then + if [[ "$slv" == *CG* ]]; then + continue # Iterative solvers are not allowed when using dense matrices. fi + fi - if [[ "$slv" == *LLT* ]] || [[ "$slv" == *LDLT* ]]; then - if [[ "$eigPb" == *nonSymPb* ]] || [[ "$genPb" == *genPb* ]]; then - continue # Skip LLT/LDLT that could fail (LLT/LDLT are meant to deal with SPD matrices). - fi - fi - - if [[ "$dsMat" == *dense* ]]; then - if [[ "$slv" == *CG* ]]; then - continue # Iterative solvers are not allowed when using dense matrices. - fi - fi + export easeCV="--nbCV 6 --maxIt 200" # Use --nbCV 6 and --maxIt 200 to ease convergence. + echo "CLI: ./arpackmm $stdPb $genPb $magOpt $mrn $shiftRI $invert $slv $rs $dsPrec $dsMat $easeCV" + echo "----------------------------------------------------------------------------------------" - # Run arpackmm: use --nbCV 6 and --maxIt 200 to ease convergence, and, --verbose 3 for debug. - export CMD="./arpackmm $eigPb $genPb $smallMag $shiftRI $invert $tol $slv $rs $dsPrec $dsMat $extraGenPb --nbCV 6 --maxIt 200 --verbose 3 --debug 3" - echo "$CMD" - eval "$CMD &> arpackmm.run.log" - echo "" - echo "========================================================================================" - echo "" + # Run arpackmm. + ./arpackmm "$stdPb" "$genPb" "$magOpt" "$mrn" "$shiftRI" "$invert" "$slv" "$rs" "$dsPrec" "$dsMat" \ + "$easeCV" --verbose 3 &> arpackmm.run.log + grep OPT arpackmm.run.log + grep OUT arpackmm.run.log + echo "----------------------------------------------------------------------------------------" - # Run arpackmm: re-run with restart. - export CMD="$CMD --restart" - echo "$CMD" - eval "$CMD &> arpackmm.run.log" - echo "" - echo "========================================================================================" - echo "" - done + # Run arpackmm: re-run with restart always with small shift to avoid the zero-starting vector problem. + ./arpackmm "$stdPb" "$genPb" "$magOpt" "$mrn" "$shiftZV" "$invert" "$slv" "$rs" "$dsPrec" "$dsMat" \ + --restart \ + "$easeCV" --verbose 3 &> arpackmm.run.log + grep OPT arpackmm.run.log + grep OUT arpackmm.run.log + echo "========================================================================================" done done done @@ -95,4 +99,4 @@ do done done -echo "OK" +echo "arpackmm: OK" diff --git a/EXAMPLES/MATRIX_MARKET/issue401.sh b/EXAMPLES/MATRIX_MARKET/issue401.sh index d089c414..dab4b0ee 100755 --- a/EXAMPLES/MATRIX_MARKET/issue401.sh +++ b/EXAMPLES/MATRIX_MARKET/issue401.sh @@ -4,7 +4,7 @@ if [ "$?" -ne "0" ]; then exit 1; fi ./arpackmm --A issue401.mtx --mag LA --nbEV 1 --nbCV 5 --restart if [ "$?" -ne "0" ]; then exit 1; fi -./arpackmm --A issue401.mtx --mag LA --nbEV 1 --nbCV 5 --restart --shift 0.9 +./arpackmm --A issue401.mtx --mag LA --nbEV 1 --nbCV 5 --restart --shiftReal 0.9 if [ "$?" -ne "0" ]; then exit 1; fi echo "OK" diff --git a/EXAMPLES/PYARPACK/pyarpack.cpp b/EXAMPLES/PYARPACK/pyarpack.cpp index a705b2c7..d91c2c56 100644 --- a/EXAMPLES/PYARPACK/pyarpack.cpp +++ b/EXAMPLES/PYARPACK/pyarpack.cpp @@ -29,8 +29,9 @@ void exportArpackSparseItr(bp::scope& pySlv, std::string const& dtype) { .def("checkEigVec", &pyarpackSparseItrSolver::checkEigVec, (bp::arg("A"), bp::arg("B") = bp::tuple(), - bp::arg("diffTol") = 1.e-3), - "check eigen vectors accuracy where A and B must be sparse and " + bp::arg("maxResNorm") = 1.e-3), + "check eigen vectors accuracy (according to max allowed residual norm) " + "where A and B must be sparse and " "provided in coo format: (dimension, row-indice array, " "column-indice array, matrice-value array) tuple") ARPACKSOLVERMEMBER(pyarpackSparseItrSolver) @@ -69,8 +70,9 @@ void exportArpackSparseDrt(bp::scope& pySlv, std::string const& dtype) { .def("checkEigVec", &pyarpackSparseDrtSolver::checkEigVec, (bp::arg("A"), bp::arg("B") = bp::tuple(), - bp::arg("diffTol") = 1.e-3), - "check eigen vectors accuracy where A and B must be sparse and " + bp::arg("maxResNorm") = 1.e-3), + "check eigen vectors accuracy (according to max allowed residual norm) " + "where A and B must be sparse and " "provided in coo format: (dimension, row-indice array, " "column-indice array, matrice-value array) tuple") ARPACKSOLVERMEMBER(pyarpackSparseDrtSolver) @@ -103,8 +105,9 @@ void exportArpackDenseDrt(bp::scope& pySlv, std::string const& dtype) { .def("checkEigVec", &pyarpackDenseDrtSolver::checkEigVec, (bp::arg("A"), bp::arg("B") = bp::tuple(), - bp::arg("diffTol") = 1.e-3), - "check eigen vectors accuracy where A and B must be dense and " + bp::arg("maxResNorm") = 1.e-3), + "check eigen vectors accuracy (according to max allowed residual norm) " + "where A and B must be dense and " "provided in raw format: (n-squared matrice-value array, row or " "column ordered boolean)") ARPACKSOLVERMEMBER(pyarpackDenseDrtSolver) diff --git a/EXAMPLES/PYARPACK/pyarpackDrtSolver.hpp b/EXAMPLES/PYARPACK/pyarpackDrtSolver.hpp index 5428bc5b..ae632104 100644 --- a/EXAMPLES/PYARPACK/pyarpackDrtSolver.hpp +++ b/EXAMPLES/PYARPACK/pyarpackDrtSolver.hpp @@ -44,7 +44,7 @@ class pyarpackSparseDrtSolver: public arpackDrtSolver { return arpackDrtSolver::solve(M, (stdPb ? NULL : &N)); }; - int checkEigVec(bp::tuple const & A, bp::tuple const B = bp::tuple(), double const diffTol = 1.e-3) { + int checkEigVec(bp::tuple const & A, bp::tuple const B = bp::tuple(), double const maxResNorm = 1.e-3) { ARPACKSOLVERDEBUGSTAT(); EM M; int rc = pyarpackServices::buildSparseMatrice(A, M, debug, "A"); @@ -55,7 +55,7 @@ class pyarpackSparseDrtSolver: public arpackDrtSolver { rc = pyarpackServices::buildSparseMatrice(B, N, debug, "B"); if (rc != 0) {pyarpackThrowError("build matrice from B KO"); return rc;} } - return arpackDrtSolver::checkEigVec(M, (stdPb ? NULL : &N), &diffTol); + return arpackDrtSolver::checkEigVec(M, (stdPb ? NULL : &N), maxResNorm); }; // Public members. @@ -101,7 +101,7 @@ class pyarpackDenseDrtSolver: public arpackDrtSolver { return arpackDrtSolver::solve(M, (stdPb ? NULL : &N)); }; - int checkEigVec(bp::tuple const & A, bp::tuple const B = bp::tuple(), double const diffTol = 1.e-3) { + int checkEigVec(bp::tuple const & A, bp::tuple const B = bp::tuple(), double const maxResNorm = 1.e-3) { ARPACKSOLVERDEBUGSTAT(); EM M; int rc = pyarpackServices::buildDenseMatrice(A, M, debug, "A"); @@ -112,7 +112,7 @@ class pyarpackDenseDrtSolver: public arpackDrtSolver { rc = pyarpackServices::buildDenseMatrice(B, N, debug, "B"); if (rc != 0) {pyarpackThrowError("build matrice from B KO"); return rc;} } - return arpackDrtSolver::checkEigVec(M, (stdPb ? NULL : &N), &diffTol); + return arpackDrtSolver::checkEigVec(M, (stdPb ? NULL : &N), maxResNorm); }; // Public members. diff --git a/EXAMPLES/PYARPACK/pyarpackItrSolver.hpp b/EXAMPLES/PYARPACK/pyarpackItrSolver.hpp index 82c5ca87..e94f9d2a 100644 --- a/EXAMPLES/PYARPACK/pyarpackItrSolver.hpp +++ b/EXAMPLES/PYARPACK/pyarpackItrSolver.hpp @@ -43,7 +43,7 @@ class pyarpackSparseItrSolver: public arpackItrSolver { return arpackItrSolver::solve(M, (stdPb ? NULL : &N)); }; - int checkEigVec(bp::tuple const & A, bp::tuple const B = bp::tuple(), double const diffTol = 1.e-3) { + int checkEigVec(bp::tuple const & A, bp::tuple const B = bp::tuple(), double const maxResNorm = 1.e-3) { ARPACKSOLVERDEBUGSTAT(); EM M; int rc = pyarpackServices::buildSparseMatrice(A, M, debug, "A"); @@ -54,7 +54,7 @@ class pyarpackSparseItrSolver: public arpackItrSolver { rc = pyarpackServices::buildSparseMatrice(B, N, debug, "B"); if (rc != 0) {pyarpackThrowError("build matrice from B KO"); return rc;} } - return arpackItrSolver::checkEigVec(M, (stdPb ? NULL : &N), &diffTol); + return arpackItrSolver::checkEigVec(M, (stdPb ? NULL : &N), maxResNorm); }; // Public members. diff --git a/cmake/tstCMakeInstall.sh.in b/cmake/tstCMakeInstall.sh.in index f22a1f24..4e13a9bf 100755 --- a/cmake/tstCMakeInstall.sh.in +++ b/cmake/tstCMakeInstall.sh.in @@ -211,7 +211,7 @@ export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:/tmp/tstCMakeInstall/local/lib" ./dnbdr1 || exit 1 ./icb_arpack_c || exit 1 ./icb_arpack_cpp || exit 1 -./arpackmm -h || exit 1 +./arpackmm -h || true # Returns 1 but not relevant: we want to check the installed arpackSolver can be used. mpirun -n 2 ./pdndrv1 || exit 1 mpirun -n 2 ./icb_parpack_c || exit 1 mpirun -n 2 ./icb_parpack_cpp || exit 1 diff --git a/pkg-config/tstAutotoolsInstall.sh.in b/pkg-config/tstAutotoolsInstall.sh.in index 31d16e69..043b2810 100755 --- a/pkg-config/tstAutotoolsInstall.sh.in +++ b/pkg-config/tstAutotoolsInstall.sh.in @@ -145,7 +145,7 @@ export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:/tmp/tstAutotoolsInstall/local/lib" ./dnbdr1 || exit 1 ./icb_arpack_c || exit 1 ./icb_arpack_cpp || exit 1 -./arpackmm -h || exit 1 +./arpackmm -h || true # Returns 1 but not relevant: we want to check the installed arpackSolver can be used. mpirun -n 2 ./pdndrv1 || exit 1 mpirun -n 2 ./icb_parpack_c || exit 1