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

arpackmm: improve command line for testing. (#470)

Browse Source
This commit is contained in:
Franck HOUSSEN
2024-09-19 18:58:07 +02:00
committed by GitHub
parent 43d3e5606b
commit f2921f232b
2 changed files with 66 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files
EXAMPLES/MATRIX_MARKET/arpackSolver.hpp
+43 -23
View File
@@ -319,28 +319,18 @@ class arpackSolver {
}
}
EigVecZ left = A.template cast<complex<double>>() * V;
EigVecZ right = stdPb ? V : B->template cast<complex<double>>() * V;
right *= lambda;
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 has norm " << left.norm() << "):" << endl;
cerr << endl << left << endl;
cerr << endl << "Error: right side (lambda*" << (stdPb ? "" : "B*") << "V has norm " << right.norm() << "):" << endl;
cerr << endl << right << endl;
cerr << endl << "Error: residual (norm " << residual.norm() << ", maxResNorm " << maxResNorm << "):" << endl;
cerr << endl << residual << endl;
double resNorm = computeResidualNorm(i, A, B);
if (resNorm > maxResNorm) {
cout << endl << "arpackSolver:" << endl;
cout << endl << rs << " value/vector " << setw(3) << i << ": check KO";
cout << ", residual (norm " << resNorm << ", maxResNorm " << maxResNorm << ")" << endl;
return 1;
}
else {
if (verbose >= 1) {
cout << endl << "arpackSolver:" << endl;
cout << endl << rs << " value/vector " << setw(3) << i << ": check OK";
cout << ", residual (norm " << residual.norm() << ", maxResNorm " << maxResNorm << ")" << endl;
cout << ", residual (norm " << resNorm << ", maxResNorm " << maxResNorm << ")" << endl;
}
}
}
@@ -348,6 +338,22 @@ class arpackSolver {
return 0;
};
double computeResidualNorm(size_t const & idx, EM const & A, EM const * B = nullptr) const {
// Check eigen value index.
if (idx >= val.size()) return -1.; // Error.
if (idx >= vec.size()) return -1.; // Error.
// Compute residual norm.
EigVecZ V = vec[idx];
complex<double> lambda = val[idx];
EigVecZ left = A.template cast<complex<double>>() * V;
EigVecZ right = stdPb ? V : B->template cast<complex<double>>() * V;
right *= lambda;
EigVecZ residual = left - right;
return residual.norm();
};
// Private methods.
private:
@@ -374,7 +380,7 @@ class arpackSolver {
string inpLine; getline(inp, inpLine); l++;
while (isspace(*inpLine.begin())) inpLine.erase(inpLine.begin()); // Suppress leading white spaces.
if (inpLine.length() == 0) continue; // Empty line.
if (inpLine[0] == '%') continue; // Comments skipped, begin reading.
if (inpLine[0] == '%' || inpLine[0] == '#') continue; // Comments skipped, begin reading.
// Read matrix market file.
@@ -384,7 +390,13 @@ class arpackSolver {
if (!inpSS) {cerr << "Error: bad header (n, m)" << endl; return 1;}
if (nnz == 0) {
inpSS >> nnz;
if (inpSS) { // OK, (optional) nnz has been provided.
if (inpSS && nnz > 0) { // OK, (optional) nnz has been provided.
i.reserve(nnz);
j.reserve(nnz);
Mij.reserve(nnz);
}
else {
nnz = n*m;
i.reserve(nnz);
j.reserve(nnz);
Mij.reserve(nnz);
@@ -406,14 +418,22 @@ class arpackSolver {
// Handle 1-based -> 0-based.
nnz = i.size(); // In case nnz was not provided.
if (nnz > 0) {
if (*max_element(begin(i), end(i)) == n || *max_element(begin(j), end(j)) == m) {
for (size_t k = 0; k < nnz; k++) i[k] -= 1;
for (size_t k = 0; k < nnz; k++) j[k] -= 1;
if (!i.empty() && !j.empty()) {
if (*max_element(begin(i), end(i)) == n && *max_element(begin(j), end(j)) == m) {
for (size_t k = 0; k < i.size(); k++) { if (i[k] > 0) i[k] -= 1; }
for (size_t k = 0; k < j.size(); k++) { if (j[k] > 0) j[k] -= 1; }
}
}
// Checking indices.
for (size_t k = 0; k < i.size(); k++) {
if (i[k] >= n) {cerr << "Error: bad index (" << fileName << ", i " << i[k] << ")" << endl; return 1;};
}
for (size_t k = 0; k < j.size(); k++) {
if (j[k] >= m) {cerr << "Error: bad index (" << fileName << ", j " << j[k] << ")" << endl; return 1;};
}
return 0;
};
EXAMPLES/MATRIX_MARKET/arpackmm.cpp
+23
View File
@@ -2,6 +2,8 @@
#include <exception>
#include <sstream>
#include <chrono>
#include <vector>
#include <tuple>
#include "arpackSolver.hpp"
#include "debug_c.hpp"
@@ -625,6 +627,7 @@ class output {
int nbIt; // Arpack number of iterations.
double imsTime; // Init mode solver time.
double rciTime; // Reverse communication interface time.
vector<tuple<double, double>> res; // Results: eigen values and associated residual.
};
template <typename RC, typename FD, typename EM, typename SLV>
@@ -712,6 +715,11 @@ int itrSolve(options& opt, output& out, double const& slvItrILUDropTol,
out.nbIt = as.nbIt;
out.imsTime = as.imsTime;
out.rciTime = as.rciTime;
for (auto idx = 0; idx < as.val.size(); idx++) {
double val = norm(as.val[idx]);
double res = as.computeResidualNorm(idx, A, opt.stdPb ? nullptr : &B);
out.res.push_back(tuple<double, double>{val, res});
}
return 0;
}
@@ -799,6 +807,11 @@ int drtSolve(options& opt, output& out) {
out.nbIt = as.nbIt;
out.imsTime = as.imsTime;
out.rciTime = as.rciTime;
for (auto idx = 0; idx < as.val.size(); idx++) {
double val = norm(as.val[idx]);
double res = as.computeResidualNorm(idx, A, opt.stdPb ? nullptr : &B);
out.res.push_back(tuple<double, double>{val, res});
}
return 0;
}
@@ -1012,6 +1025,16 @@ int run(int argc, char** argv) {
cout << "STAT: total number of restart steps "
<< nrstrt << endl;
// Output eigen values and residuals.
cout << endl;
for (auto idx = 0; idx < out.res.size(); idx++) {
tuple<double, double> res = out.res[idx];
cout << "RES: eigen value " << idx << ": norm " << get<0>(res);
cout << " (residual norm " << get<1>(res) << ")";
cout << endl;
}
return 0;
}