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

Add Ritz options (largest/smallest in real/imaginary parts) for complex routines

Browse Source
This commit is contained in:
Zhentao Wang
2022-08-20 21:30:21 +02:00
committed by Franck HOUSSEN
parent 36eaaf0791
commit 8ecf1be713
6 changed files with 432 additions and 412 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CHANGES
+1
View File
@@ -2,6 +2,7 @@ arpack-ng - 3.9.0
[ Zhentao Wang ]
* [BUG FIX] parpack.h and parpack.hpp: type of rwork should be real instead of complex.
* Allow ritz_option {"LR", "SR", "LI", "SI"} for complex eigenvalue problems in ICB.
[ Jose E. Roman ]
* Avoid using isnan() in tests, since is GNU-specific
ICB/arpack.hpp
+24
View File
@@ -16,6 +16,14 @@ enum class which : int {
largest_magnitude,
/// 'SM' - compute the NEV smallest (in magnitude) eigenvalues.
smallest_magnitude,
/// 'LR' - compute the NEV largest (real part) eigenvalues.
largest_real,
/// 'SR' - compute the NEV smallest (real part) eigenvalues.
smallest_real,
/// 'LI' - compute the NEV largest (imaginary part) eigenvalues.
largest_imaginary,
/// 'SI' - compute the NEV smallest (imaginary part) eigenvalues.
smallest_imaginary,
/// 'BE' - compute NEV eigenvalues, half from each end of the
/// spectrum. When NEV is odd, compute one more from the
/// high end than from the low end.
@@ -60,6 +68,22 @@ inline char const* convert_to_char(which const option) {
return "SM";
break;
}
case which::largest_real: {
return "LR";
break;
}
case which::smallest_real: {
return "SR";
break;
}
case which::largest_imaginary: {
return "LI";
break;
}
case which::smallest_imaginary: {
return "SI";
break;
}
case which::both_ends: {
return "BE";
break;
PARPACK/TESTS/MPI/icb_parpack_c.c
+92 -95
View File
@@ -13,149 +13,151 @@
#include <stdio.h>
#include <stdlib.h>
#include "debug_c.h" // debug parpack.
#include "mpi.h"
#include "parpack.h"
#include "stat_c.h" // arpack statistics.
#include "debug_c.h" // debug parpack.
#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(const double* x, double* y) {
int i;
for (i = 0; i < 1000; ++i) y[i] = ((double)(i + 1)) * x[i];
};
int ds() {
a_int ido = 0;
char bmat[] = "I";
a_int N = 1000;
char which[] = "LM";
a_int nev = 3;
double tol = 0.000001; // small tol => more stable checks after EV computation.
const a_int N = 1000;
const a_int nev = 9;
const a_int ncv = 2 * nev + 1;
const a_int ldv = N;
const a_int ldz = N;
const a_int lworkl = ncv * (ncv + 8);
const a_int rvec = 1; // need eigenvectors
const double tol = 0.000001; // small tol => more stable checks after EV computation.
const double sigma = 0; // not referenced in this mode
double resid[N];
a_int ncv = 2 * nev + 1;
double V[ncv * N];
a_int ldv = N;
a_int iparam[11];
a_int ipntr[14];
double V[ldv * ncv];
double z[ldz * nev];
double d[nev];
double workd[3 * N];
a_int rvec = 1;
char howmny[] = "A";
double* d = (double*)malloc((nev + 1) * sizeof(double));
double workl[lworkl];
a_int select[ncv];
int i; // C99 compliant.
for (i = 0; i < ncv; i++) select[i] = 1;
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;
a_int info = 0;
a_int iparam[11], ipntr[11];
iparam[0] = 1; // ishift
iparam[2] = 10 * N; // on input: maxit; on output: actual iteration
iparam[3] = 1; // NB, only 1 allowed
iparam[6] = 1; // mode
char bmat[] = "I";
char which[] = "LM";
char howmny[] = "A";
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
iparam[0] = 1;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[4] = 0; // number of ev found by arpack.
iparam[6] = 1;
MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
while (ido != 99) {
/* call arpack like you would have, but, use dsaupd_c instead of dsaupd_ */
a_int info = 0, ido = 0;
do {
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]));
}
if (iparam[4] != nev) {printf("Error: iparam[4] %d, nev %d\n", iparam[4], nev); return 1;} // check number of ev found by arpack.
} while (ido == 1 || ido == -1);
// check info and number of ev found by arpack.
if (info < 0 || iparam[4] < nev) {
printf("Error in saupd: iparam[4] %d, nev %d, info %d\n", iparam[4], nev, info);
return 1;
}
/* 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);
tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl, &info);
if (info < 0) {
printf("Error in seupd: info %d\n", info);
return 1;
}
int i; // C99 compliant.
for (i = 0; i < nev; ++i) {
double val = d[i];
double ref = (N-(nev-1)+i);
double eps = fabs(val - ref);
printf("rank %d : %f - %f - %f\n", rank, val, ref, eps);
/*eigen value order: smallest -> biggest*/
if(eps>1.e-05){
free(d);
return 1;
}
if (eps > 1.e-05) return 1;
}
free(d);
return 0;
}
void zMatVec(a_dcomplex* x, a_dcomplex* y) {
void zMatVec(const a_dcomplex* x, a_dcomplex* y) {
int i;
for (i = 0; i < 1000; ++i) y[i] = x[i] * CMPLXF(i + 1.0, i + 1.0);
};
int zn() {
a_int ido = 0;
char bmat[] = "I";
a_int N = 1000;
char which[] = "LM";
a_int nev = 3;
double tol = 0.000001; // small tol => more stable checks after EV computation.
const a_int N = 1000;
const a_int nev = 9;
const a_int ncv = 2 * nev + 1;
const a_int ldv = N;
const a_int ldz = N;
const a_int lworkl = ncv * (3 * ncv + 5);
const a_int rvec = 0; // eigenvectors omitted
const double tol = 0.000001; // small tol => more stable checks after EV computation.
const a_dcomplex sigma = CMPLX(0., 0.); // not referenced in this mode
a_dcomplex resid[N];
a_int ncv = 2 * nev + 1;
a_dcomplex V[ncv * N];
a_int ldv = N;
a_int iparam[11];
a_int ipntr[14];
a_dcomplex V[ldv * ncv];
a_dcomplex z[ldz * nev];
a_dcomplex d[nev];
a_dcomplex workd[3 * N];
a_int rvec = 0;
char howmny[] = "A";
a_dcomplex* d =
(a_dcomplex*)malloc((nev + 1) * sizeof(a_dcomplex));
a_int select[ncv];
int i; // C99 compliant.
for (i = 0; i < ncv; i++) select[i] = 1;
a_dcomplex z[(N + 1) * (nev + 1)];
a_int ldz = N + 1;
a_dcomplex sigma = CMPLX(0., 0.);
int k;
for (k = 0; k < 3 * N; ++k) workd[k] = CMPLX(0., 0.);
a_dcomplex workl[3 * (ncv * ncv) + 6 * ncv];
for (k = 0; k < 3 * (ncv * ncv) + 6 * ncv; ++k) workl[k] = CMPLX(0., 0.);
a_int lworkl = 3 * (ncv * ncv) + 6 * ncv;
double rwork[ncv];
a_dcomplex workl[lworkl];
a_dcomplex workev[2 * ncv];
a_int info = 0;
double rwork[ncv];
a_int select[ncv];
a_int iparam[11], ipntr[14];
iparam[0] = 1; // ishift
iparam[2] = 10 * N; // on input: maxit; on output: actual iteration
iparam[3] = 1; // NB, only 1 allowed
iparam[6] = 1; // mode
char bmat[] = "I";
char which[] = "LM";
char howmny[] = "A";
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
iparam[0] = 1;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[4] = 0; // number of ev found by arpack.
iparam[6] = 1;
MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
while (ido != 99) {
/* call arpack like you would have, but, use znaupd_c instead of znaupd_ */
a_int info = 0, ido = 0;
do {
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]));
} while (ido == 1 || ido == -1);
// check info and number of ev found by arpack.
if (info < 0 || iparam[4] < nev) {
printf("Error in naupd: iparam[4] %d, nev %d, info %d\n", iparam[4], nev, info);
return 1;
}
if (iparam[4] != nev) {printf("Error: iparam[4] %d, nev %d\n", 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);
if (info < 0) {
printf("Error in neupd: info %d\n", info);
return 1;
}
int i; // C99 compliant.
for (i = 0; i < nev; ++i) {
double rval = creal(d[i]);
double rref = (N-(nev-1)+i);
@@ -165,13 +167,8 @@ int zn() {
double ieps = fabs(ival - iref);
printf("rank %d : %f %f - %f %f - %f %f\n", rank, rval, ival, rref, iref, reps, ieps);
/*eigen value order: smallest -> biggest*/
if(reps>1.e-05 || ieps>1.e-05){
free(d);
return 1;
}
if (reps > 1.e-05 || ieps > 1.e-05) return 1;
}
free(d);
return 0;
}
PARPACK/TESTS/MPI/icb_parpack_cpp.cpp
+123 -120
View File
@@ -6,190 +6,187 @@
* Just use arpack as you would have normally done but use [ae]upd instead
* of *[ae]upd_. The main advantage is that checks of the arguments are
* performed at compile time. Note: to debug parpack, call debug_c.
* This test program solves for the 9 largest eigenvalues of
* This test program solves for the 9 eigenvalues of
* A*x = lambda*x where A is the diagonal matrix
* with entries 1000, 999, ... , 2, 1 on the diagonal.
*/
#include <array>
#include <cmath>
#include <iostream>
#include <vector>
#include "debug_c.hpp" // debug parpack.
#include "parpack.hpp"
#include "stat_c.hpp" // arpack statistics.
#include "debug_c.hpp" // debug parpack.
#include "stat_c.hpp" // arpack statistics.
void diagonal_matrix_vector_product(float const* const x, float* const y) {
template <typename Real>
void diagonal_matrix_vector_product(const Real* x, Real* y) {
for (int i = 0; i < 1000; ++i) {
y[i] = static_cast<float>(i + 1) * x[i];
}
}
void real_symmetric_runner() {
a_int N = 1000;
a_int nev = 9;
a_int ncv = 2 * nev + 1;
a_int ldz = N + 1;
a_int lworkl = 3 * (ncv * ncv) + 6 * ncv;
a_int ldv = N;
template <typename Real>
void real_symmetric_runner(double const& tol_check, arpack::which const& ritz_option) {
const a_int N = 1000;
const a_int nev = 9;
const a_int ncv = 2 * nev + 1;
const a_int ldv = N;
const a_int ldz = N;
const a_int lworkl = ncv * (ncv + 8);
const a_int rvec = 1; // need eigenvectors
a_int rvec = 1;
float tol = 0.000001; // small tol => more stable checks after EV computation.
float sigma = 0.0f;
const Real tol = 0.000001; // small tol => more stable checks after EV computation.
const Real sigma = 0.0f; // not referenced in this mode
std::array<a_int, 14> ipntr;
std::vector<Real> resid(N);
std::vector<Real> V(ldv * ncv);
std::vector<Real> z(ldz * nev);
std::vector<Real> d(nev);
std::vector<Real> workd(3 * N);
std::vector<Real> workl(lworkl);
std::vector<a_int> select(ncv); // since HOWMNY = 'A', only used as workspace here
std::vector<float> workd(3 * N, 0.0f);
std::vector<float> workl(3 * (ncv * ncv) + 6 * ncv, 0.0f);
std::vector<float> V(ncv * N);
std::vector<float> d(nev + 1);
std::vector<float> z((N + 1) * (nev + 1));
std::vector<float> resid(N);
std::vector<a_int> select(ncv);
for (int i = 0; i < ncv; i++) select[i] = 1;
a_int info = 0;
a_int iparam[11], ipntr[11];
iparam[0] = 1; // ishift
iparam[2] = 10 * N; // on input: maxit; on output: actual iteration
iparam[3] = 1; // NB, only 1 allowed
iparam[6] = 1; // mode
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
std::array<a_int, 11> iparam;
iparam[0] = 1;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[4] = 0; // number of ev found by arpack.
iparam[6] = 1;
MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
a_int ido = 0;
while (ido != 99) {
a_int info = 0, ido = 0;
do {
arpack::saupd(MCW, ido, arpack::bmat::identity, N,
arpack::which::largest_magnitude, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam.data(), ipntr.data(), workd.data(),
ritz_option, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam, ipntr, workd.data(),
workl.data(), lworkl, info);
diagonal_matrix_vector_product(&(workd[ipntr[0] - 1]),
&(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
diagonal_matrix_vector_product(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
} while (ido == 1 || ido == -1);
// check info and number of ev found by arpack.
if (info < 0 || iparam[4] < nev) { /*arpack may succeed to compute more EV than expected*/
std::cout << "ERROR in saupd: iparam[4] " << iparam[4] << ", nev " << nev
<< ", info " << info << std::endl;
throw std::domain_error("Error inside ARPACK routines");
}
arpack::seupd(MCW, rvec, arpack::howmny::ritz_vectors, select.data(),
d.data(), z.data(), ldz, sigma, arpack::bmat::identity, N,
arpack::which::largest_magnitude, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam.data(), ipntr.data(), workd.data(),
ritz_option, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam, ipntr, workd.data(),
workl.data(), lworkl, info);
if (info < 0) throw std::runtime_error("Error in seupd, info " + std::to_string(info));
for (int i = 0; i < nev; ++i) {
float val = d[i];
float ref = (N - (nev - 1) + i);
float eps = std::fabs(val - ref);
Real val = d[i];
Real ref = (N - (nev - 1) + i);
Real eps = std::fabs(val - ref);
std::cout << "rank " << rank << " : " << val << " - " << ref << " - " << eps << std::endl;
/*eigen value order: smallest -> biggest*/
if (eps > 1.) {
throw std::domain_error("Correct eigenvalues not computed");
}
if (eps > tol_check) throw std::domain_error("Correct eigenvalues not computed");
}
std::cout << "------" << std::endl;
}
void diagonal_matrix_vector_product(std::complex<float> const* const x,
std::complex<float>* const y) {
template <typename Real>
void diagonal_matrix_vector_product(const std::complex<Real>* x, std::complex<Real>* y) {
for (int i = 0; i < 1000; ++i) {
y[i] = x[i] * std::complex<float>{i + 1.0f, i + 1.0f};
// Use complex matrix (i, -i) instead of (i, i): this way "largest_magnitude"
// and "largest_imaginary" options produce different results that can be checked.
y[i] = x[i] * std::complex<Real>{Real(i + 1), -Real(i + 1)};
}
}
void complex_symmetric_runner() {
a_int N = 1000;
a_int nev = 9;
a_int ncv = 2 * nev + 1;
a_int ldv = N;
a_int ldz = N + 1;
template <typename Real>
void complex_nonsymmetric_runner(double const& tol_check, arpack::which const& ritz_option) {
const a_int N = 1000;
const a_int nev = 9;
const a_int ncv = 2 * nev + 1;
const a_int ldv = N;
const a_int ldz = N;
const a_int lworkl = ncv * (3 * ncv + 5);
const a_int rvec = 0; // eigenvectors omitted
float tol = 0.000001; // small tol => more stable checks after EV computation.
a_int rvec = 0;
std::complex<float> sigma(0.0f, 0.0f);
const Real tol = 0.000001; // small tol => more stable checks after EV computation.
const std::complex<Real> sigma(0.0f, 0.0f); // not referenced in this mode
std::vector<std::complex<float>> resid(N);
std::vector<std::complex<float>> V(ncv * N);
std::vector<std::complex<float>> workd(3 * N);
std::vector<std::complex<float>> d(nev + 1);
std::vector<std::complex<float>> z((N + 1) * (nev + 1));
std::vector<a_int> select(ncv);
for (int i = 0; i < ncv; i++) select[i] = 1;
std::vector<std::complex<Real>> resid(N);
std::vector<std::complex<Real>> V(ldv * ncv);
std::vector<std::complex<Real>> z(ldz * nev);
std::vector<std::complex<Real>> d(nev);
std::vector<std::complex<Real>> workd(3 * N);
std::vector<std::complex<Real>> workl(lworkl);
std::vector<std::complex<Real>> workev(2 * ncv);
std::vector<Real> rwork(ncv);
std::vector<a_int> select(ncv); // since HOWMNY = 'A', only used as workspace here
a_int lworkl = 3 * (ncv * ncv) + 6 * ncv;
std::vector<std::complex<float>> workl(lworkl);
std::vector<float> rwork(ncv);
std::vector<std::complex<float>> workev(2 * ncv);
a_int info = 0;
a_int iparam[11], ipntr[14];
iparam[0] = 1; // ishift
iparam[2] = 10 * N; // on input: maxit; on output: actual iteration
iparam[3] = 1; // NB, only 1 allowed
iparam[6] = 1; // mode
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
std::array<a_int, 11> iparam;
iparam[0] = 1;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[4] = 0; // number of ev found by arpack.
iparam[6] = 1;
std::array<a_int, 14> ipntr;
MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
a_int ido = 0;
while (ido != 99) {
a_int info = 0, ido = 0;
do {
arpack::naupd(MCW, ido, arpack::bmat::identity, N,
arpack::which::largest_magnitude, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam.data(), ipntr.data(), workd.data(),
ritz_option, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam, ipntr, workd.data(),
workl.data(), lworkl, rwork.data(), info);
diagonal_matrix_vector_product(&(workd[ipntr[0] - 1]),
&(workd[ipntr[1] - 1]));
}
diagonal_matrix_vector_product(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
} while (ido == 1 || ido == -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
// check info and number of ev found by arpack
if (info < 0 || iparam[4] < nev) { /*arpack may succeed to compute more EV than expected*/
std::cout << "ERROR in naupd: iparam[4] " << iparam[4] << ", nev " << nev
<< ", info " << info << std::endl;
throw std::domain_error("Error inside ARPACK routines");
}
arpack::neupd(MCW, rvec, arpack::howmny::ritz_vectors, select.data(),
d.data(), z.data(), ldz, sigma, workev.data(),
arpack::bmat::identity, N, arpack::which::largest_magnitude,
nev, tol, resid.data(), ncv, V.data(), ldv, iparam.data(),
ipntr.data(), workd.data(), workl.data(), lworkl, rwork.data(),
info);
arpack::bmat::identity, N, ritz_option,
nev, tol, resid.data(), ncv, V.data(), ldv, iparam,
ipntr, workd.data(), workl.data(), lworkl, rwork.data(), info);
if (info < 0) throw std::runtime_error("Error in neupd, info " + std::to_string(info));
for (int i = 0; i < nev; ++i) {
float rval = std::real(d[i]);
float rref = (N-(nev-1)+i);
float reps = std::fabs(rval - rref);
float ival = std::imag(d[i]);
float iref = (N-(nev-1)+i);
float ieps = std::fabs(ival - iref);
std::cout << "rank " << rank << " : " << rval << " " << ival << " - " << rref << " " << iref << " - " << reps << " " << ieps << std::endl;
if (ritz_option == arpack::which::largest_magnitude) {
for (int i = 0; i < nev; ++i) {
Real rval = std::real(d[i]);
Real rref = static_cast<Real>(N - (nev - 1) + i);
Real reps = std::fabs(rval - rref);
Real ival = std::imag(d[i]);
Real iref = -static_cast<Real>(N - (nev - 1) + i);
Real ieps = std::fabs(ival - iref);
std::cout << rval << " " << ival << " - " << rref << " " << iref << " - " << reps << " " << ieps << std::endl;
/*eigen value order: smallest -> biggest*/
if (reps > 1. || ieps > 1.) {
throw std::domain_error("Correct eigenvalues not computed");
if (reps > tol_check || ieps > tol_check) throw std::domain_error("Correct eigenvalues not computed");
}
} else if (ritz_option == arpack::which::largest_imaginary) {
for (int i = 0; i < nev; ++i) {
Real rval = std::real(d[i]);
Real rref = static_cast<Real>(nev - i);
Real reps = std::fabs(rval - rref);
Real ival = std::imag(d[i]);
Real iref = -static_cast<Real>(nev - i);
Real ieps = std::fabs(ival - iref);
std::cout << rval << " " << ival << " - " << rref << " " << iref << " - " << reps << " " << ieps << std::endl;
if (reps > tol_check || ieps > tol_check) throw std::domain_error("Correct eigenvalues not computed");
}
} else {
throw std::domain_error("The input Ritz option is not allowed in this test file.");
}
std::cout << "------" << std::endl;
}
int main() {
@@ -199,9 +196,12 @@ int main() {
try {
// parpack without debug
real_symmetric_runner();
real_symmetric_runner<float>(1., arpack::which::largest_magnitude);
real_symmetric_runner<float>(1., arpack::which::largest_algebraic);
real_symmetric_runner<double>(1.e-05, arpack::which::largest_magnitude);
real_symmetric_runner<double>(1.e-05, arpack::which::largest_algebraic);
} catch (std::domain_error& e) {
std::cout << e.what() << '\n';
std::cout << e.what() << std::endl;
MPI_Abort(MPI_COMM_WORLD, 1);
}
@@ -231,7 +231,10 @@ int main() {
1);
try {
complex_symmetric_runner();
complex_nonsymmetric_runner<float>(1., arpack::which::largest_magnitude);
complex_nonsymmetric_runner<float>(1., arpack::which::largest_imaginary);
complex_nonsymmetric_runner<double>(1.e-05, arpack::which::largest_magnitude);
complex_nonsymmetric_runner<double>(1.e-05, arpack::which::largest_imaginary);
} catch (std::domain_error& e) {
std::cout << e.what() << '\n';
MPI_Abort(MPI_COMM_WORLD, 1);
TESTS/icb_arpack_c.c
+87 -88
View File
@@ -21,59 +21,63 @@
* with entries 1000, 999, ... , 2, 1 on the diagonal.
* */
void dMatVec(double* x, double* y) {
void dMatVec(const double* x, double* y) {
int i;
for (i = 0; i < 1000; ++i) y[i] = ((double)(i + 1)) * x[i];
};
int ds() {
a_int ido = 0;
char bmat[] = "I";
a_int N = 1000;
char which[] = "LM";
a_int nev = 9;
double tol = 0.000001; // small tol => more stable checks after EV computation.
const a_int N = 1000;
const a_int nev = 9;
const a_int ncv = 2 * nev + 1;
const a_int ldv = N;
const a_int ldz = N;
const a_int lworkl = ncv * (ncv + 8);
const a_int rvec = 1; // need eigenvectors
const double tol = 0.000001; // small tol => more stable checks after EV computation.
const double sigma = 0; // not referenced in this mode
double resid[N];
a_int ncv = 2 * nev + 1;
double V[ncv * N];
a_int ldv = N;
a_int iparam[11];
a_int ipntr[14];
double V[ldv * ncv];
double z[ldz * nev];
double d[nev];
double workd[3 * N];
a_int rvec = 1;
double workl[lworkl];
a_int select[ncv]; // since HOWMNY = 'A', only used as workspace here
a_int iparam[11], ipntr[11];
iparam[0] = 1; // ishift
iparam[2] = 10 * N; // on input: maxit; on output: actual iteration
iparam[3] = 1; // NB, only 1 allowed
iparam[6] = 1; // mode
char bmat[] = "I";
char which[] = "LM";
char howmny[] = "A";
double* d = (double*)malloc((nev + 1) * sizeof(double));
a_int select[ncv];
int i; // C99 compliant.
for (i = 0; i < ncv; i++) select[i] = 1;
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;
a_int info = 0;
iparam[0] = 1;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[4] = 0; // number of ev found by arpack.
iparam[6] = 1;
while (ido != 99) {
/* call arpack like you would have, but, use dsaupd_c instead of dsaupd_ */
a_int info = 0, ido = 0;
do {
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]));
}
if (iparam[4] != nev) {printf("Error: iparam[4] %d, nev %d\n", iparam[4], nev); return 1;} // check number of ev found by arpack.
} while (ido == 1 || ido == -1);
// check info and number of ev found by arpack.
if (info < 0 || iparam[4] < nev) {
printf("Error in saupd: iparam[4] %d, nev %d, info %d\n", iparam[4], nev, info);
return 1;
}
/* 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);
if (info < 0) {
printf("Error in seupd: info %d\n", info);
return 1;
}
int i; // C99 compliant.
for (i = 0; i < nev; ++i) {
double val = d[i];
double ref = (N-(nev-1)+i);
@@ -81,72 +85,70 @@ int ds() {
printf("%f - %f - %f\n", val, ref, eps);
/*eigen value order: smallest -> biggest*/
if(eps>1.e-05){
free(d);
return 1;
}
if (eps > 1.e-05) return 1;
}
free(d);
return 0;
}
void zMatVec(a_dcomplex* x, a_dcomplex* y) {
void zMatVec(const a_dcomplex* x, a_dcomplex* y) {
int i;
for (i = 0; i < 1000; ++i) y[i] = x[i] * CMPLX(i + 1.0, i + 1.0);
};
int zn() {
a_int ido = 0;
char bmat[] = "I";
a_int N = 1000;
char which[] = "LM";
a_int nev = 9;
double tol = 0.000001; // small tol => more stable checks after EV computation.
const a_int N = 1000;
const a_int nev = 9;
const a_int ncv = 2 * nev + 1;
const a_int ldv = N;
const a_int ldz = N;
const a_int lworkl = ncv * (3 * ncv + 5);
const a_int rvec = 0; // eigenvectors omitted
const double tol = 0.000001; // small tol => more stable checks after EV computation.
const a_dcomplex sigma = CMPLX(0., 0.); // not referenced in this mode
a_dcomplex resid[N];
a_int ncv = 2 * nev + 1;
a_dcomplex V[ncv * N];
a_int ldv = N;
a_int iparam[11];
a_int ipntr[14];
a_dcomplex V[ldv * ncv];
a_dcomplex z[ldz * nev];
a_dcomplex d[nev];
a_dcomplex workd[3 * N];
a_int rvec = 0;
char howmny[] = "A";
a_dcomplex* d =
(a_dcomplex*)malloc((nev + 1) * sizeof(a_dcomplex));
a_int select[ncv];
int i; // C99 compliant.
for (i = 0; i < ncv; i++) select[i] = 1;
a_dcomplex z[(N + 1) * (nev + 1)];
a_int ldz = N + 1;
a_dcomplex sigma = CMPLX(0., 0.);
int k;
for (k = 0; k < 3 * N; ++k) workd[k] = 0;
a_dcomplex 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];
a_dcomplex workl[lworkl];
a_dcomplex workev[2 * ncv];
a_int info = 0;
double rwork[ncv];
a_int select[ncv]; // since HOWMNY = 'A', only used as workspace here
iparam[0] = 1;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[4] = 0; // number of ev found by arpack.
iparam[6] = 1;
a_int iparam[11], ipntr[14];
iparam[0] = 1; // ishift
iparam[2] = 10 * N; // on input: maxit; on output: actual iteration
iparam[3] = 1; // NB, only 1 allowed
iparam[6] = 1; // mode
while (ido != 99) {
/* call arpack like you would have, but, use znaupd_c instead of znaupd_ */
char bmat[] = "I";
char which[] = "LM";
char howmny[] = "A";
a_int info = 0, ido = 0;
do {
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]));
}
if (iparam[4] != nev) {printf("Error: iparam[4] %d, nev %d\n", iparam[4], nev); return 1;} // check number of ev found by arpack.
} while (ido == 1 || ido == -1);
// check info and number of ev found by arpack
if (info < 0 || iparam[4] < nev) {
printf("Error in naupd: iparam[4] %d, nev %d, info %d\n", iparam[4], nev, info);
return 1;
}
/* 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);
tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl, rwork, &info);
if (info < 0) {
printf("Error in neupd: info %d\n", info);
return 1;
}
int i; // C99 compliant.
for (i = 0; i < nev; ++i) {
double rval = creal(d[i]);
double rref = (N-(nev-1)+i);
@@ -157,12 +159,9 @@ int zn() {
printf("%f %f - %f %f - %f %f\n", rval, ival, rref, iref, reps, ieps);
/*eigen value order: smallest -> biggest*/
if(reps>1.e-05 || ieps>1.e-05){
free(d);
return 1;
}
if (reps > 1.e-05 || ieps > 1.e-05) return 1;
}
free(d);
return 0;
}
TESTS/icb_arpack_cpp.cpp
+105 -109
View File
@@ -5,11 +5,10 @@
* of *[ae]upd_. The main advantage is that compiler checks the argument types
* and the correct function is called based on the type (float vs double vs
* complex). Note: to debug arpack, call debug_c. This is a test program to
* solve for the 9 largest eigenvalues of A*x = lambda*x where A is the diagonal
* solve for the 9 eigenvalues of A*x = lambda*x where A is the diagonal
* matrix with entries 1000, 999, ... , 2, 1 on the diagonal.
*/
#include <array>
#include <cmath>
#include <iostream>
#include <vector>
@@ -19,73 +18,62 @@
#include "stat_c.hpp" // arpack statistics.
template <typename Real>
void diagonal_matrix_vector_product(Real const* const x, Real* const y) {
void diagonal_matrix_vector_product(const Real* x, Real* y) {
for (int i = 0; i < 1000; ++i) {
y[i] = static_cast<Real>(i + 1) * x[i];
}
}
template <typename Real>
void real_symmetric_runner(double const& tol_check) {
a_int const N = 1000;
a_int const nev = 9;
void real_symmetric_runner(double const& tol_check, arpack::which const& ritz_option) {
const a_int N = 1000;
const a_int nev = 9;
const a_int ncv = 2 * nev + 1;
const a_int ldv = N;
const a_int ldz = N;
const a_int lworkl = ncv * (ncv + 8);
const a_int rvec = 1; // need eigenvectors
a_int const ncv = 2 * nev + 1;
a_int const ldv = N;
a_int const ldz = N + 1;
a_int const lworkl = 3 * (ncv * ncv) + 6 * ncv;
Real const tol = 0.000001; // small tol => more stable checks after EV computation.
Real const sigma = 0.0;
a_int const rvec = 1;
const Real tol = 0.000001; // small tol => more stable checks after EV computation.
const Real sigma = 0.0; // not referenced in this mode
std::vector<Real> resid(N);
std::vector<Real> V(ncv * N);
std::vector<Real> workd(3 * N, 0.0);
std::vector<Real> workl(lworkl, 0.0);
std::vector<Real> d((nev + 1));
std::vector<Real> z((N + 1) * (nev + 1));
std::vector<Real> V(ldv * ncv);
std::vector<Real> z(ldz * nev);
std::vector<Real> d(nev);
std::vector<Real> workd(3 * N);
std::vector<Real> workl(lworkl);
std::vector<a_int> select(ncv); // since HOWMNY = 'A', only used as workspace here
std::array<a_int, 11> iparam{};
iparam[0] = 1;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[4] = 0; // number of ev found by arpack.
iparam[6] = 1;
std::array<a_int, 14> ipntr{};
a_int iparam[11], ipntr[11];
iparam[0] = 1; // ishift
iparam[2] = 10 * N; // on input: maxit; on output: actual iteration
iparam[3] = 1; // NB, only 1 allowed
iparam[6] = 1; // mode
a_int info = 0, ido = 0;
while (ido != 99) {
do {
arpack::saupd(ido, arpack::bmat::identity, N,
arpack::which::largest_magnitude, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam.data(), ipntr.data(), workd.data(),
ritz_option, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam, ipntr, workd.data(),
workl.data(), lworkl, info);
diagonal_matrix_vector_product(&(workd[ipntr[0] - 1]),
&(workd[ipntr[1] - 1]));
}
diagonal_matrix_vector_product(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
} while (ido == 1 || ido == -1);
// check number of ev found by arpack.
if (iparam[4] < nev) { /*arpack may succeed to compute more EV than expected*/
std::cout << "ERROR: iparam[4] " << iparam[4] << ", nev " << nev
// check info and number of ev found by arpack.
if (info < 0 || iparam[4] < nev) { /*arpack may succeed to compute more EV than expected*/
std::cout << "ERROR in saupd: iparam[4] " << iparam[4] << ", nev " << nev
<< ", info " << info << std::endl;
throw std::domain_error("Error inside ARPACK routines");
}
std::vector<a_int> select(ncv);
for (int i = 0; i < ncv; i++) select[i] = 1;
arpack::seupd(rvec, arpack::howmny::ritz_vectors, select.data(), d.data(),
z.data(), ldz, sigma, arpack::bmat::identity, N,
arpack::which::largest_magnitude, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam.data(), ipntr.data(), workd.data(),
ritz_option, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam, ipntr, workd.data(),
workl.data(), lworkl, info);
if (info < 0) throw std::runtime_error("Error in seupd, info " + std::to_string(info));
for (int i = 0; i < nev; ++i) {
Real val = d[i];
@@ -94,106 +82,111 @@ void real_symmetric_runner(double const& tol_check) {
std::cout << val << " - " << ref << " - " << eps << std::endl;
/*eigen value order: smallest -> biggest*/
if (eps > tol_check) {
throw std::domain_error("Correct eigenvalues not computed");
}
if (eps > tol_check) throw std::domain_error("Correct eigenvalues not computed");
}
std::cout << "------\n";
std::cout << "------" << std::endl;
}
template <typename Real>
void diagonal_matrix_vector_product(std::complex<Real> const* const x,
std::complex<Real>* const y) {
void diagonal_matrix_vector_product(const std::complex<Real>* x, std::complex<Real>* y) {
for (int i = 0; i < 1000; ++i) {
y[i] = x[i] * std::complex<Real>{Real(i + 1), Real(i + 1)};
// Use complex matrix (i, -i) instead of (i, i): this way "largest_magnitude"
// and "largest_imaginary" options produce different results that can be checked.
y[i] = x[i] * std::complex<Real>{Real(i + 1), -Real(i + 1)};
}
}
template <typename Real>
void complex_symmetric_runner(double const& tol_check) {
a_int const N = 1000;
a_int const nev = 9;
void complex_nonsymmetric_runner(double const& tol_check, arpack::which const& ritz_option) {
const a_int N = 1000;
const a_int nev = 9;
const a_int ncv = 2 * nev + 1;
const a_int ldv = N;
const a_int ldz = N;
const a_int lworkl = ncv * (3 * ncv + 5);
const a_int rvec = 0; // eigenvectors omitted
a_int const ncv = 2 * nev + 1;
a_int const ldv = N;
a_int const ldz = N + 1;
a_int const lworkl = 3 * (ncv * ncv) + 6 * ncv;
Real const tol = 0.000001; // small tol => more stable checks after EV computation.
std::complex<Real> const sigma(0.0, 0.0);
a_int const rvec = 0;
const Real tol = 0.000001; // small tol => more stable checks after EV computation.
const std::complex<Real> sigma(0.0, 0.0); // not referenced in this mode
std::vector<std::complex<Real>> resid(N);
std::vector<std::complex<Real>> V(ncv * N);
std::vector<std::complex<Real>> V(ldv * ncv);
std::vector<std::complex<Real>> z(ldz * nev);
std::vector<std::complex<Real>> d(nev);
std::vector<std::complex<Real>> workd(3 * N);
std::vector<std::complex<Real>> workl(lworkl);
std::vector<std::complex<Real>> d(nev + 1);
std::vector<std::complex<Real>> z((N + 1) * (nev + 1));
std::vector<Real> rwork(ncv);
std::vector<std::complex<Real>> workev(2 * ncv);
std::vector<Real> rwork(ncv);
std::vector<a_int> select(ncv); // since HOWMNY = 'A', only used as workspace here
std::array<a_int, 11> iparam{};
iparam[0] = 1;
iparam[2] = 10 * N;
iparam[3] = 1;
iparam[4] = 0; // number of ev found by arpack.
iparam[6] = 1;
std::array<a_int, 14> ipntr{};
a_int iparam[11], ipntr[14];
iparam[0] = 1; // ishift
iparam[2] = 10 * N; // on input: maxit; on output: actual iteration
iparam[3] = 1; // NB, only 1 allowed
iparam[6] = 1; // mode
a_int info = 0, ido = 0;
while (ido != 99) {
do {
arpack::naupd(ido, arpack::bmat::identity, N,
arpack::which::largest_magnitude, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam.data(), ipntr.data(), workd.data(),
ritz_option, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam, ipntr, workd.data(),
workl.data(), lworkl, rwork.data(), info);
diagonal_matrix_vector_product(&(workd[ipntr[0] - 1]),
&(workd[ipntr[1] - 1]));
}
diagonal_matrix_vector_product(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
} while (ido == 1 || ido == -1);
// check number of ev found by arpack.
if (iparam[4] < nev) { /*arpack may succeed to compute more EV than expected*/
std::cout << "ERROR: iparam[4] " << iparam[4] << ", nev " << nev
// check info and number of ev found by arpack.
if (info < 0 || iparam[4] < nev) { /*arpack may succeed to compute more EV than expected*/
std::cout << "ERROR in naupd: iparam[4] " << iparam[4] << ", nev " << nev
<< ", info " << info << std::endl;
throw std::domain_error("Error inside ARPACK routines");
}
std::vector<a_int> select(ncv);
for (int i = 0; i < ncv; i++) select[i] = 1;
arpack::neupd(rvec, arpack::howmny::ritz_vectors, select.data(), d.data(),
z.data(), ldz, sigma, workev.data(), arpack::bmat::identity, N,
arpack::which::largest_magnitude, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam.data(), ipntr.data(), workd.data(),
ritz_option, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam, ipntr, workd.data(),
workl.data(), lworkl, rwork.data(), info);
if (info < 0) throw std::runtime_error("Error in neupd, info " + std::to_string(info));
for (int i = 0; i < nev; ++i) {
Real rval = std::real(d[i]);
Real rref = static_cast<Real>(N - (nev - 1) + i);
Real reps = std::fabs(rval - rref);
Real ival = std::imag(d[i]);
Real iref = static_cast<Real>(N - (nev - 1) + i);
Real ieps = std::fabs(ival - iref);
std::cout << rval << " " << ival << " - " << rref << " " << iref << " - " << reps << " " << ieps << std::endl;
if (ritz_option == arpack::which::largest_magnitude) {
for (int i = 0; i < nev; ++i) {
Real rval = std::real(d[i]);
Real rref = static_cast<Real>(N - (nev - 1) + i);
Real reps = std::fabs(rval - rref);
Real ival = std::imag(d[i]);
Real iref = -static_cast<Real>(N - (nev - 1) + i);
Real ieps = std::fabs(ival - iref);
std::cout << rval << " " << ival << " - " << rref << " " << iref << " - " << reps << " " << ieps << std::endl;
/*eigen value order: smallest -> biggest*/
if (reps > tol_check || ieps > tol_check) {
throw std::domain_error("Correct eigenvalues not computed");
if (reps > tol_check || ieps > tol_check) throw std::domain_error("Correct eigenvalues not computed");
}
} else if (ritz_option == arpack::which::largest_imaginary) {
for (int i = 0; i < nev; ++i) {
Real rval = std::real(d[i]);
Real rref = static_cast<Real>(nev - i);
Real reps = std::fabs(rval - rref);
Real ival = std::imag(d[i]);
Real iref = -static_cast<Real>(nev - i);
Real ieps = std::fabs(ival - iref);
std::cout << rval << " " << ival << " - " << rref << " " << iref << " - " << reps << " " << ieps << std::endl;
if (reps > tol_check || ieps > tol_check) throw std::domain_error("Correct eigenvalues not computed");
}
} else {
throw std::domain_error("The input Ritz option is not allowed in this test file.");
}
std::cout << "------" << std::endl;
}
int main() {
sstats_c();
// arpack without debug
real_symmetric_runner<float>(1.);
real_symmetric_runner<double>(1.e-05);
real_symmetric_runner<float>(1., arpack::which::largest_magnitude);
real_symmetric_runner<float>(1., arpack::which::largest_algebraic);
real_symmetric_runner<double>(1.e-05, arpack::which::largest_magnitude);
real_symmetric_runner<double>(1.e-05, arpack::which::largest_algebraic);
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;
@@ -215,8 +208,11 @@ int main() {
1);
// arpack with debug
complex_symmetric_runner<float>(1.);
complex_symmetric_runner<double>(1.e-05);
complex_nonsymmetric_runner<float>(1., arpack::which::largest_magnitude);
complex_nonsymmetric_runner<float>(1., arpack::which::largest_imaginary);
complex_nonsymmetric_runner<double>(1.e-05, arpack::which::largest_magnitude);
complex_nonsymmetric_runner<double>(1.e-05, arpack::which::largest_imaginary);
return 0;
}