Skip to content

Commit 38c51c8

Browse files
phtournierphtournier
committed
- can now define non-square BIOs on a mesh X its submesh
- HMatrix transpose & multi-rhs matrix-vector product
1 parent dce0418 commit 38c51c8

4 files changed

Lines changed: 162 additions & 100 deletions

File tree

‎3rdparty/getall‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -185,10 +185,10 @@ download('hpddm','https://codeload.github.com/hpddm/hpddm/zip/7113b9a6b77fceee3f
185185
'hpddm.zip',
186186
'6910b7b974f0b60d9c247c666e7f3862');
187187

188-
download('bemtool','https://github.com/PierreMarchand20/BemTool/archive/0a884602b2fd08b72f2a14196650dd8526b840d6.zip',
188+
download('bemtool','https://github.com/PierreMarchand20/BemTool/archive/661cf0ed86998a21738c421fa14d01a355536b68.zip',
189189
'https://github.com/PierreMarchand20/BemTool',
190190
'bemtool.zip',
191-
'66d3589f3cae096267b044d163d593a5');
191+
'ee8a3f65ad4f2704b3212c2409540e9d');
192192

193193
download('Boost','https://www.ljll.fr/~tournier/boost_for_bemtool.tar.gz',
194194
'https://www.boost.org',

‎plugin/mpi/bem.cpp‎

Lines changed: 28 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -67,22 +67,33 @@ template<class Type, class K, int init>
6767
AnyType To(Stack stack,Expression emat,Expression einter)
6868
{ return To<Type, K>(stack,emat,einter,init);}
6969

70-
template<class V, class K>
70+
template<class V, class K, char N>
7171
class Prod {
7272
public:
7373
const HMatrixVirt<K>* h;
7474
const V u;
7575
Prod(HMatrixVirt<K>** v, V w) : h(*v), u(w) {}
7676

77-
void prod(V x) const {int mu = this->u->n/h->nb_cols(); h->mvprod_global(*(this->u), *x, mu);};
77+
void prod(V x) const {int mu = this->u->n/h->nb_cols(); h->mvprod_global(*(this->u), *x, 'N', mu);};
78+
void prodT(V x) const {int mu = this->u->n/h->nb_rows(); h->mvprod_global(*(this->u), *x, 'C', mu);};
7879

79-
static V mv(V Ax, Prod<V, K> A) {
80+
static V mv(V Ax, Prod<V, K, 'N'> A) {
8081
*Ax = K();
8182
A.prod(Ax);
8283
return Ax;
8384
}
84-
static V init(V Ax, Prod<V, K> A) {
85-
Ax->init(A.u->n);
85+
static V init(V Ax, Prod<V, K, 'N'> A) {
86+
Ax->init(A.h->nb_rows());
87+
return mv(Ax, A);
88+
}
89+
90+
static V mv(V Ax, Prod<V, K, 'T'> A) {
91+
*Ax = K();
92+
A.prodT(Ax);
93+
return Ax;
94+
}
95+
static V init(V Ax, Prod<V, K, 'T'> A) {
96+
Ax->init(A.h->nb_cols());
8697
return mv(Ax, A);
8798
}
8899

@@ -577,10 +588,18 @@ void addHmat() {
577588

578589
Add<HMatrixVirt<K>**>("infos",".",new OneOperator1_<std::map<std::string, std::string>*, HMatrixVirt<K>**>(get_infos));
579590

580-
Dcl_Type<Prod<KN<K>*, K>>();
581-
TheOperators->Add("*", new OneOperator2<Prod<KN<K>*, K>, HMatrixVirt<K>**, KN<K>*>(Build));
582-
TheOperators->Add("=", new OneOperator2<KN<K>*, KN<K>*, Prod<KN<K>*, K>>(Prod<KN<K>*, K>::mv));
583-
TheOperators->Add("<-", new OneOperator2<KN<K>*, KN<K>*, Prod<KN<K>*, K>>(Prod<KN<K>*, K>::init));
591+
Dcl_Type<Prod<KN<K>*, K, 'N'>>();
592+
TheOperators->Add("*", new OneOperator2<Prod<KN<K>*, K, 'N'>, HMatrixVirt<K>**, KN<K>*>(Build));
593+
TheOperators->Add("=", new OneOperator2<KN<K>*, KN<K>*, Prod<KN<K>*, K, 'N'>>(Prod<KN<K>*, K, 'N'>::mv));
594+
TheOperators->Add("<-", new OneOperator2<KN<K>*, KN<K>*, Prod<KN<K>*, K, 'N'>>(Prod<KN<K>*, K, 'N'>::init));
595+
596+
Dcl_Type<Prod<KN<K>*, K, 'T'>>();
597+
Dcl_Type<OpTrans<HMatrixVirt<K>*>>();
598+
TheOperators->Add("\'", new OneOperator1<OpTrans<HMatrixVirt<K>*>, HMatrixVirt<K>**>(Build));
599+
TheOperators->Add("*", new OneOperator2<Prod<KN<K>*, K, 'T'>, OpTrans<HMatrixVirt<K>*>, KN<K>*>(Build));
600+
TheOperators->Add("=", new OneOperator2<KN<K>*, KN<K>*, Prod<KN<K>*, K, 'T'>>(Prod<KN<K>*, K, 'T'>::mv));
601+
TheOperators->Add("<-", new OneOperator2<KN<K>*, KN<K>*, Prod<KN<K>*, K, 'T'>>(Prod<KN<K>*, K, 'T'>::init));
602+
584603

585604
SetHMatrix_Op<K>::btype = Dcl_Type<const SetHMatrix_Op<K> * >();
586605
Global.Add("set","(",new SetHMatrix<K>);
@@ -699,9 +718,6 @@ AnyType OpHMatrixtoBEMForm<R,MMesh,v_fes1,v_fes2>::Op::operator()(Stack stack)
699718
SetEnd_Data_Bem_Solver<R>(stack,ds, b->nargs,OpCall_FormBilinear_np::n_name_param); // LIST_NAME_PARM_HMAT
700719
WhereStackOfPtr2Free(stack)=new StackOfPtr2Free(stack);
701720

702-
bool samemesh = (void*)&Uh->Th == (void*)&Vh->Th; // same Fem2D::Mesh +++ pot or kernel
703-
if (VFBEM==1)
704-
ffassert (samemesh);
705721
if(init)
706722
*Hmat =0;
707723
*Hmat =0;

‎plugin/mpi/bem.hpp‎

Lines changed: 36 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -87,7 +87,7 @@ class HMatrixVirt {
8787
int state = 0;
8888
bool factinplace = 0;
8989
virtual const std::map<std::string, std::string>& get_infos() const = 0;
90-
virtual void mvprod_global(const K* const in, K* const out,const int& mu=1) const = 0;
90+
virtual void mvprod_global(const K* const in, K* const out, const char N='N', const int& mu=1) const = 0;
9191
virtual int nb_rows() const = 0;
9292
virtual int nb_cols() const = 0;
9393
virtual void cluster_to_target_permutation(const K* const in, K* const out) const = 0;
@@ -124,15 +124,20 @@ class HMatrixImpl : public HMatrixVirt<K> {
124124
std::map<std::string, std::string> infos;
125125
HMatrixImpl(htool::VirtualGenerator<K> &mat, std::shared_ptr<htool::Cluster<double>> t, std::shared_ptr<htool::Cluster<double>> s,const htool::HMatrixTreeBuilder<K,double>& hmatrix_tree_builder,string slvr,bool factinplace, MPI_Comm comm) : target_cluster(t), source_cluster(s), distributed_operator_holder(mat,*target_cluster,*source_cluster,hmatrix_tree_builder,comm),distributed_operator(distributed_operator_holder.distributed_operator), H(distributed_operator_holder.hmatrix), infos(htool::get_distributed_hmatrix_information(H,distributed_operator.get_comm())){this->solver=slvr; this->factinplace=factinplace;}
126126
const std::map<std::string, std::string>& get_infos() const { return infos; }
127-
void mvprod_global(const K* const in, K* const out,const int& mu=1) const {
127+
void mvprod_global(const K* const in, K* const out, const char N='N', const int& mu=1) const {
128128
K* work=nullptr;
129129
if (mu==1){
130-
add_distributed_operator_vector_product_global_to_global('N',K(1),distributed_operator,in,K(0),out,work);
131-
} else{
130+
add_distributed_operator_vector_product_global_to_global(N,K(1),distributed_operator,in,K(0),out,work);
131+
} else if (N=='N'){
132132
htool::MatrixView<K> out_view(distributed_operator.get_target_partition().get_global_size() ,mu,out);
133133
htool::MatrixView< const K> in_view(distributed_operator.get_source_partition().get_global_size() ,mu,in);
134-
add_distributed_operator_matrix_product_global_to_global('N',K(1),distributed_operator,in_view,K(0),out_view,work);
135-
}
134+
add_distributed_operator_matrix_product_global_to_global(N,K(1),distributed_operator,in_view,K(0),out_view,work);
135+
}
136+
else {
137+
htool::MatrixView<K> out_view(distributed_operator.get_source_partition().get_global_size() ,mu,out);
138+
htool::MatrixView< const K> in_view(distributed_operator.get_target_partition().get_global_size() ,mu,in);
139+
add_distributed_operator_matrix_product_global_to_global(N,K(1),distributed_operator,in_view,K(0),out_view,work);
140+
}
136141
}
137142
int nb_rows() const { return distributed_operator.get_target_partition().get_global_size();}
138143
int nb_cols() const { return distributed_operator.get_source_partition().get_global_size();}
@@ -463,9 +468,13 @@ void creationHMatrixtoBEMForm(const FESpace1 * Uh, const FESpace2 * Vh, const in
463468
typedef typename TMesh::RdHat TRdHat;
464469

465470
typedef typename std::conditional<SMesh::RdHat::d==1,bemtool::Mesh1D,bemtool::Mesh2D>::type MeshBemtool;
471+
typedef typename std::conditional<TMesh::RdHat::d==1,bemtool::Mesh1D,bemtool::Mesh2D>::type MeshBemtoolX;
466472
typedef typename std::conditional<SMesh::RdHat::d==1,bemtool::P0_1D,bemtool::P0_2D>::type P0;
467473
typedef typename std::conditional<SMesh::RdHat::d==1,bemtool::P1_1D,bemtool::P1_2D>::type P1;
468474
typedef typename std::conditional<SMesh::RdHat::d==1,bemtool::P2_1D,bemtool::P2_2D>::type P2;
475+
typedef typename std::conditional<TMesh::RdHat::d==1,bemtool::P0_1D,bemtool::P0_2D>::type P0X;
476+
typedef typename std::conditional<TMesh::RdHat::d==1,bemtool::P1_1D,bemtool::P1_2D>::type P1X;
477+
typedef typename std::conditional<TMesh::RdHat::d==1,bemtool::P2_1D,bemtool::P2_2D>::type P2X;
469478

470479
// size of the matrix
471480
int m=Uh->NbOfDF;
@@ -486,9 +495,10 @@ void creationHMatrixtoBEMForm(const FESpace1 * Uh, const FESpace2 * Vh, const in
486495
const TMesh & ThV =Vh->Th; // colunm
487496
bool samemesh = (void*)&Uh->Th == (void*)&Vh->Th; // same Fem2D::Mesh +++ pot or kernel
488497

489-
bemtool::Geometry node; MeshBemtool mesh;
498+
bemtool::Geometry node; bemtool::Geometry nodeX; MeshBemtool mesh; MeshBemtoolX meshX;
490499
Mesh2Bemtool(ThU, node, mesh);
491-
500+
if (!samemesh && (VFBEM==1))
501+
Mesh2Bemtool(ThV, nodeX, meshX);
492502

493503
vector<double> pt(3*n);
494504
vector<double> ps(3*m);
@@ -634,15 +644,30 @@ void creationHMatrixtoBEMForm(const FESpace1 * Uh, const FESpace2 * Vh, const in
634644

635645
if (SP0) {
636646
bemtool::Dof<P0> dof(mesh);
637-
ff_BIO_Generator<R,P0,SMesh>(generator,Ker,dof,alpha);
647+
if (samemesh)
648+
ff_BIO_Generator<R,P0,P0,SMesh,SMesh>(generator,Ker,dof,dof,alpha);
649+
else {
650+
bemtool::Dof<P0X> dofX(meshX);
651+
ff_BIO_Generator<R,P0X,P0,TMesh,SMesh>(generator,Ker,dofX,dof,alpha);
652+
}
638653
}
639654
else if (SP1) {
640655
bemtool::Dof<P1> dof(mesh,true);
641-
ff_BIO_Generator<R,P1,SMesh>(generator,Ker,dof,alpha);
656+
if (samemesh)
657+
ff_BIO_Generator<R,P1,P1,SMesh,SMesh>(generator,Ker,dof,dof,alpha);
658+
else {
659+
bemtool::Dof<P1X> dofX(meshX,true);
660+
ff_BIO_Generator<R,P1X,P1,TMesh,SMesh>(generator,Ker,dofX,dof,alpha);
661+
}
642662
}
643663
else if (SP2) {
644664
bemtool::Dof<P2> dof(mesh,true);
645-
ff_BIO_Generator<R,P2,SMesh>(generator,Ker,dof,alpha);
665+
if (samemesh)
666+
ff_BIO_Generator<R,P2,P2,SMesh,SMesh>(generator,Ker,dof,dof,alpha);
667+
else {
668+
bemtool::Dof<P2X> dofX(meshX,true);
669+
ff_BIO_Generator<R,P2X,P2,TMesh,SMesh>(generator,Ker,dofX,dof,alpha);
670+
}
646671
}
647672
else if (SRT0 && SRdHat::d == 2) {
648673
// BemKernel->typeKernel[0] == 6 :: MA_SL

0 commit comments

Comments
 (0)