diff --git a/EXAMPLES/MATRIX_MARKET/arpackmm.cpp b/EXAMPLES/MATRIX_MARKET/arpackmm.cpp index 6215e7f..8865bd8 100644 --- a/EXAMPLES/MATRIX_MARKET/arpackmm.cpp +++ b/EXAMPLES/MATRIX_MARKET/arpackmm.cpp @@ -1,339 +1,585 @@ +#include + #include "arpackSolver.hpp" #include "debug_c.hpp" #include "stat_c.hpp" -#include using namespace std; class options { - public: - options() { - fileA = "A.mtx"; - fileB = "N.A."; // Not available. - dense = false; - denseRR = true; - nbEV = 1; - nbCV = 2*nbEV + 1; - stdPb = true; // Standard or generalized (= not standard). - symPb = true; - cpxPb = false; - simplePrec = false; // Double precision. - mag = string("LM"); // Large magnitude. - shiftReal = false; shiftImag = false; - sigmaReal = 0.; sigmaImag = 0.; // Eigen value translation: look for lambda+sigma instead of lambda. - invert = false; // Eigen value invertion: look for 1./lambda instead of lambda. - tol = 1.e-06; - maxIt = 100; - schur = false; // Compute Ritz vectors. - slv = "BiCG"; - slvItrTol = 1.e-6; - slvItrMaxIt = 100; - slvItrPC = "Diag"; - slvDrtPivot = 1.e-6; - slvDrtOffset = 0.; - slvDrtScale = 1.; - check = true; - verbose = 0; - debug = 0; - restart = false; - }; + public: + options() { + fileA = "A.mtx"; + fileB = "N.A."; // Not available. + dense = false; + denseRR = true; + nbEV = 1; + nbCV = 2 * nbEV + 1; + stdPb = true; // Standard or generalized (= not standard). + symPb = true; + cpxPb = false; + simplePrec = false; // Double precision. + mag = string("LM"); // Large magnitude. + shiftReal = false; + shiftImag = false; + sigmaReal = 0.; + sigmaImag = 0.; // Eigen value translation: look for lambda+sigma instead + // of lambda. + invert = + false; // Eigen value invertion: look for 1./lambda instead of lambda. + tol = 1.e-06; + maxIt = 100; + schur = false; // Compute Ritz vectors. + slv = "BiCG"; + slvItrTol = 1.e-6; + slvItrMaxIt = 100; + slvItrPC = "Diag"; + slvDrtPivot = 1.e-6; + slvDrtOffset = 0.; + slvDrtScale = 1.; + check = true; + verbose = 0; + debug = 0; + restart = false; + }; - int readCmdLine(int argc, char ** argv) { - // Check for command line independent parameters. + int readCmdLine(int argc, char** argv) { + // Check for command line independent parameters. - for (int a = 1; argv && a < argc; a++) { - string clo = argv[a]; // Command line option. - if (clo == "--help") return usage(0); - if (clo == "--A") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - fileA = argv[a]; + for (int a = 1; argv && a < argc; a++) { + string clo = argv[a]; // Command line option. + if (clo == "--help") return usage(0); + if (clo == "--A") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); } - if (clo == "--dense") { - dense = true; - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - string rr(argv[a]); - if (rr != "true" && rr != "false") {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - denseRR = (rr == "true") ? true : false; - } - if (clo == "--nbEV") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream nEV(argv[a]); - nEV >> nbEV; if (!nEV) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - nbCV = 2*nbEV + 1; - } - if (clo == "--genPb") { - stdPb = false; - fileB = "B.mtx"; - } - if (clo == "--nonSymPb") symPb = false; - if (clo == "--cpxPb") { - symPb = false; - cpxPb = true; - } - if (clo == "--simplePrec") simplePrec = true; - if (clo == "--mag") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - mag = argv[a]; // small mag (likely poor perf) <=> large mag + invert (likely good perf). - bool ok = (mag == "LM" || mag == "SM" || mag == "LR" || mag == "SR" || mag == "LI" || mag == "SI") ? true : false; - if (!ok) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - } - if (clo == "--shiftReal") { - shiftReal = true; - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream s(argv[a]); - s >> sigmaReal; if (!s) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - } - if (clo == "--shiftImag") { - shiftImag = true; - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream s(argv[a]); - s >> sigmaImag; if (!s) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - } - if (clo == "--invert") invert = true; - if (clo == "--tol") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream t(argv[a]); - t >> tol; if (!t) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - } - if (clo == "--maxIt") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream mi(argv[a]); - mi >> maxIt; if (!mi) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - } - if (clo == "--slv") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - slv = argv[a]; - } - if (clo == "--slvItrTol") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream t(argv[a]); - double tol = 0.; - t >> slvItrTol; if (!t) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - } - if (clo == "--slvItrMaxIt") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream mi(argv[a]); - int maxIt = 0; - mi >> slvItrMaxIt; if (!mi) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - } - if (clo == "--slvItrPC") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream pc(argv[a]); - pc >> slvItrPC; if (!pc) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - } - if (clo == "--slvDrtPivot") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream pv(argv[a]); - pv >> slvDrtPivot; if (!pv) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - } - if (clo == "--slvDrtOffset") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream of(argv[a]); - of >> slvDrtOffset; if (!of) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - } - if (clo == "--slvDrtScale") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream sc(argv[a]); - sc >> slvDrtScale; if (!sc) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - } - if (clo == "--noCheck") check = false; - if (clo == "--verbose") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream vb(argv[a]); - vb >> verbose; if (!vb) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - } - if (clo == "--debug") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream dbg(argv[a]); - dbg >> debug; if (!dbg) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} - if (debug > 3) debug = 3; - debug_c(6, -6, debug, debug, debug, debug, debug, debug, debug, debug, debug, debug, debug, - debug, debug, debug, debug, debug, debug, debug, debug, debug, debug, debug); - } - if (clo == "--restart") restart = true; + fileA = argv[a]; } - - // Check for command line dependent parameters. - - for (int a = 1; argv && a < argc; a++) { - string clo = argv[a]; // Command line option. - if (clo == "--nbCV") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - stringstream nCV(argv[a]); - nCV >> nbCV; if (!nCV) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();} + if (clo == "--dense") { + dense = true; + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); } - if (clo == "--B") { - a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();} - fileB = argv[a]; + string rr(argv[a]); + if (rr != "true" && rr != "false") { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + denseRR = (rr == "true") ? true : false; + } + if (clo == "--nbEV") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream nEV(argv[a]); + nEV >> nbEV; + if (!nEV) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + nbCV = 2 * nbEV + 1; + } + if (clo == "--genPb") { + stdPb = false; + fileB = "B.mtx"; + } + if (clo == "--nonSymPb") symPb = false; + if (clo == "--cpxPb") { + symPb = false; + cpxPb = true; + } + if (clo == "--simplePrec") simplePrec = true; + if (clo == "--mag") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + mag = argv[a]; // small mag (likely poor perf) <=> large mag + invert + // (likely good perf). + bool ok = (mag == "LM" || mag == "SM" || mag == "LR" || mag == "SR" || + mag == "LI" || mag == "SI") + ? true + : false; + if (!ok) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); } } + if (clo == "--shiftReal") { + shiftReal = true; + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream s(argv[a]); + s >> sigmaReal; + if (!s) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + } + if (clo == "--shiftImag") { + shiftImag = true; + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream s(argv[a]); + s >> sigmaImag; + if (!s) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + } + if (clo == "--invert") invert = true; + if (clo == "--tol") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream t(argv[a]); + t >> tol; + if (!t) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + } + if (clo == "--maxIt") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream mi(argv[a]); + mi >> maxIt; + if (!mi) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + } + if (clo == "--slv") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + slv = argv[a]; + } + if (clo == "--slvItrTol") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream t(argv[a]); + double tol = 0.; + t >> slvItrTol; + if (!t) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + } + if (clo == "--slvItrMaxIt") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream mi(argv[a]); + int maxIt = 0; + mi >> slvItrMaxIt; + if (!mi) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + } + if (clo == "--slvItrPC") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream pc(argv[a]); + pc >> slvItrPC; + if (!pc) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + } + if (clo == "--slvDrtPivot") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream pv(argv[a]); + pv >> slvDrtPivot; + if (!pv) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + } + if (clo == "--slvDrtOffset") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream of(argv[a]); + of >> slvDrtOffset; + if (!of) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + } + if (clo == "--slvDrtScale") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream sc(argv[a]); + sc >> slvDrtScale; + if (!sc) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + } + if (clo == "--noCheck") check = false; + if (clo == "--verbose") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream vb(argv[a]); + vb >> verbose; + if (!vb) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + } + if (clo == "--debug") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream dbg(argv[a]); + dbg >> debug; + if (!dbg) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + if (debug > 3) debug = 3; + debug_c(6, -6, debug, debug, debug, debug, debug, debug, debug, debug, + debug, debug, debug, debug, debug, debug, debug, debug, debug, + debug, debug, debug, debug, debug); + } + if (clo == "--restart") restart = true; + } - return 0; - }; + // Check for command line dependent parameters. - int usage(int rc = 1) { - cout << "Usage: running arpack with matrix market files to check for eigen values/vectors." << endl; - cout << endl; - cout << " --A F: file name of matrix A such that A X = lambda X. (standard)" << endl; - cout << " the file F must be compliant with the matrix market format." << endl; - cout << " default: A.mtx" << endl; - cout << " --B F: file name of matrix B such that A X = lambda B X. (generalized)" << endl; - cout << " the file F must be compliant with the matrix market format." << endl; - cout << " default: N.A. for standard problem, or, B.mtx for generalized problem" << endl; - cout << " --dense RR: consider A and B as dense matrices." << endl; - cout << " if RR = true, use more-stable-but-slow versions of LU / QR (rank revealing)." << endl; - cout << " if RR = false, use less-stable-but-fast versions of LU / QR (depends on condition number)." << endl; - cout << " Notes:" << endl; - cout << " - only direct solvers are available when using dense matrices." << endl; - cout << " - RR does not impact the use of LLT and LDLT." << endl; - cout << " - thresholds only make sense for rank-revealing decompositions." << endl; - cout << " default: consider A and B as sparse matrices" << endl; - cout << " --nbEV: number of eigen values/vectors to compute." << endl; - cout << " default: 1" << endl; - cout << " --nbCV: number of columns of the matrix V." << endl; - cout << " default: 2*nbEV+1" << endl; - cout << " --genPb: generalized problem." << endl; - cout << " default: standard problem" << endl; - cout << " --nonSymPb: non symmetric problem (<=> use dn[ae]upd)." << endl; - cout << " default: symmetric problem (<=> use ds[ae]upd)" << endl; - cout << " --cpxPb: complex (non symmetric) problem (<=> use zn[ae]upd)." << endl; - cout << " default: false (<=> use d*[ae]upd)" << endl; - cout << " --simplePrec: use simple precision (less accurate, but, half memory footprint)." << endl; - cout << " default: false (<=> use double precision: use [dz]*upd)" << endl; - cout << " --mag M: set magnitude of eigen values to look for (LM, SM, LR, SR, LI, SI)." << endl; - cout << " default: large magnitude (LM)" << endl; - cout << " --shiftReal S: real shift where sigma = S (look for lambda+S instead of lambda)." << endl; - cout << " default: no shift, S = 0." << endl; - cout << " --shiftImag S: imaginary shift where sigma = S (look for lambda+S instead of lambda)." << endl; - cout << " default: no shift, S = 0." << endl; - cout << " --invert: invert mode (look for 1./lambda instead of lambda)." << endl; - cout << " default: no invert" << endl; - cout << " --tol T: tolerance T." << endl; - cout << " default: 1.e-06" << endl; - cout << " --maxIt M: maximum iterations M." << endl; - cout << " default: 100" << endl; - cout << " --schur: compute Schur vectors." << endl; - cout << " the Schur decomposition is such that A = Q^H x T x Q where:" << endl; - cout << " - the H superscript refers to the Hermitian transpose: Q^H = (Q^t)^*." << endl; - cout << " - Q is unitary: Q is such that Q^H x Q = I." << endl; - cout << " - T is an upper-triangular matrix whose diagonal elements are the eigenvalues of A." << endl; - cout << " every square matrix has a Schur decomposition: columns of Q are the Schur vectors." << endl; - cout << " for a general matrix A, there is no relation between Schur vectors of A and eigenvectors of A." << endl; - cout << " if q_j is the j-th Schur vector, then A x q_j is a linear combination of q_1, ..., q_j." << endl; - cout << " Schur vectors q_1, q_2, ..., q_j span an invariant subspace of A." << endl; - cout << " the Schur vectors and eigenvectors of A are the same if A is a normal matrix." << endl; - cout << " default: compute Ritz vectors (approximations of eigen vectors)" << endl; - cout << " --slv S: solver (needed if arpack mode > 1)." << endl; - cout << " BiCG: iterative method, any matrices" << endl; - cout << " CG: iterative method, sym matrices only" << endl; - cout << " LU: direct method, any matrices (pivoting needed)" << endl; - cout << " QR: direct method, any matrices (pivoting needed)" << endl; - cout << " LLT: direct method, SPD matrices only (pivoting not needed)" << endl; - cout << " LDLT: direct method, symmetric positive semi-definite matrices only (pivoting not needed)" << endl; - cout << " default: BiCG" << endl; - cout << " --slvItrTol T: solver tolerance T (for iterative solvers)." << endl; - cout << " default: 1.e-6" << endl; - cout << " --slvItrMaxIt M: solver maximum iterations M (for iterative solvers)." << endl; - cout << " default: 100" << endl; - cout << " --slvItrPC PC: solver preconditioner (for iterative solvers)." << endl; - cout << " PC preconditioner:" << endl; - cout << " Diag: eigen diagonal preconditioner (Jacobi)." << endl; - cout << " ILU#D#F: eigen ILU preconditioner." << endl; - cout << " D: drop tolerance." << endl; - cout << " F: fill factor." << endl; - cout << " default: diagonal preconditioner (Jacobi)" << endl; - cout << " --slvDrtPivot: P solver pivot P (for direct solvers)." << endl; - cout << " default: 1.e-06" << endl; - cout << " --slvDrtOffset: O solver offset O (for direct solvers)." << endl; - cout << " default: 0." << endl; - cout << " --slvDrtScale: S solver scale S (for direct solvers)." << endl; - cout << " default: 1." << endl; - cout << " --noCheck: check arpack eigen values/vectors." << endl; - cout << " check will fail if Schur vectors are computed and A is NOT a normal matrix." << endl; - cout << " default: check" << endl; - cout << " --verbose V: verbosity level (up to 3)." << endl; - cout << " default: 0" << endl; - cout << " --debug D: debug level (up to 3)." << endl; - cout << " default: 0" << endl; - cout << " --restart: restart from previous run (which had produced arpackSolver.*.out)." << endl; - cout << " restart from eigen basis approximation computed during a previous run." << endl; - cout << " default: false" << endl; - if (rc == 0) exit(0); - return rc; - }; + for (int a = 1; argv && a < argc; a++) { + string clo = argv[a]; // Command line option. + if (clo == "--nbCV") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + stringstream nCV(argv[a]); + nCV >> nbCV; + if (!nCV) { + cerr << "Error: bad " << clo << " - bad argument" << endl; + return usage(); + } + } + if (clo == "--B") { + a++; + if (a >= argc) { + cerr << "Error: bad " << clo << " - need argument" << endl; + return usage(); + } + fileB = argv[a]; + } + } - friend ostream & operator<< (ostream & ostr, options const & opt); + return 0; + }; - string fileA; - string fileB; - bool dense; - bool denseRR; - a_int nbEV; - a_int nbCV; - bool stdPb; // Standard or generalized (= not standard). - bool symPb; - bool cpxPb; - bool simplePrec; - string mag; // Magnitude <=> "which" arpack parameter. - bool shiftReal, shiftImag; - double sigmaReal, sigmaImag; // Eigen value translation: look for lambda+sigma instead of lambda. - bool invert; // Eigen value invertion: look for 1./lambda instead of lambda. - double tol; - int maxIt; - bool schur; - string slv; - double slvItrTol; - int slvItrMaxIt; - string slvItrPC; - double slvDrtPivot; - double slvDrtOffset; - double slvDrtScale; - bool check; - int verbose; - a_int debug; - bool restart; + int usage(int rc = 1) { + cout << "Usage: running arpack with matrix market files to check for eigen " + "values/vectors." + << endl; + cout << endl; + cout << " --A F: file name of matrix A such that A X = lambda " + "X. (standard)" + << endl; + cout << " the file F must be compliant with the matrix " + "market format." + << endl; + cout << " default: A.mtx" << endl; + cout << " --B F: file name of matrix B such that A X = lambda " + "B X. (generalized)" + << endl; + cout << " the file F must be compliant with the matrix " + "market format." + << endl; + cout << " default: N.A. for standard problem, or, B.mtx " + "for generalized problem" + << endl; + cout << " --dense RR: consider A and B as dense matrices." << endl; + cout << " if RR = true, use more-stable-but-slow " + "versions of LU / QR (rank revealing)." + << endl; + cout << " if RR = false, use less-stable-but-fast " + "versions of LU / QR (depends on condition number)." + << endl; + cout << " Notes:" << endl; + cout << " - only direct solvers are available when " + "using dense matrices." + << endl; + cout + << " - RR does not impact the use of LLT and LDLT." + << endl; + cout << " - thresholds only make sense for " + "rank-revealing decompositions." + << endl; + cout << " default: consider A and B as sparse matrices" + << endl; + cout << " --nbEV: number of eigen values/vectors to compute." + << endl; + cout << " default: 1" << endl; + cout << " --nbCV: number of columns of the matrix V." << endl; + cout << " default: 2*nbEV+1" << endl; + cout << " --genPb: generalized problem." << endl; + cout << " default: standard problem" << endl; + cout << " --nonSymPb: non symmetric problem (<=> use dn[ae]upd)." + << endl; + cout << " default: symmetric problem (<=> use ds[ae]upd)" + << endl; + cout << " --cpxPb: complex (non symmetric) problem (<=> use " + "zn[ae]upd)." + << endl; + cout << " default: false (<=> use d*[ae]upd)" << endl; + cout << " --simplePrec: use simple precision (less accurate, but, " + "half memory footprint)." + << endl; + cout << " default: false (<=> use double precision: use " + "[dz]*upd)" + << endl; + cout << " --mag M: set magnitude of eigen values to look for " + "(LM, SM, LR, SR, LI, SI)." + << endl; + cout << " default: large magnitude (LM)" << endl; + cout << " --shiftReal S: real shift where sigma = S (look for lambda+S " + "instead of lambda)." + << endl; + cout << " default: no shift, S = 0." << endl; + cout << " --shiftImag S: imaginary shift where sigma = S (look for " + "lambda+S instead of lambda)." + << endl; + cout << " default: no shift, S = 0." << endl; + cout << " --invert: invert mode (look for 1./lambda instead of " + "lambda)." + << endl; + cout << " default: no invert" << endl; + cout << " --tol T: tolerance T." << endl; + cout << " default: 1.e-06" << endl; + cout << " --maxIt M: maximum iterations M." << endl; + cout << " default: 100" << endl; + cout << " --schur: compute Schur vectors." << endl; + cout << " the Schur decomposition is such that A = " + "Q^H x T x Q where:" + << endl; + cout << " - the H superscript refers to the " + "Hermitian transpose: Q^H = (Q^t)^*." + << endl; + cout + << " - Q is unitary: Q is such that Q^H x Q = I." + << endl; + cout << " - T is an upper-triangular matrix whose " + "diagonal elements are the eigenvalues of A." + << endl; + cout << " every square matrix has a Schur " + "decomposition: columns of Q are the Schur vectors." + << endl; + cout << " for a general matrix A, there is no " + "relation between Schur vectors of A and eigenvectors of A." + << endl; + cout << " if q_j is the j-th Schur vector, then A x " + "q_j is a linear combination of q_1, ..., q_j." + << endl; + cout << " Schur vectors q_1, q_2, ..., q_j span an " + "invariant subspace of A." + << endl; + cout << " the Schur vectors and eigenvectors of A are " + "the same if A is a normal matrix." + << endl; + cout << " default: compute Ritz vectors (approximations " + "of eigen vectors)" + << endl; + cout << " --slv S: solver (needed if arpack mode > 1)." << endl; + cout << " BiCG: iterative method, any matrices" + << endl; + cout << " CG: iterative method, sym matrices only" + << endl; + cout << " LU: direct method, any matrices (pivoting " + "needed)" + << endl; + cout << " QR: direct method, any matrices (pivoting " + "needed)" + << endl; + cout << " LLT: direct method, SPD matrices only " + "(pivoting not needed)" + << endl; + cout << " LDLT: direct method, symmetric positive " + "semi-definite matrices only (pivoting not needed)" + << endl; + cout << " default: BiCG" << endl; + cout << " --slvItrTol T: solver tolerance T (for iterative solvers)." + << endl; + cout << " default: 1.e-6" << endl; + cout << " --slvItrMaxIt M: solver maximum iterations M (for iterative " + "solvers)." + << endl; + cout << " default: 100" << endl; + cout << " --slvItrPC PC: solver preconditioner (for iterative solvers)." + << endl; + cout << " PC preconditioner:" << endl; + cout << " Diag: eigen diagonal preconditioner " + "(Jacobi)." + << endl; + cout << " ILU#D#F: eigen ILU preconditioner." + << endl; + cout << " D: drop tolerance." << endl; + cout << " F: fill factor." << endl; + cout << " default: diagonal preconditioner (Jacobi)" + << endl; + cout << " --slvDrtPivot: P solver pivot P (for direct solvers)." << endl; + cout << " default: 1.e-06" << endl; + cout << " --slvDrtOffset: O solver offset O (for direct solvers)." << endl; + cout << " default: 0." << endl; + cout << " --slvDrtScale: S solver scale S (for direct solvers)." << endl; + cout << " default: 1." << endl; + cout << " --noCheck: check arpack eigen values/vectors." << endl; + cout << " check will fail if Schur vectors are computed " + "and A is NOT a normal matrix." + << endl; + cout << " default: check" << endl; + cout << " --verbose V: verbosity level (up to 3)." << endl; + cout << " default: 0" << endl; + cout << " --debug D: debug level (up to 3)." << endl; + cout << " default: 0" << endl; + cout << " --restart: restart from previous run (which had produced " + "arpackSolver.*.out)." + << endl; + cout << " restart from eigen basis approximation " + "computed during a previous run." + << endl; + cout << " default: false" << endl; + if (rc == 0) exit(0); + return rc; + }; + + friend ostream& operator<<(ostream& ostr, options const& opt); + + string fileA; + string fileB; + bool dense; + bool denseRR; + a_int nbEV; + a_int nbCV; + bool stdPb; // Standard or generalized (= not standard). + bool symPb; + bool cpxPb; + bool simplePrec; + string mag; // Magnitude <=> "which" arpack parameter. + bool shiftReal, shiftImag; + double sigmaReal, sigmaImag; // Eigen value translation: look for + // lambda+sigma instead of lambda. + bool invert; // Eigen value invertion: look for 1./lambda instead of lambda. + double tol; + int maxIt; + bool schur; + string slv; + double slvItrTol; + int slvItrMaxIt; + string slvItrPC; + double slvDrtPivot; + double slvDrtOffset; + double slvDrtScale; + bool check; + int verbose; + a_int debug; + bool restart; }; -ostream & operator<< (ostream & ostr, options const & opt) { +ostream& operator<<(ostream& ostr, options const& opt) { ostr << "OPT: A " << opt.fileA << ", B " << opt.fileB; - if (opt.dense && opt.denseRR) ostr << ", dense yes (RR true)"; - else if (opt.dense && !opt.denseRR) ostr << ", dense yes (RR false)"; - else ostr << ", dense no"; - ostr << ", nbEV " << opt.nbEV << ", nbCV " << opt.nbCV << ", stdPb " << (opt.stdPb ? "yes" : "no"); - ostr << ", symPb " << (opt.symPb ? "yes" : "no") << ", cpxPb " << (opt.cpxPb ? "yes" : "no"); - ostr << ", simplePrec " << (opt.simplePrec ? "yes" : "no") << ", mag " << opt.mag << endl; - ostr << "OPT: shiftReal " << (opt.shiftReal ? "yes" : "no") << ", sigmaReal " << opt.sigmaReal; - ostr << ", shiftImag " << (opt.shiftImag ? "yes" : "no") << ", sigmaImag " << opt.sigmaImag; - ostr << ", invert " << (opt.invert ? "yes" : "no") << ", tol " << opt.tol << ", maxIt " << opt.maxIt; + if (opt.dense && opt.denseRR) + ostr << ", dense yes (RR true)"; + else if (opt.dense && !opt.denseRR) + ostr << ", dense yes (RR false)"; + else + ostr << ", dense no"; + ostr << ", nbEV " << opt.nbEV << ", nbCV " << opt.nbCV << ", stdPb " + << (opt.stdPb ? "yes" : "no"); + ostr << ", symPb " << (opt.symPb ? "yes" : "no") << ", cpxPb " + << (opt.cpxPb ? "yes" : "no"); + ostr << ", simplePrec " << (opt.simplePrec ? "yes" : "no") << ", mag " + << opt.mag << endl; + ostr << "OPT: shiftReal " << (opt.shiftReal ? "yes" : "no") << ", sigmaReal " + << opt.sigmaReal; + ostr << ", shiftImag " << (opt.shiftImag ? "yes" : "no") << ", sigmaImag " + << opt.sigmaImag; + ostr << ", invert " << (opt.invert ? "yes" : "no") << ", tol " << opt.tol + << ", maxIt " << opt.maxIt; ostr << ", " << (opt.schur ? "Schur" : "Ritz") << " vectors" << endl; - ostr << "OPT: slv " << opt.slv << ", slvItrPC " << opt.slvItrPC << ", slvItrTol " << opt.slvItrTol; - ostr << ", slvItrMaxIt " << opt.slvItrMaxIt << ", slvDrtPivot " << opt.slvDrtPivot; - ostr << ", slvDrtOffset " << opt.slvDrtOffset << ", slvDrtScale " << opt.slvDrtScale << endl; - ostr << "OPT: check " << (opt.check ? "yes" : "no") << ", verbose " << opt.verbose << ", debug " << opt.debug; + ostr << "OPT: slv " << opt.slv << ", slvItrPC " << opt.slvItrPC + << ", slvItrTol " << opt.slvItrTol; + ostr << ", slvItrMaxIt " << opt.slvItrMaxIt << ", slvDrtPivot " + << opt.slvDrtPivot; + ostr << ", slvDrtOffset " << opt.slvDrtOffset << ", slvDrtScale " + << opt.slvDrtScale << endl; + ostr << "OPT: check " << (opt.check ? "yes" : "no") << ", verbose " + << opt.verbose << ", debug " << opt.debug; ostr << ", restart " << (opt.restart ? "yes" : "no") << endl; return ostr; } class output { - public: - output() { - nbVal = 0; - mode = 0; - nbIt = 0; - imsTime = 0.; - rciTime = 0.; - }; + public: + output() { + nbVal = 0; + mode = 0; + nbIt = 0; + imsTime = 0.; + rciTime = 0.; + }; - int nbVal; // Eigen values. - int mode; // Arpack mode. - int nbIt; // Arpack number of iterations. - double imsTime; // Init mode solver time. - double rciTime; // Reverse communication interface time. + int nbVal; // Eigen values. + int mode; // Arpack mode. + int nbIt; // Arpack number of iterations. + double imsTime; // Init mode solver time. + double rciTime; // Reverse communication interface time. }; -template -int itrSolve(options & opt, output & out, - double const & slvItrILUDropTol, double const & slvItrILUFillFactor) { +template +int itrSolve(options& opt, output& out, double const& slvItrILUDropTol, + double const& slvItrILUFillFactor) { // Init solver. arpackItrSolver as; @@ -359,9 +605,13 @@ int itrSolve(options & opt, output & out, EM A; auto start = chrono::high_resolution_clock::now(); int rc = as.createMatrix(opt.fileA, A); - if (rc != 0) {cerr << "Error: read A KO" << endl; return rc;} + if (rc != 0) { + cerr << "Error: read A KO" << endl; + return rc; + } auto stop = chrono::high_resolution_clock::now(); - double readATime = chrono::duration_cast(stop - start).count()/1000.; + double readATime = + chrono::duration_cast(stop - start).count() / 1000.; cout << endl; cout << "INP: create A " << readATime << " s" << endl; @@ -371,23 +621,40 @@ int itrSolve(options & opt, output & out, if (!opt.stdPb) { start = chrono::high_resolution_clock::now(); rc = as.createMatrix(opt.fileB, B); - if (rc != 0) {cerr << "Error: read B KO" << endl; return rc;} + if (rc != 0) { + cerr << "Error: read B KO" << endl; + return rc; + } stop = chrono::high_resolution_clock::now(); - double readBTime = chrono::duration_cast(stop - start).count()/1000.; + double readBTime = + chrono::duration_cast(stop - start).count() / + 1000.; cout << endl; cout << "INP: create B " << readBTime << " s" << endl; - if (A.rows() != B.rows()) {cerr << "Error: A.rows() != B.rows()" << endl; return rc;} - if (A.cols() != B.cols()) {cerr << "Error: A.cols() != B.cols()" << endl; return rc;} + if (A.rows() != B.rows()) { + cerr << "Error: A.rows() != B.rows()" << endl; + return rc; + } + if (A.cols() != B.cols()) { + cerr << "Error: A.cols() != B.cols()" << endl; + return rc; + } } // Solve. rc = as.solve(A, opt.stdPb ? nullptr : &B); - if (rc != 0) {cerr << "Error: solve KO" << endl; return rc;} + if (rc != 0) { + cerr << "Error: solve KO" << endl; + return rc; + } if (opt.check) { rc = as.checkEigVec(A, opt.stdPb ? nullptr : &B); - if (rc != 0) {cerr << "Error: check KO" << endl; return rc;} + if (rc != 0) { + cerr << "Error: check KO" << endl; + return rc; + } } // Retrieve outputs. @@ -401,8 +668,8 @@ int itrSolve(options & opt, output & out, return 0; } -template -int drtSolve(options & opt, output & out) { +template +int drtSolve(options& opt, output& out) { // Init solver. arpackDrtSolver as; @@ -427,9 +694,13 @@ int drtSolve(options & opt, output & out) { EM A; auto start = chrono::high_resolution_clock::now(); int rc = as.createMatrix(opt.fileA, A); - if (rc != 0) {cerr << "Error: read A KO" << endl; return rc;} + if (rc != 0) { + cerr << "Error: read A KO" << endl; + return rc; + } auto stop = chrono::high_resolution_clock::now(); - double readATime = chrono::duration_cast(stop - start).count()/1000.; + double readATime = + chrono::duration_cast(stop - start).count() / 1000.; cout << endl; cout << "INP: create A " << readATime << " s" << endl; @@ -439,23 +710,40 @@ int drtSolve(options & opt, output & out) { if (!opt.stdPb) { start = chrono::high_resolution_clock::now(); rc = as.createMatrix(opt.fileB, B); - if (rc != 0) {cerr << "Error: read B KO" << endl; return rc;} + if (rc != 0) { + cerr << "Error: read B KO" << endl; + return rc; + } stop = chrono::high_resolution_clock::now(); - double readBTime = chrono::duration_cast(stop - start).count()/1000.; + double readBTime = + chrono::duration_cast(stop - start).count() / + 1000.; cout << endl; cout << "INP: create B " << readBTime << " s" << endl; - if (A.rows() != B.rows()) {cerr << "Error: A.rows() != B.rows()" << endl; return rc;} - if (A.cols() != B.cols()) {cerr << "Error: A.cols() != B.cols()" << endl; return rc;} + if (A.rows() != B.rows()) { + cerr << "Error: A.rows() != B.rows()" << endl; + return rc; + } + if (A.cols() != B.cols()) { + cerr << "Error: A.cols() != B.cols()" << endl; + return rc; + } } // Solve. rc = as.solve(A, opt.stdPb ? nullptr : &B); - if (rc != 0) {cerr << "Error: solve KO" << endl; return rc;} + if (rc != 0) { + cerr << "Error: solve KO" << endl; + return rc; + } if (opt.check) { rc = as.checkEigVec(A, opt.stdPb ? nullptr : &B); - if (rc != 0) {cerr << "Error: check KO" << endl; return rc;} + if (rc != 0) { + cerr << "Error: check KO" << endl; + return rc; + } } // Retrieve outputs. @@ -469,181 +757,213 @@ int drtSolve(options & opt, output & out) { return 0; } -template -int drtSolve(options & opt, output & out) { +template +int drtSolve(options& opt, output& out) { int rc = 1; - if (opt.slv == "LU") rc = drtSolve(opt, out); - if (opt.slv == "QR") rc = drtSolve(opt, out); - if (opt.slv == "LLT") rc = drtSolve(opt, out); + if (opt.slv == "LU") rc = drtSolve(opt, out); + if (opt.slv == "QR") rc = drtSolve(opt, out); + if (opt.slv == "LLT") rc = drtSolve(opt, out); if (opt.slv == "LDLT") rc = drtSolve(opt, out); return rc; } -template -int itrSolve(options & opt, output & out) { +template +int itrSolve(options& opt, output& out) { int rc = 1; stringstream clo(opt.slvItrPC); - string slvItrPC; getline(clo, slvItrPC, '#'); + string slvItrPC; + getline(clo, slvItrPC, '#'); double slvItrILUDropTol = 1.; if (slvItrPC == "ILU") { - string dropTol; getline(clo, dropTol, '#'); - stringstream dt(dropTol); dt >> slvItrILUDropTol; + string dropTol; + getline(clo, dropTol, '#'); + stringstream dt(dropTol); + dt >> slvItrILUDropTol; } int slvItrILUFillFactor = 2; if (slvItrPC == "ILU") { - string fillFactor; getline(clo, fillFactor); - stringstream ff(fillFactor); ff >> slvItrILUFillFactor; + string fillFactor; + getline(clo, fillFactor); + stringstream ff(fillFactor); + ff >> slvItrILUFillFactor; } if (opt.slv == "BiCG") { - if (slvItrPC == "Diag") rc = itrSolve(opt, out, slvItrILUDropTol, slvItrILUFillFactor); - if (slvItrPC == "ILU") rc = itrSolve(opt, out, slvItrILUDropTol, slvItrILUFillFactor); + if (slvItrPC == "Diag") + rc = itrSolve(opt, out, slvItrILUDropTol, + slvItrILUFillFactor); + if (slvItrPC == "ILU") + rc = itrSolve(opt, out, slvItrILUDropTol, + slvItrILUFillFactor); } - if (opt.slv == "CG") { - if (slvItrPC == "Diag") rc = itrSolve(opt, out, slvItrILUDropTol, slvItrILUFillFactor); - if (slvItrPC == "ILU") rc = itrSolve(opt, out, slvItrILUDropTol, slvItrILUFillFactor); + if (opt.slv == "CG") { + if (slvItrPC == "Diag") + rc = itrSolve(opt, out, slvItrILUDropTol, + slvItrILUFillFactor); + if (slvItrPC == "ILU") + rc = itrSolve(opt, out, slvItrILUDropTol, + slvItrILUFillFactor); } return rc; } -int main(int argc, char ** argv) { +int main(int argc, char** argv) { // Check for options. options opt; int rc = opt.readCmdLine(argc, argv); - if (rc != 0) {cerr << "Error: read cmd line KO" << endl; return rc;} - cout << opt; // Print options. + if (rc != 0) { + cerr << "Error: read cmd line KO" << endl; + return rc; + } + cout << opt; // Print options. // Solve with arpack. - sstats_c(); // Reset timers. - sstatn_c(); // Reset timers. - cstatn_c(); // Reset timers. + sstats_c(); // Reset timers. + sstatn_c(); // Reset timers. + cstatn_c(); // Reset timers. - bool itrSlv = true; // Use iterative solvers. - if (opt.slv.find("LU") != string::npos || opt.slv.find("QR") != string::npos || - opt.slv.find("LLT") != string::npos || opt.slv.find("LDLT") != string::npos ) itrSlv = false; + bool itrSlv = true; // Use iterative solvers. + if (opt.slv.find("LU") != string::npos || + opt.slv.find("QR") != string::npos || + opt.slv.find("LLT") != string::npos || + opt.slv.find("LDLT") != string::npos) + itrSlv = false; output out; auto start = chrono::high_resolution_clock::now(); if (opt.dense) { if (itrSlv) { - cerr << "Error: dense matrices does not support iterative solvers" << endl; + cerr << "Error: dense matrices does not support iterative solvers" + << endl; return 1; } if (opt.simplePrec) { if (opt.cpxPb) { if (opt.denseRR) { - rc = drtSolve, float, EigDMxC, EigDFLUC, EigDFQRC, EigDLLTC, EigDLDLTC>(opt, out); + rc = drtSolve, float, EigDMxC, EigDFLUC, EigDFQRC, + EigDLLTC, EigDLDLTC>(opt, out); + } else { + rc = drtSolve, float, EigDMxC, EigDPLUC, EigDPQRC, + EigDLLTC, EigDLDLTC>(opt, out); } - else { - rc = drtSolve, float, EigDMxC, EigDPLUC, EigDPQRC, EigDLLTC, EigDLDLTC>(opt, out); - } - } - else { + } else { if (opt.denseRR) { - rc = drtSolve< float , float, EigDMxS, EigDFLUS, EigDFQRS, EigDLLTS, EigDLDLTS>(opt, out); - } - else { - rc = drtSolve< float , float, EigDMxS, EigDPLUS, EigDPQRS, EigDLLTS, EigDLDLTS>(opt, out); + rc = drtSolve(opt, out); + } else { + rc = drtSolve(opt, out); } } - } - else { + } else { if (opt.cpxPb) { if (opt.denseRR) { - rc = drtSolve, double, EigDMxZ, EigDFLUZ, EigDFQRZ, EigDLLTZ, EigDLDLTZ>(opt, out); + rc = drtSolve, double, EigDMxZ, EigDFLUZ, EigDFQRZ, + EigDLLTZ, EigDLDLTZ>(opt, out); + } else { + rc = drtSolve, double, EigDMxZ, EigDPLUZ, EigDPQRZ, + EigDLLTZ, EigDLDLTZ>(opt, out); } - else { - rc = drtSolve, double, EigDMxZ, EigDPLUZ, EigDPQRZ, EigDLLTZ, EigDLDLTZ>(opt, out); - } - } - else { + } else { if (opt.denseRR) { - rc = drtSolve< double , double, EigDMxD, EigDFLUD, EigDFQRD, EigDLLTD, EigDLDLTD>(opt, out); - } - else { - rc = drtSolve< double , double, EigDMxD, EigDPLUD, EigDPQRD, EigDLLTD, EigDLDLTD>(opt, out); + rc = drtSolve(opt, out); + } else { + rc = drtSolve(opt, out); } } } - } - else { + } else { if (opt.simplePrec) { if (opt.cpxPb) { if (itrSlv) { - rc = itrSolve, float, EigSMxC, EigSBiCGC, EigSBiCGILUC, EigSCGC, EigSCGILUC>(opt, out); + rc = itrSolve, float, EigSMxC, EigSBiCGC, EigSBiCGILUC, + EigSCGC, EigSCGILUC>(opt, out); + } else { + rc = drtSolve, float, EigSMxC, EigSLUC, EigSQRC, + EigSLLTC, EigSLDLTC>(opt, out); } - else { - rc = drtSolve, float, EigSMxC, EigSLUC, EigSQRC, EigSLLTC, EigSLDLTC >(opt, out); - } - } - else { + } else { if (itrSlv) { - rc = itrSolve< float , float, EigSMxS, EigSBiCGS, EigSBiCGILUS, EigSCGS, EigSCGILUS>(opt, out); - } - else { - rc = drtSolve< float , float, EigSMxS, EigSLUS, EigSQRS, EigSLLTS, EigSLDLTS >(opt, out); + rc = itrSolve(opt, out); + } else { + rc = drtSolve(opt, out); } } - } - else { + } else { if (opt.cpxPb) { if (itrSlv) { - rc = itrSolve, double, EigSMxZ, EigSBiCGZ, EigSBiCGILUZ, EigSCGZ, EigSCGILUZ>(opt, out); + rc = itrSolve, double, EigSMxZ, EigSBiCGZ, + EigSBiCGILUZ, EigSCGZ, EigSCGILUZ>(opt, out); + } else { + rc = drtSolve, double, EigSMxZ, EigSLUZ, EigSQRZ, + EigSLLTZ, EigSLDLTZ>(opt, out); } - else { - rc = drtSolve, double, EigSMxZ, EigSLUZ, EigSQRZ, EigSLLTZ, EigSLDLTZ >(opt, out); - } - } - else { + } else { if (itrSlv) { - rc = itrSolve< double , double, EigSMxD, EigSBiCGD, EigSBiCGILUD, EigSCGD, EigSCGILUD>(opt, out); - } - else { - rc = drtSolve< double , double, EigSMxD, EigSLUD, EigSQRD, EigSLLTD, EigSLDLTD >(opt, out); + rc = itrSolve(opt, out); + } else { + rc = drtSolve(opt, out); } } } } - if (rc != 0) {cerr << "Error: arpack solve KO" << endl; return rc;} + if (rc != 0) { + cerr << "Error: arpack solve KO" << endl; + return rc; + } // Output results and stats. auto stop = chrono::high_resolution_clock::now(); - double fullTime = chrono::duration_cast(stop - start).count()/1000.; + double fullTime = + chrono::duration_cast(stop - start).count() / 1000.; cout << endl; - cout << "OUT: mode " << out.mode << ", nb EV found " << out.nbVal << ", nb iterations " << out.nbIt << endl; - cout << "OUT: init mode solver " << out.imsTime << " s, RCI time " << out.rciTime << " s" << endl; + cout << "OUT: mode " << out.mode << ", nb EV found " << out.nbVal + << ", nb iterations " << out.nbIt << endl; + cout << "OUT: init mode solver " << out.imsTime << " s, RCI time " + << out.rciTime << " s" << endl; cout << "OUT: full time " << fullTime << " s" << endl; a_int nopx = 0, nbx = 0, nrorth = 0, nitref = 0, nrstrt = 0; - float tsaupd = 0., tsaup2 = 0., tsaitr = 0., tseigt = 0., tsgets = 0., tsapps = 0., tsconv = 0.; - float tnaupd = 0., tnaup2 = 0., tnaitr = 0., tneigt = 0., tngets = 0., tnapps = 0., tnconv = 0.; - float tcaupd = 0., tcaup2 = 0., tcaitr = 0., tceigt = 0., tcgets = 0., tcapps = 0., tcconv = 0.; + float tsaupd = 0., tsaup2 = 0., tsaitr = 0., tseigt = 0., tsgets = 0., + tsapps = 0., tsconv = 0.; + float tnaupd = 0., tnaup2 = 0., tnaitr = 0., tneigt = 0., tngets = 0., + tnapps = 0., tnconv = 0.; + float tcaupd = 0., tcaup2 = 0., tcaitr = 0., tceigt = 0., tcgets = 0., + tcapps = 0., tcconv = 0.; float tmvopx = 0., tmvbx = 0., tgetv0 = 0., titref = 0., trvec = 0.; - stat_c(nopx, nbx, nrorth, nitref, nrstrt, tsaupd, tsaup2, - tsaitr, tseigt, tsgets, tsapps, tsconv, tnaupd, tnaup2, - tnaitr, tneigt, tngets, tnapps, tnconv, tcaupd, tcaup2, - tcaitr, tceigt, tcgets, tcapps, tcconv, tmvopx, tmvbx, - tgetv0, titref, trvec); + stat_c(nopx, nbx, nrorth, nitref, nrstrt, tsaupd, tsaup2, tsaitr, tseigt, + tsgets, tsapps, tsconv, tnaupd, tnaup2, tnaitr, tneigt, tngets, tnapps, + tnconv, tcaupd, tcaup2, tcaitr, tceigt, tcgets, tcapps, tcconv, tmvopx, + tmvbx, tgetv0, titref, trvec); cout << endl; - cout << "STAT: total number of user OP*x operation " << nopx << endl; - cout << "STAT: total number of user B*x operation " << nbx << endl; - cout << "STAT: total number of reorthogonalization steps taken " << nrorth << endl; - cout << "STAT: total number of it. refinement steps in reorthogonalization " << nitref << endl; - cout << "STAT: total number of restart steps " << nrstrt << endl; + cout << "STAT: total number of user OP*x operation " + << nopx << endl; + cout << "STAT: total number of user B*x operation " + << nbx << endl; + cout << "STAT: total number of reorthogonalization steps taken " + << nrorth << endl; + cout << "STAT: total number of it. refinement steps in reorthogonalization " + << nitref << endl; + cout << "STAT: total number of restart steps " + << nrstrt << endl; return 0; } diff --git a/EXAMPLES/PYARPACK/pyarpack.cpp b/EXAMPLES/PYARPACK/pyarpack.cpp index 775eb6b..1eaf5ff 100644 --- a/EXAMPLES/PYARPACK/pyarpack.cpp +++ b/EXAMPLES/PYARPACK/pyarpack.cpp @@ -1,86 +1,124 @@ -#include -#include -#include // ostringstream. -#include +#include // PyErr_SetString. #include #include #include - -#include // PyErr_SetString. - -#include +#include #include +#include +#include // ostringstream. +#include +#include namespace bp = boost::python; namespace bn = boost::python::numpy; -template -void exportArpackSparseItr(bp::scope & pySlv, std::string const & dtype) { +template +void exportArpackSparseItr(bp::scope& pySlv, std::string const& dtype) { // Created nested namespace in module. - pySlv.attr(dtype.c_str()) = bp::class_>(dtype.c_str(), - "arpack data type (must be consistent with numpy dtype)") - .def ("solve", &pyarpackSparseItrSolver::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::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::slvTol, - "tolerance of the iterative mode solver - default: 1.e-6") - .def_readwrite("slvMaxIt", &pyarpackSparseItrSolver::slvMaxIt, - "maximum number of iterations of the iterative mode solver - default: 100") - .def_readwrite("slvILUDropTol", &pyarpackSparseItrSolver::slvILUDropTol, - "drop tolerance of the ILU preconditioner (if any) of the iterative mode solver - default: 1") - .def_readwrite("slvILUFillFactor", &pyarpackSparseItrSolver::slvILUFillFactor, - "fill factor of the ILU preconditioner (if any) of the iterative mode solver - default: 2") - ; + pySlv.attr(dtype.c_str()) = + bp::class_>( + dtype.c_str(), + "arpack data type (must be consistent with numpy dtype)") + .def("solve", &pyarpackSparseItrSolver::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::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::slvTol, + "tolerance of the iterative mode solver - default: 1.e-6") + .def_readwrite("slvMaxIt", + &pyarpackSparseItrSolver::slvMaxIt, + "maximum number of iterations of the iterative mode " + "solver - default: 100") + .def_readwrite( + "slvILUDropTol", + &pyarpackSparseItrSolver::slvILUDropTol, + "drop tolerance of the ILU preconditioner (if any) of the " + "iterative mode solver - default: 1") + .def_readwrite( + "slvILUFillFactor", + &pyarpackSparseItrSolver::slvILUFillFactor, + "fill factor of the ILU preconditioner (if any) of the iterative " + "mode solver - default: 2"); }; -template -void exportArpackSparseDrt(bp::scope & pySlv, std::string const & dtype) { +template +void exportArpackSparseDrt(bp::scope& pySlv, std::string const& dtype) { // Created nested namespace in module. - pySlv.attr(dtype.c_str()) = bp::class_>(dtype.c_str(), - "arpack data type (must be consistent with numpy dtype)") - .def ("solve", &pyarpackSparseDrtSolver::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::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::slvPvtThd, - "pivoting tolerance of the direct mode solver - default: 1.e-6") - .def_readwrite("slvOffset", &pyarpackSparseDrtSolver::slvOffset, - "cholesky offset (LLT, LDLT) of the direct mode solver - default: 0.") - .def_readwrite("slvScale", &pyarpackSparseDrtSolver::slvScale, - "cholesky scale (LLT, LDLT) of the direct mode solver - default: 1.") - ; + pySlv.attr(dtype.c_str()) = + bp::class_>( + dtype.c_str(), + "arpack data type (must be consistent with numpy dtype)") + .def("solve", &pyarpackSparseDrtSolver::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::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::slvPvtThd, + "pivoting tolerance of the direct mode solver - default: 1.e-6") + .def_readwrite("slvOffset", + &pyarpackSparseDrtSolver::slvOffset, + "cholesky offset (LLT, LDLT) of the direct mode " + "solver - default: 0.") + .def_readwrite("slvScale", + &pyarpackSparseDrtSolver::slvScale, + "cholesky scale (LLT, LDLT) of the direct mode solver " + "- default: 1."); }; -template -void exportArpackDenseDrt(bp::scope & pySlv, std::string const & dtype) { +template +void exportArpackDenseDrt(bp::scope& pySlv, std::string const& dtype) { // Created nested namespace in module. - pySlv.attr(dtype.c_str()) = bp::class_>(dtype.c_str(), - "arpack data type (must be consistent with numpy dtype)") - .def ("solve", &pyarpackDenseDrtSolver::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::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::slvPvtThd, - "pivoting tolerance of the direct mode solver - default: 1.e-6") - .def_readwrite("slvOffset", &pyarpackDenseDrtSolver::slvOffset, - "cholesky offset (LLT, LDLT) of the direct mode solver - default: 0.") - .def_readwrite("slvScale", &pyarpackDenseDrtSolver::slvScale, - "cholesky scale (LLT, LDLT) of the direct mode solver - default: 1.") - ; + pySlv.attr(dtype.c_str()) = + bp::class_>( + dtype.c_str(), + "arpack data type (must be consistent with numpy dtype)") + .def("solve", &pyarpackDenseDrtSolver::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::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::slvPvtThd, + "pivoting tolerance of the direct mode solver - default: 1.e-6") + .def_readwrite("slvOffset", + &pyarpackDenseDrtSolver::slvOffset, + "cholesky offset (LLT, LDLT) of the direct mode " + "solver - default: 0.") + .def_readwrite("slvScale", + &pyarpackDenseDrtSolver::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 EigVecZGetItem(Eigen::Matrix, Eigen::Dynamic, 1> & M, int idx) { - if (idx < 0 || idx >= M.size()) {pyarpackThrowError("index out of range"); return std::complex();} +std::complex EigVecZGetItem( + Eigen::Matrix, Eigen::Dynamic, 1>& M, int idx) { + if (idx < 0 || idx >= M.size()) { + pyarpackThrowError("index out of range"); + return std::complex(); + } 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_>>("StdVecZ") - .def(bp::vector_indexing_suite>>()) - ; + bp::class_>>("StdVecZ").def( + bp::vector_indexing_suite>>()); bp::class_, Eigen::Dynamic, 1>>("EigVecZ") - .def("__getitem__", &EigVecZGetItem) - .def("__str__", &EigVecZToString) - ; + .def("__getitem__", &EigVecZGetItem) + .def("__str__", &EigVecZToString); bp::class_>("StdVecEVZ") - .def("__iter__", bp::iterator>()) - .def(bp::vector_indexing_suite>()) - ; + .def("__iter__", bp::iterator>()) + .def(bp::vector_indexing_suite>()); // 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_(slv.c_str(), slvHelp.c_str()); - exportArpackSparseItr< float , float, EigSMxS, EigSBiCGS>(pySlvBiCGDiag, "float" ); - exportArpackSparseItr< double , double, EigSMxD, EigSBiCGD>(pySlvBiCGDiag, "double"); - exportArpackSparseItr, float, EigSMxC, EigSBiCGC>(pySlvBiCGDiag, "complexFloat" ); - exportArpackSparseItr, double, EigSMxZ, EigSBiCGZ>(pySlvBiCGDiag, "complexDouble"); + std::string slvHelp = + "arpack internal mode solver (mode > 1): BiCG with diagonal (Jacobi) " + "preconditioner"; + bp::scope pySlvBiCGDiag = + bp::class_(slv.c_str(), slvHelp.c_str()); + exportArpackSparseItr(pySlvBiCGDiag, + "float"); + exportArpackSparseItr(pySlvBiCGDiag, + "double"); + exportArpackSparseItr, float, EigSMxC, EigSBiCGC>( + pySlvBiCGDiag, "complexFloat"); + exportArpackSparseItr, 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_(slv.c_str(), slvHelp.c_str()); - exportArpackSparseItr< float , float, EigSMxS, EigSBiCGILUS>(pySlvBiCGILU, "float" ); - exportArpackSparseItr< double , double, EigSMxD, EigSBiCGILUD>(pySlvBiCGILU, "double"); - exportArpackSparseItr, float, EigSMxC, EigSBiCGILUC>(pySlvBiCGILU, "complexFloat" ); - exportArpackSparseItr, double, EigSMxZ, EigSBiCGILUZ>(pySlvBiCGILU, "complexDouble"); + std::string slvHelp = + "arpack internal mode solver (mode > 1): BiCG with ILU preconditioner"; + bp::scope pySlvBiCGILU = + bp::class_(slv.c_str(), slvHelp.c_str()); + exportArpackSparseItr(pySlvBiCGILU, + "float"); + exportArpackSparseItr(pySlvBiCGILU, + "double"); + exportArpackSparseItr, float, EigSMxC, EigSBiCGILUC>( + pySlvBiCGILU, "complexFloat"); + exportArpackSparseItr, 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_(slv.c_str(), slvHelp.c_str()); - exportArpackSparseItr< float , float, EigSMxS, EigSCGS>(pySlvCGDiag, "float" ); - exportArpackSparseItr< double , double, EigSMxD, EigSCGD>(pySlvCGDiag, "double"); - exportArpackSparseItr, float, EigSMxC, EigSCGC>(pySlvCGDiag, "complexFloat" ); - exportArpackSparseItr, double, EigSMxZ, EigSCGZ>(pySlvCGDiag, "complexDouble"); + std::string slvHelp = + "arpack internal mode solver (mode > 1): CG with diagonal (Jacobi) " + "preconditioner"; + bp::scope pySlvCGDiag = + bp::class_(slv.c_str(), slvHelp.c_str()); + exportArpackSparseItr(pySlvCGDiag, "float"); + exportArpackSparseItr(pySlvCGDiag, + "double"); + exportArpackSparseItr, float, EigSMxC, EigSCGC>( + pySlvCGDiag, "complexFloat"); + exportArpackSparseItr, 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_(slv.c_str(), slvHelp.c_str()); - exportArpackSparseItr< float , float, EigSMxS, EigSCGILUS>(pySlvCGILU, "float" ); - exportArpackSparseItr< double , double, EigSMxD, EigSCGILUD>(pySlvCGILU, "double"); - exportArpackSparseItr, float, EigSMxC, EigSCGILUC>(pySlvCGILU, "complexFloat" ); - exportArpackSparseItr, double, EigSMxZ, EigSCGILUZ>(pySlvCGILU, "complexDouble"); + std::string slvHelp = + "arpack internal mode solver (mode > 1): CG with ILU preconditioner"; + bp::scope pySlvCGILU = + bp::class_(slv.c_str(), slvHelp.c_str()); + exportArpackSparseItr(pySlvCGILU, + "float"); + exportArpackSparseItr(pySlvCGILU, + "double"); + exportArpackSparseItr, float, EigSMxC, EigSCGILUC>( + pySlvCGILU, "complexFloat"); + exportArpackSparseItr, double, EigSMxZ, EigSCGILUZ>( + pySlvCGILU, "complexDouble"); } { std::string slv = "sparseLLT"; std::string slvHelp = "arpack internal mode solver (mode > 1): LLT"; bp::scope pySlvLLT = bp::class_(slv.c_str(), slvHelp.c_str()); - exportArpackSparseDrt< float , float, EigSMxS, EigSLLTS>(pySlvLLT, "float" ); - exportArpackSparseDrt< double , double, EigSMxD, EigSLLTD>(pySlvLLT, "double"); - exportArpackSparseDrt, float, EigSMxC, EigSLLTC>(pySlvLLT, "complexFloat" ); - exportArpackSparseDrt, double, EigSMxZ, EigSLLTZ>(pySlvLLT, "complexDouble"); + exportArpackSparseDrt(pySlvLLT, "float"); + exportArpackSparseDrt(pySlvLLT, + "double"); + exportArpackSparseDrt, float, EigSMxC, EigSLLTC>( + pySlvLLT, "complexFloat"); + exportArpackSparseDrt, double, EigSMxZ, EigSLLTZ>( + pySlvLLT, "complexDouble"); } { std::string slv = "sparseLDLT"; std::string slvHelp = "arpack internal mode solver (mode > 1): LDLT"; bp::scope pySlvLDLT = bp::class_(slv.c_str(), slvHelp.c_str()); - exportArpackSparseDrt< float , float, EigSMxS, EigSLDLTS>(pySlvLDLT, "float" ); - exportArpackSparseDrt< double , double, EigSMxD, EigSLDLTD>(pySlvLDLT, "double"); - exportArpackSparseDrt, float, EigSMxC, EigSLDLTC>(pySlvLDLT, "complexFloat" ); - exportArpackSparseDrt, double, EigSMxZ, EigSLDLTZ>(pySlvLDLT, "complexDouble"); + exportArpackSparseDrt(pySlvLDLT, "float"); + exportArpackSparseDrt(pySlvLDLT, + "double"); + exportArpackSparseDrt, float, EigSMxC, EigSLDLTC>( + pySlvLDLT, "complexFloat"); + exportArpackSparseDrt, double, EigSMxZ, EigSLDLTZ>( + pySlvLDLT, "complexDouble"); } { std::string slv = "sparseLU"; std::string slvHelp = "arpack internal mode solver (mode > 1): LU"; bp::scope pySlvLU = bp::class_(slv.c_str(), slvHelp.c_str()); - exportArpackSparseDrt< float , float, EigSMxS, EigSLUS>(pySlvLU, "float" ); - exportArpackSparseDrt< double , double, EigSMxD, EigSLUD>(pySlvLU, "double"); - exportArpackSparseDrt, float, EigSMxC, EigSLUC>(pySlvLU, "complexFloat" ); - exportArpackSparseDrt, double, EigSMxZ, EigSLUZ>(pySlvLU, "complexDouble"); + exportArpackSparseDrt(pySlvLU, "float"); + exportArpackSparseDrt(pySlvLU, "double"); + exportArpackSparseDrt, float, EigSMxC, EigSLUC>( + pySlvLU, "complexFloat"); + exportArpackSparseDrt, double, EigSMxZ, EigSLUZ>( + pySlvLU, "complexDouble"); } { std::string slv = "sparseQR"; std::string slvHelp = "arpack internal mode solver (mode > 1): QR"; bp::scope pySlvQR = bp::class_(slv.c_str(), slvHelp.c_str()); - exportArpackSparseDrt< float , float, EigSMxS, EigSQRS>(pySlvQR, "float" ); - exportArpackSparseDrt< double , double, EigSMxD, EigSQRD>(pySlvQR, "double"); - exportArpackSparseDrt, float, EigSMxC, EigSQRC>(pySlvQR, "complexFloat" ); - exportArpackSparseDrt, double, EigSMxZ, EigSQRZ>(pySlvQR, "complexDouble"); + exportArpackSparseDrt(pySlvQR, "float"); + exportArpackSparseDrt(pySlvQR, "double"); + exportArpackSparseDrt, float, EigSMxC, EigSQRC>( + pySlvQR, "complexFloat"); + exportArpackSparseDrt, double, EigSMxZ, EigSQRZ>( + pySlvQR, "complexDouble"); } { std::string slv = "denseLLT"; std::string slvHelp = "arpack internal mode solver (mode > 1): LLT"; bp::scope pySlvLLT = bp::class_(slv.c_str(), slvHelp.c_str()); - exportArpackDenseDrt< float , float, EigDMxS, EigDLLTS>(pySlvLLT, "float" ); - exportArpackDenseDrt< double , double, EigDMxD, EigDLLTD>(pySlvLLT, "double"); - exportArpackDenseDrt, float, EigDMxC, EigDLLTC>(pySlvLLT, "complexFloat" ); - exportArpackDenseDrt, double, EigDMxZ, EigDLLTZ>(pySlvLLT, "complexDouble"); + exportArpackDenseDrt(pySlvLLT, "float"); + exportArpackDenseDrt(pySlvLLT, "double"); + exportArpackDenseDrt, float, EigDMxC, EigDLLTC>( + pySlvLLT, "complexFloat"); + exportArpackDenseDrt, double, EigDMxZ, EigDLLTZ>( + pySlvLLT, "complexDouble"); } { std::string slv = "denseLDLT"; std::string slvHelp = "arpack internal mode solver (mode > 1): LDLT"; bp::scope pySlvLDLT = bp::class_(slv.c_str(), slvHelp.c_str()); - exportArpackDenseDrt< float , float, EigDMxS, EigDLDLTS>(pySlvLDLT, "float" ); - exportArpackDenseDrt< double , double, EigDMxD, EigDLDLTD>(pySlvLDLT, "double"); - exportArpackDenseDrt, float, EigDMxC, EigDLDLTC>(pySlvLDLT, "complexFloat" ); - exportArpackDenseDrt, double, EigDMxZ, EigDLDLTZ>(pySlvLDLT, "complexDouble"); + exportArpackDenseDrt(pySlvLDLT, "float"); + exportArpackDenseDrt(pySlvLDLT, + "double"); + exportArpackDenseDrt, float, EigDMxC, EigDLDLTC>( + pySlvLDLT, "complexFloat"); + exportArpackDenseDrt, 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_(slv.c_str(), slvHelp.c_str()); - exportArpackDenseDrt< float , float, EigDMxS, EigDFLUS>(pySlvLURR, "float" ); - exportArpackDenseDrt< double , double, EigDMxD, EigDFLUD>(pySlvLURR, "double"); - exportArpackDenseDrt, float, EigDMxC, EigDFLUC>(pySlvLURR, "complexFloat" ); - exportArpackDenseDrt, double, EigDMxZ, EigDFLUZ>(pySlvLURR, "complexDouble"); + exportArpackDenseDrt(pySlvLURR, "float"); + exportArpackDenseDrt(pySlvLURR, + "double"); + exportArpackDenseDrt, float, EigDMxC, EigDFLUC>( + pySlvLURR, "complexFloat"); + exportArpackDenseDrt, 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_(slv.c_str(), slvHelp.c_str()); - exportArpackDenseDrt< float , float, EigDMxS, EigDFQRS>(pySlvQRRR, "float" ); - exportArpackDenseDrt< double , double, EigDMxD, EigDFQRD>(pySlvQRRR, "double"); - exportArpackDenseDrt, float, EigDMxC, EigDFQRC>(pySlvQRRR, "complexFloat" ); - exportArpackDenseDrt, double, EigDMxZ, EigDFQRZ>(pySlvQRRR, "complexDouble"); + exportArpackDenseDrt(pySlvQRRR, "float"); + exportArpackDenseDrt(pySlvQRRR, + "double"); + exportArpackDenseDrt, float, EigDMxC, EigDFQRC>( + pySlvQRRR, "complexFloat"); + exportArpackDenseDrt, 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_(slv.c_str(), slvHelp.c_str()); - exportArpackDenseDrt< float , float, EigDMxS, EigDPLUS>(pySlvLUPP, "float" ); - exportArpackDenseDrt< double , double, EigDMxD, EigDPLUD>(pySlvLUPP, "double"); - exportArpackDenseDrt, float, EigDMxC, EigDPLUC>(pySlvLUPP, "complexFloat" ); - exportArpackDenseDrt, double, EigDMxZ, EigDPLUZ>(pySlvLUPP, "complexDouble"); + exportArpackDenseDrt(pySlvLUPP, "float"); + exportArpackDenseDrt(pySlvLUPP, + "double"); + exportArpackDenseDrt, float, EigDMxC, EigDPLUC>( + pySlvLUPP, "complexFloat"); + exportArpackDenseDrt, 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_(slv.c_str(), slvHelp.c_str()); - exportArpackDenseDrt< float , float, EigDMxS, EigDPQRS>(pySlvQPPR, "float" ); - exportArpackDenseDrt< double , double, EigDMxD, EigDPQRD>(pySlvQPPR, "double"); - exportArpackDenseDrt, float, EigDMxC, EigDPQRC>(pySlvQPPR, "complexFloat" ); - exportArpackDenseDrt, double, EigDMxZ, EigDPQRZ>(pySlvQPPR, "complexDouble"); + exportArpackDenseDrt(pySlvQPPR, "float"); + exportArpackDenseDrt(pySlvQPPR, + "double"); + exportArpackDenseDrt, float, EigDMxC, EigDPQRC>( + pySlvQPPR, "complexFloat"); + exportArpackDenseDrt, double, EigDMxZ, EigDPQRZ>( + pySlvQPPR, "complexDouble"); } } diff --git a/EXAMPLES/README.CALLING-ARPACK-FROM-C-OR-CPP b/EXAMPLES/README.CALLING-ARPACK-FROM-C-OR-CPP index 33b9a36..268524b 100644 --- a/EXAMPLES/README.CALLING-ARPACK-FROM-C-OR-CPP +++ b/EXAMPLES/README.CALLING-ARPACK-FROM-C-OR-CPP @@ -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++. diff --git a/PARPACK/EXAMPLES/MPI/README.CALLING-PARPACK-FROM-C-OR-CPP b/PARPACK/EXAMPLES/MPI/README.CALLING-PARPACK-FROM-C-OR-CPP index e10e5b8..3e8e17c 100644 --- a/PARPACK/EXAMPLES/MPI/README.CALLING-PARPACK-FROM-C-OR-CPP +++ b/PARPACK/EXAMPLES/MPI/README.CALLING-PARPACK-FROM-C-OR-CPP @@ -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++. diff --git a/PARPACK/TESTS/MPI/icb_parpack_c.c b/PARPACK/TESTS/MPI/icb_parpack_c.c index a2cb077..f73e3e0 100644 --- a/PARPACK/TESTS/MPI/icb_parpack_c.c +++ b/PARPACK/TESTS/MPI/icb_parpack_c.c @@ -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 // creal, cimag. +#include #include #include -#include + +#include "debug_c.h" // debug parpack. #include "mpi.h" #include "parpack.h" -#include // 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; diff --git a/PARPACK/TESTS/MPI/icb_parpack_cpp.cpp b/PARPACK/TESTS/MPI/icb_parpack_cpp.cpp index 91ab551..ce5f367 100644 --- a/PARPACK/TESTS/MPI/icb_parpack_cpp.cpp +++ b/PARPACK/TESTS/MPI/icb_parpack_cpp.cpp @@ -11,15 +11,14 @@ * with entries 1000, 999, ... , 2, 1 on the diagonal. */ -#include "parpack.hpp" - #include #include #include #include #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(1000 - (nev - 1) + i)) > 1. || - std::abs(std::imag(d[i]) - static_cast(1000 - (nev - 1) + i)) > 1.) { + if (std::abs(std::real(d[i]) - static_cast(1000 - (nev - 1) + i)) > + 1. || + std::abs(std::imag(d[i]) - static_cast(1000 - (nev - 1) + i)) > + 1.) { throw std::domain_error("Correct eigenvalues not computed"); } } diff --git a/TESTS/bug_1315_double.c b/TESTS/bug_1315_double.c index b7d1eba..66703e8 100644 --- a/TESTS/bug_1315_double.c +++ b/TESTS/bug_1315_double.c @@ -1,8 +1,8 @@ -#include "arpackdef.h" - +#include #include #include -#include + +#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); diff --git a/TESTS/bug_1315_single.c b/TESTS/bug_1315_single.c index b265d33..c469974 100644 --- a/TESTS/bug_1315_single.c +++ b/TESTS/bug_1315_single.c @@ -1,8 +1,8 @@ -#include "arpackdef.h" - +#include #include #include -#include + +#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); diff --git a/TESTS/icb_arpack_c.c b/TESTS/icb_arpack_c.c index be7e6ad..613cf1f 100644 --- a/TESTS/icb_arpack_c.c +++ b/TESTS/icb_arpack_c.c @@ -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 // creal, cimag. +#include #include #include -#include + #include "arpack.h" -#include // 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; } diff --git a/TESTS/icb_arpack_cpp.cpp b/TESTS/icb_arpack_cpp.cpp index 0e669a6..c3c4649 100644 --- a/TESTS/icb_arpack_cpp.cpp +++ b/TESTS/icb_arpack_cpp.cpp @@ -9,25 +9,24 @@ * matrix with entries 1000, 999, ... , 2, 1 on the diagonal. */ -#include "arpack.hpp" - #include #include #include #include +#include "arpack.hpp" #include "debug_c.hpp" // debug arpack. #include "stat_c.hpp" // arpack statistics. -template +template void diagonal_matrix_vector_product(Real const* const x, Real* const y) { for (int i = 0; i < 1000; ++i) { y[i] = static_cast(i + 1) * x[i]; } } -template -void real_symmetric_runner(double const & tol_check) { +template +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 +template void diagonal_matrix_vector_product(std::complex const* const x, std::complex* const y) { for (int i = 0; i < 1000; ++i) { @@ -106,8 +107,8 @@ void diagonal_matrix_vector_product(std::complex const* const x, } } -template -void complex_symmetric_runner(double const & tol_check) { +template +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(1000 - (nev - 1) + i)) > tol_check || - std::abs(std::imag(d[i]) - static_cast(1000 - (nev - 1) + i)) > tol_check) { + if (std::abs(std::real(d[i]) - static_cast(1000 - (nev - 1) + i)) > + tol_check || + std::abs(std::imag(d[i]) - static_cast(1000 - (nev - 1) + i)) > + tol_check) { throw std::domain_error("Correct eigenvalues not computed"); } }