Parallelization

  1 //usage :
  2 //ff-mpirun [mpi parameter] MPIGMRES2d.edp [-glut ffglut] [-n N] [-k K] [-d D] [-ns] [-gmres [0|1]
  3 //arguments:
  4 //-glut ffglut : to see graphicaly the process
  5 //-n N: set the mesh cube split NxNxN
  6 //-d D: set debug flag D must be one for mpiplot
  7 //-k K: to refined by K all element
  8 //-ns: remove script dump
  9 //-gmres
 10 //0: use iterative schwarz algo.
 11 //1: Algo GMRES on residu of schwarz algo
 12 //2: DDM GMRES
 13 //3: DDM GMRES with coarse grid preconditionner (Good one)
 14 
 15 load "MPICG"
 16 load "medit"
 17 load "metis"
 18 include "getARGV.idp"
 19 include "MPIplot.idp"
 20 include "MPIGMRESmacro.idp"
 21 
 22 searchMethod = 0; //more safe seach algo (warning can be very expensive in case of lot of ouside point)
 23 assert(version >= 3.11); //need at least v3.11
 24 real[int] ttt(10);
 25 int ittt=0;
 26 macro settt {ttt[ittt++] = mpiWtime();}//
 27 
 28 // Arguments
 29 verbosity = getARGV("-vv", 0);
 30 int vdebug = getARGV("-d", 1);
 31 int ksplit = getARGV("-k", 3);
 32 int nloc = getARGV("-n", 10);
 33 string sff = getARGV("-p", "");
 34 int gmres = getARGV("-gmres", 2);
 35 bool dplot = getARGV("-dp", 0);
 36 int nC = getARGV("-N", max(nloc/10, 4));
 37 
 38 if (mpirank==0 && verbosity){
 39     cout << "ARGV: ";
 40     for (int i = 0; i < ARGV.n; ++i)
 41         cout << ARGV[i] << " ";
 42     cout << endl;
 43 }
 44 
 45 if(mpirank==0 && verbosity)
 46 cout << " vdebug: " << vdebug << ", kspilt "<< ksplit << ", nloc "<< nloc << ", sff "<< sff << "." << endl;
 47 
 48 // Parameters
 49 int withplot = 0;
 50 bool withmetis = 1;
 51 bool RAS = 1;
 52 string sPk = "P2-2gd";
 53 func Pk = P2;
 54 int sizeoverlaps = 1; //size of overlap
 55 int[int] l111 = [1, 1, 1, 1]; //mesh labels
 56 
 57 // MPI function
 58 func bool plotMPIall(mesh &Th, real[int] &u, string cm){
 59     if(vdebug)
 60         PLOTMPIALL(mesh, Pk, Th, u, {cmm=cm, nbiso=20, fill=1, dim=3, value=1});
 61     return 1;
 62 }
 63 
 64 // MPI
 65 mpiComm comm(mpiCommWorld,0,0); //trick : make a no split mpiWorld
 66 
 67 int npart = mpiSize(comm); //total number of partion
 68 int ipart = mpiRank(comm); //current partition number
 69 
 70 int njpart = 0; //Number of part with intersection (a jpart) with ipart without ipart
 71 int[int] jpart(npart); //list of jpart
 72 if(ipart==0)
 73     cout << " Final N = " << ksplit*nloc << ", nloc = " << nloc << ", split = " << ksplit << endl;
 74 settt
 75 
 76 // Mesh
 77 mesh Thg = square(nloc, nloc, label=l111);
 78 mesh ThC = square(nC, nC, label=l111);// Coarse mesh
 79 
 80 mesh Thi, Thin; //with overlap, without olverlap
 81 
 82 // Fespace
 83 fespace Phg(Thg, P0);
 84 Phg part;
 85 
 86 fespace Vhg(Thg, P1);
 87 Vhg unssd; //boolean function: 1 in the subdomain, 0 elswhere
 88 
 89 fespace VhC(ThC, P1); // of the coarse problem
 90 
 91 // Partitioning
 92 {
 93     int[int] nupart(Thg.nt);
 94     nupart = 0;
 95     if (npart > 1 && ipart == 0)
 96         metisdual(nupart, Thg, npart);
 97 
 98     broadcast(processor(0, comm), nupart);
 99     for(int i = 0; i < nupart.n; ++i)
100         part[][i] = nupart[i];
101 }
102 
103 if (withplot > 1)
104     plot(part, fill=1, cmm="dual", wait=1);
105 
106 // Overlapping partition
107 Phg suppi = abs(part-ipart) < 0.1;
108 
109 Thin = trunc(Thg, suppi>0, label=10); // non-overlapping mesh, interfaces have label 10
110 int nnn = sizeoverlaps*2;// to be sure
111 AddLayers(Thg, suppi[], nnn, unssd[]); //see above! suppi and unssd are modified
112 unssd[] *= nnn; //to put value nnn a 0
113 real nnn0 = nnn - sizeoverlaps + 0.001;
114 Thi = trunc(Thg, unssd>nnn0, label=10); //overlapping mesh, interfaces have label 10
115 
116 settt
117 
118 // Fespace
119 fespace Vhi(Thi,P1);
120 int npij = npart;
121 Vhi[int] pij(npij); //local partition of unit + pii
122 Vhi pii;
123 
124 real nnn1 = +0.001;
125 {
126     /*
127     construction of the partition of the unit,
128     let phi_i P1 FE function 1 on Thin and zero ouside of Thi and positive
129     the partition is build with
130     p_i = phi_i/ \sum phi_i
131 
132     to build the partition of one domain i
133     we nned to find all j such that supp(phi_j) \cap supp(phi_j) is not empty
134     <=> int phi_j
135     */
136     //build a local mesh of thii such that all computation of the unit partition are
137     //exact in thii
138     mesh Thii = trunc(Thg, unssd>nnn1, label=10); //overlapping mesh, interfaces have label 10
139 
140     {
141         //find all j mes (supp(p_j) cap supp(p_i)) >0
142         //compute all phi_j on Thii
143         //remark: supp p_i include in Thi
144 
145         // Fespace
146         fespace Phii(Thii, P0);
147         fespace Vhii(Thii, P1);
148         Vhi sumphi = 0;
149         Vhii phii = 0;
150 
151         jpart = 0;
152         njpart = 0;
153         int nlayer = RAS ? 1 : sizeoverlaps;
154         if (ipart == 0)
155             cout << "nlayer = " << nlayer << endl;
156         pii = max(unssd-nnn+nlayer, 0.)/nlayer;
157         if(dplot)
158             plot(pii, wait=1, cmm=" 0000");
159         sumphi[] += pii[];
160         if(dplot)
161             plot(sumphi, wait=1, cmm=" summ 0000");
162 
163         real epsmes = 1e-10*Thii.area;
164         for (int i = 0; i < npart; ++i)
165             if (i != ipart){
166             Phii suppii = abs(i-part) < 0.2;
167             if (suppii[].max > 0.5){
168                 AddLayers(Thii, suppii[], nlayer, phii[]);
169                 assert(phii[].min >= 0);
170                 real interij = int2d(Thi)(phii);
171                 if (interij > epsmes){
172                     pij[njpart] = abs(phii);
173                     if(vdebug > 2)
174                         cout << " ***** " << int2d(Thi)(real(pij[njpart])<0) << " " <<pij[njpart][].min << " " << phii[].min << endl;
175                     assert(int2d(Thi)(real(pij[njpart]) < 0) == 0);
176                     if(dplot)
177                         plot(pij[njpart], wait=1, cmm=" j = "+ i + " " + njpart);
178                     sumphi[] += pij[njpart][];
179                     if(dplot)
180                         plot(sumphi, wait=1, cmm=" sum j = "+ i + " " + njpart);
181                     jpart[njpart++] = i;
182                 }
183             }
184         }
185 
186         if(dplot)
187             plot(sumphi, wait=1, dim=3, cmm="sum ", fill=1);
188         pii[] = pii[] ./ sumphi[];
189         for (int j = 0; j < njpart; ++j)
190             pij[j][] = pij[j][] ./ sumphi[];
191         jpart.resize(njpart);
192         for (int j = 0; j < njpart; ++j)
193             assert(pij[j][].max <= 1);
194         {
195             cout << ipart << " number of jpart " << njpart << " : ";
196             for (int j = 0; j < njpart; ++j)
197                 cout << jpart[j] << " ";
198             cout << endl;
199         }
200         sumphi[] = pii[];
201         for (int j = 0; j < njpart; ++j)
202             sumphi[] += pij[j][];
203         if(vdebug > 2)
204             cout << "sum min " << sumphi[].min << " " << sumphi[].max << endl;
205         assert(sumphi[].min > 1.-1e-6 && sumphi[].max < 1.+1e-6);
206     }
207 } //Thii is remove here
208 // end of the construction of the local partition of the unity ...
209 // on Thi
210 if (ipart == 0)
211     cout << "End build partition" << endl;
212 
213 // Computation of number of intersection
214 //here pii and the pij is the local partition of the unit on
215 //Thi (mesh with overlap)
216 if ( dplot){
217     plot(Thi, wait=1);
218     for(int j = 0; j < njpart; ++j)
219         plot(pij[j], cmm=" j="+j, wait=1);
220 }
221 
222 //Partition of the unity on Thi
223 //computation of message
224 //all j > we have to receive
225 //data on intersection of the support of pij[0] and pij[j]
226 settt
227 
228 if(vdebug)
229     plotMPIall(Thi, pii[], "pi_i");
230 
231 mesh[int] aThij(njpart);
232 matrix Pii;
233 matrix[int] sMj(njpart); //M of send to j
234 matrix[int] rMj(njpart); //M to recv from j
235 fespace Whi(Thi, Pk);
236 mesh Thij = Thi;
237 fespace Whij(Thij, Pk);//
238 
239 //construction of the mesh intersect i,j part
240 for(int jp = 0; jp < njpart; ++jp)
241     aThij[jp] = trunc(Thi, pij[jp] > 1e-6, label=10); //mesh of the supp of pij
242 
243 for(int jp = 0; jp < njpart; ++jp)
244     aThij[jp] = trunc(aThij[jp], 1, split=ksplit);
245 
246 Thi = trunc(Thi, 1, split=ksplit);
247 
248 settt
249 
250 if (ipart == 0)
251     cout << "End build mesh intersection" << endl;
252 
253 // Construction of transfert matrix
254 {
255     Whi wpii = pii;
256     Pii = wpii[];
257     for(int jp = 0; jp < njpart; ++jp){
258         int j = jpart[jp];
259         Thij = aThij[jp];
260         matrix I = interpolate(Whij, Whi); //Whji <- Whi
261         sMj[jp] = I*Pii; //Whi -> s Whij
262         rMj[jp] = interpolate(Whij, Whi, t=1); //Whji -> Whi
263         if(vdebug > 10){
264             {Whi uuu=1; Whij vvv=-1; vvv[]+=I*uuu[]; cout << jp << " %%% " << vvv[].linfty << endl; assert(vvv[].linfty < 1e-6);}
265             {Whi uuu=1; Whij vvv=-1; vvv[]+=rMj[jp]'*uuu[]; cout << jp << " ### " << vvv[].linfty << endl; assert(vvv[].linfty < 1e-6);}
266         }
267     }
268 }
269 if (ipart == 0)
270     cout << "End build transfert matrix" << endl;
271 
272 // Allocate array of send and recv data
273 InitU(njpart, Whij, Thij, aThij, Usend) //initU(n, Vh, Th, aTh, U)
274 InitU(njpart, Whij, Thij, aThij, Vrecv)
275 if (ipart == 0)
276     cout << "End init data for send/revc" << endl;
277 
278 Whi ui, vi;
279 
280 func bool Update(real[int] &ui, real[int] &vi){
281     for(int j = 0; j < njpart; ++j)
282         Usend[j][] = sMj[j]*ui;
283     SendRecvUV(comm, jpart, Usend, Vrecv)
284     vi = Pii*ui;
285     for(int j = 0; j < njpart; ++j)
286         vi += rMj[j]*Vrecv[j][];
287     return true;
288 }
289 
290 // Definition of the Problem
291 func G = x*0.1;
292 func F = 1.;
293 macro grad(u) [dx(u),dy(u)] //
294 varf vBC (U, V) = on(1, U=G);
295 varf vPb (U, V) = int2d(Thi)(grad(U)'*grad(V)) + int2d(Thi)(F*V) + on(10, U=0) + on(1, U=G);
296 varf vPbC (U, V) = int2d(ThC)(grad(U)'*grad(V)) + on(1, U=0);
297 varf vPbon (U, V) = on(10, U=1) + on(1, U=1);
298 varf vPbon10only (U, V) = on(10, U=1) + on(1, U=0);
299 //remark the order is important we want 0 part on 10 and 1
300 
301 matrix Ai = vPb(Whi, Whi, solver=sparsesolver);
302 matrix AC, Rci, Pci;
303 
304 if (mpiRank(comm) == 0)
305     AC = vPbC(VhC, VhC, solver=sparsesolver);
306 
307 Pci = interpolate(Whi, VhC);
308 Rci = Pci'*Pii;
309 
310 real[int] onG10 = vPbon10only(0, Whi);
311 real[int] onG = vPbon(0, Whi);
312 real[int] Bi=vPb(0, Whi);
313 
314 int kiter = -1;
315 
316 func bool CoarseSolve(real[int] &V, real[int] &U, mpiComm &comm){
317     //solving the coarse probleme
318     real[int] Uc(Rci.n), Bc(Uc.n);
319     Uc = Rci*U;
320     mpiReduce(Uc, Bc, processor(0, comm), mpiSUM);
321     if (mpiRank(comm) == 0)
322         Uc = AC^-1*Bc;
323     broadcast(processor(0, comm), Uc);
324     V = Pci*Uc;
325 }
326 
327 func real[int] DJ (real[int] &U){
328     ++kiter;
329     real[int] V(U.n);
330     V = Ai*U;
331     V = onG10 ? 0.: V; //remove internal boundary
332     return V;
333 }
334 
335 func real[int] PDJ (real[int] &U){
336     real[int] V(U.n);
337 
338     real[int] b = onG10 ? 0. : U;
339     V = Ai^-1*b;
340     Update(V, U);
341     return U;
342 }
343 
344 func real[int] PDJC (real[int] &U){
345     real[int] V(U.n);
346     CoarseSolve(V, U, comm);
347     V = -V; //-C2*Uo
348     U += Ai*V; //U = (I-A C2) Uo
349     real[int] b = onG10 ? 0. : U;
350     U = Ai^-1*b; // (C1( I -A C2) Uo
351     V = U -V;
352     Update(V, U);
353     return U;
354 }
355 
356 func real[int] DJ0(real[int] &U){
357     ++kiter;
358     real[int] V(U.n);
359     real[int] b = onG .* U;
360     b = onG ? b : Bi ;
361     V = Ai^-1*b;
362     Update(V, U);
363     V -= U;
364     return V;
365 }
366 
367 Whi u = 0, v;
368 { //verification
369     Whi u = 1, v;
370     Update(u[], v[]);
371     u[] -= v[];
372     assert(u[].linfty < 1e-6);
373 }
374 
375 settt
376 u[] = vBC(0, Whi, tgv=1); //set u with tgv BC value
377 
378 real epss = 1e-6;
379 int rgmres = 0;
380 if (gmres == 1){
381     rgmres = MPIAffineGMRES(DJ0, u[], veps=epss, nbiter=300, comm=comm, dimKrylov=100, verbosity=ipart ? 0: 50);
382     real[int] b = onG .* u[];
383     b = onG ? b : Bi;
384     v[] = Ai^-1*b;
385     Update(v[], u[]);
386 }
387 else if (gmres == 2)
388     rgmres = MPILinearGMRES(DJ, precon=PDJ, u[], Bi, veps=epss, nbiter=300, comm=comm, dimKrylov=100, verbosity=ipart ? 0: 50);
389 else if (gmres == 3)
390     rgmres = MPILinearGMRES(DJ, precon=PDJC, u[], Bi, veps=epss, nbiter=300, comm=comm, dimKrylov=100, verbosity=ipart ? 0: 50);
391 else //algo Shwarz for demo
392     for(int iter = 0; iter < 10; ++iter){
393         real[int] b = onG .* u[];
394         b = onG ? b : Bi ;
395         v[] = Ai^-1*b;
396 
397         Update(v[], u[]);
398         if(vdebug)
399             plotMPIall(Thi, u[], "u-"+iter);
400         v[] -= u[];
401 
402         real err = v[].linfty;
403         real umax = u[].max;
404         real[int] aa = [err, umax], bb(2);
405         mpiAllReduce(aa, bb, comm, mpiMAX);
406         real errg = bb[0];
407         real umaxg = bb[1];
408 
409         if (ipart == 0)
410             cout << ipart << " err = " << errg << " u. max " << umaxg << endl;
411         if (errg < 1e-5) break;
412     }
413 
414 if (vdebug)
415     plotMPIall(Thi, u[], "u-final");
416 
417 settt
418 
419 real errg = 1, umaxg;
420 {
421     real umax = u[].max, umaxg;
422     real[int] aa = [umax], bb(1);
423     mpiAllReduce(aa, bb, comm, mpiMAX);
424     errg = bb[0];
425     if (ipart == 0)
426         cout << "umax global = " << bb[0] << " Wtime = " << (ttt[ittt-1]-ttt[ittt-2]) << " s " << " " << kiter << endl;
427 }
428 
429 if (sff != ""){
430     ofstream ff(sff+".txt", append);
431     cout << " ++++ ";
432     cout << mpirank << "/" << mpisize << " k=" << ksplit << " n= " << nloc << " " << sizeoverlaps << " it= " << kiter;
433     for (int i = 1; i < ittt; ++i)
434         cout << " " << ttt[i]-ttt[i-1] << " ";
435     cout << epss << " " << Ai.nbcoef << " " << Ai.n << endl;
436 
437     /*
438     1 mpirank
439     2 mpisize
440     3 ksplit
441     4 nloc
442     5 sizeoverlaps
443     6 kiter
444     7 mesh & part build
445     8 build the partion
446     9 build mesh, transfere , and the fine mesh ..
447     10 build the matrix, the trans matrix, factorizatioon
448     11 GMRES
449     */
450 
451     ff << mpirank << " " << mpisize << " " << sPk << " ";
452     ff << ksplit << " " << nloc << " " << sizeoverlaps << " " << kiter;
453     for (int i = 1; i < ittt; ++i)
454         ff << " " << ttt[i]-ttt[i-1] << " ";
455     ff << epss << " " << Ai.nbcoef << " " << Ai.n << " " << gmres << endl;
456 }

 1 load "MUMPS_FreeFem"
 2 //default solver: real-> MUMPS, complex -> MUMPS
 3 load "real_SuperLU_DIST_FreeFem"
 4 default solver: real-> SuperLU_DIST, complex -> MUMPS
 5 load "real_pastix_FreeFem"
 6 //default solver: real-> pastix, complex -> MUMPS
 7 
 8 // Solving with pastix
 9 {
10     matrix A =
11         [[1, 2, 2, 1, 1],
12         [ 2, 12, 0, 10 , 10],
13         [ 2, 0, 1, 0, 2],
14         [ 1, 10, 0, 22, 0.],
15         [ 1, 10, 2, 0., 22]];
16 
17     real[int] xx = [1, 32, 45, 7, 2], x(5), b(5), di(5);
18     b = A*xx;
19     cout << "b = " << b << endl;
20     cout << "xx = " << xx << endl;
21 
22     set(A, solver=sparsesolver, datafilename="ffpastix_iparm_dparm.txt");
23     cout << "solve" << endl;
24     x = A^-1*b;
25     cout << "b = " << b << endl;
26     cout << "x = " << endl;
27     cout << x << endl;
28     di = xx - x;
29     if (mpirank == 0){
30         cout << "x-xx = " << endl;
31         cout << "Linf = " << di.linfty << ", L2 = " << di.l2 << endl;
32     }
33 }
34 
35 // Solving with SuperLU_DIST
36 realdefaulttoSuperLUdist();
37 //default solver: real-> SuperLU_DIST, complex -> MUMPS
38 {
39     matrix A =
40         [[1, 2, 2, 1, 1],
41         [ 2, 12, 0, 10 , 10],
42         [ 2, 0, 1, 0, 2],
43         [ 1, 10, 0, 22, 0.],
44         [ 1, 10, 2, 0., 22]];
45 
46     real[int] xx = [1, 32, 45, 7, 2], x(5), b(5), di(5);
47     b = A*xx;
48     cout << "b = " << b << endl;
49     cout << "xx = " << xx << endl;
50 
51     set(A, solver=sparsesolver, datafilename="ffsuperlu_dist_fileparam.txt");
52     cout << "solve" << endl;
53     x = A^-1*b;
54     cout << "b = " << b << endl;
55     cout << "x = " << endl;
56     cout << x << endl;
57     di = xx - x;
58     if (mpirank == 0){
59         cout << "x-xx = " << endl;
60         cout << "Linf = " << di.linfty << ", L2 = " << di.l2 << endl;
61     }
62 }
63 
64 // Solving with MUMPS
65 defaulttoMUMPS();
66 //default solver: real-> MUMPS, complex -> MUMPS
67 {
68     matrix A =
69         [[1, 2, 2, 1, 1],
70         [ 2, 12, 0, 10 , 10],
71         [ 2, 0, 1, 0, 2],
72         [ 1, 10, 0, 22, 0.],
73         [ 1, 10, 2, 0., 22]];
74 
75     real[int] xx = [1, 32, 45, 7, 2], x(5), b(5), di(5);
76     b = A*xx;
77     cout << "b = " << b << endl;
78     cout << "xx = " << xx << endl;
79 
80     set(A, solver=sparsesolver, datafilename="ffmumps_fileparam.txt");
81     cout << "solving solution" << endl;
82     x = A^-1*b;
83     cout << "b = " << b << endl;
84     cout << "x = " << endl;
85     cout << x << endl;
86     di = xx - x;
87     if (mpirank == 0){
88         cout << "x-xx = " << endl;
89         cout << "Linf = " << di.linfty << ", L2 " << di.l2 << endl;
90     }
91 }

  1 load "MUMPS_FreeFem"
  2 
  3 // Parameters
  4 int[int] ICNTL(40); //declaration of ICNTL parameter for MUMPS
  5 
  6 //get value of ICNTL from file
  7 if (mpirank == 0){
  8     ifstream ff("ffmumps_fileparam.txt");
  9     string line;
 10     getline(ff, line);
 11     getline(ff, line);
 12     for (int iii = 0; iii < 40; iii++){
 13         ff >> ICNTL[iii];
 14         getline(ff, line);
 15     }
 16 }
 17 
 18 broadcast(processor(0), ICNTL);
 19 
 20 // Given data of MUMPS solver in array lparams(SYM, PAR, ICNTL)
 21 // There is no symmetric storage for a matrix associated with a sparse solver.
 22 // Therefore, the matrix will be considered unsymmetric for parallel sparse solver even if symmetric.
 23 {
 24     // Problem
 25     int SYM = 0;
 26     int PAR = 1;
 27     matrix A =
 28         [
 29             [40, 0, 45, 0, 0],
 30             [0, 12, 0, 0, 0],
 31             [0, 0, 40, 0, 0],
 32             [12, 0, 0, 22, 0],
 33             [0, 0, 20, 0, 22]
 34         ];
 35 
 36     // Construction of integer parameter for MUMPS
 37     int[int] MumpsLParams(42);
 38     MumpsLParams[0] = SYM;
 39     MumpsLParams[1] = PAR;
 40     for (int ii = 0; ii < 40; ii++)
 41         MumpsLParams[ii+2] = ICNTL[ii]; //ICNTL begin with index 0 here
 42 
 43     real[int] xx = [1, 32, 45, 7, 2], x(5), b(5), di(5);
 44     b = A*xx;
 45     if (mpirank == 0)
 46         cout << "xx = " << xx << endl;
 47 
 48     set(A, solver=sparsesolver, lparams=MumpsLParams); //we take the default value for CNTL MUMPS parameter
 49 
 50     // Solve
 51     if (mpirank == 0)
 52         cout << "Solve" << endl;
 53     x = A^-1*b;
 54     if (mpirank == 0)
 55         cout << "b = " << b << endl;
 56     if (mpirank == 0)
 57         cout << "x = " << endl; cout << x << endl;
 58     di = xx-x;
 59     if (mpirank == 0){
 60         cout << "x-xx = " << endl;
 61         cout << "Linf = " << di.linfty << ", L2 = " << di.l2 << endl;
 62     }
 63 }
 64 
 65 // Read parameter of MUMPS solver in file ffmumps_fileparam.txt
 66 {
 67     // Problem
 68     matrix A =
 69         [
 70             [40, 0, 45, 0, 0],
 71             [0, 12, 0, 0 , 0],
 72             [0, 0, 40, 0, 0],
 73             [12, 0, 0, 22, 0],
 74             [0, 0, 20, 0, 22]
 75         ];
 76 
 77     real[int] xx = [1, 32, 45, 7000, 2], x(5), b(5), di(5);
 78     b = A*xx;
 79     if (mpirank == 0){
 80         cout << "b = " << b << endl;
 81         cout << "xx = " << xx << endl;
 82     }
 83 
 84     set(A, solver=sparsesolver, datafilename="ffmumps_fileparam.txt");
 85 
 86     // Solve
 87     if (mpirank == 0)
 88         cout << "Solve" << endl;
 89     x = A^-1*b;
 90 
 91     if (mpirank == 0){
 92         cout << "b = " << b << endl;
 93         cout << "x = " << x << endl;
 94     }
 95     di = xx-x;
 96     if (mpirank == 0){
 97         cout << "x-xx = " << endl;
 98         cout << "Linf = " << di.linfty << ", L2 = " << di.l2 << endl;
 99     }
100 }