Skip to content

Commit a5f4584

Browse files
prj-prj-
committed
Add support for KNM solves with Mat
1 parent 08464d6 commit a5f4584

File tree

2 files changed

+114
-28
lines changed

2 files changed

+114
-28
lines changed

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

Lines changed: 19 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,8 @@ func real wedge(real a, real b) {
2424
return 3;
2525
}
2626
real omega = 2 * pi * 5;
27-
func f = 80 * 100 * s * exp(-20 * 100 * s * ((x-0.5)^2 + (y-0.25)^2));
27+
real xi = 0, yi = 0;
28+
func f = 80 * 100 * s * exp(-20 * 100 * s * ((x-0.5-xi)^2 + (y-0.25-yi)^2));
2829
complex[int] rhs(Wh.ndof); // local right-hand side
2930
matrix<complex> Loc; // local operator
3031
{ // local weak form
@@ -47,3 +48,20 @@ u[] = A^-1 * rhs;
4748

4849
macro def(u)u//
4950
plotMPI(Th, u, Pk, def, complex, cmm = "Global solution");
51+
52+
int M = 4;
53+
complex[int, int] RHS(rhs.n, M);
54+
RHS(:, 0) = rhs;
55+
for(int i = 1; i < M; ++i) {
56+
xi += 0.1;
57+
yi += 0.1;
58+
varf vPb(u, v) = int2d(Th)(f * v) + on(1, u = 0.0);
59+
RHS(:, i) = vPb(0, Wh, tgv = -1);
60+
}
61+
set(A, sparams = "-ksp_type hpddm -ksp_converged_reason -ksp_hpddm_type bgmres");
62+
complex[int, int] B = A^-1 * RHS;
63+
for(int i = 0; i < M; ++i) {
64+
u[] = B(:, i);
65+
macro params()cmm = "Global solution #" + i, wait = 1, fill = 1, value = 1, dim = 3//
66+
plotMPI(Th, u, Pk, def, complex, params);
67+
}

‎plugin/mpi/PETSc-code.hpp‎

Lines changed: 95 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -2239,6 +2239,16 @@ namespace PETSc {
22392239
!std::is_same< T, KN< PetscScalar > >::value && !std::is_same< T, KNM< PetscScalar > >::value &&
22402240
!std::is_same< T, KN< long > >::value && !std::is_same< T, KNM< long > >::value>::type* = nullptr >
22412241
void resize(T* v, int n, int m) {}
2242+
template< class T, typename std::enable_if< std::is_same< T, KN< PetscScalar > >::value || std::is_same< T, KN< long > >::value >::type* =
2243+
nullptr >
2244+
void init_KN_or_KNM(T* v, T* w) {
2245+
v->init(w->n);
2246+
}
2247+
template< class T, typename std::enable_if< std::is_same< T, KNM< PetscScalar > >::value || std::is_same< T, KNM< long > >::value >::type* =
2248+
nullptr >
2249+
void init_KN_or_KNM(T* v, T* w) {
2250+
v->init(w->N(), w->M());
2251+
}
22422252
template< class T, class U >
22432253
void changeNumbering_func(PetscInt* const num, PetscInt first, PetscInt last,
22442254
PetscInt m, PetscInt n, PetscInt bs, T* ptIn, U* ptOut, bool inverse) {
@@ -5015,6 +5025,7 @@ namespace PETSc {
50155025
InvPETSc(T v, U w) : t(v), u(w) {}
50165026
void solve(U out) const {
50175027
if ((*t)._petsc) {
5028+
Mat X, Y;
50185029
Vec x, y;
50195030
double timing = MPI_Wtime( );
50205031
MatCreateVecs((*t)._petsc, &x, &y);
@@ -5026,6 +5037,7 @@ namespace PETSc {
50265037
PetscStrcmp(type, MATNEST, &isType);
50275038
PetscScalar* p = reinterpret_cast<PetscScalar*>(u->operator upscaled_type<PetscScalar>*());
50285039
if ((*t)._vector_global) {
5040+
ffassert(out->M() == 1);
50295041
for(int i = 0; i < u->n; ++i)
50305042
p[i] = u->operator[](i);
50315043
MPI_Allreduce(MPI_IN_PLACE, p, u->n, HPDDM::Wrapper<K>::mpi_type(), MPI_SUM, PetscObjectComm((PetscObject)(*t)._petsc));
@@ -5070,26 +5082,53 @@ namespace PETSc {
50705082
}
50715083
else if (std::is_same< typename std::remove_reference< decltype(*t.A->_A) >::type,
50725084
HpSchwarz< PetscScalar > >::value) {
5073-
VecGetArray(x, &ptr);
5074-
for(int i = 0; i < u->n; ++i)
5075-
p[i] = u->operator[](i);
5076-
if(isType) {
5077-
ffassert((std::is_same<PetscReal, upscaled_type<PetscReal>>::value));
5078-
loopDistributedVec<0, 'N'>((*t)._petsc, (*t)._exchange, u, ptr);
5085+
#if !defined(PETSC_USE_REAL_DOUBLE)
5086+
for(int j = 0; j < out->M(); ++j)
5087+
for(int i = 0; i < u->N(); ++i)
5088+
p[i + j * u->N()] = u->operator[](i, j);
5089+
#else
5090+
static_assert(std::is_same<PetscReal, upscaled_type<PetscReal>>::value, "Wrong types");
5091+
#endif
5092+
if (out->M() == 1) {
5093+
VecGetArray(x, &ptr);
5094+
if(isType) {
5095+
ffassert((std::is_same<PetscReal, upscaled_type<PetscReal>>::value));
5096+
loopDistributedVec<0, 'N'>((*t)._petsc, (*t)._exchange, u, ptr);
5097+
}
5098+
else
5099+
HPDDM::Subdomain< K >::template distributedVec< 0 >((*t)._num, (*t)._first, (*t)._last,
5100+
p, ptr,
5101+
static_cast<PetscInt>(u->N()), 1);
5102+
VecRestoreArray(x, &ptr);
5103+
}
5104+
else {
5105+
ffassert(!isType);
5106+
PetscInt m, M;
5107+
MatGetLocalSize((*t)._petsc, &m, nullptr);
5108+
MatGetSize((*t)._petsc, &M, nullptr);
5109+
MatCreateDense(PetscObjectComm((PetscObject)(*t)._petsc), m, PETSC_DECIDE, M, out->M(), nullptr, &X);
5110+
MatCreateDense(PetscObjectComm((PetscObject)(*t)._petsc), m, PETSC_DECIDE, M, out->M(), nullptr, &Y);
5111+
MatDenseGetArrayWrite(X, &ptr);
5112+
for(int j = 0; j < out->M(); ++j) {
5113+
PetscScalar* in = p + j * u->N();
5114+
PetscScalar* out = ptr + j * m;
5115+
HPDDM::Subdomain< K >::template distributedVec< 0 >((*t)._num, (*t)._first, (*t)._last,
5116+
in, out,
5117+
static_cast<PetscInt>(u->N()), 1);
5118+
}
5119+
MatDenseRestoreArrayWrite(X, &ptr);
50795120
}
5080-
else
5081-
HPDDM::Subdomain< K >::template distributedVec< 0 >((*t)._num, (*t)._first, (*t)._last,
5082-
p, ptr,
5083-
static_cast<PetscInt>(u->n), 1);
5084-
VecRestoreArray(x, &ptr);
50855121
if ((*t)._ksp) {
50865122
PetscBool nonZero;
50875123
KSPGetInitialGuessNonzero((*t)._ksp, &nonZero);
50885124
if (nonZero) {
5125+
ffassert(out->M() == 1);
50895126
VecGetArray(y, &ptr);
50905127
p = reinterpret_cast<PetscScalar*>(out->operator upscaled_type<PetscScalar>*());
5128+
#if !defined(PETSC_USE_REAL_DOUBLE)
50915129
for(int i = 0; i < out->n; ++i)
50925130
p[i] = out->operator[](i);
5131+
#endif
50935132
if(isType)
50945133
loopDistributedVec<0, 'N'>((*t)._petsc, (*t)._exchange, out, ptr);
50955134
else
@@ -5101,6 +5140,7 @@ namespace PETSc {
51015140
}
51025141
*out = K();
51035142
} else {
5143+
ffassert(out->M() == 1);
51045144
VecSet(x, PetscScalar( ));
51055145
Vec isVec;
51065146
VecCreateMPIWithArray(PETSC_COMM_SELF, 1, (*t)._A->getMatrix( )->HPDDM_n,
@@ -5129,12 +5169,20 @@ namespace PETSc {
51295169
isType = PETSC_TRUE;
51305170
}
51315171
}
5132-
if (trans == 'N')
5133-
KSPSolve((*t)._ksp, x, y);
5172+
if (trans == 'N') {
5173+
if (out->M() == 1) KSPSolve((*t)._ksp, x, y);
5174+
else KSPMatSolve((*t)._ksp, X, Y);
5175+
}
51345176
else {
5135-
if (!std::is_same< PetscScalar, PetscReal >::value && t.conjugate) VecConjugate(x);
5136-
KSPSolveTranspose((*t)._ksp, x, y);
5137-
if (!std::is_same< PetscScalar, PetscReal >::value && t.conjugate) VecConjugate(y);
5177+
if (out->M() == 1) {
5178+
if (!std::is_same< PetscScalar, PetscReal >::value && t.conjugate) VecConjugate(x);
5179+
KSPSolveTranspose((*t)._ksp, x, y);
5180+
if (!std::is_same< PetscScalar, PetscReal >::value && t.conjugate) VecConjugate(y);
5181+
} else {
5182+
if (!std::is_same< PetscScalar, PetscReal >::value && t.conjugate) MatConjugate(X);
5183+
KSPMatSolveTranspose((*t)._ksp, X, Y);
5184+
if (!std::is_same< PetscScalar, PetscReal >::value && t.conjugate) MatConjugate(Y);
5185+
}
51385186
}
51395187
if (verbosity > 0 && mpirank == 0)
51405188
cout << " --- system solved with PETSc (in " << MPI_Wtime( ) - timing << ")" << endl;
@@ -5184,13 +5232,28 @@ namespace PETSc {
51845232
}
51855233
else if (std::is_same< typename std::remove_reference< decltype(*t.A->_A) >::type,
51865234
HpSchwarz< PetscScalar > >::value) {
5187-
VecGetArray(y, &ptr);
5188-
if(isType)
5189-
loopDistributedVec<1, 'N'>((*t)._petsc, (*t)._exchange, out, ptr);
5190-
else
5191-
HPDDM::Subdomain< K >::template distributedVec< 1 >((*t)._num, (*t)._first, (*t)._last, p,
5192-
ptr, static_cast<PetscInt>(out->n), 1);
5193-
VecRestoreArray(y, &ptr);
5235+
if (out->M() == 1) {
5236+
VecGetArray(y, &ptr);
5237+
if(isType)
5238+
loopDistributedVec<1, 'N'>((*t)._petsc, (*t)._exchange, out, ptr);
5239+
else
5240+
HPDDM::Subdomain< K >::template distributedVec< 1 >((*t)._num, (*t)._first, (*t)._last, p,
5241+
ptr, static_cast<PetscInt>(out->N()), 1);
5242+
VecRestoreArray(y, &ptr);
5243+
} else {
5244+
PetscInt m;
5245+
MatGetLocalSize((*t)._petsc, &m, nullptr);
5246+
MatDenseGetArrayWrite(Y, &ptr);
5247+
for(int j = 0; j < out->M(); ++j) {
5248+
PetscScalar* in = p + j * u->N();
5249+
PetscScalar* out = ptr + j * m;
5250+
HPDDM::Subdomain< K >::template distributedVec< 1 >((*t)._num, (*t)._first, (*t)._last, in,
5251+
out, static_cast<PetscInt>(u->N()), 1);
5252+
}
5253+
MatDenseRestoreArrayWrite(Y, &ptr);
5254+
MatDestroy(&X);
5255+
MatDestroy(&Y);
5256+
}
51945257
} else {
51955258
Vec isVec;
51965259
VecCreateMPIWithArray(PETSC_COMM_SELF, 1, (*t)._A->getMatrix( )->HPDDM_n,
@@ -5204,12 +5267,13 @@ namespace PETSc {
52045267
if (std::is_same< typename std::remove_reference< decltype(*t.A->_A) >::type,
52055268
HpSchwarz< PetscScalar > >::value &&
52065269
t.A->_A)
5207-
(*t)._A->exchange(p);
5270+
(*t)._A->exchange(p, out->M());
52085271
if(!std::is_same<PetscReal, upscaled_type<PetscReal>>::value) {
5209-
for(int i = out->n - 1; i >= 0; --i)
5272+
ffassert(out->M() == 1);
5273+
for(int i = out->N() - 1; i >= 0; --i)
52105274
out->operator[](i) = p[i];
52115275
p = reinterpret_cast<PetscScalar*>(u->operator upscaled_type<PetscScalar>*());
5212-
for(int i = u->n - 1; i >= 0; --i)
5276+
for(int i = u->N() - 1; i >= 0; --i)
52135277
u->operator[](i) = p[i];
52145278
}
52155279
}
@@ -5222,7 +5286,7 @@ namespace PETSc {
52225286
PetscInt n, m;
52235287
MatGetSize(A.t.A->_petsc, &n, &m);
52245288
ffassert(n == m);
5225-
Ax->init(A.u->n);
5289+
init_KN_or_KNM(Ax, A.u);
52265290
return inv(Ax, A);
52275291
}
52285292
};
@@ -6717,6 +6781,10 @@ static void Init_PETSc( ) {
67176781
addProd< Dmat, PETSc::ProdPETSc, KN< upscaled_type<PetscScalar> >, PetscScalar, 'T' >( );
67186782
addInv< Dmat, PETSc::InvPETSc, KN< upscaled_type<PetscScalar> >, PetscScalar, 'N' >( );
67196783
addInv< Dmat, PETSc::InvPETSc, KN< upscaled_type<PetscScalar> >, PetscScalar, 'T' >( );
6784+
Dcl_Type<PETSc::InvPETSc<pwr<Dmat>, KNM< upscaled_type<PetscScalar> >*, PetscScalar, 'N'>>();
6785+
TheOperators->Add("*", new OneBinaryOperator_st<assign<PETSc::InvPETSc<pwr<Dmat>, KNM< upscaled_type<PetscScalar> >*, PetscScalar, 'N'>, pwr<Dmat>, KNM< upscaled_type<PetscScalar> >*>>);
6786+
TheOperators->Add("=", new OneOperator2<KNM< upscaled_type<PetscScalar> >*, KNM< upscaled_type<PetscScalar> >*, PETSc::InvPETSc<pwr<Dmat>, KNM< upscaled_type<PetscScalar> >*, PetscScalar, 'N'>>(PETSc::InvPETSc<pwr<Dmat>, KNM< upscaled_type<PetscScalar> >*, PetscScalar, 'N'>::inv));
6787+
TheOperators->Add("<-", new OneOperator2<KNM< upscaled_type<PetscScalar> >*, KNM< upscaled_type<PetscScalar> >*, PETSc::InvPETSc<pwr<Dmat>, KNM< upscaled_type<PetscScalar> >*, PetscScalar, 'N'>>(PETSc::InvPETSc<pwr<Dmat>, KNM< upscaled_type<PetscScalar> >*, PetscScalar, 'N'>::init));
67206788
Global.Add("set", "(", new PETSc::setOptions< Dmat >( ));
67216789
Global.Add("set", "(", new PETSc::setOptions< Dmat >(1));
67226790
Global.Add("set", "(", new PETSc::setOptions< Dmat >(1, 1));

0 commit comments

Comments
 (0)