1188 lines
45 KiB
C++
1188 lines
45 KiB
C++
#ifndef __ARPACKSOLVER_HPP__
|
|
#define __ARPACKSOLVER_HPP__
|
|
|
|
#include "arpack.h"
|
|
|
|
#include <iostream>
|
|
#include <string>
|
|
#include <sstream> // stringstream.
|
|
#include <fstream> // [io]fstream.
|
|
#include <chrono>
|
|
#include <iomanip> // setw.
|
|
#include <cassert>
|
|
#include <vector>
|
|
#include <type_traits> // is_same.
|
|
#include <cmath> // abs
|
|
#include <complex>
|
|
#include <limits> // epsilon
|
|
|
|
#include <Eigen/Sparse>
|
|
#include <Eigen/IterativeLinearSolvers>
|
|
#include <Eigen/SparseLU>
|
|
#include <Eigen/SparseQR>
|
|
#include <Eigen/SparseCholesky>
|
|
#include <Eigen/Dense>
|
|
#include <Eigen/LU>
|
|
#include <Eigen/QR>
|
|
#include <Eigen/Cholesky>
|
|
|
|
using namespace std;
|
|
|
|
// Sparse matrix related types.
|
|
|
|
typedef Eigen::SparseMatrix< float> EigSMxS; // Real.
|
|
typedef Eigen::SparseMatrix< double> EigSMxD; // Real.
|
|
typedef Eigen::SparseMatrix<complex< float>> EigSMxC; // Complex.
|
|
typedef Eigen::SparseMatrix<complex<double>> EigSMxZ; // Complex.
|
|
|
|
// Iterative solvers for sparse matrices.
|
|
|
|
typedef Eigen::BiCGSTAB <EigSMxS> EigSBiCGS; // Real.
|
|
typedef Eigen::BiCGSTAB <EigSMxD> EigSBiCGD; // Real.
|
|
typedef Eigen::BiCGSTAB <EigSMxC> EigSBiCGC; // Complex.
|
|
typedef Eigen::BiCGSTAB <EigSMxZ> EigSBiCGZ; // Complex.
|
|
|
|
typedef Eigen::ConjugateGradient<EigSMxS> EigSCGS; // Real.
|
|
typedef Eigen::ConjugateGradient<EigSMxD> EigSCGD; // Real.
|
|
typedef Eigen::ConjugateGradient<EigSMxC> EigSCGC; // Complex.
|
|
typedef Eigen::ConjugateGradient<EigSMxZ> EigSCGZ; // Complex.
|
|
|
|
typedef Eigen::IncompleteLUT< float> EigILUS; // Real.
|
|
typedef Eigen::IncompleteLUT< double> EigILUD; // Real.
|
|
typedef Eigen::IncompleteLUT<complex< float>> EigILUC; // Complex.
|
|
typedef Eigen::IncompleteLUT<complex<double>> EigILUZ; // Complex.
|
|
|
|
typedef Eigen::BiCGSTAB <EigSMxS, EigILUS> EigSBiCGILUS; // Real.
|
|
typedef Eigen::BiCGSTAB <EigSMxD, EigILUD> EigSBiCGILUD; // Real.
|
|
typedef Eigen::BiCGSTAB <EigSMxC, EigILUC> EigSBiCGILUC; // Complex.
|
|
typedef Eigen::BiCGSTAB <EigSMxZ, EigILUZ> EigSBiCGILUZ; // Complex.
|
|
|
|
typedef Eigen::ConjugateGradient<EigSMxS, Eigen::Lower|Eigen::Upper, EigILUS> EigSCGILUS; // Real.
|
|
typedef Eigen::ConjugateGradient<EigSMxD, Eigen::Lower|Eigen::Upper, EigILUD> EigSCGILUD; // Real.
|
|
typedef Eigen::ConjugateGradient<EigSMxC, Eigen::Lower|Eigen::Upper, EigILUC> EigSCGILUC; // Complex.
|
|
typedef Eigen::ConjugateGradient<EigSMxZ, Eigen::Lower|Eigen::Upper, EigILUZ> EigSCGILUZ; // Complex.
|
|
|
|
// Direct solvers for sparse matrices.
|
|
|
|
typedef Eigen::SimplicialLLT <EigSMxS, Eigen::Lower, Eigen::COLAMDOrdering<int>> EigSLLTS; // Real.
|
|
typedef Eigen::SimplicialLLT <EigSMxD, Eigen::Lower, Eigen::COLAMDOrdering<int>> EigSLLTD; // Real.
|
|
typedef Eigen::SimplicialLLT <EigSMxC, Eigen::Lower, Eigen::COLAMDOrdering<int>> EigSLLTC; // Complex.
|
|
typedef Eigen::SimplicialLLT <EigSMxZ, Eigen::Lower, Eigen::COLAMDOrdering<int>> EigSLLTZ; // Complex.
|
|
|
|
typedef Eigen::SimplicialLDLT<EigSMxS, Eigen::Lower, Eigen::COLAMDOrdering<int>> EigSLDLTS; // Real.
|
|
typedef Eigen::SimplicialLDLT<EigSMxD, Eigen::Lower, Eigen::COLAMDOrdering<int>> EigSLDLTD; // Real.
|
|
typedef Eigen::SimplicialLDLT<EigSMxC, Eigen::Lower, Eigen::COLAMDOrdering<int>> EigSLDLTC; // Complex.
|
|
typedef Eigen::SimplicialLDLT<EigSMxZ, Eigen::Lower, Eigen::COLAMDOrdering<int>> EigSLDLTZ; // Complex.
|
|
|
|
typedef Eigen::SparseLU<EigSMxS, Eigen::COLAMDOrdering<int>> EigSLUS; // Real.
|
|
typedef Eigen::SparseLU<EigSMxD, Eigen::COLAMDOrdering<int>> EigSLUD; // Real.
|
|
typedef Eigen::SparseLU<EigSMxC, Eigen::COLAMDOrdering<int>> EigSLUC; // Complex.
|
|
typedef Eigen::SparseLU<EigSMxZ, Eigen::COLAMDOrdering<int>> EigSLUZ; // Complex.
|
|
|
|
typedef Eigen::SparseQR<EigSMxS, Eigen::COLAMDOrdering<int>> EigSQRS; // Real.
|
|
typedef Eigen::SparseQR<EigSMxD, Eigen::COLAMDOrdering<int>> EigSQRD; // Real.
|
|
typedef Eigen::SparseQR<EigSMxC, Eigen::COLAMDOrdering<int>> EigSQRC; // Complex.
|
|
typedef Eigen::SparseQR<EigSMxZ, Eigen::COLAMDOrdering<int>> EigSQRZ; // Complex.
|
|
|
|
// Dense matrix related types.
|
|
|
|
typedef Eigen::Matrix< float, Eigen::Dynamic, Eigen::Dynamic> EigDMxS; // Real.
|
|
typedef Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic> EigDMxD; // Real.
|
|
typedef Eigen::Matrix<complex< float>, Eigen::Dynamic, Eigen::Dynamic> EigDMxC; // Complex.
|
|
typedef Eigen::Matrix<complex<double>, Eigen::Dynamic, Eigen::Dynamic> EigDMxZ; // Complex.
|
|
|
|
// Direct solvers for dense matrices.
|
|
|
|
typedef Eigen::LLT<EigDMxS> EigDLLTS; // Real.
|
|
typedef Eigen::LLT<EigDMxD> EigDLLTD; // Real.
|
|
typedef Eigen::LLT<EigDMxC> EigDLLTC; // Complex.
|
|
typedef Eigen::LLT<EigDMxZ> EigDLLTZ; // Complex.
|
|
|
|
typedef Eigen::LDLT<EigDMxS> EigDLDLTS; // Real.
|
|
typedef Eigen::LDLT<EigDMxD> EigDLDLTD; // Real.
|
|
typedef Eigen::LDLT<EigDMxC> EigDLDLTC; // Complex.
|
|
typedef Eigen::LDLT<EigDMxZ> EigDLDLTZ; // Complex.
|
|
|
|
typedef Eigen::FullPivLU<EigDMxS> EigDFLUS; // Real.
|
|
typedef Eigen::FullPivLU<EigDMxD> EigDFLUD; // Real.
|
|
typedef Eigen::FullPivLU<EigDMxC> EigDFLUC; // Complex.
|
|
typedef Eigen::FullPivLU<EigDMxZ> EigDFLUZ; // Complex.
|
|
|
|
typedef Eigen::FullPivHouseholderQR<EigDMxS> EigDFQRS; // Real.
|
|
typedef Eigen::FullPivHouseholderQR<EigDMxD> EigDFQRD; // Real.
|
|
typedef Eigen::FullPivHouseholderQR<EigDMxC> EigDFQRC; // Complex.
|
|
typedef Eigen::FullPivHouseholderQR<EigDMxZ> EigDFQRZ; // Complex.
|
|
|
|
typedef Eigen::PartialPivLU<EigDMxS> EigDPLUS; // Real.
|
|
typedef Eigen::PartialPivLU<EigDMxD> EigDPLUD; // Real.
|
|
typedef Eigen::PartialPivLU<EigDMxC> EigDPLUC; // Complex.
|
|
typedef Eigen::PartialPivLU<EigDMxZ> EigDPLUZ; // Complex.
|
|
|
|
typedef Eigen::HouseholderQR<EigDMxS> EigDPQRS; // Real.
|
|
typedef Eigen::HouseholderQR<EigDMxD> EigDPQRD; // Real.
|
|
typedef Eigen::HouseholderQR<EigDMxC> EigDPQRC; // Complex.
|
|
typedef Eigen::HouseholderQR<EigDMxZ> EigDPQRZ; // Complex.
|
|
|
|
// Definition of arpackSolver class.
|
|
|
|
typedef Eigen::Matrix<complex<double>, Eigen::Dynamic, 1> EigVecZ;
|
|
|
|
typedef vector<complex<double>> StdVecZ;
|
|
typedef vector<EigVecZ> StdVecEVZ;
|
|
|
|
// RC: Real or Complex.
|
|
// FD: Float or Double.
|
|
// EM: Eigen Matrix (sparse or dense).
|
|
// SLV: Solver.
|
|
template<typename RC, typename FD, typename EM, typename SLV>
|
|
class arpackSolver {
|
|
// Nested typedef.
|
|
|
|
typedef Eigen::Map<Eigen::Matrix<RC, Eigen::Dynamic, 1>> EV;
|
|
|
|
// Public methods.
|
|
|
|
public:
|
|
|
|
arpackSolver() {
|
|
symPb = true;
|
|
nbEV = 1;
|
|
nbCV = 2*nbEV+1;
|
|
tol = 1.e-6;
|
|
sigmaReal = sigmaImag = 0.;
|
|
dumpToFile = false;
|
|
restartFromFile = false;
|
|
mag = "LM";
|
|
maxIt = 100;
|
|
schur = false;
|
|
verbose = 0;
|
|
|
|
stdPb = true;
|
|
mode = 1;
|
|
nbIt = 0;
|
|
imsTime = 0.;
|
|
rciTime = 0.;
|
|
|
|
nbDim = 0;
|
|
resid = nullptr;
|
|
v = nullptr;
|
|
};
|
|
|
|
~arpackSolver() {
|
|
if (v) {delete [] v; v = nullptr;}
|
|
if (resid) {delete [] resid; resid = nullptr;}
|
|
};
|
|
|
|
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;
|
|
int rc = readMatrixMarket(fileName, n, m, i, j, Mij);
|
|
if (rc != 0) {cerr << "Error: read matrix market file KO" << endl; return rc;}
|
|
|
|
// Create matrix from file.
|
|
|
|
M = Eigen::SparseMatrix<RC>(n, m); // Set matrice dimensions.
|
|
vector<Eigen::Triplet<RC>> triplets;
|
|
a_uint nnz = Mij.size();
|
|
triplets.reserve(nnz);
|
|
for (size_t k = 0; k < nnz; k++) triplets.emplace_back(i[k], j[k], Mij[k]);
|
|
M.setFromTriplets(triplets.begin(), triplets.end()); // Set all (i, j, Mij).
|
|
|
|
return 0;
|
|
};
|
|
|
|
int createMatrix(string const & fileName, Eigen::Matrix<RC, Eigen::Dynamic, Eigen::Dynamic> & M) {
|
|
// Read matrix from file.
|
|
|
|
a_uint n = 0, m = 0;
|
|
vector<a_uint> i, j;
|
|
vector<RC> Mij;
|
|
int rc = readMatrixMarket(fileName, n, m, i, j, Mij);
|
|
if (rc != 0) {cerr << "Error: read matrix market file KO" << endl; return rc;}
|
|
|
|
// Create matrix from file.
|
|
|
|
M = Eigen::Matrix<RC, Eigen::Dynamic, Eigen::Dynamic>(n, m); // Set matrice dimensions.
|
|
M.setZero(n, m); // Avoid spurious/random values which may break solves (LU, QR, ...).
|
|
a_uint nnz = Mij.size();
|
|
for (size_t k = 0; k < nnz; k++) M(i[k], j[k]) = Mij[k];
|
|
|
|
return 0;
|
|
};
|
|
|
|
void dumpParameters() {
|
|
if (verbose >= 1) {
|
|
cout << endl << "arpackSolver:" << endl;
|
|
cout << endl << "symPb: " << symPb << endl;
|
|
cout << endl << "nbEV: " << nbEV << endl;
|
|
cout << endl << "nbCV: " << nbCV << endl;
|
|
cout << endl << "tol: " << tol << endl;
|
|
cout << endl << "sigmaReal: " << sigmaReal << endl;
|
|
cout << endl << "sigmaImag: " << sigmaImag << endl;
|
|
cout << endl << "dumpToFile: " << dumpToFile << endl;
|
|
cout << endl << "restartFromFile: " << restartFromFile << endl;
|
|
cout << endl << "mag: " << mag << endl;
|
|
cout << endl << "maxIt: " << maxIt << endl;
|
|
cout << endl << "schur: " << schur << endl;
|
|
}
|
|
};
|
|
|
|
int solve(EM & A, EM const * B = nullptr) {
|
|
stdPb = !B ? true : false;
|
|
if (nbCV > A.cols()) nbCV = A.cols(); // Cut-off arpack workspace dim.
|
|
|
|
dumpParameters();
|
|
if (verbose == 3) {
|
|
cout << endl << "arpackSolver:" << endl;
|
|
cout << endl << "A:" << endl;
|
|
cout << endl << A << endl;
|
|
if (B) {
|
|
cout << endl << "B:" << endl;
|
|
cout << endl << *B << endl;
|
|
}
|
|
}
|
|
nbDim = A.rows();
|
|
|
|
// If needed, transform the initial problem into a new one that arpack can handle.
|
|
|
|
auto eps = numeric_limits<FD>::epsilon();
|
|
bool shiftReal = (fabs(sigmaReal) > eps) ? true : false;
|
|
bool shiftImag = (fabs(sigmaImag) > eps) ? true : false;
|
|
bool backTransform = false;
|
|
mode = 0;
|
|
if (stdPb) {
|
|
mode = 1;
|
|
// 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);
|
|
A -= sigma*I;
|
|
backTransform = true;
|
|
}
|
|
}
|
|
else {
|
|
// CAUTION: back transform must NOT be done if mode > 1.
|
|
mode = 2;
|
|
if (shiftReal || shiftImag) mode = 3;
|
|
}
|
|
if (verbose >= 1) {
|
|
cout << endl << "arpackSolver:" << endl;
|
|
cout << endl << "mode " << mode << ", backTransform " << (backTransform ? "yes" : "no") << endl;
|
|
}
|
|
|
|
// Solve with arpack.
|
|
// Note: when initializing Eigen solvers, API differ depending on solvers.
|
|
|
|
SLV solver;
|
|
int rc = initSolver(solver);
|
|
if (rc != 0) {cerr << "Error: initialize solver KO" << endl; return rc;}
|
|
rc = solve(A, B, solver);
|
|
if (rc != 0) {cerr << "Error: arpack solve KO" << endl; return rc;}
|
|
|
|
// If needed, transform back the arpack problem into the initial problem.
|
|
|
|
if (backTransform) {
|
|
EM I(A.rows(), A.cols());
|
|
I.setIdentity();
|
|
RC sigma; makeSigma(sigma);
|
|
A += sigma*I;
|
|
for (size_t i = 0; i < val.size(); i++) val[i] += sigma;
|
|
}
|
|
|
|
return 0;
|
|
};
|
|
|
|
int checkEigVec(EM const & A, EM const * B = nullptr, double const maxResNorm = 1.e-3) const {
|
|
// Check eigen vectors.
|
|
|
|
string rs = schur ? "Schur" : "Ritz";
|
|
|
|
if (vec.size() == 0) {
|
|
cerr << "Error: no " << rs << " value / vector found" << endl;
|
|
return 1;
|
|
}
|
|
|
|
for (size_t i = 0; i < vec.size(); i++) {
|
|
EigVecZ V = vec[i];
|
|
complex<double> lambda = val[i];
|
|
if (verbose >= 1) {
|
|
cout << endl << "arpackSolver:" << endl;
|
|
cout << endl << rs << " value " << setw(3) << i << ": " << lambda << endl;
|
|
if (verbose >= 2) {
|
|
cout << endl << rs << " vector " << setw(3) << i << " (norm " << V.norm() << "): " << endl;
|
|
cout << endl << V << 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 " << 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:
|
|
|
|
void makeConstant( float & cst, float const & val) {cst = val;};
|
|
|
|
void makeConstant( double & cst, double const & val) {cst = val;};
|
|
|
|
void makeConstant(complex< float> & cst, float const & val) {cst = complex<float>(val, val);};
|
|
|
|
void makeConstant(complex<double> & cst, double const & val) {cst = complex<double>(val, val);};
|
|
|
|
int readMatrixMarket(string const & fileName,
|
|
a_uint & n, a_uint & m, vector<a_uint> & i, vector<a_uint> & j, vector<RC> & Mij) {
|
|
ifstream inp(fileName);
|
|
if (!inp) {cerr << "Error: can not open " << fileName << endl; return 1;}
|
|
|
|
// Read matrix from file.
|
|
|
|
a_uint l = 0, nnz = 0;
|
|
do {
|
|
// Skip comments.
|
|
|
|
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] == '%' || inpLine[0] == '#') continue; // Comments skipped, begin reading.
|
|
|
|
// Read matrix market file.
|
|
|
|
stringstream inpSS(inpLine);
|
|
if (n == 0 && m == 0) { // Header.
|
|
inpSS >> n >> m;
|
|
if (!inpSS) {cerr << "Error: bad header (n, m)" << endl; return 1;}
|
|
if (nnz == 0) {
|
|
inpSS >> nnz;
|
|
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);
|
|
}
|
|
}
|
|
}
|
|
else { // Body.
|
|
a_uint k = 0, l = 0;
|
|
RC zero; makeConstant(zero, 0.);
|
|
RC Mkl = zero;
|
|
inpSS >> k >> l >> Mkl;
|
|
if (!inpSS) {cerr << "Error: bad line (" << fileName << ", line " << l << ")" << endl; return 1;}
|
|
i.push_back(k);
|
|
j.push_back(l);
|
|
Mij.push_back(Mkl);
|
|
}
|
|
}
|
|
while (inp);
|
|
|
|
// Handle 1-based -> 0-based.
|
|
|
|
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;
|
|
};
|
|
|
|
void aupd(a_int * ido, char const * bMat, a_int nbDim, char const * which, float * resid, float * v,
|
|
a_int ldv, a_int * iparam, a_int * ipntr, float * workd, float * workl, a_int lworkl, float * & rwork,
|
|
a_int * info) {
|
|
assert(rwork == nullptr);
|
|
if (symPb) {
|
|
ssaupd_c(ido, bMat, nbDim, which, nbEV, tol, resid, nbCV, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
|
|
}
|
|
else {
|
|
snaupd_c(ido, bMat, nbDim, which, nbEV, tol, resid, nbCV, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
|
|
}
|
|
};
|
|
|
|
void aupd(a_int * ido, char const * bMat, a_int nbDim, char const * which, double * resid, double * v,
|
|
a_int ldv, a_int * iparam, a_int * ipntr, double * workd, double * workl, a_int lworkl, double * & rwork,
|
|
a_int * info) {
|
|
assert(rwork == nullptr);
|
|
if (symPb) {
|
|
dsaupd_c(ido, bMat, nbDim, which, nbEV, tol, resid, nbCV, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
|
|
}
|
|
else {
|
|
dnaupd_c(ido, bMat, nbDim, which, nbEV, tol, resid, nbCV, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
|
|
}
|
|
};
|
|
|
|
void aupd(a_int * ido, char const * bMat, a_int nbDim, char const * which, complex<float> * resid, complex<float> * v,
|
|
a_int ldv, a_int * iparam, a_int * ipntr, complex<float> * workd, complex<float> * workl, a_int lworkl, float * & rwork,
|
|
a_int * info) {
|
|
if (!rwork) rwork = new float[nbCV];
|
|
cnaupd_c(ido, bMat, nbDim, which, nbEV, tol, reinterpret_cast<a_fcomplex*>(resid), nbCV,
|
|
reinterpret_cast<a_fcomplex*>(v), ldv, iparam, ipntr, reinterpret_cast<a_fcomplex*>(workd),
|
|
reinterpret_cast<a_fcomplex*>(workl), lworkl, rwork, info);
|
|
};
|
|
|
|
void aupd(a_int * ido, char const * bMat, a_int nbDim, char const * which, complex<double> * resid, complex<double> * v,
|
|
a_int ldv, a_int * iparam, a_int * ipntr, complex<double> * workd, complex<double> * workl, a_int lworkl, double * & rwork,
|
|
a_int * info) {
|
|
if (!rwork) rwork = new double[nbCV];
|
|
znaupd_c(ido, bMat, nbDim, which, nbEV, tol, reinterpret_cast<a_dcomplex*>(resid), nbCV,
|
|
reinterpret_cast<a_dcomplex*>(v), ldv, iparam, ipntr, reinterpret_cast<a_dcomplex*>(workd),
|
|
reinterpret_cast<a_dcomplex*>(workl), lworkl, rwork, info);
|
|
};
|
|
|
|
void spectrum(RC * d, RC * z, a_int nbDim, a_int * iparam) {
|
|
// Arpack compute the whole spectrum.
|
|
|
|
a_int nbConv = iparam[4];
|
|
val.reserve(nbConv);
|
|
for (a_int i = 0; d && i < nbConv; i++) {
|
|
complex<double> lambda(d[i]);
|
|
val.push_back(lambda);
|
|
if (val.size() == (size_t) nbEV) break; // If more converged than requested, likely not accurate (check KO).
|
|
}
|
|
|
|
vec.reserve(nbConv);
|
|
for (a_int i = 0; z && i < nbConv; i++) {
|
|
EV V = EV(z + i*nbDim, nbDim);
|
|
vec.push_back(V.template cast<complex<double>>());
|
|
if (vec.size() == (size_t) nbEV) break; // If more converged than requested, likely not accurate (check KO).
|
|
}
|
|
};
|
|
|
|
void halfSpectrum(RC * dr, RC * di, RC * z, a_int nbDim, a_int * iparam) {
|
|
// Arpack compute only half of the spectrum.
|
|
|
|
a_int nbConv = iparam[4];
|
|
val.reserve(nbConv);
|
|
for (a_int i = 0; dr && di && i <= nbConv/2; i++) { // Scan first half of the spectrum.
|
|
// Get first half of the spectrum.
|
|
|
|
complex<double> lambda(dr[i], di[i]);
|
|
val.push_back(lambda);
|
|
if (val.size() == (size_t) nbEV) break; // If more converged than requested, likely not accurate (check KO).
|
|
|
|
// Deduce second half of the spectrum.
|
|
|
|
val.push_back(complex<double>(lambda.real(), -1.*lambda.imag()));
|
|
if (val.size() == (size_t) nbEV) break; // If more converged than requested, likely not accurate (check KO).
|
|
}
|
|
|
|
vec.reserve(nbConv);
|
|
for (a_int i = 0; z && i <= nbConv/2; i++) { // Scan half spectrum.
|
|
// Get first half of the spectrum.
|
|
|
|
EV Vr = EV(z + (2*i+0)*nbDim, nbDim); // Real part.
|
|
EV Vi = EV(z + (2*i+1)*nbDim, nbDim); // Imaginary part.
|
|
complex<double> imag(0., 1.);
|
|
EigVecZ V = Vr.template cast<complex<double>>() + imag * Vi.template cast<complex<double>>();
|
|
vec.push_back(V);
|
|
if (vec.size() == (size_t) nbEV) break; // If more converged than requested, likely not accurate (check KO).
|
|
|
|
// Deduce second half of the spectrum.
|
|
|
|
V = Vr.template cast<complex<double>>() - imag * Vi.template cast<complex<double>>();
|
|
vec.push_back(V);
|
|
if (vec.size() == (size_t) nbEV) break; // If more converged than requested, likely not accurate (check KO).
|
|
}
|
|
};
|
|
|
|
int eupd(a_int rvec, char const * howmny, a_int const * select, float * z,
|
|
a_int ldz, char const * bMat, a_int nbDim, char const * which, float * resid, float * v,
|
|
a_int ldv, a_int * iparam, a_int * ipntr, float * workd, float * workl, a_int lworkl, float * rwork,
|
|
a_int & info) {
|
|
assert(rwork == nullptr);
|
|
if (symPb) {
|
|
float * d = new float[nbEV]; for (a_int k = 0; k < nbEV; k++) d[k] = 0.;
|
|
|
|
sseupd_c(rvec, howmny, select, d, z, ldz, sigmaReal,
|
|
bMat, nbDim, which, nbEV, tol, resid, nbCV, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
|
|
if (info == -14) cerr << "Error: sseupd - KO: ssaupd did not find any eigenvalues to sufficient accuracy" << endl;
|
|
if (info < 0 && info != -14 /*-14: don't break*/) {cerr << "Error: sseupd - KO with info " << info << endl; return 1;}
|
|
|
|
spectrum(d, z, nbDim, iparam);
|
|
|
|
if (d) {delete [] d; d = nullptr;}
|
|
}
|
|
else {
|
|
float * dr = new float[nbEV+1]; for (a_int k = 0; k < nbEV+1; k++) dr[k] = 0.;
|
|
float * di = new float[nbEV+1]; for (a_int k = 0; k < nbEV+1; k++) di[k] = 0.;
|
|
float * workev = new float[3*nbCV];
|
|
|
|
sneupd_c(rvec, howmny, select, dr, di, z, ldz, sigmaReal, sigmaImag, workev,
|
|
bMat, nbDim, which, nbEV, tol, resid, nbCV, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
|
|
if (info == -14) cerr << "Error: sneupd - KO: snaupd did not find any eigenvalues to sufficient accuracy" << endl;
|
|
if (info < 0 && info != -14 /*-14: don't break*/) {cerr << "Error: sneupd - KO with info " << info << endl; return 1;}
|
|
|
|
halfSpectrum(dr, di, z, nbDim, iparam);
|
|
|
|
if (workev) {delete [] workev; workev = nullptr;}
|
|
if (dr) {delete [] dr; dr = nullptr;}
|
|
if (di) {delete [] di; di = nullptr;}
|
|
}
|
|
|
|
return 0;
|
|
};
|
|
|
|
int eupd(a_int rvec, char const * howmny, a_int const * select, double * z,
|
|
a_int ldz, char const * bMat, a_int nbDim, char const * which, double * resid, double * v,
|
|
a_int ldv, a_int * iparam, a_int * ipntr, double * workd, double * workl, a_int lworkl, double * rwork,
|
|
a_int & info) {
|
|
assert(rwork == nullptr);
|
|
if (symPb) {
|
|
double * d = new double[nbEV]; for (a_int k = 0; k < nbEV; k++) d[k] = 0.;
|
|
|
|
dseupd_c(rvec, howmny, select, d, z, ldz, sigmaReal,
|
|
bMat, nbDim, which, nbEV, tol, resid, nbCV, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
|
|
if (info == -14) cerr << "Error: dseupd - KO: dsaupd did not find any eigenvalues to sufficient accuracy" << endl;
|
|
if (info < 0 && info != -14 /*-14: don't break*/) {cerr << "Error: dseupd - KO with info " << info << endl; return 1;}
|
|
|
|
spectrum(d, z, nbDim, iparam);
|
|
|
|
if (d) {delete [] d; d = nullptr;}
|
|
}
|
|
else {
|
|
double * dr = new double[nbEV+1]; for (a_int k = 0; k < nbEV+1; k++) dr[k] = 0.;
|
|
double * di = new double[nbEV+1]; for (a_int k = 0; k < nbEV+1; k++) di[k] = 0.;
|
|
double * workev = new double[3*nbCV];
|
|
|
|
dneupd_c(rvec, howmny, select, dr, di, z, ldz, sigmaReal, sigmaImag, workev,
|
|
bMat, nbDim, which, nbEV, tol, resid, nbCV, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
|
|
if (info == -14) cerr << "Error: dneupd - KO: dnaupd did not find any eigenvalues to sufficient accuracy" << endl;
|
|
if (info < 0 && info != -14 /*-14: don't break*/) {cerr << "Error: dneupd - KO with info " << info << endl; return 1;}
|
|
|
|
halfSpectrum(dr, di, z, nbDim, iparam);
|
|
|
|
if (workev) {delete [] workev; workev = nullptr;}
|
|
if (dr) {delete [] dr; dr = nullptr;}
|
|
if (di) {delete [] di; di = nullptr;}
|
|
}
|
|
|
|
return 0;
|
|
};
|
|
|
|
int eupd(a_int rvec, char const * howmny, a_int const * select, complex<float> * z,
|
|
a_int ldz, char const * bMat, a_int nbDim, char const * which, complex<float> * resid, complex<float> * v,
|
|
a_int ldv, a_int * iparam, a_int * ipntr, complex<float> * workd, complex<float> * workl, a_int lworkl, float * rwork,
|
|
a_int & info) {
|
|
complex<float> * d = new complex<float>[nbEV+1]; for (a_int k = 0; k < nbEV+1; k++) d[k] = complex<float>(0., 0.);
|
|
complex<float> * workev = new complex<float>[2*nbCV];
|
|
|
|
complex<float> sigma = complex<float>((float) sigmaReal, (float) sigmaImag);
|
|
cneupd_c(rvec, howmny, select, reinterpret_cast<a_fcomplex*>(d), reinterpret_cast<a_fcomplex*>(z), ldz,
|
|
reinterpret_cast<a_fcomplex &>(sigma), reinterpret_cast<a_fcomplex*>(workev),
|
|
bMat, nbDim, which, nbEV, tol, reinterpret_cast<a_fcomplex*>(resid), nbCV,
|
|
reinterpret_cast<a_fcomplex*>(v), ldv, iparam, ipntr,
|
|
reinterpret_cast<a_fcomplex*>(workd), reinterpret_cast<a_fcomplex*>(workl), lworkl, rwork, &info);
|
|
if (info == -14) cerr << "Error: cneupd - KO: cnaupd did not find any eigenvalues to sufficient accuracy" << endl;
|
|
if (info < 0 && info != -14 /*-14: don't break*/) {cerr << "Error: cneupd - KO with info " << info << endl; return 1;}
|
|
|
|
spectrum(d, z, nbDim, iparam);
|
|
|
|
if (workev) {delete [] workev; workev = nullptr;}
|
|
if (d) {delete [] d; d = nullptr;}
|
|
|
|
return 0;
|
|
};
|
|
|
|
int eupd(a_int rvec, char const * howmny, a_int const * select, complex<double> * z,
|
|
a_int ldz, char const * bMat, a_int nbDim, char const * which, complex<double> * resid, complex<double> * v,
|
|
a_int ldv, a_int * iparam, a_int * ipntr, complex<double> * workd, complex<double> * workl, a_int lworkl, double * rwork,
|
|
a_int & info) {
|
|
complex<double> * d = new complex<double>[nbEV+1]; for (a_int k = 0; k < nbEV+1; k++) d[k] = complex<double>(0., 0.);
|
|
complex<double> * workev = new complex<double>[2*nbCV];
|
|
|
|
complex<double> sigma = complex<double>(sigmaReal, sigmaImag);
|
|
zneupd_c(rvec, howmny, select, reinterpret_cast<a_dcomplex*>(d), reinterpret_cast<a_dcomplex*>(z), ldz,
|
|
reinterpret_cast<a_dcomplex &>(sigma), reinterpret_cast<a_dcomplex*>(workev),
|
|
bMat, nbDim, which, nbEV, tol, reinterpret_cast<a_dcomplex*>(resid), nbCV,
|
|
reinterpret_cast<a_dcomplex*>(v), ldv, iparam, ipntr,
|
|
reinterpret_cast<a_dcomplex*>(workd), reinterpret_cast<a_dcomplex*>(workl), lworkl, rwork, &info);
|
|
if (info == -14) cerr << "Error: zneupd - KO: znaupd did not find any eigenvalues to sufficient accuracy" << endl;
|
|
if (info < 0 && info != -14 /*-14: don't break*/) {cerr << "Error: zneupd - KO with info " << info << endl; return 1;}
|
|
|
|
spectrum(d, z, nbDim, iparam);
|
|
|
|
if (workev) {delete [] workev; workev = nullptr;}
|
|
if (d) {delete [] d; d = nullptr;}
|
|
|
|
return 0;
|
|
};
|
|
|
|
int setMode(int const mode, EM const & A, EM const * B, SLV & solver) {
|
|
int rc = 1;
|
|
|
|
if (mode == 1) {
|
|
rc = 0;
|
|
}
|
|
else if (mode == 2 || mode == 3) {
|
|
if (!stdPb && !B) {cerr << "Error: generalized problem without B" << endl; return 1;}
|
|
|
|
if (mode == 2) { // Invert mode.
|
|
solver.compute(*B);
|
|
}
|
|
else { // Shift invert mode.
|
|
RC sigma; makeSigma(sigma);
|
|
auto S = A - sigma*(*B);
|
|
solver.compute(S);
|
|
}
|
|
|
|
rc = 0;
|
|
}
|
|
else {cerr << "Error: arpack mode must be 1, 2 or 3 - KO" << endl; rc = 1;}
|
|
|
|
return rc;
|
|
};
|
|
|
|
template <typename T>
|
|
int saveSolve(string const & fileName,
|
|
a_int const & nbDim, T * rv) {
|
|
ofstream ofs(fileName.c_str(), ofstream::trunc);
|
|
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;
|
|
};
|
|
|
|
template <typename T>
|
|
int restartSolve(string const & fileName,
|
|
a_int const & nbDim, T * rv,
|
|
bool allowZero = true) {
|
|
ifstream ifs(fileName.c_str());
|
|
if (ifs.is_open()) {
|
|
a_int nDim = 0;
|
|
ifs >> nDim;
|
|
if (nDim != nbDim) {cerr << "Error: bad dim - restart KO" << endl; return 1;}
|
|
for (a_int n = 0; rv && n < nbDim; n++) {
|
|
RC val; makeConstant(val, 0.);
|
|
ifs >> val;
|
|
if (abs(val) < 1.e-6 && !allowZero) {
|
|
// Do NOT let residual be zero: this stops arpack to iterate (info = -9).
|
|
auto eps = numeric_limits<FD>::epsilon();
|
|
RC epsilon; makeConstant(epsilon, eps);
|
|
val = epsilon;
|
|
}
|
|
rv[n] = val;
|
|
}
|
|
|
|
if (verbose >= 1) {
|
|
cout << endl << "arpackSolver:" << endl;
|
|
cout << endl << fileName << ": restart OK" << endl;
|
|
if (verbose >= 2) {
|
|
for (a_int n = 0; rv && n < nbDim; n++) cout << rv[n] << endl;
|
|
}
|
|
}
|
|
}
|
|
return 0;
|
|
};
|
|
|
|
int initPointerSize(a_int & iparamSz, a_int & ipntrSz, string const & aeupd) {
|
|
iparamSz = 11; ipntrSz = 14;
|
|
if (aeupd == "aupd") {
|
|
if (is_same<RC, double>::value && symPb) ipntrSz = 11;
|
|
if (is_same<RC, float>::value && symPb) ipntrSz = 11;
|
|
return 0;
|
|
}
|
|
if (aeupd == "eupd") {
|
|
if (is_same<RC, double>::value && symPb) {iparamSz = 7; ipntrSz = 11;}
|
|
if (is_same<RC, float>::value && symPb) {iparamSz = 7; ipntrSz = 11;}
|
|
return 0;
|
|
}
|
|
return 1;
|
|
};
|
|
|
|
int solve(EM const & A, EM const * B, SLV & solver) {
|
|
if (!stdPb && !B) {cerr << "Error: generalized problem without B" << endl; return 1;}
|
|
|
|
// Arpack set up.
|
|
|
|
// Note: some in/out parameters passed to d[sn][ae]upd are set to 0. before use.
|
|
// d[sn][ae]upd uses dgetv0 to generate a random starting vector (when info is initialized to 0).
|
|
// dgetv0 rely on resid/v: resid/v should be initialized to 0.0 to avoid "bad" starting random vectors.
|
|
|
|
char const * which = mag.c_str();
|
|
a_int ido = 0; // First call to arpack.
|
|
char const * iMat = "I";
|
|
char const * gMat = "G";
|
|
char const * bMat = (mode == 1) ? iMat : gMat;
|
|
a_int nbDim = A.rows();
|
|
auto eps = numeric_limits<FD>::epsilon();
|
|
RC epsilon; makeConstant(epsilon, eps);
|
|
if (!resid) {
|
|
resid = new RC[nbDim];
|
|
// Do NOT let residual be zero: this stops arpack to iterate (info = -9).
|
|
for (a_int n = 0; n < nbDim; n++) resid[n] = epsilon;
|
|
};
|
|
a_int ldv = nbDim;
|
|
RC cst; makeConstant(cst, 10.);
|
|
RC v0 = cst * epsilon; // Start with something close to zero but *not* exactly zero.
|
|
if (!v) {
|
|
v = new RC[ldv*nbCV];
|
|
for (a_int n = 0; n < ldv*nbCV; n++) v[n] = v0;
|
|
};
|
|
a_int iparamSz = 0, ipntrSz = 0;
|
|
int rc = initPointerSize(iparamSz, ipntrSz, "aupd");
|
|
if (rc != 0) {cerr << "Error: bad iparam/ipntr initialization for aupd" << endl; return rc;}
|
|
vector<a_int> iparamAupd(iparamSz, 0);
|
|
iparamAupd[0] = 1; // Use exact shifts (=> we'll never have ido == 3).
|
|
iparamAupd[2] = maxIt; // Maximum number of iterations.
|
|
iparamAupd[3] = 1; // Block size.
|
|
iparamAupd[4] = 0; // Number of ev found by arpack.
|
|
iparamAupd[6] = mode;
|
|
RC zero; makeConstant(zero, 0.);
|
|
vector<a_int> ipntrAupd(ipntrSz, 0);
|
|
RC * workd = new RC[3*nbDim]; for (a_int n = 0; n < 3*nbDim; n++) workd[n] = zero; // Avoid "bad" X/Y vector.
|
|
a_int lworkl = symPb ? nbCV*nbCV + 8*nbCV : 3*nbCV*nbCV + 6*nbCV;
|
|
RC * workl = new RC[lworkl];
|
|
a_int info = 0; // Use random initial residual vector.
|
|
|
|
// Handling restart.
|
|
|
|
if (restartFromFile) {
|
|
info = 1; // Restart.
|
|
rc = restartSolve("arpackSolver.resid.out", nbDim, resid, false);
|
|
if (rc != 0) {cerr << "Error: bad restart (resid)" << endl; return rc;}
|
|
rc = restartSolve("arpackSolver.v.out", ldv*nbCV, v);
|
|
if (rc != 0) {cerr << "Error: bad restart (v)" << endl; return rc;}
|
|
}
|
|
|
|
// Initialize solver.
|
|
|
|
auto start = chrono::high_resolution_clock::now();
|
|
rc = setMode(mode, A, B, solver);
|
|
if (rc != 0) {cerr << "Error: bad arpack mode" << endl; return rc;}
|
|
auto stop = chrono::high_resolution_clock::now();
|
|
imsTime = chrono::duration_cast<chrono::milliseconds>(stop - start).count()/1000.;
|
|
|
|
// Arpack solve.
|
|
|
|
FD * rwork = nullptr;
|
|
do {
|
|
// Call arpack.
|
|
|
|
aupd(&ido, bMat, nbDim, which, resid, v, ldv, iparamAupd.data(), ipntrAupd.data(), workd, workl, lworkl, rwork, &info);
|
|
if (info == 1) cerr << "Error: [dz][sn]aupd - KO: maximum number of iterations taken. Increase --maxIt..." << endl;
|
|
if (info == 3) cerr << "Error: [dz][sn]aupd - KO: no shifts could be applied. Increase --nbCV..." << endl;
|
|
if (info == -9) cerr << "Error: [dz][sn]aupd - KO: starting vector is zero. Retry: play with shift..." << endl;
|
|
if (info < 0) {cerr << "Error: [dz][sn]aupd - KO with info " << info << ", nbIt " << iparamAupd[2] << endl; return 1;}
|
|
|
|
// Reverse Communication Interface: perform actions according to arpack.
|
|
|
|
start = chrono::high_resolution_clock::now();
|
|
|
|
a_int xIdx = ipntrAupd[0] - 1; // 0-based (Fortran is 1-based).
|
|
a_int yIdx = ipntrAupd[1] - 1; // 0-based (Fortran is 1-based).
|
|
|
|
EV X(workd + xIdx, nbDim); // Arpack provides X.
|
|
EV Y(workd + yIdx, nbDim); // Arpack provides Y.
|
|
|
|
if (ido == -1) {
|
|
if (iparamAupd[6] == 1) {
|
|
Y = A * X;
|
|
}
|
|
else if (iparamAupd[6] == 2) {
|
|
Y = A * X;
|
|
auto YY = Y; // Use copy of Y (not Y) for solve (avoid potential memory overwrite as Y is both in/out).
|
|
Y = solver.solve(YY); // Y = B^-1 * A * X.
|
|
}
|
|
else if (iparamAupd[6] == 3) {
|
|
auto Z = (*B) * X; // Z = B * X.
|
|
Y = solver.solve(Z); // Y = (A - sigma * B)^-1 * B * X.
|
|
}
|
|
}
|
|
else if (ido == 1) {
|
|
if (iparamAupd[6] == 1) {
|
|
Y = A * X;
|
|
}
|
|
else if (iparamAupd[6] == 2) {
|
|
Y = A * X;
|
|
if (symPb) X = Y; // Remark 5 in dsaupd documentation.
|
|
auto YY = Y; // Use copy of Y (not Y) for solve (avoid potential memory overwrite as Y is both in/out).
|
|
Y = solver.solve(YY); // Y = B^-1 * A * X.
|
|
}
|
|
else if (iparamAupd[6] == 3) {
|
|
a_int zIdx = ipntrAupd[2] - 1; // 0-based (Fortran is 1-based).
|
|
EV Z(workd + zIdx, nbDim); // Arpack provides Z.
|
|
Y = solver.solve(Z); // Y = (A - sigma * B)^-1 * B * X.
|
|
}
|
|
}
|
|
else if (ido == 2) {
|
|
if (iparamAupd[6] == 1) Y = X; // Y = I * X.
|
|
else if (iparamAupd[6] == 2) Y = (*B) * X; // Y = B * X.
|
|
else if (iparamAupd[6] == 3) Y = (*B) * X; // Y = B * X.
|
|
}
|
|
else if (ido != 99) {cerr << "Error: unexpected ido " << ido << " - KO" << endl; return 1;}
|
|
|
|
stop = chrono::high_resolution_clock::now();
|
|
rciTime += chrono::duration_cast<chrono::milliseconds>(stop - start).count()/1000.;
|
|
|
|
} while (ido != 99);
|
|
|
|
// Get arpack results (computed eigen values and vectors).
|
|
|
|
nbIt = iparamAupd[2]; // Actual number of iterations.
|
|
a_int rvec = 1;
|
|
char const * howmnyA = "A"; // Ritz vectors.
|
|
char const * howmnyP = "P"; // Schur vectors.
|
|
char const * howmny = schur ? howmnyP : howmnyA;
|
|
a_int * select = new a_int[nbCV]; for (a_int n = 0; n < nbCV; n++) select[n] = 1;
|
|
a_int const nbZ = nbDim*(nbEV+1); // Caution: nbEV+1 for dneupd.
|
|
RC * z = new RC[nbZ]; for (a_int n = 0; n < nbZ; n++) z[n] = zero;
|
|
a_int ldz = nbDim;
|
|
rc = initPointerSize(iparamSz, ipntrSz, "eupd");
|
|
if (rc != 0) {cerr << "Error: bad iparam/ipntr initialization for eupd" << endl; return rc;}
|
|
vector<a_int> iparamEupd(iparamSz, 0);
|
|
for (a_int p = 0; p < iparamSz; p++) iparamEupd[p] = iparamAupd[p]; // Initialize eupd parameters with aupd ones.
|
|
vector<a_int> ipntrEupd(ipntrSz, 0);
|
|
for (a_int p = 0; p < ipntrSz; p++) ipntrEupd[p] = ipntrAupd[p]; // Initialize eupd parameters with aupd ones.
|
|
rc = eupd(rvec, howmny, select, z, ldz, bMat, nbDim, which, resid, v, ldv, iparamEupd.data(), ipntrEupd.data(), workd, workl, lworkl, rwork, info);
|
|
if (rc != 0) {cerr << "Error: bad arpack eupd" << endl; return rc;}
|
|
|
|
// Save solve (dump results to file) to allow later restart.
|
|
|
|
if (dumpToFile) {
|
|
saveSolve("arpackSolver.resid.out", nbDim, resid);
|
|
saveSolve("arpackSolver.v.out", ldv*nbCV, v);
|
|
}
|
|
|
|
// Clean.
|
|
|
|
if (rwork) {delete [] rwork; rwork = nullptr;}
|
|
if (z) {delete [] z; z = nullptr;}
|
|
if (select) {delete [] select; select = nullptr;}
|
|
if (workl) {delete [] workl; workl = nullptr;}
|
|
if (workd) {delete [] workd; workd = nullptr;}
|
|
|
|
return 0;
|
|
};
|
|
|
|
void makeSigma( float & sigma) {sigma = (float) sigmaReal;};
|
|
|
|
void makeSigma( double & sigma) {sigma = sigmaReal;};
|
|
|
|
void makeSigma(complex< float> & sigma) {sigma = complex< float>((float) sigmaReal, (float) sigmaImag);};
|
|
|
|
void makeSigma(complex<double> & sigma) {sigma = complex<double>(sigmaReal, sigmaImag);};
|
|
|
|
virtual int initSolver(SLV & solver) = 0;
|
|
|
|
virtual void dumpAllParameters() = 0;
|
|
|
|
// Public members.
|
|
|
|
public:
|
|
|
|
// Arpack parameters.
|
|
|
|
bool symPb; // Symmetric problem.
|
|
a_int nbEV;
|
|
a_int nbCV;
|
|
double tol;
|
|
double sigmaReal, sigmaImag; // Eigen value translation: look for lambda+sigma instead of lambda.
|
|
bool dumpToFile; // Dump resid and v to arpackSolver.*.out files after solve.
|
|
bool restartFromFile; // Restart solve with resid and v values provided in arpackSolver.*.out files.
|
|
string mag; // Magnitude <=> "which" arpack parameter.
|
|
int maxIt;
|
|
bool schur;
|
|
int verbose;
|
|
|
|
// Arpack outputs.
|
|
|
|
bool stdPb; // Standard or generalized (= not standard).
|
|
StdVecZ val; // Eigen values.
|
|
StdVecEVZ vec; // Eigen vectors.
|
|
int mode;
|
|
int nbIt;
|
|
double imsTime; // Init mode solver time.
|
|
double rciTime; // Reverse communication interface time.
|
|
|
|
// Protected members.
|
|
|
|
protected:
|
|
|
|
a_int nbDim;
|
|
|
|
// Private members.
|
|
|
|
private:
|
|
|
|
RC * resid; // Saved: enable restart from previous solve.
|
|
RC * v; // Saved: enable restart from previous solve.
|
|
};
|
|
|
|
// Definition of arpackItrSolver class: specialization of arpackSolver using iterative solvers.
|
|
// Note: Eigen provides iterative solvers only for sparse matrices.
|
|
|
|
// RC: Real or Complex.
|
|
// FD: Float or Double.
|
|
// EM: Eigen Matrix (sparse only, not dense).
|
|
// SLV: Solver.
|
|
template<typename RC, typename FD, typename EM, typename SLV>
|
|
class arpackItrSolver: public arpackSolver<RC, FD, EM, SLV> {
|
|
// Public methods.
|
|
|
|
public:
|
|
|
|
arpackItrSolver(): arpackSolver<RC, FD, EM, SLV>() {
|
|
slvTol = 1.e-6;
|
|
slvMaxIt = 100;
|
|
slvILUDropTol = 1.;
|
|
slvILUFillFactor = 2;
|
|
};
|
|
|
|
void dumpAllParameters() {
|
|
this->dumpParameters();
|
|
if (this->verbose >= 1) {
|
|
cout << endl << "arpackItrSolver:" << endl;
|
|
cout << endl << "slvTol: " << slvTol << endl;
|
|
cout << endl << "slvMaxIt: " << slvMaxIt << endl;
|
|
cout << endl << "slvILUDropTol: " << slvILUDropTol << endl;
|
|
cout << endl << "slvILUFillFactor: " << slvILUFillFactor << endl;
|
|
}
|
|
};
|
|
|
|
virtual int initSolver(Eigen::BiCGSTAB<EM> & solver) {
|
|
// Solve with arpack using sparse matrices and iterative solvers.
|
|
|
|
solver.setTolerance(slvTol);
|
|
solver.setMaxIterations(slvMaxIt);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::ConjugateGradient<EM> & solver) {
|
|
// Solve with arpack using sparse matrices and iterative solvers.
|
|
|
|
solver.setTolerance(slvTol);
|
|
solver.setMaxIterations(slvMaxIt);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::BiCGSTAB<EM, Eigen::IncompleteLUT<RC>> & solver) {
|
|
// Solve with arpack using sparse matrices and iterative solvers.
|
|
|
|
solver.setTolerance(slvTol);
|
|
solver.setMaxIterations(slvMaxIt);
|
|
solver.preconditioner().setDroptol(slvILUDropTol);
|
|
solver.preconditioner().setFillfactor(slvILUFillFactor);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::ConjugateGradient<EM, Eigen::Lower|Eigen::Upper, Eigen::IncompleteLUT<RC>> & solver) {
|
|
// Solve with arpack using sparse matrices and iterative solvers.
|
|
|
|
solver.setTolerance(slvTol);
|
|
solver.setMaxIterations(slvMaxIt);
|
|
solver.preconditioner().setDroptol(slvILUDropTol);
|
|
solver.preconditioner().setFillfactor(slvILUFillFactor);
|
|
|
|
return 0;
|
|
};
|
|
|
|
// Public members.
|
|
|
|
public:
|
|
|
|
// Iterative solvers parameters.
|
|
|
|
double slvTol; // Tolerance of the iterative mode solver.
|
|
int slvMaxIt; // Maximum number of iterations of the iterative mode solver.
|
|
double slvILUDropTol; // Drop tolerance of the ILU preconditioner (if any) of the iterative mode solver.
|
|
int slvILUFillFactor; // Fill factor of the ILU preconditioner (if any) of the iterative mode solver.
|
|
};
|
|
|
|
// Definition of arpackDrtSolver class: specialization of arpackSolver using direct solvers.
|
|
// Note: Eigen provides direct solvers for both sparse and dense matrices.
|
|
|
|
// RC: Real or Complex.
|
|
// FD: Float or Double.
|
|
// EM: Eigen Matrix (sparse or dense).
|
|
// SLV: Solver.
|
|
template<typename RC, typename FD, typename EM, typename SLV>
|
|
class arpackDrtSolver: public arpackSolver<RC, FD, EM, SLV> {
|
|
// Public methods.
|
|
|
|
public:
|
|
|
|
arpackDrtSolver(): arpackSolver<RC, FD, EM, SLV>() {
|
|
slvPvtThd = 1.e-6;
|
|
slvOffset = 0.;
|
|
slvScale = 1.;
|
|
};
|
|
|
|
void dumpAllParameters() {
|
|
this->dumpParameters();
|
|
if (this->verbose >= 1) {
|
|
cout << endl << "arpackDrtSolver:" << endl;
|
|
cout << endl << "slvPvtThd: " << slvPvtThd << endl;
|
|
cout << endl << "slvOffset: " << slvOffset << endl;
|
|
cout << endl << "slvScale: " << slvScale << endl;
|
|
}
|
|
};
|
|
|
|
virtual int initSolver(Eigen::SparseLU<EM, Eigen::COLAMDOrdering<int>> & solver) {
|
|
// Solve with arpack using sparse matrices and direct solvers.
|
|
|
|
solver.setPivotThreshold(slvPvtThd);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::SparseQR<EM, Eigen::COLAMDOrdering<int>> & solver) {
|
|
// Solve with arpack using sparse matrices and direct solvers.
|
|
|
|
solver.setPivotThreshold(slvPvtThd);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::SimplicialLLT<EM, Eigen::Lower, Eigen::COLAMDOrdering<int>> & solver) {
|
|
// Solve with arpack using sparse matrices and direct solvers.
|
|
|
|
solver.setShift(slvOffset, slvScale);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::SimplicialLDLT<EM, Eigen::Lower, Eigen::COLAMDOrdering<int>> & solver) {
|
|
// Solve with arpack using sparse matrices and direct solvers.
|
|
|
|
solver.setShift(slvOffset, slvScale);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::LLT<EM> & solver) {
|
|
// Solve with arpack using dense matrices and direct solvers.
|
|
|
|
if (this->mode == 1) return 0;
|
|
|
|
solver = Eigen::LLT<EM>(this->nbDim);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::LDLT<EM> & solver) {
|
|
// Solve with arpack using dense matrices and direct solvers.
|
|
|
|
if (this->mode == 1) return 0;
|
|
|
|
solver = Eigen::LDLT<EM>(this->nbDim);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::FullPivLU<EM> & solver) {
|
|
// Solve with arpack using dense matrices and direct solvers.
|
|
|
|
if (this->mode == 1) return 0;
|
|
|
|
solver = Eigen::FullPivLU<EM>(this->nbDim, this->nbDim);
|
|
solver.setThreshold(slvPvtThd);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::FullPivHouseholderQR<EM> & solver) {
|
|
// Solve with arpack using dense matrices and direct solvers.
|
|
|
|
if (this->mode == 1) return 0;
|
|
|
|
solver = Eigen::FullPivHouseholderQR<EM>(this->nbDim, this->nbDim);
|
|
solver.setThreshold(slvPvtThd);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::PartialPivLU<EM> & solver) {
|
|
// Solve with arpack using dense matrices and direct solvers.
|
|
|
|
if (this->mode == 1) return 0;
|
|
|
|
solver = Eigen::PartialPivLU<EM>(this->nbDim);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::HouseholderQR<EM> & solver) {
|
|
// Solve with arpack using dense matrices and direct solvers.
|
|
|
|
if (this->mode == 1) return 0;
|
|
|
|
solver = Eigen::HouseholderQR<EM>(this->nbDim, this->nbDim);
|
|
|
|
return 0;
|
|
};
|
|
|
|
// Public members.
|
|
|
|
public:
|
|
|
|
// Direct solvers parameters.
|
|
|
|
double slvPvtThd; // Pivoting tolerance of the direct mode solver.
|
|
double slvOffset; // Cholesky offset (LLT, LDLT) of the direct mode solver.
|
|
double slvScale; // Cholesky scale (LLT, LDLT) of the direct mode solver.
|
|
};
|
|
|
|
#endif
|
|
|
|
// Local Variables:
|
|
// mode: c++
|
|
// c-file-style:"stroustrup"
|
|
// show-trailing-whitespace: t
|
|
// End:
|
|
/* vim: set sw=2 ts=2 et smartindent :*/
|