Parallelization

True parallel Schwarz example

Thank you to F. Nataf

This is a explanation of the two examples MPI-GMRES 2D and MPI-GMRES 3D, a Schwarz parallel with a complexity almost independent of the number of process (with a coarse grid preconditioner).

To solve the following Poisson problem on domain \(\Omega\) with boundary \(\Gamma\) in \(L^2(\Omega)\) :

where \(f\) and \(g\) are two given functions of \(L^2(\Omega)\) and of \(H^{\frac12}(\Gamma)\),

Lets introduce \((\pi_i)_{i=1,.., N_p}\) a regular partition of the unity of \(\Omega\), q-e-d:

Denote \(\Omega_i\) the sub domain which is the support of \(\pi_i\) function and also denote \(\Gamma_i\) the boundary of \(\Omega_i\).

The parallel Schwarz method is:

Let \(\ell=0\) the iterator and a initial guest \(u^0\) respecting the boundary condition (i.e. \(u^0_{|\Gamma} = g\)).

After discretization with the Lagrange finite element method, with a compatible mesh \({\mathcal{T}_h}_i\) of \(\Omega_i\), i. e., the exist a global mesh \({\mathcal{T}_h}\) such that \({\mathcal{T}_h}_i\) is include in \({\mathcal{T}_h}\).

Let us denote:

• \({V_h}_i\) the finite element space corresponding to domain \(\Omega_i\),

• \({\mathcal{N}_h}_i\) is the set of the degree of freedom \(\sigma_i^k\),

• \({\mathcal{N}^{\Gamma_i}_{hi}}\) is the set of the degree of freedom of \({V_h}_i\) on the boundary \(\Gamma_i\) of \(\Omega_i\),

• \(\sigma_i^k({v_h})\) is the value the degree of freedom \(k\),

• \({V_{0h}}_i= \{ {v_h} \in {V_h}_i :\forall k \in {\mathcal{N}^{\Gamma_i}_{hi}}, \quad \sigma_i^k({v_h})=0 \}\),

