phynexis-deps/arpack-ng
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fix arpack solver (#468)

Browse Source 
* 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.
This commit is contained in:
Franck HOUSSEN
2024-09-18 21:55:38 +02:00
committed by GitHub
parent e55d3b2dcb
commit 43d3e5606b
13 changed files with 245 additions and 189 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.github/workflows/emulated.yml
+2 -2
View File
@@ -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
+2
View File
@@ -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
+5
View File
@@ -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
+3 -3
View File
@@ -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
+16 -13
View File
@@ -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;
};
EXAMPLES/MATRIX_MARKET/arpackmm.cpp
+127 -88
View File
@@ -1,3 +1,6 @@
#include <iostream>
#include <exception>
#include <sstream>
#include <chrono>
#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) {
// Convert command line to stream.
stringstream ssIP; // Independent parameters.
for (int a = 1; argv && a < argc; a++) {
ssIP << argv[a] << " ";
}
// Check for command line 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);
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 == "--maxIt") {
a++;
if (a >= argc) {
if (clo == "--maxResNorm") {
ssIP >> clo;
if (!ssIP) {
cerr << "Error: bad " << clo << " - need argument" << endl;
return usage();
}
stringstream mi(argv[a]);
stringstream mrn(clo);
mrn >> maxResNorm;
if (!mrn) {
cerr << "Error: bad " << clo << " - bad argument" << endl;
return usage();
}
}
if (clo == "--maxIt") {
ssIP >> clo;
if (!ssIP) {
cerr << "Error: bad " << clo << " - need argument" << endl;
return usage();
}
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; }
}
// Convert command line to stream.
stringstream ssDP; // Dependent parameters.
for (int a = 1; argv && a < argc; a++) {
ssDP << argv[a] << " ";
}
// Check for command line dependent parameters.
for (int a = 1; argv && a < argc; a++) {
string clo = argv[a]; // Command line option.
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"
EXAMPLES/MATRIX_MARKET/arpackmm.sh
+72 -68
View File
@@ -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
if [[ "$genPb" == *genPb* ]]; then
continue # Skip to ensure stable test: tricky to convergence.
else
export shiftOpt="--shiftReal 100.0"
fi
export shiftZV="--shiftReal 1.0"
fi
for shiftRI in "" "$shiftOpt"
# 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
export shiftEV="--shiftReal 380.0"
fi
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.
# 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" == *" CG "* ]]; then
if [[ "$eigPb" == *nonSymPb* ]]; then
continue # Skip CG that could fail (CG is meant to deal with sym 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 "----------------------------------------------------------------------------------------"
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
# 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 "----------------------------------------------------------------------------------------"
if [[ "$dsMat" == *dense* ]]; then
if [[ "$slv" == *CG* ]]; then
continue # Iterative solvers are not allowed when using dense matrices.
fi
fi
# 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: 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"
EXAMPLES/MATRIX_MARKET/issue401.sh
+1 -1
View File
@@ -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"
EXAMPLES/PYARPACK/pyarpack.cpp
+9 -6
View File
@@ -29,8 +29,9 @@ void exportArpackSparseItr(bp::scope& pySlv, std::string const& dtype) {
.def("checkEigVec",
&pyarpackSparseItrSolver<RC, FD, EM, SLV>::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<RC, FD, EM, SLV>::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<RC, FD, EM, SLV>::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)
EXAMPLES/PYARPACK/pyarpackDrtSolver.hpp
+4 -4
View File
@@ -44,7 +44,7 @@ class pyarpackSparseDrtSolver: public arpackDrtSolver<RC, FD, EM, SLV> {
return arpackDrtSolver<RC, FD, EM, SLV>::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<RC, EM>::buildSparseMatrice(A, M, debug, "A");
@@ -55,7 +55,7 @@ class pyarpackSparseDrtSolver: public arpackDrtSolver<RC, FD, EM, SLV> {
rc = pyarpackServices<RC, EM>::buildSparseMatrice(B, N, debug, "B");
if (rc != 0) {pyarpackThrowError("build matrice from B KO"); return rc;}
}
return arpackDrtSolver<RC, FD, EM, SLV>::checkEigVec(M, (stdPb ? NULL : &N), &diffTol);
return arpackDrtSolver<RC, FD, EM, SLV>::checkEigVec(M, (stdPb ? NULL : &N), maxResNorm);
};
// Public members.
@@ -101,7 +101,7 @@ class pyarpackDenseDrtSolver: public arpackDrtSolver<RC, FD, EM, SLV> {
return arpackDrtSolver<RC, FD, EM, SLV>::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<RC, EM>::buildDenseMatrice(A, M, debug, "A");
@@ -112,7 +112,7 @@ class pyarpackDenseDrtSolver: public arpackDrtSolver<RC, FD, EM, SLV> {
rc = pyarpackServices<RC, EM>::buildDenseMatrice(B, N, debug, "B");
if (rc != 0) {pyarpackThrowError("build matrice from B KO"); return rc;}
}
return arpackDrtSolver<RC, FD, EM, SLV>::checkEigVec(M, (stdPb ? NULL : &N), &diffTol);
return arpackDrtSolver<RC, FD, EM, SLV>::checkEigVec(M, (stdPb ? NULL : &N), maxResNorm);
};
// Public members.
EXAMPLES/PYARPACK/pyarpackItrSolver.hpp
+2 -2
View File
@@ -43,7 +43,7 @@ class pyarpackSparseItrSolver: public arpackItrSolver<RC, FD, EM, SLV> {
return arpackItrSolver<RC, FD, EM, SLV>::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<RC, EM>::buildSparseMatrice(A, M, debug, "A");
@@ -54,7 +54,7 @@ class pyarpackSparseItrSolver: public arpackItrSolver<RC, FD, EM, SLV> {
rc = pyarpackServices<RC, EM>::buildSparseMatrice(B, N, debug, "B");
if (rc != 0) {pyarpackThrowError("build matrice from B KO"); return rc;}
}
return arpackItrSolver<RC, FD, EM, SLV>::checkEigVec(M, (stdPb ? NULL : &N), &diffTol);
return arpackItrSolver<RC, FD, EM, SLV>::checkEigVec(M, (stdPb ? NULL : &N), maxResNorm);
};
// Public members.
cmake/tstCMakeInstall.sh.in
+1 -1
View File
@@ -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
pkg-config/tstAutotoolsInstall.sh.in
+1 -1
View File
@@ -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