Skip to content

Commit

Permalink
Fix icb_parpack_[c|cpp] to split the problem across MPI processes
Browse files Browse the repository at this point in the history
  • Loading branch information
@FabienPean
FabienPean committed Aug 2, 2023
1 parent 7cec434 commit d53f32a
Show file tree
Hide file tree
Showing 3 changed files with 69 additions and 33 deletions.
3 changes: 3 additions & 0 deletions CHANGES
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,8 @@
arpack-ng - 3.9.1

[ Fabien Péan ]
* [BUG FIX] Tests for PARPACK with C/C++ bindings icb_parpack_c and icb_parpack_cpp are now really parallel and split the problem across MPI processes.

[ Szabolcs Horvát ]
* [BUG FIX] Ensure that LAPACK RNG state is propagated (regression in 3.9.0).
* [BUG FIX] Ensure that separate random seeds are used on different parallel thread in D and S versions of functions (issue from original ARPACK).
Expand Down
48 changes: 32 additions & 16 deletions PARPACK/TESTS/MPI/icb_parpack_c.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,9 +22,9 @@
* with entries 1000, 999, ... , 2, 1 on the diagonal.
* */

void dMatVec(const double* x, double* y) {
void dMatVec(int start, int N, const double* x, double* y) {
int i;
for (i = 0; i < 1000; ++i) y[i] = ((double)(i + 1)) * x[i];
for (i = 0; i < N; ++i) y[i] = ((double)(start + i + 1)) * x[i];
};

int ds() {
Expand Down Expand Up @@ -57,16 +57,24 @@ int ds() {
char which[] = "LM";
char howmny[] = "A";

int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);

/// Split problem across each process/////////////////////
int rank, nprocs;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
a_int N_local = N / nprocs;
if (rank < N % nprocs) // spread the remaining on each process
N_local = N_local + 1;
printf("rank: %d, size: %d\n", rank, N_local);
int start = rank*(N / nprocs) + ((rank < N % nprocs) ? rank : N % nprocs);
printf("rank: %d, start: %d\n", rank, start);
//////////////////////////////////////////////////////////
a_int info = 0, ido = 0;
do {
pdsaupd_c(MCW, &ido, bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam,
pdsaupd_c(MCW, &ido, bmat, N_local, which, nev, tol, resid, ncv, V, ldv, iparam,
ipntr, workd, workl, lworkl, &info);

dMatVec(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
dMatVec(start, N_local, &(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
} while (ido == 1 || ido == -1);

// check info and number of ev found by arpack.
Expand All @@ -75,7 +83,7 @@ int ds() {
return 1;
}

pdseupd_c(MCW, rvec, howmny, select, d, z, ldz, sigma, bmat, N, which, nev,
pdseupd_c(MCW, rvec, howmny, select, d, z, ldz, sigma, bmat, N_local, which, nev,
tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl, &info);
if (info < 0) {
printf("Error in seupd: info %d\n", info);
Expand All @@ -94,9 +102,9 @@ int ds() {
return 0;
}

void zMatVec(const a_dcomplex* x, a_dcomplex* y) {
void zMatVec(int start, int N, 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);
for (i = 0; i < N; ++i) y[i] = x[i] * CMPLXF(start + i + 1.0, start + i + 1.0);
};

int zn() {
Expand Down Expand Up @@ -131,16 +139,24 @@ int zn() {
char which[] = "LM";
char howmny[] = "A";

int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);

/// Split problem across each process/////////////////////
int rank, nprocs;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
a_int N_local = N / nprocs;
if (rank < N % nprocs) // spread the remaining on each process
N_local = N_local + 1;
printf("rank: %d, size: %d\n", rank, N_local);
int start = rank*(N / nprocs) + ((rank < N % nprocs) ? rank : N % nprocs);
printf("rank: %d, start: %d\n", rank, start);
//////////////////////////////////////////////////////////
a_int info = 0, ido = 0;
do {
pznaupd_c(MCW, &ido, bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam,
pznaupd_c(MCW, &ido, bmat, N_local, which, nev, tol, resid, ncv, V, ldv, iparam,
ipntr, workd, workl, lworkl, rwork, &info);

zMatVec(&(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
zMatVec(start, N_local, &(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
} while (ido == 1 || ido == -1);

// check info and number of ev found by arpack.
Expand All @@ -149,7 +165,7 @@ int zn() {
return 1;
}

pzneupd_c(MCW, rvec, howmny, select, d, z, ldz, sigma, workev, bmat, N, which,
pzneupd_c(MCW, rvec, howmny, select, d, z, ldz, sigma, workev, bmat, N_local, which,
nev, tol, resid, ncv, V, ldv, iparam, ipntr, workd, workl, lworkl,
rwork, &info);
if (info < 0) {
Expand Down
51 changes: 34 additions & 17 deletions PARPACK/TESTS/MPI/icb_parpack_cpp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,9 +20,9 @@
#include "stat_c.hpp" // arpack statistics.

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 diagonal_matrix_vector_product(int start, int nloc, const Real* x, Real* y) {
for (int i = 0; i < nloc; ++i) {
y[i] = static_cast<Real>(start + i + 1) * x[i];
}
}

Expand Down Expand Up @@ -53,18 +53,27 @@ void real_symmetric_runner(double const& tol_check, arpack::which const& ritz_op
iparam[3] = 1; // NB, only 1 allowed
iparam[6] = 1; // mode

int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
/// Split problem across each process/////////////////////
int rank, nprocs;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
a_int N_local = N / nprocs;
if (rank < N % nprocs) // spread the remaining on each process
N_local = N_local + 1;
std::cout<< "rank: " << rank << ", size :" << N_local << "\n";
int start = rank*(N / nprocs) + std::min(rank, N % nprocs);
std::cout<< "rank: " <<rank<< ", start:" << start<<"\n";
//////////////////////////////////////////////////////////

a_int info = 0, ido = 0;
do {
arpack::saupd(MCW, ido, arpack::bmat::identity, N,
arpack::saupd(MCW, ido, arpack::bmat::identity, N_local,
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(start, N_local, &(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
} while (ido == 1 || ido == -1);

// check info and number of ev found by arpack.
Expand All @@ -75,7 +84,7 @@ void real_symmetric_runner(double const& tol_check, arpack::which const& ritz_op
}

arpack::seupd(MCW, rvec, arpack::howmny::ritz_vectors, select.data(),
d.data(), z.data(), ldz, sigma, arpack::bmat::identity, N,
d.data(), z.data(), ldz, sigma, arpack::bmat::identity, N_local,
ritz_option, nev, tol, resid.data(), ncv,
V.data(), ldv, iparam, ipntr, workd.data(),
workl.data(), lworkl, info);
Expand All @@ -94,11 +103,11 @@ void real_symmetric_runner(double const& tol_check, arpack::which const& ritz_op
}

template <typename Real>
void diagonal_matrix_vector_product(const std::complex<Real>* x, std::complex<Real>* y) {
for (int i = 0; i < 1000; ++i) {
void diagonal_matrix_vector_product(int start, int N_local, const std::complex<Real>* x, std::complex<Real>* y) {
for (int i = 0; i < N_local; ++i) {
// 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)};
y[i] = x[i] * std::complex<Real>{Real(start + i + 1), -Real(start + i + 1)};
}
}

Expand Down Expand Up @@ -131,18 +140,26 @@ void complex_nonsymmetric_runner(double const& tol_check, arpack::which const& r
iparam[3] = 1; // NB, only 1 allowed
iparam[6] = 1; // mode

int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);

/// Split problem across each process/////////////////////
int rank, nprocs;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
a_int N_local = N / nprocs;
if (rank < N % nprocs) // spread the remaining on each process
N_local = N_local + 1;
std::cout<< "rank: "<<rank<<", size :"<<N_local<<"\n";
int start = rank * (N / nprocs) + std::min(rank, N % nprocs);
std::cout<< "rank: " << rank << ", start:" << start << "\n";
//////////////////////////////////////////////////////////
a_int info = 0, ido = 0;
do {
arpack::naupd(MCW, ido, arpack::bmat::identity, N,
arpack::naupd(MCW, ido, arpack::bmat::identity, N_local,
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(start, N_local, &(workd[ipntr[0] - 1]), &(workd[ipntr[1] - 1]));
} while (ido == 1 || ido == -1);

// check info and number of ev found by arpack
Expand All @@ -154,7 +171,7 @@ void complex_nonsymmetric_runner(double const& tol_check, arpack::which const& r

arpack::neupd(MCW, rvec, arpack::howmny::ritz_vectors, select.data(),
d.data(), z.data(), ldz, sigma, workev.data(),
arpack::bmat::identity, N, ritz_option,
arpack::bmat::identity, N_local, 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));
Expand Down

0 comments on commit d53f32a

Please sign in to comment.