• The conditional expression \(a\;?\;b:c\) is defined like in :c`C` of C++ language by

  \[\begin{split}a?b: c \equiv \left\{ \begin{array}{l} \mbox{if } a \mbox{ is true then return b}\\ \mbox{else return } c\\ \end{array} \right..\end{split}\]

We have to defined to operator to build the previous algorithm:

We denote \({u_h^{\ell}}_{|i}\) the restriction of \(u_h^\ell\) on \({V_h}_i\), so the discrete problem on \(\Omega_i\) of problem (37) is find \({u_h^{\ell}}_{i}\in {V_h}_i\) such that:

where \(g_i^k\) is the value of \(g\) associated to the degree of freedom \(k\in {\mathcal{N}^{\Gamma_i}_{hi}}\).

In FreeFEM, it can be written has with U is the vector corresponding to \({u_h^{\ell}}_{|i}\) and the vector U1 is the vector corresponding to \({u_h^{\ell}}_{i}\) is the solution of:

1real[int] U1(Ui.n);
2real[int] b = onG .* U;
3b = onG ? b : Bi ;
4U1 = Ai^-1*b;

where \(\mathtt{onG}[i] =(i \in \Gamma_i\setminus\Gamma) ? 1 : 0\), and \(\mathtt{Bi}\) the right of side of the problem, are defined by

 1// Fespace
 2fespace Whi(Thi, P2);
 3
 4// Problem
 5varf vPb (U, V)
 6    = int3d(Thi)(
 7          grad(U)'*grad(V)
 8    )
 9    + int3d(Thi)(
10          F*V
11    )
12    + on(1, U=g)
13    + on(10, U=G)
14    ;
15
16varf vPbon (U, V) = on(10, U=1) + on(1, U=0);
17
18matrix Ai = vPb (Whi, Whi, solver=sparsesolver);
19real[int] onG = vPbon(0, Whi);
20real[int] Bi=vPb(0, Whi);

where the FreeFEM label of \(\Gamma\) is 1 and the label of \(\Gamma_i\setminus \Gamma\) is \(10\).

To build the transfer/update part corresponding to (38) equation on process \(i\), let us call njpart the number the neighborhood of domain of \(\Omega_i\) (i.e: \(\pi_j\) is none \(0\) of \(\Omega_i\)), we store in an array jpart of size njpart all this neighborhood.

Let us introduce two array of matrix, Smj[j] to defined the vector to send from \(i\) to \(j\) a neighborhood process, and the matrix \(rMj[j]\) to after to reduce owith neighborhood \(j\) domain.

So the tranfert and update part compute \(v_i= \pi_i u_i + \sum_{j\in J_i} \pi_j u_j\) and can be write the FreeFEM function Update:

 1func bool Update (real[int] &ui, real[int] &vi){
 2    int n = jpart.n;
 3    for (int j = 0; j < njpart; ++j) Usend[j][] = sMj[j]*ui;
 4    mpiRequest[int] rq(n*2);
 5    for (int j = 0; j < n; ++j) Irecv(processor(jpart[j], comm,rq[j]), Ri[j][]);
 6    for (int j = 0; j < n; ++j) Isend(processor(jpart[j], comm, rq[j+n]), Si[j][]);
 7    for (int j = 0; j < n*2; ++j) int k = mpiWaitAny(rq);
 8    // apply the unity local partition
 9    vi = Pii*ui; //set to pi_i u_i
10    for (int j = 0; j < njpart; ++j) vi += rMj[j]*Vrecv[j][]; //add pi_j u_j
11    return true;
12}

where the buffer are defined by:

1InitU(njpart, Whij, Thij, aThij, Usend) //defined the send buffer
2InitU(njpart, Whij, Thij, aThij, Vrecv) //defined the revc buffer

with the following macro definition:

1macro InitU(n, Vh, Th, aTh, U) Vh[int] U(n); for (int j = 0; j < n; ++j){Th = aTh[j]; U[j] = 0;}

First GMRES algorithm: you can easily accelerate the fixed point algorithm by using a parallel GMRES algorithm after the introduction the following affine \(\mathcal{A}_i\) operator sub domain \(\Omega_i\).

1func real[int] DJ0 (real[int]& U){
2    real[int] V(U.n), b = onG .* U;
3    b = onG ? b : Bi ;
4    V = Ai^-1*b;
5    Update(V, U);
6    V -= U;
7    return V;
8}

Where the parallel MPIGMRES or MPICG algorithm is just a simple way to solve in parallel the following \(A_i x_i = b_i, i = 1, .., N_p\) by just changing the dot product by reduce the local dot product of all process with the following MPI code:

1template<class R> R ReduceSum1(R s, MPI_Comm *comm){
2    R r = 0;
3    MPI_Allreduce(&s, &r, 1, MPI_TYPE<R>::TYPE(), MPI_SUM, *comm );
4    return r;
5}

This is done in MPIGC dynamics library tool.

Second GMRES algorithm: Use scharwz algorithm as a preconditioner of basic GMRES method to solving the parallel problem.

 1func real[int] DJ (real[int]& U){ //the original problem
 2    ++kiter;
 3    real[int] V(U.n);
 4    V = Ai*U;
 5    V = onGi ? 0.: V; //remove boundary term
 6    return V;
 7}
 8
 9func real[int] PDJ (real[int]& U){ //the preconditioner
10    real[int] V(U.n);
11    real[int] b = onG ? 0. : U;
12    V = Ai^-1*b;
13    Update(V, U);
14    return U;
15}

Third GMRES algorithm: Add a coarse solver to the previous algorithm

First build a coarse grid on processor 0, and the

 1matrix AC, Rci, Pci;
 2if (mpiRank(comm) == 0)
 3    AC = vPbC(VhC, VhC, solver=sparsesolver); //the coarse problem
 4
 5Pci = interpolate(Whi, VhC); //the projection on coarse grid
 6Rci = Pci'*Pii; //the restriction on Process i grid with the partition pi_i
 7
 8func bool CoarseSolve (real[int]& V, real[int]& U, mpiComm& comm){
 9    // solving the coarse problem
10    real[int] Uc(Rci.n), Bc(Uc.n);
11    Uc = Rci*U;
12    mpiReduce(Uc, Bc, processor(0, comm), mpiSUM);
13    if (mpiRank(comm) == 0)
14    Uc = AC^-1*Bc;
15    broadcast(processor(0, comm), Uc);
16    V = Pci*Uc;
17}

The New preconditionner

 1func real[int] PDJC (real[int]& U){
 2    // Idea: F. Nataf.
 3    // 0 ~ (I C1A)(I-C2A) => I ~ - C1AC2A +C1A +C2A
 4    // New Prec P= C1+C2 - C1AC2 = C1(I- A C2) +C2
 5    // ( C1(I- A C2) +C2 ) Uo
 6    // V = - C2*Uo
 7    // ....
 8    real[int] V(U.n);
 9    CoarseSolve(V, U, comm);
10    V = -V; //-C2*Uo
11    U += Ai*V; //U = (I-A C2) Uo
12    real[int] b = onG ? 0. : U;
13    U = Ai^-1*b; //C1( I -A C2) Uo
14    V = U - V;
15    Update(V, U);
16    return U;
17}

The code of the 4 algorithms:

 1real epss = 1e-6;
 2int rgmres = 0;
 3if (gmres == 1){
 4    rgmres = MPIAffineGMRES(DJ0, u[], veps=epss, nbiter=300,
 5        comm=comm, dimKrylov=100, verbosity=ipart?0: 50);
 6    real[int] b = onG .* u[];
 7    b = onG ? b : Bi ;
 8    v[] = Ai^-1*b;
 9    Update(v[], u[]);
10}
11else if (gmres == 2)
12    rgmres = MPILinearGMRES(DJ, precon=PDJ, u[], Bi, veps=epss,
13        nbiter=300, comm=comm, dimKrylov=100, verbosity=ipart?0: 50);
14else if (gmres == 3)
15    rgmres = MPILinearGMRES(DJ, precon=PDJC, u[], Bi, veps=epss,
16        nbiter=300, comm=comm, dimKrylov=100, verbosity=ipart?0: 50);
17else //algo Shwarz for demo
18    for(int iter = 0; iter < 10; ++iter)
19        ...

We have all ingredient to solve in parallel if we have et the partitions of the unity. To build this partition we do:

The initial step on process \(1\) to build a coarse mesh, \({\mathcal{T}_h}^*\) of the full domain, and build the partition \(\pi\) function constant equal to \(i\) on each sub domain \(\mathcal{O}_i, i =1 ,.., N_p\), of the grid with the metis graph partitioner [KARYPIS1995] and on each process \(i\) in \(1..,N_p\) do

1. Broadcast from process \(1\), the mesh \({\mathcal{T}_h}^*\) (call Thii in FreeFEM script), and \(\pi\) function,

2. remark that the characteristic function \(\mathrm{1\!\!I}_{\mathcal{O}_i}\) of domain \(\mathcal{O}_i\), is defined by \((\pi=i)?1:0\),

3. Let us call \(\Pi^2_P\) (resp. \(\Pi^2_V\)) the \(L^2\) on \(P_h^*\) the space of the constant finite element function per element on \({\mathcal{T}_h}^*\) (resp. \(V_h^*\) the space of the affine continuous finite element per element on \({\mathcal{T}_h}^*\)) and build in parallel the \(\pi_i\) and \(\Omega_i\), such that \(\mathcal{O}_i\ \subset \Omega_i\) where \(\mathcal{O}_i= supp ((\Pi^2_V \Pi^2_C)^m \mathrm{1\!\!I}_{O_i})\), and \(m\) is a the overlaps size on the coarse mesh (generally one), (this is done in function AddLayers(Thii,suppii[],nlayer,phii[]); We choose a function \(\pi^*_i = (\Pi^2_1 \Pi^2_0)^m \mathrm{1\!\!I}_{\mathcal{O}_i}\) so the partition of the unity is simply defined by

   \[\pi_i = \frac{\pi_i^*}{\sum_{j=1}^{N_p} \pi_j^*}\]

   The set \(J_i\) of neighborhood of the domain \(\Omega_i\), and the local version on \(V_{hi}\) can be defined the array jpart and njpart with:

    1Vhi pii = piistar;
    2Vhi[int] pij(npij); //local partition of 1 = pii + sum_j pij[j]
    3int[int] jpart(npart);
    4int njpart = 0;
    5Vhi sumphi = piistar;
    6for (int i = 0; i < npart; ++i)
    7    if (i != ipart){
    8        if (int3d(Thi)(pijstar,j) > 0){
    9            pij[njpart] = pijstar;
   10            sumphi[] += pij[njpart][];
   11            jpart[njpart++] = i;
   12        }
   13    }
   14pii[] = pii[] ./ sumphi[];
   15for (int j = 0; j < njpart; ++j)
   16pij[j][] = pij[j][] ./ sumphi[];
   17jpart.resize(njpart);
   

4. We call \({\mathcal{T}_h}^*_{ij}\) the sub mesh part of \({\mathcal{T}_h}_i\) where \(\pi_j\) are none zero. And thanks to the function trunc to build this array,

   1for(int jp = 0; jp < njpart; ++jp)
   2    aThij[jp] = trunc(Thi, pij[jp] > 1e-10, label=10);
   

5. At this step we have all on the coarse mesh, so we can build the fine final mesh by splitting all meshes: Thi, Thij[j], Thij[j] with FreeFEM trunc mesh function which do restriction and slipping.

6. The construction of the send/recv matrices sMj and freefem:`rMj: can done with this code:

    1mesh3 Thij = Thi;
    2fespace Whij(Thij, Pk);
    3matrix Pii; Whi wpii = pii; Pii = wpii[]; //Diagonal matrix corresponding X pi_i
    4matrix[int] sMj(njpart), rMj(njpart); //M send/recive case
    5for (int jp = 0; jp < njpart; ++jp){
    6    int j = jpart[jp];
    7    Thij = aThij[jp]; //change mesh to change Whij, Whij
    8    matrix I = interpolate(Whij, Whi); //Whij <- Whi
    9    sMj[jp] = I*Pii; //Whi -> s Whij
   10    rMj[jp] = interpolate(Whij, Whi, t=1); //Whij -> Whi
   11}
   

To buil a not too bad application, all variables come from parameters value with the following code

1include "getARGV.idp"
2verbosity = getARGV("-vv", 0);
3int vdebug = getARGV("-d", 1);
4int ksplit = getARGV("-k", 10);
5int nloc = getARGV("-n", 25);
6string sff = getARGV("-p, ", "");
7int gmres = getARGV("-gmres", 3);
8bool dplot = getARGV("-dp", 0);
9int nC = getARGV("-N", max(nloc/10, 4));

And small include to make graphic in parallel of distribute solution of vector \(u\) on mesh \(T_h\) with the following interface:

1include "MPIplot.idp"
2func bool plotMPIall(mesh &Th, real[int] &u, string cm){
3    PLOTMPIALL(mesh, Pk, Th, u, {cmm=cm, nbiso=20, fill=1, dim=3, value=1});
4    return 1;
5}

MUMPS solver

MUltifrontal Massively Parallel Solver (MUMPS) is an open-source library.

This package solves linear system of the form \(A \: x = b\) where \(A\) is a square sparse matrix with a direct method. The square matrix considered in MUMPS can be either unsymmetric, symmetric positive definite or general symmetric.

The method implemented in MUMPS is a direct method based on a multifrontal approach. It constructs a direct factorization \(A \:= \: L\:U\), \(A\: = \: L^t \: D \: L\) depending of the symmetry of the matrix \(A\).

MUMPS uses the following libraries :

• BLAS,

• BLACS,

• ScaLAPACK.

MUMPS parameters:

There are four input parameters in MUMPS. Two integers SYM and PAR, a vector of integer of size 40 INCTL and a vector of real of size 15 CNTL.

The first parameter gives the type of the matrix: 0 for unsymmetric matrix, 1 for symmetric positive matrix and 2 for general symmetric.

The second parameter defined if the host processor work during the factorization and solves steps : PAR=1 host processor working and PAR=0 host processor not working.

The parameter INCTL and CNTL is the control parameter of MUMPS. The vectors ICNTL and CNTL in MUMPS becomes with index 1 like vector in Fortran. For more details see the MUMPS user’s guide.

We describe now some elements of the main parameters of ICNTL for MUMPS.

• Input matrix parameter The input matrix is controlled by parameters ICNTL(5) and ICNTL(18).

  The matrix format (resp. matrix pattern and matrix entries) are controlled by INCTL(5) (resp. INCTL(18)).

  The different values of ICNTL(5) are 0 for assembled format and 1 for element format. In the current release of FreeFEM, we consider that FE matrix or matrix is storage in assembled format. Therefore, INCTL(5) is treated as 0 value.

  The main option for ICNTL(18): INCLTL(18)=0 centrally on the host processor, ICNTL(18)=3 distributed the input matrix pattern and the entries (recommended option for distributed matrix by developer of MUMPS). For other values of ICNTL(18) see the MUMPS user’s guide. These values can be used also in FreeFEM.

  The default option implemented in FreeFEM are ICNTL(5)=0 and ICNTL(18)=0.

• Preprocessing parameter The preprocessed matrix \(A_{p}\) that will be effectively factored is defined by

  \[A_{p} = P \: D_r \: A \: Q_c \ D_c P^t\]

  where \(P\) is the permutation matrix, \(Q_c\) is the column permutation, \(D_r\) and \(D_c\) are diagonal matrix for respectively row and column scaling.

  The ordering strategy to obtain \(P\) is controlled by parameter ICNTL(7). The permutation of zero free diagonal \(Q_c\) is controlled by parameter ICNTL(6). The row and column scaling is controlled by parameter ICNTL(18). These option are connected and also strongly related with ICNTL(12) (see the MUMPS user’s guide for more details).

  The parameters permr, scaler, and scalec in FreeFEM allow to give permutation matrix(\(P\)), row scaling (\(D_r\)) and column scaling (\(D_c\)) of the user respectively.

Calling MUMPS in FreeFEM

To call MUMPS in FreeFEM, we need to load the dynamic library MUMPS_freefem.dylib (MacOSX), MUMPS_freefem.so (Unix) or MUMPS_freefem.dll (Windows).

This is done in typing load "MUMPS_FreeFem" in the .edp file. We give now the two methods to give the option of MUMPS solver in FreeFEM.

• Solver parameters is defined in .edp file: In this method, we need to give the parameters lparams and dparams.

  These parameters are defined for MUMPS by :

  • lparams[0] = SYM, lparams[1] = PAR,

  • \(\forall i\) = 1,…,40, lparams[i+1] = ICNTL(i)

  • \(\forall i\) = 1,…,15, dparams[i-1] = CNTL(i)

• Reading solver parameters on a file:

  The structure of data file for MUMPS in FreeFEM is : first line parameter SYM and second line parameter PAR and in the following line the different value of vectors ICNTL and CNTL. An example of this parameter file is given in ffmumpsfileparam.txt.

   10 /* SYM :: 0 for non symmetric matrix, 1 for symmetric definite positive matrix and 2 general symmetric matrix*/
   21 /* PAR :: 0 host not working during factorization and solves steps, 1 host working during factorization and solves steps*/
   3-1 /* ICNTL(1) :: output stream for error message */
   4-1 /* ICNTL(2) :: output for diagnostic printing, statics and warning message */
   5-1 /* ICNTL(3) :: for global information */
   60 /* ICNTL(4) :: Level of printing for error, warning and diagnostic message */
   70 /* ICNTL(5) :: matrix format : 0 assembled format, 1 elemental format. */
   87 /* ICNTL(6) :: control option for permuting and/or scaling the matrix in analysis phase */
   93 /* ICNTL(7) :: pivot order strategy : AMD, AMF, metis, pord scotch*/
  1077 /* ICNTL(8) :: Row and Column scaling strategy */
  111 /* ICNTL(9) :: 0 solve Ax = b, 1 solve the transposed system A^t x = b : parameter is not considered in the current release of FreeFEM*/
  120 /* ICNTL(10) :: number of steps of iterative refinement */
  130 /* ICNTL(11) :: statics related to linear system depending on ICNTL(9) */
  141 /* ICNTL(12) :: constrained ordering strategy for general symmetric matrix */
  150 /* ICNTL(13) :: method to control splitting of the root frontal matrix */
  1620 /* ICNTL(14) :: percentage increase in the estimated working space (default 20\%)*/
  170 /* ICNTL(15) :: not used in this release of MUMPS */
  180 /* ICNTL(16) :: not used in this release of MUMPS */
  190 /* ICNTL(17) :: not used in this release of MUMPS */
  203 /* ICNTL(18) :: method for given : matrix pattern and matrix entries : */
  210 /* ICNTL(19) :: method to return the Schur complement matrix */
  220 /* ICNTL(20) :: right hand side form ( 0 dense form, 1 sparse form) : parameter will be set to 0 for FreeFEM */
  230 /* ICNTL(21) :: 0, 1 kept distributed solution : parameter is not considered in the current release of FreeFEM */
  240 /* ICNTL(22) :: controls the in-core/out-of-core (OOC) facility */
  250 /* ICNTL(23) :: maximum size of the working memory in Megabyte than MUMPS can allocate per working processor */
  260 /* ICNTL(24) :: control the detection of null pivot */
  270 /* ICNTL(25) :: control the computation of a null space basis */
  280 /* ICNTL(26) :: This parameter is only significant with Schur option (ICNTL(19) not zero). : parameter is not considered in the current release of FreeFEM */
  29-8 /* ICNTL(27) (Experimental parameter subject to change in next release of MUMPS) :: control the blocking factor for multiple righthand side during the solution phase : parameter is not considered in the current release of FreeFEM */
  300 /* ICNTL(28) :: not used in this release of MUMPS*/
  310 /* ICNTL(29) :: not used in this release of MUMPS*/
  320 /* ICNTL(30) :: not used in this release of MUMPS*/
  330 /* ICNTL(31) :: not used in this release of MUMPS*/
  340 /* ICNTL(32) :: not used in this release of MUMPS*/
  350 /* ICNTL(33) :: not used in this release of MUMPS*/
  360 /* ICNTL(34) :: not used in this release of MUMPS*/
  370 /* ICNTL(35) :: not used in this release of MUMPS*/
  380 /* ICNTL(36) :: not used in this release of MUMPS*/
  390 /* ICNTL(37) :: not used in this release of MUMPS*/
  400 /* ICNTL(38) :: not used in this release of MUMPS*/
  411 /* ICNTL(39) :: not used in this release of MUMPS*/
  420 /* ICNTL(40) :: not used in this release of MUMPS*/
  430.01 /* CNTL(1) :: relative threshold for numerical pivoting */
  441e-8 /* CNTL(2) :: stopping criteria for iterative refinement */
  45-1 /* CNTL(3) :: threshold for null pivot detection */
  46-1 /* CNTL(4) :: determine the threshold for partial pivoting */
  470.0 /* CNTL(5) :: fixation for null pivots */
  480 /* CNTL(6) :: not used in this release of MUMPS */
  490 /* CNTL(7) :: not used in this release of MUMPS */
  500 /* CNTL(8) :: not used in this release of MUMPS */
  510 /* CNTL(9) :: not used in this release of MUMPS */
  520 /* CNTL(10) :: not used in this release of MUMPS */
  530 /* CNTL(11) :: not used in this release of MUMPS */
  540 /* CNTL(12) :: not used in this release of MUMPS */
  550 /* CNTL(13) :: not used in this release of MUMPS */
  560 /* CNTL(14) :: not used in this release of MUMPS */
  570 /* CNTL(15) :: not used in this release of MUMPS */
  

  If no solver parameter is given, we used default option of MUMPS solver.

SuperLU distributed solver

The package SuperLU_DIST solves linear systems using LU factorization. It is a free scientific library

This library provides functions to handle square or rectangular matrix in real and complex arithmetics. The method implemented in SuperLU_DIST is a supernodal method. New release of this package includes a parallel symbolic factorization. This scientific library is written in C and MPI for communications.

SuperLU_DIST parameters:

We describe now some parameters of SuperLU_DIST. The SuperLU_DIST library use a 2D-logical process group. This process grid is specified by \(nprow\) (process row) and \(npcol\) (process column) such that \(N_{p} = nprow \: npcol\) where \(N_{p}\) is the number of all process allocated for SuperLU_DIST.

The input matrix parameters is controlled by “matrix=” in sparams for internal parameter or in the third line of parameters file. The different value are

• matrix=assembled global matrix are available on all process

• matrix=distributedglobal The global matrix is distributed among all the process

• matrix=distributed The input matrix is distributed (not yet implemented)

The option arguments of SuperLU_DIST are described in the section Users-callable routine of the SuperLU users’ guide.

The parameter Fact and TRANS are specified in FreeFEM interfaces to SuperLU_DIST during the different steps. For this reason, the value given by the user for this option is not considered.

The factorization LU is calculated in SuperLU_DIST on the matrix \(A_p\).

where \(P_c\) and \(P_r\) is the row and column permutation matrix respectively, \(D_r\) and \(D_c\) are diagonal matrix for respectively row and column scaling.

The option argument RowPerm (resp. ColPerm) control the row (resp. column) permutation matrix. \(D_r\) and \(D_c\) is controlled by the parameter DiagScale.

The parameter permr, permc, scaler, and scalec in FreeFEM is provided to give row permutation, column permutation, row scaling and column scaling of the user respectively.

The other parameters for LU factorization are ParSymFact and ReplaceTinyPivot. The parallel symbolic factorization works only on a power of two processes and need the ParMetis ordering. The default option argument of SuperLU_DIST are given in the file ffsuperlu_dist_fileparam.txt.

Calling SuperLU_DIST in FreeFEM

To call SuperLU_DIST in FreeFEM, we need to load the library dynamic correspond to interface. This done by the following line load "real_superlu _DIST_FreeFem" (resp. load "complex_superlu_DIST_FreeFem") for real (resp. complex) arithmetics in the file .edp.

Solver parameters is defined in .edp file:

To call SuperLU_DIST with internal parameter, we used the parameters sparams. The value of parameters of SuperLU_DIST in sparams are defined by :

• nprow=1,

• npcol=1,

• matrix= distributedgloba,

• Fact= DOFACT,

• Equil=NO,

• ParSymbFact=NO,

• ColPerm= MMD_AT_PLUS_A,

• RowPerm= LargeDiag,

• DiagPivotThresh=1.0,

• IterRefine=DOUBLE,

• Trans=NOTRANS,

• ReplaceTinyPivot=NO,

• SolveInitialized=NO,

• PrintStat=NO,

• DiagScale=NOEQUIL

This value correspond to the parameter in the file ffsuperlu_dist_fileparam.txt. If one parameter is not specified by the user, we take the default value of SuperLU_DIST.

Reading solver parameters on a file: The structure of data file for SuperLU_DIST in FreeFEM is given in the file ffsuperlu_dist_fileparam.txt (default value of the FreeFEM interface).

 11 /* nprow : integer value */
 21 /* npcol : integer value */
 3distributedglobal /* matrix input : assembled, distributedglobal, distributed */
 4DOFACT /* Fact : DOFACT, SamePattern, SamePattern_SameRowPerm, FACTORED */
 5NO /* Equil : NO, YES */
 6NO /* ParSymbFact : NO, YES */
 7MMD_AT_PLUS_A /* ColPerm : NATURAL, MMD_AT_PLUS_A, MMD_ATA, METIS_AT_PLUS_A, PARMETIS, MY_PERMC */
 8LargeDiag /* RowPerm : NOROWPERM, LargeDiag, MY_PERMR */
 91.0 /* DiagPivotThresh : real value */
10DOUBLE /* IterRefine : NOREFINE, SINGLE, DOUBLE, EXTRA */
11NOTRANS /* Trans : NOTRANS, TRANS, CONJ*/
12NO /* ReplaceTinyPivot : NO, YES*/
13NO /* SolveInitialized : NO, YES*/
14NO /* RefineInitialized : NO, YES*/
15NO /* PrintStat : NO, YES*/
16NOEQUIL /* DiagScale : NOEQUIL, ROW, COL, BOTH*/

If no solver parameter is given, we used default option of SuperLU_DIST solver.

PaStiX solver

PaStiX (Parallel Sparse matrix package) is a free scientific library under CECILL-C license. This package solves sparse linear system with a direct and block ILU(k) iterative methods. his solver can be applied to a real or complex matrix with a symmetric pattern.

PaStiX parameters:

The input matrix parameter of FreeFEM depend on PaStiX interface. matrix = assembled for non distributed matrix. It is the same parameter for SuperLU_DIST.

There are four parameters in PaStiX : iparm, dparm, perm and invp. These parameters are respectively the integer parameters (vector of size 64), real parameters (vector of size 64), permutation matrix and inverse permutation matrix respectively. iparm and dparm vectors are described in PaStiX RefCard.

The parameters permr and permc in FreeFEM are provided to give permutation matrix and inverse permutation matrix of the user respectively.

Solver parameters defined in .edp file:

To call PaStiX in FreeFEM in this case, we need to specify the parameters lparams and dparams. These parameters are defined by :

\(\forall i\) = 0,… ,63, lparams[i] = iparm[i].

\(\forall i\) = 0,… ,63, dparams[i] = dparm[i].

Reading solver parameters on a file:

The structure of data file for PaStiX parameters in FreeFEM is: first line structure parameters of the matrix and in the following line the value of vectors iparm and dparm in this order.

 1assembled /* matrix input :: assembled, distributed global and distributed */
 2iparm[0]
 3iparm[1]
 4...
 5...
 6iparm[63]
 7dparm[0]
 8dparm[1]
 9...
10...
11dparm[63]

An example of this file parameter is given in ffpastix_iparm_dparm.txt with a description of these parameters. This file is obtained with the example file iparm.txt and dparm.txt including in the PaStiX package.

If no solver parameter is given, we use the default option of PaStiX solver.

In Table 4, we recall the different matrix considering in the different direct solvers.

pARMS solver

pARMS (parallel Algebraic Multilevel Solver) is a software developed by Youssef Saad and al at University of Minnesota.

This software is specialized in the resolution of large sparse non symmetric linear systems of equation. Solvers developed in pARMS are of type “Krylov’s subspace”.

It consists of variants of GMRES like FGMRES (Flexible GMRES), DGMRES (Deflated GMRES) [SAAD2003] and BICGSTAB. pARMS also implements parallel preconditioner like RAS (Restricted Additive Schwarz) [CAI1989] and Schur Complement type preconditioner.

All these parallel preconditioners are based on the principle of domain decomposition. Thus, the matrix \(A\) is partitioned into sub matrices \(A_i\)(\(i=1,...,p\)) where p represents the number of partitions one needs. The union of \(A_i\) forms the original matrix. The solution of the overall system is obtained by solving the local systems on \(A_i\) (see [SMITH1996]). Therefore, a distinction is made between iterations on \(A\) and the local iterations on \(A_i\).

To solve the local problem on \(A_i\) there are several preconditioners as ilut (Incomplete LU with threshold), iluk (Incomplete LU with level of fill in) and ARMS (Algebraic Recursive Multilevel Solver).

 1load "parms_FreeFem" //Tell FreeFem that you will use pARMS
 2
 3// Mesh
 4border C(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}
 5mesh Th = buildmesh (C(50));
 6
 7// Fespace
 8fespace Vh(Th, P2); Vh u, v;
 9
10// Function
11func f= x*y;
12
13// Problem
14problem Poisson (u, v, solver=sparsesolver)
15    = int2d(Th)(
16          dx(u)*dx(v)
17        + dy(u)*dy(v) )
18    + int2d(Th)(
19        - f*v
20    )
21    + on(1, u=0) ;
22
23// Solve
24real cpu = clock();
25Poisson;
26cout << " CPU time = " << clock()-cpu << endl;
27
28// Plot
29plot(u);

Here are some default parameters:

• solver=FGMRES,

• Krylov dimension=30,

• Maximum of Krylov=1000,

• Tolerance for convergence=1e-08 (see book [SAAD2003] to understand all this parameters),

• preconditionner=Restricted Additif Schwarz [CAI1989],

• Inner Krylov dimension=5,

• Maximum of inner Krylov dimension=5,

• Inner preconditionner=ILUK.

To specify the parameters to apply to the solver, the user can either give an integer vector for integer parameters and real vectors for real parameters or provide a file which contains those parameters.

 1load "parms_FreeFem"
 2
 3// Parameters
 4real nu = 1.;
 5int[int] iparm(16);
 6real[int] dparm(6);
 7for (int ii = 0; ii < 16; ii++)
 8    iparm[ii] = -1;
 9for (int ii = 0; ii < 6; ii++)
10    dparm[ii] = -1.0; iparm[0]=0;
11
12// Mesh
13mesh Th = square(10, 10);
14int[int] wall = [1, 3];
15int inlet = 4;
16
17// Fespace
18fespace Vh(Th, [P2, P2, P1]);
19
20// Function
21func uc = 1.;
22
23// Problem
24varf Stokes ([u, v, p], [ush, vsh, psh], solver=sparsesolver)
25    = int2d(Th)(
26          nu*(
27            dx(u)*dx(ush)
28            + dy(u)*dy(ush)
29            + dx(v)*dx(vsh)
30            + dy(v)*dy(vsh)
31        )
32        - p*psh*(1.e-6)
33        - p*(dx(ush) + dy(vsh))
34        - (dx(u) + dy(v))*psh
35    )
36    + on(wall, wall, u=0., v=0.)
37    + on(inlet, u=uc, v=0) ;
38
39matrix AA = Stokes(Vh, Vh);
40set(AA, solver=sparsesolver, lparams=iparm, dparams=dparm); //set pARMS as linear solver
41real[int] bb = Stokes(0, Vh);
42real[int] sol(AA.n);
43sol = AA^-1 * bb;

Interfacing with HIPS

HIPS (Hierarchical Iterative Parallel Solver) is a scientific library that provides an efficient parallel iterative solver for very large sparse linear systems. HIPS is available as free software under the CeCILL-C licence.

HIPS implements two solver classes which are the iteratives class (GMRES, PCG) and the Direct class. Concerning preconditionners, HIPS implements a type of multilevel ILU. For further informations on those preconditionners see the HIPS documentation.

 1load "msh3"
 2load "hips_FreeFem" //load Hips library
 3
 4// Parameters
 5int nn = 10;
 6real zmin = 0, zmax = 1;
 7int[int] iparm(14);
 8real[int] dparm(6);
 9for (int iii = 0; iii < 14; iii++)
10    iparm[iii] = -1;
11for (int iii = 0; iii < 6; iii++)
12    dparm[iii] = -1;
13iparm[0] = 0; //use iterative solver
14iparm[1] = 1; //PCG as Krylov method
15iparm[4] = 0; //Matrix are symmetric
16iparm[5] = 1; //Pattern are also symmetric
17iparm[9] = 1; //Scale matrix
18dparm[0] = 1e-13; //Tolerance to convergence
19dparm[1] = 5e-4; //Threshold in ILUT
20dparm[2] = 5e-4; //Threshold for Schur preconditionner
21
22// Functions
23func ue = 2*x*x + 3*y*y + 4*z*z + 5*x*y + 6*x*z + 1;
24func uex = 4*x + 5*y + 6*z;
25func uey = 6*y + 5*x;
26func uez = 8*z + 6*x;
27func f = -18.;
28
29// Mesh
30mesh Th2 = square(nn, nn);
31
32int[int] rup = [0,2], rdown=[0, 1];
33int[int] rmid=[1, 1, 2, 1, 3, 1, 4, 1];
34
35mesh3 Th=buildlayers(Th2, nn, zbound=[zmin, zmax], reffacemid=rmid,
36    reffaceup = rup, reffacelow = rdown);
37
38// Fespace
39fespace Vh2(Th2, P2);
40Vh2 ux, uz, p2;
41
42fespace Vh(Th, P2);
43Vh uhe = ue;
44cout << "uhe min =" << uhe[].min << ", max =" << uhe[].max << endl;
45Vh u, v;
46Vh F;
47
48// Macro
49macro Grad3(u) [dx(u), dy(u), dz(u)] //
50
51// Problem
52varf va (u, v)
53    = int3d(Th)(
54          Grad3(v)' * Grad3(u)
55    )
56    + int2d(Th, 2)(
57          u*v
58    )
59    - int3d(Th)(
60          f*v
61    )
62    - int2d(Th, 2)(
63          ue*v + (uex*N.x + uey*N.y + uez*N.z)*v
64    )
65    + on(1, u=ue);
66
67varf l (unused, v) = int3d(Th)(f*v);
68
69real cpu=clock();
70matrix Aa = va(Vh, Vh);
71
72F[] = va(0, Vh);
73
74if (mpirank == 0){
75    cout << "Size of A =" << Aa.n << endl;
76    cout << "Non zero coefficients =" << Aa.nbcoef << endl;
77    cout << "CPU TIME FOR FORMING MATRIX =" << clock()-cpu << endl;
78}
79
80set(Aa, solver=sparsesolver, dparams=dparm, lparams=iparm); //Set hips as linear solver
81
82// Solve
83u[] = Aa^-1*F[];
84
85// Plot
86plot(u);