Franck HOUSSEN
6e2893f17f
* TravisCI: fix ILP64 job. No need for security.debian.org. * Extract arpackSolver.hpp from arpackmm.cpp. arpackSolver.hpp is added in the install directory.
1040 lines
40 KiB
C++
1040 lines
40 KiB
C++
#ifndef __ARPACKSOLVER_HPP__
|
|
#define __ARPACKSOLVER_HPP__
|
|
|
|
#include "arpack.h"
|
|
|
|
#include <iostream>
|
|
#include <string>
|
|
#include <sstream> // stringstream.
|
|
#include <fstream> // [io]fstream.
|
|
#include <memory> // unique_ptr.
|
|
#include <chrono>
|
|
#include <iomanip> // setw.
|
|
#include <cassert>
|
|
|
|
#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;
|
|
|
|
// Vector related types.
|
|
|
|
typedef Eigen::Matrix< float, Eigen::Dynamic, 1> EigVecS; // Real.
|
|
typedef Eigen::Matrix< double, Eigen::Dynamic, 1> EigVecD; // Real.
|
|
typedef Eigen::Matrix<complex< float>, Eigen::Dynamic, 1> EigVecC; // Complex.
|
|
typedef Eigen::Matrix<complex<double>, Eigen::Dynamic, 1> EigVecZ; // Complex.
|
|
|
|
typedef Eigen::Map<EigVecS> EigMpVS; // Real.
|
|
typedef Eigen::Map<EigVecD> EigMpVD; // Real.
|
|
typedef Eigen::Map<EigVecC> EigMpVC; // Complex.
|
|
typedef Eigen::Map<EigVecZ> EigMpVZ; // Complex.
|
|
|
|
// Sparse matrix related types.
|
|
|
|
typedef Eigen::SparseMatrix< float> EigSMxS; // Real.
|
|
typedef Eigen::Triplet < float> EigSMTS; // Real.
|
|
typedef Eigen::SparseMatrix< double> EigSMxD; // Real.
|
|
typedef Eigen::Triplet < double> EigSMTD; // Real.
|
|
typedef Eigen::SparseMatrix<complex< float>> EigSMxC; // Complex.
|
|
typedef Eigen::Triplet <complex< float>> EigSMTC; // Complex.
|
|
typedef Eigen::SparseMatrix<complex<double>> EigSMxZ; // Complex.
|
|
typedef Eigen::Triplet <complex<double>> EigSMTZ; // 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.
|
|
|
|
template<typename RC, typename FD,
|
|
typename EM, typename ET, typename EV,
|
|
typename SLV>
|
|
class arpackSolver {
|
|
// Public methods.
|
|
|
|
public:
|
|
|
|
arpackSolver() {
|
|
stdPb = true;
|
|
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;
|
|
|
|
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.
|
|
|
|
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<ET> 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;
|
|
};
|
|
|
|
int solve(EM & A, EM & B) {
|
|
if (verbose == 3) {
|
|
cout << endl << "arpackSolver:" << endl;
|
|
cout << endl << "A:" << endl;
|
|
cout << endl << A << endl;
|
|
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<double>::epsilon();
|
|
bool shiftReal = (shiftReal && fabs(sigmaReal) > eps) ? true : false;
|
|
bool shiftImag = (shiftImag && fabs(sigmaImag) > eps) ? true : false;
|
|
bool backTransform = false;
|
|
mode = 0;
|
|
if (stdPb) {
|
|
mode = 1;
|
|
if (shiftReal && !shiftImag) {
|
|
EM I(A.rows(), A.cols());
|
|
I.setIdentity();
|
|
RC sigma; makeSigma(sigma);
|
|
A -= sigma*I;
|
|
backTransform = true;
|
|
}
|
|
}
|
|
else {
|
|
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) {
|
|
// 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;
|
|
}
|
|
}
|
|
|
|
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() > sqrt(tol)) {
|
|
cerr << endl << "Error: bad vector " << setw(3) << i << " (norm " << V.norm() << "):" << endl;
|
|
cerr << endl << V << endl;
|
|
cerr << endl << "Error: left side (A*V - norm " << left.norm() << "):" << endl;
|
|
cerr << endl << left << endl;
|
|
cerr << endl << "Error: right side (lambda*" << (stdPb ? "" : "B*") << "V - norm " << right.norm() << "):" << endl;
|
|
cerr << endl << right << endl;
|
|
cerr << endl << "Error: diff (norm " << diff.norm() << ", sqrt(tol) " << sqrt(tol) << "):" << endl;
|
|
cerr << endl << diff << 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() << ", sqrt(tol) " << sqrt(tol) << ")" << endl;
|
|
}
|
|
}
|
|
}
|
|
|
|
return 0;
|
|
};
|
|
|
|
// Private methods.
|
|
|
|
private:
|
|
|
|
void makeZero( float & zero) {zero = 0.f;};
|
|
|
|
void makeZero( double & zero) {zero = 0.;};
|
|
|
|
void makeZero(complex< float> & zero) {zero = complex<double>(0.f, 0.f);};
|
|
|
|
void makeZero(complex<double> & zero) {zero = complex<double>(0., 0.);};
|
|
|
|
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] == '%') 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) { // OK, (optional) nnz has been provided.
|
|
i.reserve(nnz);
|
|
j.reserve(nnz);
|
|
Mij.reserve(nnz);
|
|
}
|
|
}
|
|
}
|
|
else { // Body.
|
|
a_uint k = 0, l = 0;
|
|
RC zero; makeZero(zero);
|
|
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.
|
|
|
|
nnz = i.size(); // In case nnz was not provided.
|
|
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;
|
|
}
|
|
|
|
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<_Complex float*>(resid), nbCV,
|
|
reinterpret_cast<_Complex float*>(v), ldv, iparam, ipntr, reinterpret_cast<_Complex float*>(workd),
|
|
reinterpret_cast<_Complex float*>(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<_Complex double*>(resid), nbCV,
|
|
reinterpret_cast<_Complex double*>(v), ldv, iparam, ipntr, reinterpret_cast<_Complex double*>(workd),
|
|
reinterpret_cast<_Complex double*>(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(bool 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: 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 {
|
|
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: dneupd - KO: [dz]naupd 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(bool 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: [dz]naupd 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(bool 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<_Complex float*>(d), reinterpret_cast<_Complex float*>(z), ldz,
|
|
reinterpret_cast<_Complex float &>(sigma), reinterpret_cast<_Complex float*>(workev),
|
|
bMat, nbDim, which, nbEV, tol, reinterpret_cast<_Complex float*>(resid), nbCV,
|
|
reinterpret_cast<_Complex float*>(v), ldv, iparam, ipntr,
|
|
reinterpret_cast<_Complex float*>(workd), reinterpret_cast<_Complex float*>(workl), lworkl, rwork, &info);
|
|
if (info == -14) cerr << "Error: zneupd - KO: dsaupd 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 eupd(bool 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<_Complex double*>(d), reinterpret_cast<_Complex double*>(z), ldz,
|
|
reinterpret_cast<_Complex double &>(sigma), reinterpret_cast<_Complex double*>(workev),
|
|
bMat, nbDim, which, nbEV, tol, reinterpret_cast<_Complex double*>(resid), nbCV,
|
|
reinterpret_cast<_Complex double*>(v), ldv, iparam, ipntr,
|
|
reinterpret_cast<_Complex double*>(workd), reinterpret_cast<_Complex double*>(workl), lworkl, rwork, &info);
|
|
if (info == -14) cerr << "Error: zneupd - KO: dsaupd 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 (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;
|
|
};
|
|
|
|
void restartSolve(string const & fileName,
|
|
a_int const & nbDim, RC * rv, bool readNbCV) {
|
|
if (restartFromFile) {
|
|
ifstream ifs(fileName.c_str());
|
|
if (ifs.is_open()) {
|
|
a_int nbCV = 1; if (readNbCV) {ifs >> nbCV; if (nbCV < nbCV) nbCV = nbCV;}
|
|
for (a_int n = 0; rv && n < nbDim*nbCV; n++) ifs >> rv[n];
|
|
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;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
};
|
|
|
|
int solve(EM const & A, EM const & B, SLV & solver) {
|
|
// 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();
|
|
RC zero; makeZero(zero);
|
|
if (!resid) {
|
|
resid = new RC[nbDim];
|
|
for (a_int n = 0; n < nbDim; n++) resid[n] = zero; // Avoid "bad" starting vector.
|
|
};
|
|
restartSolve("arpackSolver.resid.out", nbDim, resid, false);
|
|
a_int ldv = nbDim;
|
|
if (!v) {
|
|
v = new RC[ldv*nbCV];
|
|
for (a_int n = 0; n < ldv*nbCV; n++) v[n] = zero; // Avoid "bad" starting vector.
|
|
};
|
|
restartSolve("arpackSolver.v.out", ldv, v, true);
|
|
a_int iparam[11];
|
|
iparam[0] = 1; // Use exact shifts (=> we'll never have ido == 3).
|
|
iparam[2] = maxIt; // Maximum number of iterations.
|
|
iparam[3] = 1; // Block size.
|
|
iparam[4] = 0; // Number of ev found by arpack.
|
|
iparam[6] = mode;
|
|
a_int ipntr[14];
|
|
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.
|
|
if (restartFromFile) info = 1;
|
|
|
|
// Initialize solver.
|
|
|
|
auto start = chrono::high_resolution_clock::now();
|
|
int 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, iparam, ipntr, 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 " << iparam[2] << endl; return 1;}
|
|
|
|
// Reverse Communication Interface: perform actions according to arpack.
|
|
|
|
start = chrono::high_resolution_clock::now();
|
|
|
|
a_int xIdx = ipntr[0] - 1; // 0-based (Fortran is 1-based).
|
|
a_int yIdx = ipntr[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 (iparam[6] == 1) {
|
|
Y = A * X;
|
|
}
|
|
else if (iparam[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 (iparam[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 (iparam[6] == 1) {
|
|
Y = A * X;
|
|
}
|
|
else if (iparam[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 (iparam[6] == 3) {
|
|
a_int zIdx = ipntr[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 (iparam[6] == 1) Y = X; // Y = I * X.
|
|
else if (iparam[6] == 2) Y = B * X; // Y = B * X.
|
|
else if (iparam[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 = iparam[2]; // Actual number of iterations.
|
|
bool rvec = true;
|
|
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 = eupd(rvec, howmny, select, z, ldz, bMat, nbDim, which, resid, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info);
|
|
if (rc != 0) {cerr << "Error: bad arpack eupd" << endl; return rc;}
|
|
|
|
if (dumpToFile) {
|
|
ofstream rfs("arpackSolver.resid.out"); for (a_int n = 0; n < nbDim; n++) rfs << resid[n] << endl;
|
|
ofstream vfs("arpackSolver.v.out"); vfs << nbCV << endl; for (a_int n = 0; n < ldv*nbCV; n++) vfs << v[n] << endl;
|
|
}
|
|
|
|
// 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;
|
|
|
|
// Public members.
|
|
|
|
public:
|
|
|
|
// Arpack parameters.
|
|
|
|
bool stdPb; // Standard or generalized (= not standard).
|
|
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.
|
|
|
|
vector<complex<double>> val; // Eigen values.
|
|
vector<EigVecZ> 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.
|
|
};
|
|
|
|
template<typename RC, typename FD,
|
|
typename EM, typename ET, typename EV,
|
|
typename SLV>
|
|
class arpackItrSolver: public arpackSolver<RC, FD, EM, ET, EV, SLV> {
|
|
// Public methods.
|
|
|
|
public:
|
|
|
|
arpackItrSolver(): arpackSolver<RC, FD, EM, ET, EV, SLV>() {
|
|
slvItrTol = unique_ptr<double>(new double); *slvItrTol = 1.e-6;
|
|
slvItrMaxIt = unique_ptr<int>(new int); *slvItrMaxIt = 100;
|
|
slvItrILUDropTol = unique_ptr<double>(new double); *slvItrILUDropTol = 1.;
|
|
slvItrILUFillFactor = unique_ptr<int>(new int); *slvItrILUFillFactor = 2;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::BiCGSTAB<EM> & solver) {
|
|
// Solve with arpack using sparse matrices and iterative solvers.
|
|
|
|
if (slvItrTol) solver.setTolerance(*slvItrTol);
|
|
if (slvItrMaxIt) solver.setMaxIterations(*slvItrMaxIt);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::ConjugateGradient<EM> & solver) {
|
|
// Solve with arpack using sparse matrices and iterative solvers.
|
|
|
|
if (slvItrTol) solver.setTolerance(*slvItrTol);
|
|
if (slvItrMaxIt) solver.setMaxIterations(*slvItrMaxIt);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::BiCGSTAB<EM, Eigen::IncompleteLUT<RC>> & solver) {
|
|
// Solve with arpack using sparse matrices and iterative solvers.
|
|
|
|
if (slvItrTol) solver.setTolerance(*slvItrTol);
|
|
if (slvItrMaxIt) solver.setMaxIterations(*slvItrMaxIt);
|
|
if (slvItrILUDropTol) solver.preconditioner().setDroptol(*slvItrILUDropTol);
|
|
if (slvItrILUFillFactor) solver.preconditioner().setFillfactor(*slvItrILUFillFactor);
|
|
|
|
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.
|
|
|
|
if (slvItrTol) solver.setTolerance(*slvItrTol);
|
|
if (slvItrMaxIt) solver.setMaxIterations(*slvItrMaxIt);
|
|
if (slvItrILUDropTol) solver.preconditioner().setDroptol(*slvItrILUDropTol);
|
|
if (slvItrILUFillFactor) solver.preconditioner().setFillfactor(*slvItrILUFillFactor);
|
|
|
|
return 0;
|
|
};
|
|
|
|
// Public members.
|
|
|
|
public:
|
|
|
|
// Iterative solvers parameters.
|
|
|
|
unique_ptr<double> slvItrTol;
|
|
unique_ptr<int> slvItrMaxIt;
|
|
unique_ptr<double> slvItrILUDropTol;
|
|
unique_ptr<int> slvItrILUFillFactor;
|
|
};
|
|
|
|
template<typename RC, typename FD,
|
|
typename EM, typename ET, typename EV,
|
|
typename SLV>
|
|
class arpackDrtSolver: public arpackSolver<RC, FD, EM, ET, EV, SLV> {
|
|
// Public methods.
|
|
|
|
public:
|
|
|
|
arpackDrtSolver(): arpackSolver<RC, FD, EM, ET, EV, SLV>() {
|
|
slvDrtPvtThd = unique_ptr<double>(new double); *slvDrtPvtThd = 1.e-6;
|
|
slvDrtOffset = unique_ptr<double>(new double); *slvDrtOffset = 0.;
|
|
slvDrtScale = unique_ptr<double>(new double); *slvDrtScale = 1.;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::SparseLU<EM, Eigen::COLAMDOrdering<int>> & solver) {
|
|
// Solve with arpack using sparse matrices and direct solvers.
|
|
|
|
if (slvDrtPvtThd) solver.setPivotThreshold(*slvDrtPvtThd);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::SparseQR<EM, Eigen::COLAMDOrdering<int>> & solver) {
|
|
// Solve with arpack using sparse matrices and direct solvers.
|
|
|
|
if (slvDrtPvtThd) solver.setPivotThreshold(*slvDrtPvtThd);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::SimplicialLLT<EM, Eigen::Lower, Eigen::COLAMDOrdering<int>> & solver) {
|
|
// Solve with arpack using sparse matrices and direct solvers.
|
|
|
|
if (slvDrtOffset && slvDrtScale) solver.setShift(*slvDrtOffset, *slvDrtScale);
|
|
|
|
return 0;
|
|
};
|
|
|
|
virtual int initSolver(Eigen::SimplicialLDLT<EM, Eigen::Lower, Eigen::COLAMDOrdering<int>> & solver) {
|
|
// Solve with arpack using sparse matrices and direct solvers.
|
|
|
|
if (slvDrtOffset && slvDrtScale) solver.setShift(*slvDrtOffset, *slvDrtScale);
|
|
|
|
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(Eigen::Default);
|
|
if (slvDrtPvtThd) solver.setThreshold(*slvDrtPvtThd);
|
|
|
|
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(Eigen::Default);
|
|
if (slvDrtPvtThd) solver.setThreshold(*slvDrtPvtThd);
|
|
|
|
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.
|
|
|
|
unique_ptr<double> slvDrtPvtThd;
|
|
unique_ptr<double> slvDrtOffset;
|
|
unique_ptr<double> slvDrtScale;
|
|
};
|
|
|
|
#endif
|
|
|
|
// Local Variables:
|
|
// mode: c++
|
|
// c-file-style:"stroustrup"
|
|
// show-trailing-whitespace: t
|
|
// End:
|
|
/* vim: set sw=2 ts=2 et smartindent :*/
|