Skip to content

Commit 553b4bd

Browse files
prj-prj-
committed
Support for block solves with .asarray
1 parent d9046be commit 553b4bd

File tree

3 files changed

+32
-19
lines changed

3 files changed

+32
-19
lines changed

‎3rdparty/ff-petsc/Makefile‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -73,7 +73,7 @@ ifeq ($(FF_generic_petsc), yes)
7373
FLAGS_MTUNE := -mtune=generic
7474
endif
7575
unexport MAKEFLAGS
76-
COMMON_FLAGS := --with-debugging=0 COPTFLAGS='-O3 $(FLAGS_MTUNE)' CXXOPTFLAGS='-O3 $(FLAGS_MTUNE)' FOPTFLAGS='-O3 $(FLAGS_MTUNE)' --with-ssl=0 --with-x=0 --with-fortran-bindings=0 --with-hg=0 --with-c2html=0 --with-coverage=0 --with-syclc=0 --with-hipc=0 --with-cudac=0 --with-tau-perfstubs=0 \
76+
COMMON_FLAGS := --with-debugging=0 COPTFLAGS='-O3 $(FLAGS_MTUNE)' CXXOPTFLAGS='-O3 $(FLAGS_MTUNE)' FOPTFLAGS='-O3 $(FLAGS_MTUNE)' --with-ssl=0 --with-x=0 --with-fortran-bindings=0 --with-hg=0 --with-c2html=0 --with-coverage=0 --with-syclc=0 --with-hipc=0 --with-cudac=0 --with-tau-perfstubs=0 --with-python=0 \
7777
--download-slepc --download-hpddm
7878

7979
ifneq ($(wildcard ~/.petsc_pkg),)

‎examples/hpddm/helmholtz-2d-PETSc-complex.edp‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -65,3 +65,11 @@ for(int i = 0; i < M; ++i) {
6565
macro params()cmm = "Global solution #" + i, wait = 1, fill = 1, value = 1, dim = 3//
6666
plotMPI(Th, u, Pk, def, complex, params);
6767
}
68+
complex[int] C(B.n * B.m);
69+
C = A^-1 * RHS.asarray;
70+
for(int i = 0; i < M; ++i) {
71+
u[] = B(:, i);
72+
u[] -= C(i * Wh.ndof:(i + 1) * Wh.ndof - 1);
73+
macro params()cmm = "Error #" + i, wait = 1, fill = 1, value = 1, dim = 3//
74+
plotMPI(Th, u, Pk, def, complex, params);
75+
}

‎plugin/mpi/PETSc-code.hpp‎

