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

Update the coding style (#256)

Browse Source 
* clang-format the C++ codebase

* clang-format the C code
This commit is contained in:
Sylvestre Ledru
2020-03-21 18:47:59 +01:00
committed by GitHub
parent 4a12554877
commit 6dcecc8eef
10 changed files with 1247 additions and 817 deletions
Download Patch File Download Diff File Expand all files Collapse all files
EXAMPLES/MATRIX_MARKET/arpackmm.cpp
+738 -418
View File
File diff suppressed because it is too large Load Diff
EXAMPLES/PYARPACK/pyarpack.cpp
+246 -147
View File
@@ -1,86 +1,124 @@
#include <string>
#include <vector>
#include <sstream> // ostringstream.
#include <complex>
#include <Python.h> // PyErr_SetString.
#include <boost/python.hpp>
#include <boost/python/numpy.hpp>
#include <boost/python/suite/indexing/vector_indexing_suite.hpp>
#include <Python.h> // PyErr_SetString.
#include <pyarpackItrSolver.hpp>
#include <complex>
#include <pyarpackDrtSolver.hpp>
#include <pyarpackItrSolver.hpp>
#include <sstream> // ostringstream.
#include <string>
#include <vector>
namespace bp = boost::python;
namespace bn = boost::python::numpy;
template<typename RC, typename FD, typename EM, typename SLV>
void exportArpackSparseItr(bp::scope & pySlv, std::string const & dtype) {
template <typename RC, typename FD, typename EM, typename SLV>
void exportArpackSparseItr(bp::scope& pySlv, std::string const& dtype) {
// Created nested namespace in module.
pySlv.attr(dtype.c_str()) = bp::class_<pyarpackSparseItrSolver<RC, FD, EM, SLV>>(dtype.c_str(),
"arpack data type (must be consistent with numpy dtype)")
.def ("solve", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::solve,
(bp::arg("A"), bp::arg("B") = bp::tuple()),
"solve standard or generalised eigen problem where A and B must be sparse and provided in coo format: (dimension, row-indice array, column-indice array, matrice-value array) tuple")
.def ("checkEigVec", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::checkEigVec,
(bp::arg("A"), bp::arg("B") = bp::tuple(), bp::arg("diffTol") = 1.e-3),
"check eigen vectors accuracy where A and B must be sparse and provided in coo format: (dimension, row-indice array, column-indice array, matrice-value array) tuple")
ARPACKSOLVERMEMBER(pyarpackSparseItrSolver)
.def_readwrite("slvTol", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::slvTol,
"tolerance of the iterative mode solver - default: 1.e-6")
.def_readwrite("slvMaxIt", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::slvMaxIt,
"maximum number of iterations of the iterative mode solver - default: 100")
.def_readwrite("slvILUDropTol", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::slvILUDropTol,
"drop tolerance of the ILU preconditioner (if any) of the iterative mode solver - default: 1")
.def_readwrite("slvILUFillFactor", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::slvILUFillFactor,
"fill factor of the ILU preconditioner (if any) of the iterative mode solver - default: 2")
;
pySlv.attr(dtype.c_str()) =
bp::class_<pyarpackSparseItrSolver<RC, FD, EM, SLV>>(
dtype.c_str(),
"arpack data type (must be consistent with numpy dtype)")
.def("solve", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::solve,
(bp::arg("A"), bp::arg("B") = bp::tuple()),
"solve standard or generalised eigen problem where A and B must "
"be sparse and provided in coo format: (dimension, row-indice "
"array, column-indice array, matrice-value array) tuple")
.def("checkEigVec",
&pyarpackSparseItrSolver<RC, FD, EM, SLV>::checkEigVec,
(bp::arg("A"), bp::arg("B") = bp::tuple(),
bp::arg("diffTol") = 1.e-3),
"check eigen vectors accuracy where A and B must be sparse and "
"provided in coo format: (dimension, row-indice array, "
"column-indice array, matrice-value array) tuple")
ARPACKSOLVERMEMBER(pyarpackSparseItrSolver)
.def_readwrite(
"slvTol", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::slvTol,
"tolerance of the iterative mode solver - default: 1.e-6")
.def_readwrite("slvMaxIt",
&pyarpackSparseItrSolver<RC, FD, EM, SLV>::slvMaxIt,
"maximum number of iterations of the iterative mode "
"solver - default: 100")
.def_readwrite(
"slvILUDropTol",
&pyarpackSparseItrSolver<RC, FD, EM, SLV>::slvILUDropTol,
"drop tolerance of the ILU preconditioner (if any) of the "
"iterative mode solver - default: 1")
.def_readwrite(
"slvILUFillFactor",
&pyarpackSparseItrSolver<RC, FD, EM, SLV>::slvILUFillFactor,
"fill factor of the ILU preconditioner (if any) of the iterative "
"mode solver - default: 2");
};
template<typename RC, typename FD, typename EM, typename SLV>
void exportArpackSparseDrt(bp::scope & pySlv, std::string const & dtype) {
template <typename RC, typename FD, typename EM, typename SLV>
void exportArpackSparseDrt(bp::scope& pySlv, std::string const& dtype) {
// Created nested namespace in module.
pySlv.attr(dtype.c_str()) = bp::class_<pyarpackSparseDrtSolver<RC, FD, EM, SLV>>(dtype.c_str(),
"arpack data type (must be consistent with numpy dtype)")
.def ("solve", &pyarpackSparseDrtSolver<RC, FD, EM, SLV>::solve,
(bp::arg("A"), bp::arg("B") = bp::tuple()),
"solve standard or generalised eigen problem where A and B must be sparse and provided in coo format: (dimension, row-indice array, column-indice array, matrice-value array) tuple")
.def ("checkEigVec", &pyarpackSparseDrtSolver<RC, FD, EM, SLV>::checkEigVec,
(bp::arg("A"), bp::arg("B") = bp::tuple(), bp::arg("diffTol") = 1.e-3),
"check eigen vectors accuracy where A and B must be sparse and provided in coo format: (dimension, row-indice array, column-indice array, matrice-value array) tuple")
ARPACKSOLVERMEMBER(pyarpackSparseDrtSolver)
.def_readwrite("slvPvtThd", &pyarpackSparseDrtSolver<RC, FD, EM, SLV>::slvPvtThd,
"pivoting tolerance of the direct mode solver - default: 1.e-6")
.def_readwrite("slvOffset", &pyarpackSparseDrtSolver<RC, FD, EM, SLV>::slvOffset,
"cholesky offset (LLT, LDLT) of the direct mode solver - default: 0.")
.def_readwrite("slvScale", &pyarpackSparseDrtSolver<RC, FD, EM, SLV>::slvScale,
"cholesky scale (LLT, LDLT) of the direct mode solver - default: 1.")
;
pySlv.attr(dtype.c_str()) =
bp::class_<pyarpackSparseDrtSolver<RC, FD, EM, SLV>>(
dtype.c_str(),
"arpack data type (must be consistent with numpy dtype)")
.def("solve", &pyarpackSparseDrtSolver<RC, FD, EM, SLV>::solve,
(bp::arg("A"), bp::arg("B") = bp::tuple()),
"solve standard or generalised eigen problem where A and B must "
"be sparse and provided in coo format: (dimension, row-indice "
"array, column-indice array, matrice-value array) tuple")
.def("checkEigVec",
&pyarpackSparseDrtSolver<RC, FD, EM, SLV>::checkEigVec,
(bp::arg("A"), bp::arg("B") = bp::tuple(),
bp::arg("diffTol") = 1.e-3),
"check eigen vectors accuracy where A and B must be sparse and "
"provided in coo format: (dimension, row-indice array, "
"column-indice array, matrice-value array) tuple")
ARPACKSOLVERMEMBER(pyarpackSparseDrtSolver)
.def_readwrite(
"slvPvtThd", &pyarpackSparseDrtSolver<RC, FD, EM, SLV>::slvPvtThd,
"pivoting tolerance of the direct mode solver - default: 1.e-6")
.def_readwrite("slvOffset",
&pyarpackSparseDrtSolver<RC, FD, EM, SLV>::slvOffset,
"cholesky offset (LLT, LDLT) of the direct mode "
"solver - default: 0.")
.def_readwrite("slvScale",
&pyarpackSparseDrtSolver<RC, FD, EM, SLV>::slvScale,
"cholesky scale (LLT, LDLT) of the direct mode solver "
"- default: 1.");
};
template<typename RC, typename FD, typename EM, typename SLV>
void exportArpackDenseDrt(bp::scope & pySlv, std::string const & dtype) {
template <typename RC, typename FD, typename EM, typename SLV>
void exportArpackDenseDrt(bp::scope& pySlv, std::string const& dtype) {
// Created nested namespace in module.
pySlv.attr(dtype.c_str()) = bp::class_<pyarpackDenseDrtSolver<RC, FD, EM, SLV>>(dtype.c_str(),
"arpack data type (must be consistent with numpy dtype)")
.def ("solve", &pyarpackDenseDrtSolver<RC, FD, EM, SLV>::solve,
(bp::arg("A"), bp::arg("B") = bp::tuple()),
"solve standard or generalised eigen problem where A and B must be dense and provided in raw format: (n-squared matrice-value array, row or column ordered boolean)")
.def ("checkEigVec", &pyarpackDenseDrtSolver<RC, FD, EM, SLV>::checkEigVec,
(bp::arg("A"), bp::arg("B") = bp::tuple(), bp::arg("diffTol") = 1.e-3),
"check eigen vectors accuracy where A and B must be dense and provided in raw format: (n-squared matrice-value array, row or column ordered boolean)")
ARPACKSOLVERMEMBER(pyarpackDenseDrtSolver)
.def_readwrite("slvPvtThd", &pyarpackDenseDrtSolver<RC, FD, EM, SLV>::slvPvtThd,
"pivoting tolerance of the direct mode solver - default: 1.e-6")
.def_readwrite("slvOffset", &pyarpackDenseDrtSolver<RC, FD, EM, SLV>::slvOffset,
"cholesky offset (LLT, LDLT) of the direct mode solver - default: 0.")
.def_readwrite("slvScale", &pyarpackDenseDrtSolver<RC, FD, EM, SLV>::slvScale,
"cholesky scale (LLT, LDLT) of the direct mode solver - default: 1.")
;
pySlv.attr(dtype.c_str()) =
bp::class_<pyarpackDenseDrtSolver<RC, FD, EM, SLV>>(
dtype.c_str(),
"arpack data type (must be consistent with numpy dtype)")
.def("solve", &pyarpackDenseDrtSolver<RC, FD, EM, SLV>::solve,
(bp::arg("A"), bp::arg("B") = bp::tuple()),
"solve standard or generalised eigen problem where A and B must "
"be dense and provided in raw format: (n-squared matrice-value "
"array, row or column ordered boolean)")
.def("checkEigVec",
&pyarpackDenseDrtSolver<RC, FD, EM, SLV>::checkEigVec,
(bp::arg("A"), bp::arg("B") = bp::tuple(),
bp::arg("diffTol") = 1.e-3),
"check eigen vectors accuracy where A and B must be dense and "
"provided in raw format: (n-squared matrice-value array, row or "
"column ordered boolean)")
ARPACKSOLVERMEMBER(pyarpackDenseDrtSolver)
.def_readwrite(
"slvPvtThd", &pyarpackDenseDrtSolver<RC, FD, EM, SLV>::slvPvtThd,
"pivoting tolerance of the direct mode solver - default: 1.e-6")
.def_readwrite("slvOffset",
&pyarpackDenseDrtSolver<RC, FD, EM, SLV>::slvOffset,
"cholesky offset (LLT, LDLT) of the direct mode "
"solver - default: 0.")
.def_readwrite("slvScale",
&pyarpackDenseDrtSolver<RC, FD, EM, SLV>::slvScale,
"cholesky scale (LLT, LDLT) of the direct mode solver "
"- default: 1.");
};
class sparseBiCGDiag {};
@@ -99,36 +137,36 @@ class denseQRRR {};
class denseLUPP {};
class denseQRPP {};
std::complex<double> EigVecZGetItem(Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> & M, int idx) {
if (idx < 0 || idx >= M.size()) {pyarpackThrowError("index out of range"); return std::complex<double>();}
std::complex<double> EigVecZGetItem(
Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1>& M, int idx) {
if (idx < 0 || idx >= M.size()) {
pyarpackThrowError("index out of range");
return std::complex<double>();
}
return M[idx];
};
std::string EigVecZToString(EigVecZ const & vec) {
std::string EigVecZToString(EigVecZ const& vec) {
std::ostringstream s;
s << vec;
return s.str();
};
BOOST_PYTHON_MODULE(pyarpack)
{
BOOST_PYTHON_MODULE(pyarpack) {
// Initialize.
bn::initialize();
bp::class_<std::vector<std::complex<double>>>("StdVecZ")
.def(bp::vector_indexing_suite<std::vector<std::complex<double>>>())
;
bp::class_<std::vector<std::complex<double>>>("StdVecZ").def(
bp::vector_indexing_suite<std::vector<std::complex<double>>>());
bp::class_<Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1>>("EigVecZ")
.def("__getitem__", &EigVecZGetItem)
.def("__str__", &EigVecZToString)
;
.def("__getitem__", &EigVecZGetItem)
.def("__str__", &EigVecZToString);
bp::class_<std::vector<EigVecZ>>("StdVecEVZ")
.def("__iter__", bp::iterator<std::vector<EigVecZ>>())
.def(bp::vector_indexing_suite<std::vector<EigVecZ>>())
;
.def("__iter__", bp::iterator<std::vector<EigVecZ>>())
.def(bp::vector_indexing_suite<std::vector<EigVecZ>>());
// Documentation of the python module.
@@ -179,148 +217,209 @@ BOOST_PYTHON_MODULE(pyarpack)
// Create python module.
std::string module = "pyarpack";
bp::object pyModule(bp::handle<>(bp::borrowed(PyImport_AddModule(module.c_str()))));
bp::object pyModule(
bp::handle<>(bp::borrowed(PyImport_AddModule(module.c_str()))));
// Create modules.
{
std::string slv = "sparseBiCGDiag";
std::string slvHelp = "arpack internal mode solver (mode > 1): BiCG with diagonal (Jacobi) preconditioner";
bp::scope pySlvBiCGDiag = bp::class_<sparseBiCGDiag>(slv.c_str(), slvHelp.c_str());
exportArpackSparseItr< float , float, EigSMxS, EigSBiCGS>(pySlvBiCGDiag, "float" );
exportArpackSparseItr< double , double, EigSMxD, EigSBiCGD>(pySlvBiCGDiag, "double");
exportArpackSparseItr<std::complex< float>, float, EigSMxC, EigSBiCGC>(pySlvBiCGDiag, "complexFloat" );
exportArpackSparseItr<std::complex<double>, double, EigSMxZ, EigSBiCGZ>(pySlvBiCGDiag, "complexDouble");
std::string slvHelp =
"arpack internal mode solver (mode > 1): BiCG with diagonal (Jacobi) "
"preconditioner";
bp::scope pySlvBiCGDiag =
bp::class_<sparseBiCGDiag>(slv.c_str(), slvHelp.c_str());
exportArpackSparseItr<float, float, EigSMxS, EigSBiCGS>(pySlvBiCGDiag,
"float");
exportArpackSparseItr<double, double, EigSMxD, EigSBiCGD>(pySlvBiCGDiag,
"double");
exportArpackSparseItr<std::complex<float>, float, EigSMxC, EigSBiCGC>(
pySlvBiCGDiag, "complexFloat");
exportArpackSparseItr<std::complex<double>, double, EigSMxZ, EigSBiCGZ>(
pySlvBiCGDiag, "complexDouble");
}
{
std::string slv = "sparseBiCGILU";
std::string slvHelp = "arpack internal mode solver (mode > 1): BiCG with ILU preconditioner";
bp::scope pySlvBiCGILU = bp::class_<sparseBiCGILU>(slv.c_str(), slvHelp.c_str());
exportArpackSparseItr< float , float, EigSMxS, EigSBiCGILUS>(pySlvBiCGILU, "float" );
exportArpackSparseItr< double , double, EigSMxD, EigSBiCGILUD>(pySlvBiCGILU, "double");
exportArpackSparseItr<std::complex< float>, float, EigSMxC, EigSBiCGILUC>(pySlvBiCGILU, "complexFloat" );
exportArpackSparseItr<std::complex<double>, double, EigSMxZ, EigSBiCGILUZ>(pySlvBiCGILU, "complexDouble");
std::string slvHelp =
"arpack internal mode solver (mode > 1): BiCG with ILU preconditioner";
bp::scope pySlvBiCGILU =
bp::class_<sparseBiCGILU>(slv.c_str(), slvHelp.c_str());
exportArpackSparseItr<float, float, EigSMxS, EigSBiCGILUS>(pySlvBiCGILU,
"float");
exportArpackSparseItr<double, double, EigSMxD, EigSBiCGILUD>(pySlvBiCGILU,
"double");
exportArpackSparseItr<std::complex<float>, float, EigSMxC, EigSBiCGILUC>(
pySlvBiCGILU, "complexFloat");
exportArpackSparseItr<std::complex<double>, double, EigSMxZ, EigSBiCGILUZ>(
pySlvBiCGILU, "complexDouble");
}
{
std::string slv = "sparseCGDiag";
std::string slvHelp = "arpack internal mode solver (mode > 1): CG with diagonal (Jacobi) preconditioner";
bp::scope pySlvCGDiag = bp::class_<sparseCGDiag>(slv.c_str(), slvHelp.c_str());
exportArpackSparseItr< float , float, EigSMxS, EigSCGS>(pySlvCGDiag, "float" );
exportArpackSparseItr< double , double, EigSMxD, EigSCGD>(pySlvCGDiag, "double");
exportArpackSparseItr<std::complex< float>, float, EigSMxC, EigSCGC>(pySlvCGDiag, "complexFloat" );
exportArpackSparseItr<std::complex<double>, double, EigSMxZ, EigSCGZ>(pySlvCGDiag, "complexDouble");
std::string slvHelp =
"arpack internal mode solver (mode > 1): CG with diagonal (Jacobi) "
"preconditioner";
bp::scope pySlvCGDiag =
bp::class_<sparseCGDiag>(slv.c_str(), slvHelp.c_str());
exportArpackSparseItr<float, float, EigSMxS, EigSCGS>(pySlvCGDiag, "float");
exportArpackSparseItr<double, double, EigSMxD, EigSCGD>(pySlvCGDiag,
"double");
exportArpackSparseItr<std::complex<float>, float, EigSMxC, EigSCGC>(
pySlvCGDiag, "complexFloat");
exportArpackSparseItr<std::complex<double>, double, EigSMxZ, EigSCGZ>(
pySlvCGDiag, "complexDouble");
}
{
std::string slv = "sparseCGILU";
std::string slvHelp = "arpack internal mode solver (mode > 1): CG with ILU preconditioner";
bp::scope pySlvCGILU = bp::class_<sparseCGILU>(slv.c_str(), slvHelp.c_str());
exportArpackSparseItr< float , float, EigSMxS, EigSCGILUS>(pySlvCGILU, "float" );
exportArpackSparseItr< double , double, EigSMxD, EigSCGILUD>(pySlvCGILU, "double");
exportArpackSparseItr<std::complex< float>, float, EigSMxC, EigSCGILUC>(pySlvCGILU, "complexFloat" );
exportArpackSparseItr<std::complex<double>, double, EigSMxZ, EigSCGILUZ>(pySlvCGILU, "complexDouble");
std::string slvHelp =
"arpack internal mode solver (mode > 1): CG with ILU preconditioner";
bp::scope pySlvCGILU =
bp::class_<sparseCGILU>(slv.c_str(), slvHelp.c_str());
exportArpackSparseItr<float, float, EigSMxS, EigSCGILUS>(pySlvCGILU,
"float");
exportArpackSparseItr<double, double, EigSMxD, EigSCGILUD>(pySlvCGILU,
"double");
exportArpackSparseItr<std::complex<float>, float, EigSMxC, EigSCGILUC>(
pySlvCGILU, "complexFloat");
exportArpackSparseItr<std::complex<double>, double, EigSMxZ, EigSCGILUZ>(
pySlvCGILU, "complexDouble");
}
{
std::string slv = "sparseLLT";
std::string slvHelp = "arpack internal mode solver (mode > 1): LLT";
bp::scope pySlvLLT = bp::class_<sparseLLT>(slv.c_str(), slvHelp.c_str());
exportArpackSparseDrt< float , float, EigSMxS, EigSLLTS>(pySlvLLT, "float" );
exportArpackSparseDrt< double , double, EigSMxD, EigSLLTD>(pySlvLLT, "double");
exportArpackSparseDrt<std::complex< float>, float, EigSMxC, EigSLLTC>(pySlvLLT, "complexFloat" );
exportArpackSparseDrt<std::complex<double>, double, EigSMxZ, EigSLLTZ>(pySlvLLT, "complexDouble");
exportArpackSparseDrt<float, float, EigSMxS, EigSLLTS>(pySlvLLT, "float");
exportArpackSparseDrt<double, double, EigSMxD, EigSLLTD>(pySlvLLT,
"double");
exportArpackSparseDrt<std::complex<float>, float, EigSMxC, EigSLLTC>(
pySlvLLT, "complexFloat");
exportArpackSparseDrt<std::complex<double>, double, EigSMxZ, EigSLLTZ>(
pySlvLLT, "complexDouble");
}
{
std::string slv = "sparseLDLT";
std::string slvHelp = "arpack internal mode solver (mode > 1): LDLT";
bp::scope pySlvLDLT = bp::class_<sparseLDLT>(slv.c_str(), slvHelp.c_str());
exportArpackSparseDrt< float , float, EigSMxS, EigSLDLTS>(pySlvLDLT, "float" );
exportArpackSparseDrt< double , double, EigSMxD, EigSLDLTD>(pySlvLDLT, "double");
exportArpackSparseDrt<std::complex< float>, float, EigSMxC, EigSLDLTC>(pySlvLDLT, "complexFloat" );
exportArpackSparseDrt<std::complex<double>, double, EigSMxZ, EigSLDLTZ>(pySlvLDLT, "complexDouble");
exportArpackSparseDrt<float, float, EigSMxS, EigSLDLTS>(pySlvLDLT, "float");
exportArpackSparseDrt<double, double, EigSMxD, EigSLDLTD>(pySlvLDLT,
"double");
exportArpackSparseDrt<std::complex<float>, float, EigSMxC, EigSLDLTC>(
pySlvLDLT, "complexFloat");
exportArpackSparseDrt<std::complex<double>, double, EigSMxZ, EigSLDLTZ>(
pySlvLDLT, "complexDouble");
}
{
std::string slv = "sparseLU";
std::string slvHelp = "arpack internal mode solver (mode > 1): LU";
bp::scope pySlvLU = bp::class_<sparseLU>(slv.c_str(), slvHelp.c_str());
exportArpackSparseDrt< float , float, EigSMxS, EigSLUS>(pySlvLU, "float" );
exportArpackSparseDrt< double , double, EigSMxD, EigSLUD>(pySlvLU, "double");
exportArpackSparseDrt<std::complex< float>, float, EigSMxC, EigSLUC>(pySlvLU, "complexFloat" );
exportArpackSparseDrt<std::complex<double>, double, EigSMxZ, EigSLUZ>(pySlvLU, "complexDouble");
exportArpackSparseDrt<float, float, EigSMxS, EigSLUS>(pySlvLU, "float");
exportArpackSparseDrt<double, double, EigSMxD, EigSLUD>(pySlvLU, "double");
exportArpackSparseDrt<std::complex<float>, float, EigSMxC, EigSLUC>(
pySlvLU, "complexFloat");
exportArpackSparseDrt<std::complex<double>, double, EigSMxZ, EigSLUZ>(
pySlvLU, "complexDouble");
}
{
std::string slv = "sparseQR";
std::string slvHelp = "arpack internal mode solver (mode > 1): QR";
bp::scope pySlvQR = bp::class_<sparseQR>(slv.c_str(), slvHelp.c_str());
exportArpackSparseDrt< float , float, EigSMxS, EigSQRS>(pySlvQR, "float" );
exportArpackSparseDrt< double , double, EigSMxD, EigSQRD>(pySlvQR, "double");
exportArpackSparseDrt<std::complex< float>, float, EigSMxC, EigSQRC>(pySlvQR, "complexFloat" );
exportArpackSparseDrt<std::complex<double>, double, EigSMxZ, EigSQRZ>(pySlvQR, "complexDouble");
exportArpackSparseDrt<float, float, EigSMxS, EigSQRS>(pySlvQR, "float");
exportArpackSparseDrt<double, double, EigSMxD, EigSQRD>(pySlvQR, "double");
exportArpackSparseDrt<std::complex<float>, float, EigSMxC, EigSQRC>(
pySlvQR, "complexFloat");
exportArpackSparseDrt<std::complex<double>, double, EigSMxZ, EigSQRZ>(
pySlvQR, "complexDouble");
}
{
std::string slv = "denseLLT";
std::string slvHelp = "arpack internal mode solver (mode > 1): LLT";
bp::scope pySlvLLT = bp::class_<denseLLT>(slv.c_str(), slvHelp.c_str());
exportArpackDenseDrt< float , float, EigDMxS, EigDLLTS>(pySlvLLT, "float" );
exportArpackDenseDrt< double , double, EigDMxD, EigDLLTD>(pySlvLLT, "double");
exportArpackDenseDrt<std::complex< float>, float, EigDMxC, EigDLLTC>(pySlvLLT, "complexFloat" );
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDLLTZ>(pySlvLLT, "complexDouble");
exportArpackDenseDrt<float, float, EigDMxS, EigDLLTS>(pySlvLLT, "float");
exportArpackDenseDrt<double, double, EigDMxD, EigDLLTD>(pySlvLLT, "double");
exportArpackDenseDrt<std::complex<float>, float, EigDMxC, EigDLLTC>(
pySlvLLT, "complexFloat");
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDLLTZ>(
pySlvLLT, "complexDouble");
}
{
std::string slv = "denseLDLT";
std::string slvHelp = "arpack internal mode solver (mode > 1): LDLT";
bp::scope pySlvLDLT = bp::class_<denseLDLT>(slv.c_str(), slvHelp.c_str());
exportArpackDenseDrt< float , float, EigDMxS, EigDLDLTS>(pySlvLDLT, "float" );
exportArpackDenseDrt< double , double, EigDMxD, EigDLDLTD>(pySlvLDLT, "double");
exportArpackDenseDrt<std::complex< float>, float, EigDMxC, EigDLDLTC>(pySlvLDLT, "complexFloat" );
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDLDLTZ>(pySlvLDLT, "complexDouble");
exportArpackDenseDrt<float, float, EigDMxS, EigDLDLTS>(pySlvLDLT, "float");
exportArpackDenseDrt<double, double, EigDMxD, EigDLDLTD>(pySlvLDLT,
"double");
exportArpackDenseDrt<std::complex<float>, float, EigDMxC, EigDLDLTC>(
pySlvLDLT, "complexFloat");
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDLDLTZ>(
pySlvLDLT, "complexDouble");
}
{
std::string slv = "denseLURR";
std::string slvHelp = "arpack internal mode solver (mode > 1): LU Rank Revealing (slower, more stable)";
std::string slvHelp =
"arpack internal mode solver (mode > 1): LU Rank Revealing (slower, "
"more stable)";
bp::scope pySlvLURR = bp::class_<denseLURR>(slv.c_str(), slvHelp.c_str());
exportArpackDenseDrt< float , float, EigDMxS, EigDFLUS>(pySlvLURR, "float" );
exportArpackDenseDrt< double , double, EigDMxD, EigDFLUD>(pySlvLURR, "double");
exportArpackDenseDrt<std::complex< float>, float, EigDMxC, EigDFLUC>(pySlvLURR, "complexFloat" );
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDFLUZ>(pySlvLURR, "complexDouble");
exportArpackDenseDrt<float, float, EigDMxS, EigDFLUS>(pySlvLURR, "float");
exportArpackDenseDrt<double, double, EigDMxD, EigDFLUD>(pySlvLURR,
"double");
exportArpackDenseDrt<std::complex<float>, float, EigDMxC, EigDFLUC>(
pySlvLURR, "complexFloat");
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDFLUZ>(
pySlvLURR, "complexDouble");
}
{
std::string slv = "denseQRRR";
std::string slvHelp = "arpack internal mode solver (mode > 1): QR Rank Revealing (slower, more stable)";
std::string slvHelp =
"arpack internal mode solver (mode > 1): QR Rank Revealing (slower, "
"more stable)";
bp::scope pySlvQRRR = bp::class_<denseQRRR>(slv.c_str(), slvHelp.c_str());
exportArpackDenseDrt< float , float, EigDMxS, EigDFQRS>(pySlvQRRR, "float" );
exportArpackDenseDrt< double , double, EigDMxD, EigDFQRD>(pySlvQRRR, "double");
exportArpackDenseDrt<std::complex< float>, float, EigDMxC, EigDFQRC>(pySlvQRRR, "complexFloat" );
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDFQRZ>(pySlvQRRR, "complexDouble");
exportArpackDenseDrt<float, float, EigDMxS, EigDFQRS>(pySlvQRRR, "float");
exportArpackDenseDrt<double, double, EigDMxD, EigDFQRD>(pySlvQRRR,
"double");
exportArpackDenseDrt<std::complex<float>, float, EigDMxC, EigDFQRC>(
pySlvQRRR, "complexFloat");
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDFQRZ>(
pySlvQRRR, "complexDouble");
}
{
std::string slv = "denseLUPP";
std::string slvHelp = "arpack internal mode solver (mode > 1): LU Partial Pivoting (faster, less stable)";
std::string slvHelp =
"arpack internal mode solver (mode > 1): LU Partial Pivoting (faster, "
"less stable)";
bp::scope pySlvLUPP = bp::class_<denseLUPP>(slv.c_str(), slvHelp.c_str());
exportArpackDenseDrt< float , float, EigDMxS, EigDPLUS>(pySlvLUPP, "float" );
exportArpackDenseDrt< double , double, EigDMxD, EigDPLUD>(pySlvLUPP, "double");
exportArpackDenseDrt<std::complex< float>, float, EigDMxC, EigDPLUC>(pySlvLUPP, "complexFloat" );
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDPLUZ>(pySlvLUPP, "complexDouble");
exportArpackDenseDrt<float, float, EigDMxS, EigDPLUS>(pySlvLUPP, "float");
exportArpackDenseDrt<double, double, EigDMxD, EigDPLUD>(pySlvLUPP,
"double");
exportArpackDenseDrt<std::complex<float>, float, EigDMxC, EigDPLUC>(
pySlvLUPP, "complexFloat");
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDPLUZ>(
pySlvLUPP, "complexDouble");
}
{
std::string slv = "denseQRPP";
std::string slvHelp = "arpack internal mode solver (mode > 1): QR Partial Pivoting (faster, less stable)";
std::string slvHelp =
"arpack internal mode solver (mode > 1): QR Partial Pivoting (faster, "
"less stable)";
bp::scope pySlvQPPR = bp::class_<denseQRPP>(slv.c_str(), slvHelp.c_str());
exportArpackDenseDrt< float , float, EigDMxS, EigDPQRS>(pySlvQPPR, "float" );
exportArpackDenseDrt< double , double, EigDMxD, EigDPQRD>(pySlvQPPR, "double");
exportArpackDenseDrt<std::complex< float>, float, EigDMxC, EigDPQRC>(pySlvQPPR, "complexFloat" );
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDPQRZ>(pySlvQPPR, "complexDouble");
exportArpackDenseDrt<float, float, EigDMxS, EigDPQRS>(pySlvQPPR, "float");
exportArpackDenseDrt<double, double, EigDMxD, EigDPQRD>(pySlvQPPR,
"double");
exportArpackDenseDrt<std::complex<float>, float, EigDMxC, EigDPQRC>(
pySlvQPPR, "complexFloat");
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDPQRZ>(
pySlvQPPR, "complexDouble");
}
}
EXAMPLES/README.CALLING-ARPACK-FROM-C-OR-CPP
+4 -2
View File
@@ -1,2 +1,4 @@
In ../TESTS, the file icb_arpack_c.c is an example of how to call arpack from C.
In ../TESTS, the file icb_arpack_cpp.cpp is an example of how to call arpack from C++.
In../ TESTS,
the file icb_arpack_c.c is an example of how to call arpack from C.In../
TESTS,
the file icb_arpack_cpp.cpp is an example of how to call arpack from C++.
PARPACK/EXAMPLES/MPI/README.CALLING-PARPACK-FROM-C-OR-CPP
+5 -2
View File
@@ -1,2 +1,5 @@
In PARPACK/TESTS/MPI, the file icb_parpack_c.c is an example of how to call parpack from C.
In PARPACK/TESTS/MPI, the file icb_parpack_cpp.cpp is an example of how to call parpack from C++.
In PARPACK / TESTS / MPI,
the file icb_parpack_c.c is an example of how to call parpack from
C.In PARPACK /
TESTS / MPI,
the file icb_parpack_cpp.cpp is an example of how to call parpack from C++.
PARPACK/TESTS/MPI/icb_parpack_c.c
+79 -75
View File
@@ -1,30 +1,31 @@
/*
* This example demonstrates the use of ISO_C_BINDING to call arpack (portability).
* IMPORTANT: MPI communicators MUST be passed from C to Fortran using MPI_Comm_c2f.
* This example demonstrates the use of ISO_C_BINDING to call arpack
* (portability). IMPORTANT: MPI communicators MUST be passed from C to Fortran
* using MPI_Comm_c2f.
*
* Just use arpack as you would have normally done, but, use *[ae]upd_c instead of *[ae]upd_.
* The main advantage is that compiler checks (arguments) are performed at build time.
* Note: to debug parpack, call debug_c.
* Just use arpack as you would have normally done, but, use *[ae]upd_c instead
* of *[ae]upd_. The main advantage is that compiler checks (arguments) are
* performed at build time. Note: to debug parpack, call debug_c.
*/
#include <complex.h> // creal, cimag.
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "debug_c.h" // debug parpack.
#include "mpi.h"
#include "parpack.h"
#include <complex.h> // creal, cimag.
#include "debug_c.h" // debug parpack.
#include "stat_c.h" // arpack statistics.
#include "stat_c.h" // arpack statistics.
/* test program to solve for the 9 largest eigenvalues of
* A*x = lambda*x where A is the diagonal matrix
* with entries 1000, 999, ... , 2, 1 on the diagonal.
* */
void dMatVec(double * x, double * y) {
void dMatVec(double* x, double* y) {
int i;
for ( i = 0; i < 1000; ++i)
y[i] = ((double) (i+1))*x[i];
for (i = 0; i < 1000; ++i) y[i] = ((double)(i + 1)) * x[i];
};
int ds() {
@@ -35,55 +36,54 @@ int ds() {
a_int nev = 3;
double tol = 0;
double resid[N];
a_int ncv = 2*nev+1;
double V[ncv*N];
a_int ncv = 2 * nev + 1;
double V[ncv * N];
a_int ldv = N;
a_int iparam[11];
a_int ipntr[14];
double workd[3*N];
double workd[3 * N];
a_int rvec = 1;
char howmny[] = "A";
double* d = (double*) malloc((nev+1)*sizeof(double));
double* d = (double*)malloc((nev + 1) * sizeof(double));
a_int select[ncv];
for (int i = 0; i < ncv; i++) select[i] = 1;
double z[(N+1)*(nev+1)];
a_int ldz = N+1;
double sigma=0;
double z[(N + 1) * (nev + 1)];
a_int ldz = N + 1;
double sigma = 0;
int k;
for (k=0; k < 3*N; ++k )
workd[k] = 0;
double workl[3*(ncv*ncv) + 6*ncv];
for (k=0; k < 3*(ncv*ncv) + 6*ncv; ++k )
workl[k] = 0;
a_int lworkl = 3*(ncv*ncv) + 6*ncv;
for (k = 0; k < 3 * N; ++k) workd[k] = 0;
double workl[3 * (ncv * ncv) + 6 * ncv];
for (k = 0; k < 3 * (ncv * ncv) + 6 * ncv; ++k) workl[k] = 0;
a_int lworkl = 3 * (ncv * ncv) + 6 * ncv;
a_int info = 0;
int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
iparam[0] = 1;
iparam[2] = 10*N;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[4] = 0; // number of ev found by arpack.
iparam[4] = 0; // number of ev found by arpack.
iparam[6] = 1;
MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
while(ido != 99) {
while (ido != 99) {
/* call arpack like you would have, but, use dsaupd_c instead of dsaupd_ */
pdsaupd_c(MCW, &ido, bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr,
workd, workl, lworkl, &info);
pdsaupd_c(MCW, &ido, bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam,
ipntr, workd, workl, lworkl, &info);
dMatVec(&(workd[ipntr[0]-1]), &(workd[ipntr[1]-1]));
dMatVec(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
}
if (iparam[4] != nev) return 1; // check number of ev found by arpack.
if (iparam[4] != nev) return 1; // check number of ev found by arpack.
/* call arpack like you would have, but, use dseupd_c instead of dseupd_ */
pdseupd_c(MCW, rvec, howmny, select, d, z, ldz, sigma,
bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr,
workd, workl, lworkl, &info);
pdseupd_c(MCW, rvec, howmny, select, d, z, ldz, sigma, bmat, N, which, nev,
tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl,
&info);
int i;
for (i = 0; i < nev; ++i) {
printf("rank %d - %f\n", rank, d[i]);
/*eigen value order: smallest -> biggest*/
if(fabs(d[i] - (double)(1000-(nev-1)+i))>1e-6){
if (fabs(d[i] - (double)(1000 - (nev - 1) + i)) > 1e-6) {
free(d);
return 1;
}
@@ -92,10 +92,9 @@ int ds() {
return 0;
}
void zMatVec(double _Complex * x, double _Complex * y) {
void zMatVec(double _Complex* x, double _Complex* y) {
int i;
for (i = 0; i < 1000; ++i)
y[i] = x[i] * (i+1.0 + _Complex_I * (i+1.0));
for (i = 0; i < 1000; ++i) y[i] = x[i] * (i + 1.0 + _Complex_I * (i + 1.0));
};
int zn() {
@@ -106,57 +105,58 @@ int zn() {
a_int nev = 1;
double tol = 0;
double _Complex resid[N];
a_int ncv = 2*nev+1;
double _Complex V[ncv*N];
a_int ncv = 2 * nev + 1;
double _Complex V[ncv * N];
a_int ldv = N;
a_int iparam[11];
a_int ipntr[14];
double _Complex workd[3*N];
double _Complex workd[3 * N];
a_int rvec = 0;
char howmny[] = "A";
double _Complex* d = (double _Complex*) malloc((nev+1)*sizeof(double _Complex));
double _Complex* d =
(double _Complex*)malloc((nev + 1) * sizeof(double _Complex));
a_int select[ncv];
for (int i = 0; i < ncv; i++) select[i] = 1;
double _Complex z[(N+1)*(nev+1)];
a_int ldz = N+1;
double _Complex sigma=0. + I*0.;
double _Complex z[(N + 1) * (nev + 1)];
a_int ldz = N + 1;
double _Complex sigma = 0. + I * 0.;
int k;
for (k=0; k < 3*N; ++k )
workd[k] = 0. + I * 0.;
double _Complex workl[3*(ncv*ncv) + 6*ncv];
for (k=0; k < 3*(ncv*ncv) + 6*ncv; ++k )
workl[k] = 0. + I * 0.;
a_int lworkl = 3*(ncv*ncv) + 6*ncv;
for (k = 0; k < 3 * N; ++k) workd[k] = 0. + I * 0.;
double _Complex workl[3 * (ncv * ncv) + 6 * ncv];
for (k = 0; k < 3 * (ncv * ncv) + 6 * ncv; ++k) workl[k] = 0. + I * 0.;
a_int lworkl = 3 * (ncv * ncv) + 6 * ncv;
double _Complex rwork[ncv];
double _Complex workev[2*ncv];
double _Complex workev[2 * ncv];
a_int info = 0;
int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
iparam[0] = 1;
iparam[2] = 10*N;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[4] = 0; // number of ev found by arpack.
iparam[4] = 0; // number of ev found by arpack.
iparam[6] = 1;
MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
while(ido != 99) {
while (ido != 99) {
/* call arpack like you would have, but, use znaupd_c instead of znaupd_ */
pznaupd_c(MCW, &ido, bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr,
workd, workl, lworkl, rwork, &info);
pznaupd_c(MCW, &ido, bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam,
ipntr, workd, workl, lworkl, rwork, &info);
zMatVec(&(workd[ipntr[0]-1]), &(workd[ipntr[1]-1]));
zMatVec(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
}
if (iparam[4] != nev) return 1; // check number of ev found by arpack.
if (iparam[4] != nev) return 1; // check number of ev found by arpack.
/* call arpack like you would have, but, use zneupd_c instead of zneupd_ */
pzneupd_c(MCW, rvec, howmny, select, d, z, ldz, sigma, workev,
bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr,
workd, workl, lworkl, rwork, &info);
pzneupd_c(MCW, rvec, howmny, select, d, z, ldz, sigma, workev, bmat, N, which,
nev, tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl,
rwork, &info);
int i;
for (i = 0; i < nev; ++i) {
printf("rank %d - %f %f\n", rank, creal(d[i]), cimag(d[i]));
/*eigen value order: smallest -> biggest*/
if(fabs(creal(d[i]) - (double)(1000-(nev-1)+i))>1e-6 || fabs(cimag(d[i]) - (double)(1000-(nev-1)+i))>1e-6){
if (fabs(creal(d[i]) - (double)(1000 - (nev - 1) + i)) > 1e-6 ||
fabs(cimag(d[i]) - (double)(1000 - (nev - 1) + i)) > 1e-6) {
free(d);
return 1;
}
@@ -169,7 +169,7 @@ int main() {
MPI_Init(NULL, NULL);
sstats_c();
int rc = ds(); // parpack without debug.
int rc = ds(); // parpack without debug.
fflush(stdout);
MPI_Barrier(MPI_COMM_WORLD);
if (rc != 0) return rc;
@@ -178,21 +178,25 @@ int main() {
float tnaupd_c, tnaup2_c, tnaitr_c, tneigt_c, tngets_c, tnapps_c, tnconv_c;
float tcaupd_c, tcaup2_c, tcaitr_c, tceigt_c, tcgets_c, tcapps_c, tcconv_c;
float tmvopx_c, tmvbx_c, tgetv0_c, titref_c, trvec_c;
stat_c( &nopx_c, &nbx_c, &nrorth_c, &nitref_c, &nrstrt_c,
&tsaupd_c, &tsaup2_c, &tsaitr_c, &tseigt_c, &tsgets_c, &tsapps_c, &tsconv_c,
&tnaupd_c, &tnaup2_c, &tnaitr_c, &tneigt_c, &tngets_c, &tnapps_c, &tnconv_c,
&tcaupd_c, &tcaup2_c, &tcaitr_c, &tceigt_c, &tcgets_c, &tcapps_c, &tcconv_c,
&tmvopx_c, &tmvbx_c, &tgetv0_c, &titref_c, &trvec_c);
printf("Timers : nopx %d, tmvopx %f - nbx %d, tmvbx %f\n", nopx_c, tmvopx_c, nbx_c, tmvbx_c);
stat_c(&nopx_c, &nbx_c, &nrorth_c, &nitref_c, &nrstrt_c, &tsaupd_c, &tsaup2_c,
&tsaitr_c, &tseigt_c, &tsgets_c, &tsapps_c, &tsconv_c, &tnaupd_c,
&tnaup2_c, &tnaitr_c, &tneigt_c, &tngets_c, &tnapps_c, &tnconv_c,
&tcaupd_c, &tcaup2_c, &tcaitr_c, &tceigt_c, &tcgets_c, &tcapps_c,
&tcconv_c, &tmvopx_c, &tmvbx_c, &tgetv0_c, &titref_c, &trvec_c);
printf("Timers : nopx %d, tmvopx %f - nbx %d, tmvbx %f\n", nopx_c, tmvopx_c,
nbx_c, tmvbx_c);
int rank = 0; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
int rank = 0;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
if (rank == 0) printf("------\n");
// clang-format off
debug_c(6, -6, 1,
1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1); // set debug flags.
rc = zn(); // parpack with debug.
// clang-format on
rc = zn(); // parpack with debug.
fflush(stdout);
MPI_Barrier(MPI_COMM_WORLD);
if (rc != 0) return rc;
PARPACK/TESTS/MPI/icb_parpack_cpp.cpp
+14 -9
View File
@@ -11,15 +11,14 @@
* with entries 1000, 999, ... , 2, 1 on the diagonal.
*/
#include "parpack.hpp"
#include <array>
#include <cmath>
#include <iostream>
#include <vector>
#include "debug_c.hpp" // debug parpack.
#include "stat_c.hpp" // arpack statistics.
#include "parpack.hpp"
#include "stat_c.hpp" // arpack statistics.
void diagonal_matrix_vector_product(float const* const x, float* const y) {
for (int i = 0; i < 1000; ++i) {
@@ -76,8 +75,10 @@ void real_symmetric_runner() {
&(workd[ipntr[1] - 1]));
}
// check number of ev found by arpack.
if (iparam[4] < nev /*arpack may succeed to compute more EV than expected*/ || info != 0) {
std::cout << "ERROR: iparam[4] " << iparam[4] << ", nev " << nev << ", info " << info << std::endl;
if (iparam[4] < nev /*arpack may succeed to compute more EV than expected*/ ||
info != 0) {
std::cout << "ERROR: iparam[4] " << iparam[4] << ", nev " << nev
<< ", info " << info << std::endl;
throw std::domain_error("Error inside ARPACK routines");
}
@@ -157,8 +158,10 @@ void complex_symmetric_runner() {
}
// check number of ev found by arpack
if (iparam[4] < nev /*arpack may succeed to compute more EV than expected*/ || info != 0) {
std::cout << "ERROR: iparam[4] " << iparam[4] << ", nev " << nev << ", info " << info << std::endl;
if (iparam[4] < nev /*arpack may succeed to compute more EV than expected*/ ||
info != 0) {
std::cout << "ERROR: iparam[4] " << iparam[4] << ", nev " << nev
<< ", info " << info << std::endl;
throw std::domain_error("Error inside ARPACK routines");
}
@@ -173,8 +176,10 @@ void complex_symmetric_runner() {
std::cout << "rank " << rank << " - " << std::real(d[i]) << " "
<< std::imag(d[i]) << '\n';
/*eigen value order: smallest -> biggest*/
if (std::abs(std::real(d[i]) - static_cast<float>(1000 - (nev - 1) + i)) > 1. ||
std::abs(std::imag(d[i]) - static_cast<float>(1000 - (nev - 1) + i)) > 1.) {
if (std::abs(std::real(d[i]) - static_cast<float>(1000 - (nev - 1) + i)) >
1. ||
std::abs(std::imag(d[i]) - static_cast<float>(1000 - (nev - 1) + i)) >
1.) {
throw std::domain_error("Correct eigenvalues not computed");
}
}
TESTS/bug_1315_double.c
+38 -41
View File
@@ -1,8 +1,8 @@
#include "arpackdef.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "arpackdef.h"
// This test calls fortran from C the old-fashion cumbersome way.
// Note: icb_arpack_c tests the same kind of things using ICB.
@@ -19,19 +19,19 @@
* symmetric but is done to exhibit the bug.
* */
extern void dnaupd(a_int *, char *, a_int *, char *, a_int *,
double *, double *, a_int *, double *,
a_int *, a_int *, a_int *, double *,
double *, a_int *, a_int *);
extern void dnaupd(a_int *, char *, a_int *, char *, a_int *, double *,
double *, a_int *, double *, a_int *, a_int *, a_int *,
double *, double *, a_int *, a_int *);
extern void dneupd( a_int*, char*, a_int *, double *, double *, double *, a_int*, double *,
double *, double *, char *, a_int *, char *, a_int *, double *, double *, a_int *,
double *, a_int *, a_int *, a_int *, double *, double *, a_int *, a_int * );
extern void dneupd(a_int *, char *, a_int *, double *, double *, double *,
a_int *, double *, double *, double *, char *, a_int *,
char *, a_int *, double *, double *, a_int *, double *,
a_int *, a_int *, a_int *, double *, double *, a_int *,
a_int *);
void matVec(double * x, double * y) {
void matVec(double *x, double *y) {
int i;
for ( i = 0; i < 1000; ++i)
y[i] = ((double) (i+1))*x[i];
for (i = 0; i < 1000; ++i) y[i] = ((double)(i + 1)) * x[i];
};
int main() {
@@ -42,54 +42,51 @@ int main() {
a_int nev = 9;
double tol = 0;
double resid[N];
a_int ncv = 2*nev+1;
double V[ncv*N];
a_int ncv = 2 * nev + 1;
double V[ncv * N];
a_int ldv = N;
a_int iparam[11];
a_int ipntr[14];
double workd[3*N];
double workd[3 * N];
a_int rvec = 1;
char howmny[] = "A";
double* dr = (double*) malloc((nev+1)*sizeof(double));
double* di = (double*) malloc((nev+1)*sizeof(double));
a_int select[3*ncv];
double z[(N+1)*(nev+1)];
a_int ldz = N+1;
double sigmar=0;
double sigmai=0;
double workev[3*ncv];
double *dr = (double *)malloc((nev + 1) * sizeof(double));
double *di = (double *)malloc((nev + 1) * sizeof(double));
a_int select[3 * ncv];
double z[(N + 1) * (nev + 1)];
a_int ldz = N + 1;
double sigmar = 0;
double sigmai = 0;
double workev[3 * ncv];
int k;
for (k=0; k < 3*N; ++k )
workd[k] = 0;
double workl[3*(ncv*ncv) + 6*ncv];
for (k=0; k < 3*(ncv*ncv) + 6*ncv; ++k )
workl[k] = 0;
a_int lworkl = 3*(ncv*ncv) + 6*ncv;
for (k = 0; k < 3 * N; ++k) workd[k] = 0;
double workl[3 * (ncv * ncv) + 6 * ncv];
for (k = 0; k < 3 * (ncv * ncv) + 6 * ncv; ++k) workl[k] = 0;
a_int lworkl = 3 * (ncv * ncv) + 6 * ncv;
a_int info = 0;
iparam[0] = 1;
iparam[2] = 10*N;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[6] = 1;
dnaupd(&ido, bmat, &N, which, &nev, &tol, resid, &ncv, V, &ldv, iparam, ipntr,
workd, workl, &lworkl, &info);
workd, workl, &lworkl, &info);
while(ido == -1 || ido == 1) {
while (ido == -1 || ido == 1) {
matVec(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
matVec(&(workd[ipntr[0]-1]), &(workd[ipntr[1]-1]));
dnaupd(&ido, bmat, &N, which, &nev, &tol, resid, &ncv, V, &ldv, iparam, ipntr,
workd, workl, &lworkl, &info);
dnaupd(&ido, bmat, &N, which, &nev, &tol, resid, &ncv, V, &ldv, iparam,
ipntr, workd, workl, &lworkl, &info);
}
dneupd( &rvec, howmny, select, dr,di, z, &ldz, &sigmar, &sigmai,workev,
bmat, &N, which, &nev, &tol, resid, &ncv, V, &ldv, iparam, ipntr,
workd, workl, &lworkl, &info);
dneupd(&rvec, howmny, select, dr, di, z, &ldz, &sigmar, &sigmai, workev, bmat,
&N, which, &nev, &tol, resid, &ncv, V, &ldv, iparam, ipntr, workd,
workl, &lworkl, &info);
int i;
for (i = 0; i < nev; ++i) {
printf("%f\n", dr[i]);
if(fabs(dr[i] - (double)(1000-i))>1e-6){
if (fabs(dr[i] - (double)(1000 - i)) > 1e-6) {
free(dr);
free(di);
exit(EXIT_FAILURE);
TESTS/bug_1315_single.c
+37 -42
View File
@@ -1,8 +1,8 @@
#include "arpackdef.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "arpackdef.h"
// This test calls fortran from C the old-fashion cumbersome way.
// Note: icb_arpack_c tests the same kind of things using ICB.
@@ -19,20 +19,18 @@
* symmetric but is done to exhibit the bug.
*/
extern void snaupd(a_int *, char *, a_int *, char *, a_int *,
float *, float *, a_int *, float *,
a_int *, a_int *, a_int *, float *,
float *, a_int *, a_int *);
extern void snaupd(a_int *, char *, a_int *, char *, a_int *, float *, float *,
a_int *, float *, a_int *, a_int *, a_int *, float *,
float *, a_int *, a_int *);
extern void sneupd( a_int*, char*, a_int *, float *, float *, float *, a_int*, float *,
float *, float *, char *, a_int *, char *, a_int *, float *, float *, a_int *,
float *, a_int *, a_int *, a_int *, float *, float *, a_int *, a_int * );
extern void sneupd(a_int *, char *, a_int *, float *, float *, float *, a_int *,
float *, float *, float *, char *, a_int *, char *, a_int *,
float *, float *, a_int *, float *, a_int *, a_int *,
a_int *, float *, float *, a_int *, a_int *);
void matVec(float * x, float * y) {
void matVec(float *x, float *y) {
int i;
for ( i = 0; i < 1000; ++i)
y[i] = ((float) (i+1))*x[i];
for (i = 0; i < 1000; ++i) y[i] = ((float)(i + 1)) * x[i];
};
int main() {
@@ -43,54 +41,51 @@ int main() {
a_int nev = 9;
float tol = 0;
float resid[N];
a_int ncv = 2*nev+1;
float V[ncv*N];
a_int ncv = 2 * nev + 1;
float V[ncv * N];
a_int ldv = N;
a_int iparam[11];
a_int ipntr[14];
float workd[3*N];
float workd[3 * N];
a_int rvec = 1;
char howmny[] = "A";
float* dr = (float*) malloc((nev+1)*sizeof(float));
float* di = (float*) malloc((nev+1)*sizeof(float));
a_int select[3*ncv];
float z[(N+1)*(nev+1)];
a_int ldz = N+1;
float sigmar=0;
float sigmai=0;
float workev[3*ncv];
float *dr = (float *)malloc((nev + 1) * sizeof(float));
float *di = (float *)malloc((nev + 1) * sizeof(float));
a_int select[3 * ncv];
float z[(N + 1) * (nev + 1)];
a_int ldz = N + 1;
float sigmar = 0;
float sigmai = 0;
float workev[3 * ncv];
int k;
for (k=0; k < 3*N; ++k )
workd[k] = 0;
float workl[3*(ncv*ncv) + 6*ncv];
for (k=0; k < 3*(ncv*ncv) + 6*ncv; ++k )
workl[k] = 0;
a_int lworkl = 3*(ncv*ncv) + 6*ncv;
for (k = 0; k < 3 * N; ++k) workd[k] = 0;
float workl[3 * (ncv * ncv) + 6 * ncv];
for (k = 0; k < 3 * (ncv * ncv) + 6 * ncv; ++k) workl[k] = 0;
a_int lworkl = 3 * (ncv * ncv) + 6 * ncv;
a_int info = 0;
iparam[0] = 1;
iparam[2] = 10*N;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[6] = 1;
snaupd(&ido, bmat, &N, which, &nev, &tol, resid, &ncv, V, &ldv, iparam, ipntr,
workd, workl, &lworkl, &info);
workd, workl, &lworkl, &info);
while(ido == -1 || ido == 1) {
while (ido == -1 || ido == 1) {
matVec(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
matVec(&(workd[ipntr[0]-1]), &(workd[ipntr[1]-1]));
snaupd(&ido, bmat, &N, which, &nev, &tol, resid, &ncv, V, &ldv, iparam, ipntr,
workd, workl, &lworkl, &info);
snaupd(&ido, bmat, &N, which, &nev, &tol, resid, &ncv, V, &ldv, iparam,
ipntr, workd, workl, &lworkl, &info);
}
sneupd( &rvec, howmny, select, dr,di, z, &ldz, &sigmar, &sigmai,workev,
bmat, &N, which, &nev, &tol, resid, &ncv, V, &ldv, iparam, ipntr,
workd, workl, &lworkl, &info);
sneupd(&rvec, howmny, select, dr, di, z, &ldz, &sigmar, &sigmai, workev, bmat,
&N, which, &nev, &tol, resid, &ncv, V, &ldv, iparam, ipntr, workd,
workl, &lworkl, &info);
int i;
for (i = 0; i < nev; ++i) {
printf("%f\n", dr[i]);
if(fabs(dr[i] - (float)(1000-i))>1e-2){
if (fabs(dr[i] - (float)(1000 - i)) > 1e-2) {
free(dr);
free(di);
exit(EXIT_FAILURE);
TESTS/icb_arpack_c.c
+67 -67
View File
@@ -1,28 +1,29 @@
/*
* This example demonstrates the use of ISO_C_BINDING to call arpack (portability).
* This example demonstrates the use of ISO_C_BINDING to call arpack
* (portability).
*
* Just use arpack as you would have normally done, but, use *[ae]upd_c instead of *[ae]upd_.
* The main advantage is that compiler checks (arguments) are performed at build time.
* Note: to debug arpack, call debug_c.
* Just use arpack as you would have normally done, but, use *[ae]upd_c instead
* of *[ae]upd_. The main advantage is that compiler checks (arguments) are
* performed at build time. Note: to debug arpack, call debug_c.
*/
#include <complex.h> // creal, cimag.
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "arpack.h"
#include <complex.h> // creal, cimag.
#include "debug_c.h" // debug arpack.
#include "stat_c.h" // arpack statistics.
#include "debug_c.h" // debug arpack.
#include "stat_c.h" // arpack statistics.
/* test program to solve for the 9 largest eigenvalues of
* A*x = lambda*x where A is the diagonal matrix
* with entries 1000, 999, ... , 2, 1 on the diagonal.
* */
void dMatVec(double * x, double * y) {
void dMatVec(double* x, double* y) {
int i;
for ( i = 0; i < 1000; ++i)
y[i] = ((double) (i+1))*x[i];
for (i = 0; i < 1000; ++i) y[i] = ((double)(i + 1)) * x[i];
};
int ds() {
@@ -33,53 +34,50 @@ int ds() {
a_int nev = 9;
double tol = 0;
double resid[N];
a_int ncv = 2*nev+1;
double V[ncv*N];
a_int ncv = 2 * nev + 1;
double V[ncv * N];
a_int ldv = N;
a_int iparam[11];
a_int ipntr[14];
double workd[3*N];
double workd[3 * N];
a_int rvec = 1;
char howmny[] = "A";
double* d = (double*) malloc((nev+1)*sizeof(double));
double* d = (double*)malloc((nev + 1) * sizeof(double));
a_int select[ncv];
for (int i = 0; i < ncv; i++) select[i] = 1;
double z[(N+1)*(nev+1)];
a_int ldz = N+1;
double sigma=0;
double z[(N + 1) * (nev + 1)];
a_int ldz = N + 1;
double sigma = 0;
int k;
for (k=0; k < 3*N; ++k )
workd[k] = 0;
double workl[3*(ncv*ncv) + 6*ncv];
for (k=0; k < 3*(ncv*ncv) + 6*ncv; ++k )
workl[k] = 0;
a_int lworkl = 3*(ncv*ncv) + 6*ncv;
for (k = 0; k < 3 * N; ++k) workd[k] = 0;
double workl[3 * (ncv * ncv) + 6 * ncv];
for (k = 0; k < 3 * (ncv * ncv) + 6 * ncv; ++k) workl[k] = 0;
a_int lworkl = 3 * (ncv * ncv) + 6 * ncv;
a_int info = 0;
iparam[0] = 1;
iparam[2] = 10*N;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[4] = 0; // number of ev found by arpack.
iparam[4] = 0; // number of ev found by arpack.
iparam[6] = 1;
while(ido != 99) {
while (ido != 99) {
/* call arpack like you would have, but, use dsaupd_c instead of dsaupd_ */
dsaupd_c(&ido, bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr,
workd, workl, lworkl, &info);
dMatVec(&(workd[ipntr[0]-1]), &(workd[ipntr[1]-1]));
dMatVec(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
}
if (iparam[4] != nev) return 1; // check number of ev found by arpack.
if (iparam[4] != nev) return 1; // check number of ev found by arpack.
/* call arpack like you would have, but, use dseupd_c instead of dseupd_ */
dseupd_c(rvec, howmny, select, d, z, ldz, sigma,
bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr,
workd, workl, lworkl, &info);
dseupd_c(rvec, howmny, select, d, z, ldz, sigma, bmat, N, which, nev, tol,
resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl, &info);
int i;
for (i = 0; i < nev; ++i) {
printf("%f\n", d[i]);
/*eigen value order: smallest -> biggest*/
if(fabs(d[i] - (double)(1000-(nev-1)+i))>1e-6){
if (fabs(d[i] - (double)(1000 - (nev - 1) + i)) > 1e-6) {
free(d);
return 1;
}
@@ -88,10 +86,9 @@ int ds() {
return 0;
}
void zMatVec(double _Complex * x, double _Complex * y) {
void zMatVec(double _Complex* x, double _Complex* y) {
int i;
for (i = 0; i < 1000; ++i)
y[i] = x[i] * (i+1.0 + _Complex_I * (i+1.0));
for (i = 0; i < 1000; ++i) y[i] = x[i] * (i + 1.0 + _Complex_I * (i + 1.0));
};
int zn() {
@@ -102,55 +99,55 @@ int zn() {
a_int nev = 9;
double tol = 0;
double _Complex resid[N];
a_int ncv = 2*nev+1;
double _Complex V[ncv*N];
a_int ncv = 2 * nev + 1;
double _Complex V[ncv * N];
a_int ldv = N;
a_int iparam[11];
a_int ipntr[14];
double _Complex workd[3*N];
double _Complex workd[3 * N];
a_int rvec = 0;
char howmny[] = "A";
double _Complex* d = (double _Complex*) malloc((nev+1)*sizeof(double _Complex));
double _Complex* d =
(double _Complex*)malloc((nev + 1) * sizeof(double _Complex));
a_int select[ncv];
for (int i = 0; i < ncv; i++) select[i] = 1;
double _Complex z[(N+1)*(nev+1)];
a_int ldz = N+1;
double _Complex sigma=0. + I*0.;
double _Complex z[(N + 1) * (nev + 1)];
a_int ldz = N + 1;
double _Complex sigma = 0. + I * 0.;
int k;
for (k=0; k < 3*N; ++k )
workd[k] = 0;
double _Complex workl[3*(ncv*ncv) + 6*ncv];
for (k=0; k < 3*(ncv*ncv) + 6*ncv; ++k )
workl[k] = 0;
a_int lworkl = 3*(ncv*ncv) + 6*ncv;
for (k = 0; k < 3 * N; ++k) workd[k] = 0;
double _Complex workl[3 * (ncv * ncv) + 6 * ncv];
for (k = 0; k < 3 * (ncv * ncv) + 6 * ncv; ++k) workl[k] = 0;
a_int lworkl = 3 * (ncv * ncv) + 6 * ncv;
double rwork[ncv];
double _Complex workev[2*ncv];
double _Complex workev[2 * ncv];
a_int info = 0;
iparam[0] = 1;
iparam[2] = 10*N;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[4] = 0; // number of ev found by arpack.
iparam[4] = 0; // number of ev found by arpack.
iparam[6] = 1;
while(ido != 99) {
while (ido != 99) {
/* call arpack like you would have, but, use znaupd_c instead of znaupd_ */
znaupd_c(&ido, bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr,
workd, workl, lworkl, rwork, &info);
zMatVec(&(workd[ipntr[0]-1]), &(workd[ipntr[1]-1]));
zMatVec(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
}
if (iparam[4] != nev) return 1; // check number of ev found by arpack.
if (iparam[4] != nev) return 1; // check number of ev found by arpack.
/* call arpack like you would have, but, use zneupd_c instead of zneupd_ */
zneupd_c(rvec, howmny, select, d, z, ldz, sigma, workev,
bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr,
workd, workl, lworkl, rwork, &info);
zneupd_c(rvec, howmny, select, d, z, ldz, sigma, workev, bmat, N, which, nev,
tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl, rwork,
&info);
int i;
for (i = 0; i < nev; ++i) {
printf("%f %f\n", creal(d[i]), cimag(d[i]));
/*eigen value order: smallest -> biggest*/
if(fabs(creal(d[i]) - (double)(1000-(nev-1)+i))>1e-6 || fabs(cimag(d[i]) - (double)(1000-(nev-1)+i))>1e-6){
if (fabs(creal(d[i]) - (double)(1000 - (nev - 1) + i)) > 1e-6 ||
fabs(cimag(d[i]) - (double)(1000 - (nev - 1) + i)) > 1e-6) {
free(d);
return 1;
}
@@ -161,27 +158,30 @@ int zn() {
int main() {
sstats_c();
int rc = ds(); // arpack without debug.
int rc = ds(); // arpack without debug.
if (rc != 0) return rc;
a_int nopx_c, nbx_c, nrorth_c, nitref_c, nrstrt_c;
float tsaupd_c, tsaup2_c, tsaitr_c, tseigt_c, tsgets_c, tsapps_c, tsconv_c;
float tnaupd_c, tnaup2_c, tnaitr_c, tneigt_c, tngets_c, tnapps_c, tnconv_c;
float tcaupd_c, tcaup2_c, tcaitr_c, tceigt_c, tcgets_c, tcapps_c, tcconv_c;
float tmvopx_c, tmvbx_c, tgetv0_c, titref_c, trvec_c;
stat_c( &nopx_c, &nbx_c, &nrorth_c, &nitref_c, &nrstrt_c,
&tsaupd_c, &tsaup2_c, &tsaitr_c, &tseigt_c, &tsgets_c, &tsapps_c, &tsconv_c,
&tnaupd_c, &tnaup2_c, &tnaitr_c, &tneigt_c, &tngets_c, &tnapps_c, &tnconv_c,
&tcaupd_c, &tcaup2_c, &tcaitr_c, &tceigt_c, &tcgets_c, &tcapps_c, &tcconv_c,
&tmvopx_c, &tmvbx_c, &tgetv0_c, &titref_c, &trvec_c);
printf("Timers : nopx %d, tmvopx %f - nbx %d, tmvbx %f\n", nopx_c, tmvopx_c, nbx_c, tmvbx_c);
stat_c(&nopx_c, &nbx_c, &nrorth_c, &nitref_c, &nrstrt_c, &tsaupd_c, &tsaup2_c,
&tsaitr_c, &tseigt_c, &tsgets_c, &tsapps_c, &tsconv_c, &tnaupd_c,
&tnaup2_c, &tnaitr_c, &tneigt_c, &tngets_c, &tnapps_c, &tnconv_c,
&tcaupd_c, &tcaup2_c, &tcaitr_c, &tceigt_c, &tcgets_c, &tcapps_c,
&tcconv_c, &tmvopx_c, &tmvbx_c, &tgetv0_c, &titref_c, &trvec_c);
printf("Timers : nopx %d, tmvopx %f - nbx %d, tmvbx %f\n", nopx_c, tmvopx_c,
nbx_c, tmvbx_c);
printf("------\n");
// clang-format off
debug_c(6, -6, 1,
1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1); // set debug flags.
rc = zn(); // arpack with debug.
// clang-format on
rc = zn(); // arpack with debug.
return rc;
}
TESTS/icb_arpack_cpp.cpp
+19 -14
View File
@@ -9,25 +9,24 @@
* matrix with entries 1000, 999, ... , 2, 1 on the diagonal.
*/
#include "arpack.hpp"
#include <array>
#include <cmath>
#include <iostream>
#include <vector>
#include "arpack.hpp"
#include "debug_c.hpp" // debug arpack.
#include "stat_c.hpp" // arpack statistics.
template<typename Real>
template <typename Real>
void diagonal_matrix_vector_product(Real const* const x, Real* const y) {
for (int i = 0; i < 1000; ++i) {
y[i] = static_cast<Real>(i + 1) * x[i];
}
}
template<typename Real>
void real_symmetric_runner(double const & tol_check) {
template <typename Real>
void real_symmetric_runner(double const& tol_check) {
a_int const N = 1000;
a_int const nev = 9;
@@ -73,8 +72,10 @@ void real_symmetric_runner(double const & tol_check) {
}
// check number of ev found by arpack.
if (iparam[4] < nev /*arpack may succeed to compute more EV than expected*/ || info != 0) {
std::cout << "ERROR: iparam[4] " << iparam[4] << ", nev " << nev << ", info " << info << std::endl;
if (iparam[4] < nev /*arpack may succeed to compute more EV than expected*/ ||
info != 0) {
std::cout << "ERROR: iparam[4] " << iparam[4] << ", nev " << nev
<< ", info " << info << std::endl;
throw std::domain_error("Error inside ARPACK routines");
}
@@ -98,7 +99,7 @@ void real_symmetric_runner(double const & tol_check) {
std::cout << "------\n";
}
template<typename Real>
template <typename Real>
void diagonal_matrix_vector_product(std::complex<Real> const* const x,
std::complex<Real>* const y) {
for (int i = 0; i < 1000; ++i) {
@@ -106,8 +107,8 @@ void diagonal_matrix_vector_product(std::complex<Real> const* const x,
}
}
template<typename Real>
void complex_symmetric_runner(double const & tol_check) {
template <typename Real>
void complex_symmetric_runner(double const& tol_check) {
a_int const N = 1000;
a_int const nev = 9;
@@ -154,8 +155,10 @@ void complex_symmetric_runner(double const & tol_check) {
}
// check number of ev found by arpack.
if (iparam[4] < nev /*arpack may succeed to compute more EV than expected*/ || info != 0) {
std::cout << "ERROR: iparam[4] " << iparam[4] << ", nev " << nev << ", info " << info << std::endl;
if (iparam[4] < nev /*arpack may succeed to compute more EV than expected*/ ||
info != 0) {
std::cout << "ERROR: iparam[4] " << iparam[4] << ", nev " << nev
<< ", info " << info << std::endl;
throw std::domain_error("Error inside ARPACK routines");
}
@@ -172,8 +175,10 @@ void complex_symmetric_runner(double const & tol_check) {
std::cout << d[i] << "\n";
/*eigen value order: smallest -> biggest*/
if (std::abs(std::real(d[i]) - static_cast<Real>(1000 - (nev - 1) + i)) > tol_check ||
std::abs(std::imag(d[i]) - static_cast<Real>(1000 - (nev - 1) + i)) > tol_check) {
if (std::abs(std::real(d[i]) - static_cast<Real>(1000 - (nev - 1) + i)) >
tol_check ||
std::abs(std::imag(d[i]) - static_cast<Real>(1000 - (nev - 1) + i)) >
tol_check) {
throw std::domain_error("Correct eigenvalues not computed");
}
}