Skip to content

Commit

Permalink
arpackmm: improve command line for testing. (#470)
Browse files Browse the repository at this point in the history
  • Loading branch information
@fghoussen
fghoussen authored Sep 19, 2024
1 parent 43d3e56 commit f2921f2
Show file tree
Hide file tree
Showing 2 changed files with 66 additions and 23 deletions.
66 changes: 43 additions & 23 deletions EXAMPLES/MATRIX_MARKET/arpackSolver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -319,35 +319,41 @@ 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;
}
}
}

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:
Expand All @@ -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.

Expand All @@ -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);
Expand All @@ -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;
};

Expand Down
23 changes: 23 additions & 0 deletions EXAMPLES/MATRIX_MARKET/arpackmm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
#include <exception>
#include <sstream>
#include <chrono>
#include <vector>
#include <tuple>

#include "arpackSolver.hpp"
#include "debug_c.hpp"
Expand Down Expand Up @@ -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>
Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -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;
}

Expand Down

0 comments on commit f2921f2

Please sign in to comment.