Lines changed: 23 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -5031,14 +5031,19 @@ namespace PETSc {
50315031
double timing = MPI_Wtime( );
50325032
MatCreateVecs((*t)._petsc, &x, &y);
50335033
ffassert(out->N() == u->N() && out->M() == u->M());
5034+
int mu = 1;
5035+
if ((*t)._A) {
5036+
ffassert(out->n % (*t)._A->getDof() == 0);
5037+
mu = out->n / (*t)._A->getDof();
5038+
} else ffassert(out->M() == 1);
50345039
PetscScalar* ptr;
50355040
MatType type;
50365041
PetscBool isType;
50375042
MatGetType((*t)._petsc, &type);
50385043
PetscStrcmp(type, MATNEST, &isType);
50395044
PetscScalar* p = reinterpret_cast<PetscScalar*>(u->operator upscaled_type<PetscScalar>*());
50405045
if ((*t)._vector_global) {
5041-
ffassert(out->M() == 1);
5046+
ffassert(mu == 1);
50425047
for(int i = 0; i < u->n; ++i)
50435048
p[i] = u->operator[](i);
50445049
MPI_Allreduce(MPI_IN_PLACE, p, u->n, HPDDM::Wrapper<K>::mpi_type(), MPI_SUM, PetscObjectComm((PetscObject)(*t)._petsc));
@@ -5084,13 +5089,13 @@ namespace PETSc {
50845089
else if (std::is_same< typename std::remove_reference< decltype(*t.A->_A) >::type,
50855090
HpSchwarz< PetscScalar > >::value) {
50865091
#if !defined(PETSC_USE_REAL_DOUBLE)
5087-
ffassert(u->M() == 1);
5092+
ffassert(mu == 1);
50885093
for(int i = 0; i < u->N(); ++i)
50895094
p[i] = u->operator[](i);
50905095
#else
50915096
static_assert(std::is_same<PetscReal, upscaled_type<PetscReal>>::value, "Wrong types");
50925097
#endif
5093-
if (out->M() == 1) {
5098+
if (mu == 1) {
50945099
VecGetArray(x, &ptr);
50955100
if(isType) {
50965101
ffassert((std::is_same<PetscReal, upscaled_type<PetscReal>>::value));
@@ -5107,23 +5112,23 @@ namespace PETSc {
51075112
PetscInt m, M;
51085113
MatGetLocalSize((*t)._petsc, &m, nullptr);
51095114
MatGetSize((*t)._petsc, &M, nullptr);
5110-
MatCreateDense(PetscObjectComm((PetscObject)(*t)._petsc), m, PETSC_DECIDE, M, out->M(), nullptr, &X);
5111-
MatCreateDense(PetscObjectComm((PetscObject)(*t)._petsc), m, PETSC_DECIDE, M, out->M(), nullptr, &Y);
5115+
MatCreateDense(PetscObjectComm((PetscObject)(*t)._petsc), m, PETSC_DECIDE, M, mu, nullptr, &X);
5116+
MatCreateDense(PetscObjectComm((PetscObject)(*t)._petsc), m, PETSC_DECIDE, M, mu, nullptr, &Y);
51125117
MatDenseGetArrayWrite(X, &ptr);
5113-
for(int j = 0; j < out->M(); ++j) {
5114-
PetscScalar* in = p + j * u->N();
5118+
for(int j = 0; j < mu; ++j) {
5119+
PetscScalar* in = p + j * (*t)._A->getDof();
51155120
PetscScalar* out = ptr + j * m;
51165121
HPDDM::Subdomain< K >::template distributedVec< 0 >((*t)._num, (*t)._first, (*t)._last,
51175122
in, out,
5118-
static_cast<PetscInt>(u->N()), 1);
5123+
static_cast<PetscInt>((*t)._A->getDof()), 1);
51195124
}
51205125
MatDenseRestoreArrayWrite(X, &ptr);
51215126
}
51225127
if ((*t)._ksp) {
51235128
PetscBool nonZero;
51245129
KSPGetInitialGuessNonzero((*t)._ksp, &nonZero);
51255130
if (nonZero) {
5126-
ffassert(out->M() == 1);
5131+
ffassert(mu == 1);
51275132
VecGetArray(y, &ptr);
51285133
p = reinterpret_cast<PetscScalar*>(out->operator upscaled_type<PetscScalar>*());
51295134
#if !defined(PETSC_USE_REAL_DOUBLE)
@@ -5141,7 +5146,7 @@ namespace PETSc {
51415146
}
51425147
*out = K();
51435148
} else {
5144-
ffassert(out->M() == 1);
5149+
ffassert(mu == 1);
51455150
VecSet(x, PetscScalar( ));
51465151
Vec isVec;
51475152
VecCreateMPIWithArray(PETSC_COMM_SELF, 1, (*t)._A->getMatrix( )->HPDDM_n,
@@ -5171,11 +5176,11 @@ namespace PETSc {
51715176
}
51725177
}
51735178
if (trans == 'N') {
5174-
if (out->M() == 1) KSPSolve((*t)._ksp, x, y);
5179+
if (mu == 1) KSPSolve((*t)._ksp, x, y);
51755180
else KSPMatSolve((*t)._ksp, X, Y);
51765181
}
51775182
else {
5178-
if (out->M() == 1) {
5183+
if (mu == 1) {
51795184
if (!std::is_same< PetscScalar, PetscReal >::value && t.conjugate) VecConjugate(x);
51805185
KSPSolveTranspose((*t)._ksp, x, y);
51815186
if (!std::is_same< PetscScalar, PetscReal >::value && t.conjugate) VecConjugate(y);
@@ -5233,7 +5238,7 @@ namespace PETSc {
52335238
}
52345239
else if (std::is_same< typename std::remove_reference< decltype(*t.A->_A) >::type,
52355240
HpSchwarz< PetscScalar > >::value) {
5236-
if (out->M() == 1) {
5241+
if (mu == 1) {
52375242
VecGetArray(y, &ptr);
52385243
if(isType)
52395244
loopDistributedVec<1, 'N'>((*t)._petsc, (*t)._exchange, out, ptr);
@@ -5245,11 +5250,11 @@ namespace PETSc {
52455250
PetscInt m;
52465251
MatGetLocalSize((*t)._petsc, &m, nullptr);
52475252
MatDenseGetArrayWrite(Y, &ptr);
5248-
for(int j = 0; j < out->M(); ++j) {
5249-
PetscScalar* in = p + j * u->N();
5253+
for(int j = 0; j < mu; ++j) {
5254+
PetscScalar* in = p + j * (*t)._A->getDof();
52505255
PetscScalar* out = ptr + j * m;
52515256
HPDDM::Subdomain< K >::template distributedVec< 1 >((*t)._num, (*t)._first, (*t)._last, in,
5252-
out, static_cast<PetscInt>(u->N()), 1);
5257+
out, static_cast<PetscInt>((*t)._A->getDof()), 1);
52535258
}
52545259
MatDenseRestoreArrayWrite(Y, &ptr);
52555260
MatDestroy(&X);
@@ -5268,9 +5273,9 @@ namespace PETSc {
52685273
if (std::is_same< typename std::remove_reference< decltype(*t.A->_A) >::type,
52695274
HpSchwarz< PetscScalar > >::value &&
52705275
t.A->_A)
5271-
(*t)._A->exchange(p, out->M());
5276+
(*t)._A->exchange(p, mu);
52725277
if(!std::is_same<PetscReal, upscaled_type<PetscReal>>::value) {
5273-
ffassert(out->M() == 1);
5278+
ffassert(mu == 1);
52745279
for(int i = out->N() - 1; i >= 0; --i)
52755280
out->operator[](i) = p[i];
52765281
p = reinterpret_cast<PetscScalar*>(u->operator upscaled_type<PetscScalar>*());

0 commit comments

Comments
 (0)