xref: /petsc/src/mat/tests/ex192.c (revision bef158480efac06de457f7a665168877ab3c2fd7) 
1 
2 static char help[] = "Tests MatSolve() and MatMatSolve() with MUMPS or MKL_PARDISO sequential solvers in Schur complement mode.\n\
3 Example: mpiexec -n 1 ./ex192 -f <matrix binary file> -nrhs 4 -symmetric_solve -hermitian_solve -schur_ratio 0.3\n\n";
4 
5 #include <petscmat.h>
6 
7 int main(int argc,char **args)
8 {
9   Mat            A,RHS,C,F,X,S;
10   Vec            u,x,b;
11   Vec            xschur,bschur,uschur;
12   IS             is_schur;
13   PetscErrorCode ierr;
14   PetscMPIInt    size;
15   PetscInt       isolver=0,size_schur,m,n,nfact,nsolve,nrhs;
16   PetscReal      norm,tol=PETSC_SQRT_MACHINE_EPSILON;
17   PetscRandom    rand;
18   PetscBool      data_provided,herm,symm,use_lu,cuda = PETSC_FALSE;
19   PetscReal      sratio = 5.1/12.;
20   PetscViewer    fd;              /* viewer */
21   char           solver[256];
22   char           file[PETSC_MAX_PATH_LEN]; /* input file name */
23 
24   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
25   ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
26   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor test");
27   /* Determine which type of solver we want to test for */
28   herm = PETSC_FALSE;
29   symm = PETSC_FALSE;
30   ierr = PetscOptionsGetBool(NULL,NULL,"-symmetric_solve",&symm,NULL);CHKERRQ(ierr);
31   ierr = PetscOptionsGetBool(NULL,NULL,"-hermitian_solve",&herm,NULL);CHKERRQ(ierr);
32   if (herm) symm = PETSC_TRUE;
33   ierr = PetscOptionsGetBool(NULL,NULL,"-cuda_solve",&cuda,NULL);CHKERRQ(ierr);
34   ierr = PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL);CHKERRQ(ierr);
35 
36   /* Determine file from which we read the matrix A */
37   ierr = PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&data_provided);CHKERRQ(ierr);
38   if (!data_provided) { /* get matrices from PETSc distribution */
39     ierr = PetscStrncpy(file,"${PETSC_DIR}/share/petsc/datafiles/matrices/",sizeof(file));CHKERRQ(ierr);
40     if (symm) {
41 #if defined (PETSC_USE_COMPLEX)
42       ierr = PetscStrlcat(file,"hpd-complex-",sizeof(file));CHKERRQ(ierr);
43 #else
44       ierr = PetscStrlcat(file,"spd-real-",sizeof(file));CHKERRQ(ierr);
45 #endif
46     } else {
47 #if defined (PETSC_USE_COMPLEX)
48       ierr = PetscStrlcat(file,"nh-complex-",sizeof(file));CHKERRQ(ierr);
49 #else
50       ierr = PetscStrlcat(file,"ns-real-",sizeof(file));CHKERRQ(ierr);
51 #endif
52     }
53 #if defined(PETSC_USE_64BIT_INDICES)
54     ierr = PetscStrlcat(file,"int64-",sizeof(file));CHKERRQ(ierr);
55 #else
56     ierr = PetscStrlcat(file,"int32-",sizeof(file));CHKERRQ(ierr);
57 #endif
58 #if defined (PETSC_USE_REAL_SINGLE)
59     ierr = PetscStrlcat(file,"float32",sizeof(file));CHKERRQ(ierr);
60 #else
61     ierr = PetscStrlcat(file,"float64",sizeof(file));CHKERRQ(ierr);
62 #endif
63   }
64   /* Load matrix A */
65   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);
66   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
67   ierr = MatLoad(A,fd);CHKERRQ(ierr);
68   ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
69   ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);
70   if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
71 
72   /* Create dense matrix C and X; C holds true solution with identical colums */
73   nrhs = 2;
74   ierr = PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);CHKERRQ(ierr);
75   ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);
76   ierr = MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);CHKERRQ(ierr);
77   ierr = MatSetType(C,MATDENSE);CHKERRQ(ierr);
78   ierr = MatSetFromOptions(C);CHKERRQ(ierr);
79   ierr = MatSetUp(C);CHKERRQ(ierr);
80 
81   ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand);CHKERRQ(ierr);
82   ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr);
83   ierr = MatSetRandom(C,rand);CHKERRQ(ierr);
84   ierr = MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);CHKERRQ(ierr);
85 
86   /* Create vectors */
87   ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
88   ierr = VecSetSizes(x,n,PETSC_DECIDE);CHKERRQ(ierr);
89   ierr = VecSetFromOptions(x);CHKERRQ(ierr);
90   ierr = VecDuplicate(x,&b);CHKERRQ(ierr);
91   ierr = VecDuplicate(x,&u);CHKERRQ(ierr); /* save the true solution */
92 
93   ierr = PetscOptionsGetInt(NULL,NULL,"-solver",&isolver,NULL);CHKERRQ(ierr);
94   switch (isolver) {
95 #if defined(PETSC_HAVE_MUMPS)
96     case 0:
97       ierr = PetscStrcpy(solver,MATSOLVERMUMPS);CHKERRQ(ierr);
98       break;
99 #endif
100 #if defined(PETSC_HAVE_MKL_PARDISO)
101     case 1:
102       ierr = PetscStrcpy(solver,MATSOLVERMKL_PARDISO);CHKERRQ(ierr);
103       break;
104 #endif
105     default:
106       ierr = PetscStrcpy(solver,MATSOLVERPETSC);CHKERRQ(ierr);
107       break;
108   }
109 
110 #if defined (PETSC_USE_COMPLEX)
111   if (isolver == 0 && symm && !data_provided) { /* MUMPS (5.0.0) does not have support for hermitian matrices, so make them symmetric */
112     PetscScalar im = PetscSqrtScalar((PetscScalar)-1.);
113     PetscScalar val = -1.0;
114     val = val + im;
115     ierr = MatSetValue(A,1,0,val,INSERT_VALUES);CHKERRQ(ierr);
116     ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
117     ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
118   }
119 #endif
120 
121   ierr = PetscOptionsGetReal(NULL,NULL,"-schur_ratio",&sratio,NULL);CHKERRQ(ierr);
122   if (sratio < 0. || sratio > 1.) {
123     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Invalid ratio for schur degrees of freedom %f", sratio);
124   }
125   size_schur = (PetscInt)(sratio*m);
126 
127   ierr = PetscPrintf(PETSC_COMM_SELF,"Solving with %s: nrhs %D, sym %d, herm %d, size schur %D, size mat %D\n",solver,nrhs,symm,herm,size_schur,m);CHKERRQ(ierr);
128 
129   /* Test LU/Cholesky Factorization */
130   use_lu = PETSC_FALSE;
131   if (!symm) use_lu = PETSC_TRUE;
132 #if defined (PETSC_USE_COMPLEX)
133   if (isolver == 1) use_lu = PETSC_TRUE;
134 #endif
135   if (cuda && symm && !herm) use_lu = PETSC_TRUE;
136 
137   if (herm && !use_lu) { /* test also conversion routines inside the solver packages */
138     ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
139     ierr = MatConvert(A,MATSEQSBAIJ,MAT_INPLACE_MATRIX,&A);CHKERRQ(ierr);
140   }
141 
142   if (use_lu) {
143     ierr = MatGetFactor(A,solver,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
144   } else {
145     if (herm) {
146       ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
147       ierr = MatSetOption(A,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr);
148     } else {
149       ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
150       ierr = MatSetOption(A,MAT_SPD,PETSC_FALSE);CHKERRQ(ierr);
151     }
152     ierr = MatGetFactor(A,solver,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr);
153   }
154   ierr = ISCreateStride(PETSC_COMM_SELF,size_schur,m-size_schur,1,&is_schur);CHKERRQ(ierr);
155   ierr = MatFactorSetSchurIS(F,is_schur);CHKERRQ(ierr);
156 
157   ierr = ISDestroy(&is_schur);CHKERRQ(ierr);
158   if (use_lu) {
159     ierr = MatLUFactorSymbolic(F,A,NULL,NULL,NULL);CHKERRQ(ierr);
160   } else {
161     ierr = MatCholeskyFactorSymbolic(F,A,NULL,NULL);CHKERRQ(ierr);
162   }
163 
164   for (nfact = 0; nfact < 3; nfact++) {
165     Mat AD;
166 
167     if (!nfact) {
168       ierr = VecSetRandom(x,rand);CHKERRQ(ierr);
169       if (symm && herm) {
170         ierr = VecAbs(x);CHKERRQ(ierr);
171       }
172       ierr = MatDiagonalSet(A,x,ADD_VALUES);CHKERRQ(ierr);
173     }
174     if (use_lu) {
175       ierr = MatLUFactorNumeric(F,A,NULL);CHKERRQ(ierr);
176     } else {
177       ierr = MatCholeskyFactorNumeric(F,A,NULL);CHKERRQ(ierr);
178     }
179     if (cuda) {
180       ierr = MatFactorGetSchurComplement(F,&S,NULL);CHKERRQ(ierr);
181       ierr = MatSetType(S,MATSEQDENSECUDA);CHKERRQ(ierr);
182       ierr = MatCreateVecs(S,&xschur,&bschur);CHKERRQ(ierr);
183       ierr = MatFactorRestoreSchurComplement(F,&S,MAT_FACTOR_SCHUR_UNFACTORED);CHKERRQ(ierr);
184     }
185     ierr = MatFactorCreateSchurComplement(F,&S,NULL);CHKERRQ(ierr);
186     if (!cuda) {
187       ierr = MatCreateVecs(S,&xschur,&bschur);CHKERRQ(ierr);
188     }
189     ierr = VecDuplicate(xschur,&uschur);CHKERRQ(ierr);
190     if (nfact == 1 && (!cuda || (herm && symm))) {
191       ierr = MatFactorInvertSchurComplement(F);CHKERRQ(ierr);
192     }
193     for (nsolve = 0; nsolve < 2; nsolve++) {
194       ierr = VecSetRandom(x,rand);CHKERRQ(ierr);
195       ierr = VecCopy(x,u);CHKERRQ(ierr);
196 
197       if (nsolve) {
198         ierr = MatMult(A,x,b);CHKERRQ(ierr);
199         ierr = MatSolve(F,b,x);CHKERRQ(ierr);
200       } else {
201         ierr = MatMultTranspose(A,x,b);CHKERRQ(ierr);
202         ierr = MatSolveTranspose(F,b,x);CHKERRQ(ierr);
203       }
204       /* Check the error */
205       ierr = VecAXPY(u,-1.0,x);CHKERRQ(ierr);  /* u <- (-1.0)x + u */
206       ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr);
207       if (norm > tol) {
208         PetscReal resi;
209         if (nsolve) {
210           ierr = MatMult(A,x,u);CHKERRQ(ierr); /* u = A*x */
211         } else {
212           ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr); /* u = A*x */
213         }
214         ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr);  /* u <- (-1.0)b + u */
215         ierr = VecNorm(u,NORM_2,&resi);CHKERRQ(ierr);
216         if (nsolve) {
217           ierr = PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatSolve error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);CHKERRQ(ierr);
218         } else {
219           ierr = PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatSolveTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);CHKERRQ(ierr);
220         }
221       }
222       ierr = VecSetRandom(xschur,rand);CHKERRQ(ierr);
223       ierr = VecCopy(xschur,uschur);CHKERRQ(ierr);
224       if (nsolve) {
225         ierr = MatMult(S,xschur,bschur);CHKERRQ(ierr);
226         ierr = MatFactorSolveSchurComplement(F,bschur,xschur);CHKERRQ(ierr);
227       } else {
228         ierr = MatMultTranspose(S,xschur,bschur);CHKERRQ(ierr);
229         ierr = MatFactorSolveSchurComplementTranspose(F,bschur,xschur);CHKERRQ(ierr);
230       }
231       /* Check the error */
232       ierr = VecAXPY(uschur,-1.0,xschur);CHKERRQ(ierr);  /* u <- (-1.0)x + u */
233       ierr = VecNorm(uschur,NORM_2,&norm);CHKERRQ(ierr);
234       if (norm > tol) {
235         PetscReal resi;
236         if (nsolve) {
237           ierr = MatMult(S,xschur,uschur);CHKERRQ(ierr); /* u = A*x */
238         } else {
239           ierr = MatMultTranspose(S,xschur,uschur);CHKERRQ(ierr); /* u = A*x */
240         }
241         ierr = VecAXPY(uschur,-1.0,bschur);CHKERRQ(ierr);  /* u <- (-1.0)b + u */
242         ierr = VecNorm(uschur,NORM_2,&resi);CHKERRQ(ierr);
243         if (nsolve) {
244           ierr = PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatFactorSolveSchurComplement error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);CHKERRQ(ierr);
245         } else {
246           ierr = PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatFactorSolveSchurComplementTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);CHKERRQ(ierr);
247         }
248       }
249     }
250     ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&AD);CHKERRQ(ierr);
251     if (!nfact) {
252       ierr = MatMatMult(AD,C,MAT_INITIAL_MATRIX,2.0,&RHS);CHKERRQ(ierr);
253     } else {
254       ierr = MatMatMult(AD,C,MAT_REUSE_MATRIX,2.0,&RHS);CHKERRQ(ierr);
255     }
256     ierr = MatDestroy(&AD);CHKERRQ(ierr);
257     for (nsolve = 0; nsolve < 2; nsolve++) {
258       ierr = MatMatSolve(F,RHS,X);CHKERRQ(ierr);
259 
260       /* Check the error */
261       ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
262       ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr);
263       if (norm > tol) {
264         ierr = PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatMatSolve: Norm of error %g\n",nfact,nsolve,norm);CHKERRQ(ierr);
265       }
266     }
267     if (isolver == 0) {
268       Mat spRHS,spRHST,RHST;
269 
270       ierr = MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);CHKERRQ(ierr);
271       ierr = MatConvert(RHST,MATSEQAIJ,MAT_INITIAL_MATRIX,&spRHST);CHKERRQ(ierr);
272       ierr = MatCreateTranspose(spRHST,&spRHS);CHKERRQ(ierr);
273       for (nsolve = 0; nsolve < 2; nsolve++) {
274         ierr = MatMatSolve(F,spRHS,X);CHKERRQ(ierr);
275 
276         /* Check the error */
277         ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
278         ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr);
279         if (norm > tol) {
280           ierr = PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) sparse MatMatSolve: Norm of error %g\n",nfact,nsolve,norm);CHKERRQ(ierr);
281         }
282       }
283       ierr = MatDestroy(&spRHST);CHKERRQ(ierr);
284       ierr = MatDestroy(&spRHS);CHKERRQ(ierr);
285       ierr = MatDestroy(&RHST);CHKERRQ(ierr);
286     }
287     ierr = MatDestroy(&S);CHKERRQ(ierr);
288     ierr = VecDestroy(&xschur);CHKERRQ(ierr);
289     ierr = VecDestroy(&bschur);CHKERRQ(ierr);
290     ierr = VecDestroy(&uschur);CHKERRQ(ierr);
291   }
292   /* Free data structures */
293   ierr = MatDestroy(&A);CHKERRQ(ierr);
294   ierr = MatDestroy(&C);CHKERRQ(ierr);
295   ierr = MatDestroy(&F);CHKERRQ(ierr);
296   ierr = MatDestroy(&X);CHKERRQ(ierr);
297   ierr = MatDestroy(&RHS);CHKERRQ(ierr);
298   ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr);
299   ierr = VecDestroy(&x);CHKERRQ(ierr);
300   ierr = VecDestroy(&b);CHKERRQ(ierr);
301   ierr = VecDestroy(&u);CHKERRQ(ierr);
302   ierr = PetscFinalize();
303   return ierr;
304 }
305 
306 
307 /*TEST
308 
309    testset:
310      requires: mkl_pardiso double !complex !define(PETSC_USE_64BIT_INDICES)
311      args: -solver 1
312 
313      test:
314        suffix: mkl_pardiso
315      test:
316        requires: cuda
317        suffix: mkl_pardiso_cuda
318        args: -cuda_solve
319        output_file: output/ex192_mkl_pardiso.out
320      test:
321        suffix: mkl_pardiso_1
322        args: -symmetric_solve
323        output_file: output/ex192_mkl_pardiso_1.out
324      test:
325        requires: cuda
326        suffix: mkl_pardiso_cuda_1
327        args: -symmetric_solve -cuda_solve
328        output_file: output/ex192_mkl_pardiso_1.out
329      test:
330        suffix: mkl_pardiso_3
331        args: -symmetric_solve -hermitian_solve
332        output_file: output/ex192_mkl_pardiso_3.out
333      test:
334        requires: cuda
335        suffix: mkl_pardiso_cuda_3
336        args: -symmetric_solve -hermitian_solve -cuda_solve
337        output_file: output/ex192_mkl_pardiso_3.out
338 
339    testset:
340      requires: mumps double !complex
341      args: -solver 0
342 
343      test:
344        suffix: mumps
345      test:
346        requires: cuda
347        suffix: mumps_cuda
348        args: -cuda_solve
349        output_file: output/ex192_mumps.out
350      test:
351        suffix: mumps_2
352        args: -symmetric_solve
353        output_file: output/ex192_mumps_2.out
354      test:
355        requires: cuda
356        suffix: mumps_cuda_2
357        args: -symmetric_solve -cuda_solve
358        output_file: output/ex192_mumps_2.out
359      test:
360        suffix: mumps_3
361        args: -symmetric_solve -hermitian_solve
362        output_file: output/ex192_mumps_3.out
363      test:
364        requires: cuda
365        suffix: mumps_cuda_3
366        args: -symmetric_solve -hermitian_solve -cuda_solve
367        output_file: output/ex192_mumps_3.out
368 
369 TEST*/
370