1 static char help[] = "Tests MatSolve() and MatMatSolve() (interface to superlu_dist, mumps and mkl_pardiso).\n\
2 Example: mpiexec -n <np> ./ex125 -f <matrix binary file> -nrhs 4 -mat_solver_type <>\n\n";
3
4 /*
5 -mat_solver_type:
6 superlu
7 superlu_dist
8 mumps
9 mkl_pardiso
10 cusparse
11 petsc
12 */
13
14 #include <petscmat.h>
15
CreateRandom(PetscInt n,PetscInt m,Mat * A)16 PetscErrorCode CreateRandom(PetscInt n, PetscInt m, Mat *A)
17 {
18 PetscFunctionBeginUser;
19 PetscCall(MatCreate(PETSC_COMM_WORLD, A));
20 PetscCall(MatSetType(*A, MATAIJ));
21 PetscCall(MatSetFromOptions(*A));
22 PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, n, m));
23 PetscCall(MatSeqAIJSetPreallocation(*A, 5, NULL));
24 PetscCall(MatMPIAIJSetPreallocation(*A, 5, NULL, 5, NULL));
25 PetscCall(MatSetRandom(*A, NULL));
26 PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
27 PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
28 PetscFunctionReturn(PETSC_SUCCESS);
29 }
30
CreateIdentity(PetscInt n,Mat * A)31 PetscErrorCode CreateIdentity(PetscInt n, Mat *A)
32 {
33 PetscFunctionBeginUser;
34 PetscCall(MatCreate(PETSC_COMM_WORLD, A));
35 PetscCall(MatSetType(*A, MATAIJ));
36 PetscCall(MatSetFromOptions(*A));
37 PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, n, n));
38 PetscCall(MatSetUp(*A));
39 PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
40 PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
41 PetscCall(MatShift(*A, 1.0));
42 PetscFunctionReturn(PETSC_SUCCESS);
43 }
44
main(int argc,char ** args)45 int main(int argc, char **args)
46 {
47 Mat A, Ae, RHS = NULL, RHS1 = NULL, C, F, X;
48 Vec u, x, b;
49 PetscMPIInt size;
50 PetscInt m, n, nfact, nsolve, nrhs, ipack = 5;
51 PetscReal norm, tol = 10 * PETSC_SQRT_MACHINE_EPSILON;
52 IS perm = NULL, iperm = NULL;
53 MatFactorInfo info;
54 PetscRandom rand;
55 PetscBool flg, symm, testMatSolve = PETSC_TRUE, testMatMatSolve = PETSC_TRUE, testMatMatSolveTranspose = PETSC_TRUE, testMatSolveTranspose = PETSC_TRUE, match = PETSC_FALSE;
56 PetscBool chol = PETSC_FALSE, view = PETSC_FALSE, matsolvexx = PETSC_FALSE, test_inertia;
57 #if defined(PETSC_HAVE_MUMPS)
58 PetscBool test_mumps_opts = PETSC_FALSE;
59 #endif
60 PetscViewer fd; /* viewer */
61 char file[PETSC_MAX_PATH_LEN]; /* input file name */
62 char pack[PETSC_MAX_PATH_LEN];
63
64 PetscFunctionBeginUser;
65 PetscCall(PetscInitialize(&argc, &args, NULL, help));
66 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
67
68 /* Determine file from which we read the matrix A */
69 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
70 if (flg) { /* Load matrix A */
71 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
72 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
73 PetscCall(MatSetFromOptions(A));
74 PetscCall(MatLoad(A, fd));
75 PetscCall(PetscViewerDestroy(&fd));
76 } else {
77 n = 13;
78 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
79 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
80 PetscCall(MatSetType(A, MATAIJ));
81 PetscCall(MatSetFromOptions(A));
82 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
83 PetscCall(MatSetUp(A));
84 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
85 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
86 PetscCall(MatShift(A, 1.0));
87 }
88
89 /* if A is symmetric, set its flag -- required by MatGetInertia() */
90 PetscCall(MatIsSymmetric(A, 0.0, &symm));
91 PetscCall(MatSetOption(A, MAT_SYMMETRIC, symm));
92
93 test_inertia = symm;
94 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_inertia", &test_inertia, NULL));
95
96 PetscCall(PetscOptionsGetBool(NULL, NULL, "-cholesky", &chol, NULL));
97 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL));
98
99 /* test MATNEST support */
100 flg = PETSC_FALSE;
101 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nest", &flg, NULL));
102 if (flg) {
103 Mat B;
104
105 flg = PETSC_FALSE;
106 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nest_bordered", &flg, NULL));
107 if (!flg) {
108 Mat mats[9] = {NULL, NULL, A, NULL, A, NULL, A, NULL, NULL};
109
110 /* Create a nested matrix representing
111 | 0 0 A |
112 | 0 A 0 |
113 | A 0 0 |
114 */
115 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 3, NULL, 3, NULL, mats, &B));
116 flg = PETSC_TRUE;
117 } else {
118 Mat mats[4];
119
120 /* Create a nested matrix representing
121 | Id R |
122 | R^t A |
123 */
124 PetscCall(MatGetSize(A, NULL, &n));
125 m = n + 12;
126 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
127 PetscCall(CreateIdentity(m, &mats[0]));
128 PetscCall(CreateRandom(m, n, &mats[1]));
129 mats[3] = A;
130
131 /* use CreateTranspose/CreateHermitianTranspose or explicit matrix for debugging purposes */
132 flg = PETSC_FALSE;
133 PetscCall(PetscOptionsGetBool(NULL, NULL, "-expl", &flg, NULL));
134 #if PetscDefined(USE_COMPLEX)
135 if (chol) { /* Hermitian transpose not supported by MUMPS Cholesky factor */
136 if (!flg) PetscCall(MatCreateTranspose(mats[1], &mats[2]));
137 else PetscCall(MatTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2]));
138 flg = PETSC_TRUE;
139 } else {
140 if (!flg) {
141 Mat B;
142
143 PetscCall(MatDuplicate(mats[1], MAT_COPY_VALUES, &B));
144 PetscCall(MatCreateHermitianTranspose(B, &mats[2]));
145 PetscCall(MatDestroy(&B));
146 if (n == m) {
147 PetscCall(MatScale(mats[2], PetscCMPLX(4.0, -2.0)));
148 PetscCall(MatShift(mats[2], PetscCMPLX(-2.0, 1.0))); // mats[2] = (4 - 2i) B* - (2 - i) I
149 PetscCall(MatCreateHermitianTranspose(mats[2], &B));
150 PetscCall(MatDestroy(mats + 2));
151 PetscCall(MatScale(B, 0.5));
152 PetscCall(MatShift(B, PetscCMPLX(1.0, 0.5))); // B = 0.5 mats[2]* - (1 - 0.5i) I = (2 + i) B - (1 + 0.5i) I + (1 + 0.5i) I = (2 + i) B
153 PetscCall(MatCreateHermitianTranspose(B, &mats[2])); // mats[2] = B* = (2 - i) B*
154 PetscCall(MatDestroy(&B));
155 PetscCall(MatScale(mats[1], PetscCMPLX(2.0, 1.0))); // mats[1] = (2 + i) B = mats[2]*
156 } else flg = PETSC_TRUE;
157 } else PetscCall(MatHermitianTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2]));
158 }
159 #else
160 if (!flg) {
161 Mat B;
162
163 PetscCall(MatDuplicate(mats[1], MAT_COPY_VALUES, &B));
164 PetscCall(MatCreateTranspose(B, &mats[2]));
165 PetscCall(MatDestroy(&B));
166 if (n == m) {
167 PetscCall(MatScale(mats[2], 4.0));
168 PetscCall(MatShift(mats[2], -2.0)); // mats[2] = 4 B' - 2 I
169 PetscCall(MatCreateTranspose(mats[2], &B));
170 PetscCall(MatDestroy(mats + 2));
171 PetscCall(MatScale(B, 0.5));
172 PetscCall(MatShift(B, 1.0)); // B = 0.5 mats[2]' + I = 0.5 (4 B' - 2 I)' + I = 2 B
173 PetscCall(MatCreateTranspose(B, &mats[2])); // mats[2] = B' = 2 B'
174 PetscCall(MatDestroy(&B));
175 PetscCall(MatScale(mats[1], 2.0)); // mats[1] = 2 B = mats[2]'
176 } else flg = PETSC_TRUE;
177 } else PetscCall(MatTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2]));
178 #endif
179 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, mats, &B));
180 PetscCall(MatDestroy(&mats[0]));
181 PetscCall(MatDestroy(&mats[1]));
182 PetscCall(MatDestroy(&mats[2]));
183 }
184 PetscCall(MatDestroy(&A));
185 A = B;
186 PetscCall(MatSetOption(A, MAT_SYMMETRIC, symm));
187
188 /* not all the combinations of MatMat operations are supported by MATNEST. */
189 PetscCall(MatComputeOperator(A, MATAIJ, &Ae));
190 } else {
191 PetscCall(PetscObjectReference((PetscObject)A));
192 Ae = A;
193 flg = PETSC_TRUE;
194 }
195 PetscCall(MatGetLocalSize(A, &m, &n));
196 PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n);
197
198 PetscCall(MatViewFromOptions(A, NULL, "-A_view"));
199 PetscCall(MatViewFromOptions(Ae, NULL, "-A_view_expl"));
200
201 /* Create dense matrix C and X; C holds true solution with identical columns */
202 nrhs = 2;
203 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
204 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex125: nrhs %" PetscInt_FMT "\n", nrhs));
205 PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
206 PetscCall(MatSetOptionsPrefix(C, "rhs_"));
207 PetscCall(MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs));
208 PetscCall(MatSetType(C, MATDENSE));
209 PetscCall(MatSetFromOptions(C));
210 PetscCall(MatSetUp(C));
211
212 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_factor", &view, NULL));
213 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatsolve", &testMatMatSolve, NULL));
214 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatsolvetranspose", &testMatMatSolveTranspose, NULL));
215 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matsolvetranspose", &testMatSolveTranspose, NULL));
216 #if defined(PETSC_HAVE_MUMPS)
217 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_mumps_opts", &test_mumps_opts, NULL));
218 #endif
219
220 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
221 PetscCall(PetscRandomSetFromOptions(rand));
222 PetscCall(MatSetRandom(C, rand));
223 PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X));
224
225 /* Create vectors */
226 PetscCall(MatCreateVecs(A, &x, &b));
227 PetscCall(VecDuplicate(x, &u)); /* save the true solution */
228
229 /* Test Factorization */
230 if (flg) PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm)); // TODO FIXME: MatConvert_Nest_AIJ() does not support chained MatCreate[Hermitian]Transpose()
231
232 PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", pack, sizeof(pack), NULL));
233 #if defined(PETSC_HAVE_SUPERLU)
234 PetscCall(PetscStrcmp(MATSOLVERSUPERLU, pack, &match));
235 if (match) {
236 PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!");
237 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU LU:\n"));
238 PetscCall(MatGetFactor(A, MATSOLVERSUPERLU, MAT_FACTOR_LU, &F));
239 matsolvexx = PETSC_FALSE; /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix, need further work */
240 ipack = 0;
241 goto skipoptions;
242 }
243 #endif
244 #if defined(PETSC_HAVE_SUPERLU_DIST)
245 PetscCall(PetscStrcmp(MATSOLVERSUPERLU_DIST, pack, &match));
246 if (match) {
247 PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!");
248 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU_DIST LU:\n"));
249 PetscCall(MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_LU, &F));
250 matsolvexx = PETSC_TRUE;
251 if (symm) { /* A is symmetric */
252 testMatMatSolveTranspose = PETSC_TRUE;
253 testMatSolveTranspose = PETSC_TRUE;
254 } else { /* superlu_dist does not support solving A^t x = rhs yet */
255 testMatMatSolveTranspose = PETSC_FALSE;
256 testMatSolveTranspose = PETSC_FALSE;
257 }
258 ipack = 1;
259 goto skipoptions;
260 }
261 #endif
262 #if defined(PETSC_HAVE_MUMPS)
263 PetscCall(PetscStrcmp(MATSOLVERMUMPS, pack, &match));
264 if (match) {
265 if (chol) {
266 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS CHOLESKY:\n"));
267 PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_CHOLESKY, &F));
268 } else {
269 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS LU:\n"));
270 PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_LU, &F));
271 }
272 matsolvexx = PETSC_TRUE;
273 if (test_mumps_opts) {
274 /* test mumps options */
275 PetscInt icntl;
276 PetscReal cntl;
277
278 icntl = 2; /* sequential matrix ordering */
279 PetscCall(MatMumpsSetIcntl(F, 7, icntl));
280
281 cntl = 1.e-6; /* threshold for row pivot detection */
282 PetscCall(MatMumpsSetIcntl(F, 24, 1));
283 PetscCall(MatMumpsSetCntl(F, 3, cntl));
284 }
285 ipack = 2;
286 goto skipoptions;
287 }
288 #endif
289 #if defined(PETSC_HAVE_MKL_PARDISO)
290 PetscCall(PetscStrcmp(MATSOLVERMKL_PARDISO, pack, &match));
291 if (match) {
292 if (chol) {
293 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO CHOLESKY:\n"));
294 PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_CHOLESKY, &F));
295 } else {
296 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO LU:\n"));
297 PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_LU, &F));
298 }
299 ipack = 3;
300 goto skipoptions;
301 }
302 #endif
303 #if defined(PETSC_HAVE_CUDA)
304 PetscCall(PetscStrcmp(MATSOLVERCUSPARSE, pack, &match));
305 if (match) {
306 if (chol) {
307 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE CHOLESKY:\n"));
308 PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_CHOLESKY, &F));
309 } else {
310 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE LU:\n"));
311 PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_LU, &F));
312 }
313 testMatSolveTranspose = PETSC_FALSE;
314 testMatMatSolveTranspose = PETSC_FALSE;
315 ipack = 4;
316 goto skipoptions;
317 }
318 #endif
319 /* PETSc */
320 match = PETSC_TRUE;
321 if (match) {
322 if (chol) {
323 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC CHOLESKY:\n"));
324 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &F));
325 } else {
326 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC LU:\n"));
327 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F));
328 }
329 matsolvexx = PETSC_TRUE;
330 ipack = 5;
331 goto skipoptions;
332 }
333
334 skipoptions:
335 PetscCall(MatFactorInfoInitialize(&info));
336 info.fill = 5.0;
337 info.shifttype = (PetscReal)MAT_SHIFT_NONE;
338 if (chol) {
339 PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
340 } else {
341 PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
342 }
343
344 for (nfact = 0; nfact < 2; nfact++) {
345 if (chol) {
346 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the CHOLESKY numfactorization \n", nfact));
347 PetscCall(MatCholeskyFactorNumeric(F, A, &info));
348 } else {
349 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the LU numfactorization \n", nfact));
350 PetscCall(MatLUFactorNumeric(F, A, &info));
351 }
352 if (view) {
353 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO));
354 PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD));
355 PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
356 view = PETSC_FALSE;
357 }
358
359 #if defined(PETSC_HAVE_SUPERLU_DIST)
360 if (ipack == 1) { /* Test MatSuperluDistGetDiagU()
361 -- input: matrix factor F; output: main diagonal of matrix U on all processes */
362 PetscInt M;
363 PetscScalar *diag;
364 #if !defined(PETSC_USE_COMPLEX)
365 PetscInt nneg, nzero, npos;
366 #endif
367
368 PetscCall(MatGetSize(F, &M, NULL));
369 PetscCall(PetscMalloc1(M, &diag));
370 PetscCall(MatSuperluDistGetDiagU(F, diag));
371 PetscCall(PetscFree(diag));
372
373 #if !defined(PETSC_USE_COMPLEX)
374 /* Test MatGetInertia() */
375 if (test_inertia) { /* A is symmetric */
376 PetscCall(MatGetInertia(F, &nneg, &nzero, &npos));
377 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos));
378 }
379 #endif
380 }
381 #endif
382
383 #if defined(PETSC_HAVE_MUMPS)
384 /* mumps interface allows repeated call of MatCholeskyFactorSymbolic(), while the succession calls do nothing */
385 if (ipack == 2) {
386 if (chol) {
387 PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
388 PetscCall(MatCholeskyFactorNumeric(F, A, &info));
389 } else {
390 PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
391 PetscCall(MatLUFactorNumeric(F, A, &info));
392 }
393 }
394 #endif
395
396 /* Test MatMatSolve(), A X = B, where B can be dense or sparse */
397 if (testMatMatSolve) {
398 if (!nfact) {
399 PetscCall(MatMatMult(Ae, C, MAT_INITIAL_MATRIX, 2.0, &RHS));
400 } else {
401 PetscCall(MatMatMult(Ae, C, MAT_REUSE_MATRIX, 2.0, &RHS));
402 }
403 for (nsolve = 0; nsolve < 2; nsolve++) {
404 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatMatSolve \n", nsolve));
405 PetscCall(MatMatSolve(F, RHS, X));
406
407 /* Check the error */
408 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
409 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
410 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
411 }
412
413 if (matsolvexx) {
414 /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */
415 PetscCall(MatCopy(RHS, X, SAME_NONZERO_PATTERN));
416 PetscCall(MatMatSolve(F, X, X));
417 /* Check the error */
418 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
419 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
420 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve(F,RHS,RHS): Norm of error %g\n", (double)norm));
421 }
422
423 if (ipack == 2 && size == 1) {
424 Mat spRHS, spRHST, RHST;
425
426 PetscCall(MatTranspose(RHS, MAT_INITIAL_MATRIX, &RHST));
427 PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST));
428 PetscCall(MatCreateTranspose(spRHST, &spRHS));
429 for (nsolve = 0; nsolve < 2; nsolve++) {
430 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the sparse MatMatSolve \n", nsolve));
431 PetscCall(MatMatSolve(F, spRHS, X));
432
433 /* Check the error */
434 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
435 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
436 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the sparse MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
437 }
438 PetscCall(MatDestroy(&spRHST));
439 PetscCall(MatDestroy(&spRHS));
440 PetscCall(MatDestroy(&RHST));
441 }
442 }
443
444 /* Test testMatMatSolveTranspose(), A^T X = B, where B can be dense or sparse */
445 if (testMatMatSolveTranspose) {
446 if (!nfact) {
447 PetscCall(MatTransposeMatMult(Ae, C, MAT_INITIAL_MATRIX, 2.0, &RHS1));
448 } else {
449 PetscCall(MatTransposeMatMult(Ae, C, MAT_REUSE_MATRIX, 2.0, &RHS1));
450 }
451
452 for (nsolve = 0; nsolve < 2; nsolve++) {
453 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatMatSolveTranspose\n", nsolve));
454 PetscCall(MatMatSolveTranspose(F, RHS1, X));
455
456 /* Check the error */
457 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
458 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
459 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
460 }
461
462 if (ipack == 2 && size == 1) {
463 Mat spRHS, spRHST, RHST;
464
465 PetscCall(MatTranspose(RHS1, MAT_INITIAL_MATRIX, &RHST));
466 PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST));
467 PetscCall(MatCreateTranspose(spRHST, &spRHS));
468 for (nsolve = 0; nsolve < 2; nsolve++) {
469 PetscCall(MatMatSolveTranspose(F, spRHS, X));
470
471 /* Check the error */
472 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
473 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
474 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the sparse MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
475 }
476 PetscCall(MatDestroy(&spRHST));
477 PetscCall(MatDestroy(&spRHS));
478 PetscCall(MatDestroy(&RHST));
479 }
480 }
481
482 /* Test MatSolve() */
483 if (testMatSolve) {
484 for (nsolve = 0; nsolve < 2; nsolve++) {
485 PetscCall(VecSetRandom(x, rand));
486 PetscCall(VecCopy(x, u));
487 PetscCall(MatMult(Ae, x, b));
488
489 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatSolve \n", nsolve));
490 PetscCall(MatSolve(F, b, x));
491
492 /* Check the error */
493 PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
494 PetscCall(VecNorm(u, NORM_2, &norm));
495 if (norm > tol) {
496 PetscReal resi;
497 PetscCall(MatMult(Ae, x, u)); /* u = A*x */
498 PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
499 PetscCall(VecNorm(u, NORM_2, &resi));
500 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolve: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact));
501 }
502 }
503 }
504
505 /* Test MatSolveTranspose() */
506 if (testMatSolveTranspose) {
507 for (nsolve = 0; nsolve < 2; nsolve++) {
508 PetscCall(VecSetRandom(x, rand));
509 PetscCall(VecCopy(x, u));
510 PetscCall(MatMultTranspose(Ae, x, b));
511
512 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatSolveTranspose\n", nsolve));
513 PetscCall(MatSolveTranspose(F, b, x));
514
515 /* Check the error */
516 PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
517 PetscCall(VecNorm(u, NORM_2, &norm));
518 if (norm > tol) {
519 PetscReal resi;
520 PetscCall(MatMultTranspose(Ae, x, u)); /* u = A*x */
521 PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
522 PetscCall(VecNorm(u, NORM_2, &resi));
523 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveTranspose: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact));
524 }
525 }
526 }
527 }
528
529 /* Free data structures */
530 PetscCall(MatDestroy(&Ae));
531 PetscCall(MatDestroy(&A));
532 PetscCall(MatDestroy(&C));
533 PetscCall(MatDestroy(&F));
534 PetscCall(MatDestroy(&X));
535 PetscCall(MatDestroy(&RHS));
536 PetscCall(MatDestroy(&RHS1));
537
538 PetscCall(PetscRandomDestroy(&rand));
539 PetscCall(ISDestroy(&perm));
540 PetscCall(ISDestroy(&iperm));
541 PetscCall(VecDestroy(&x));
542 PetscCall(VecDestroy(&b));
543 PetscCall(VecDestroy(&u));
544 PetscCall(PetscFinalize());
545 return 0;
546 }
547
548 /*TEST
549
550 test:
551 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
552 args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type petsc
553 output_file: output/ex125.out
554
555 test:
556 suffix: 2
557 args: -mat_solver_type petsc
558 output_file: output/ex125.out
559
560 test:
561 suffix: mkl_pardiso
562 requires: mkl_pardiso datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
563 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mkl_pardiso
564
565 test:
566 suffix: mkl_pardiso_2
567 requires: mkl_pardiso
568 args: -mat_solver_type mkl_pardiso
569 output_file: output/ex125_mkl_pardiso.out
570
571 testset:
572 requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
573 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps
574 output_file: output/ex125_mumps_seq.out
575
576 test:
577 requires: defined(PETSC_HAVE_MUMPS_MIXED_PRECISION)
578 suffix: mumps_single
579 args: -pc_precision single -tol 1e-5
580 test:
581 suffix: mumps_double
582 args: -pc_precision double
583
584 testset:
585 requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
586 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
587 output_file: output/ex125_mumps_seq.out
588
589 test:
590 requires: defined(PETSC_HAVE_MUMPS_MIXED_PRECISION)
591 suffix: mumps_nest_single
592 args: -pc_precision single -tol 1e-4
593 test:
594 suffix: mumps_nest_double
595 args: -pc_precision double
596
597 testset:
598 nsize: 3
599 requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
600 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps
601 output_file: output/ex125_mumps_par.out
602
603 test:
604 requires: defined(PETSC_HAVE_MUMPS_MIXED_PRECISION)
605 suffix: mumps_2_single
606 args: -pc_precision single -tol 1e-5
607 test:
608 suffix: mumps_2_double
609 args: -pc_precision double
610
611 test:
612 suffix: mumps_2_nest
613 nsize: 3
614 requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
615 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
616 output_file: output/ex125_mumps_par.out
617
618 test:
619 suffix: mumps_3
620 requires: mumps
621 args: -mat_solver_type mumps
622 output_file: output/ex125_mumps_seq.out
623
624 testset:
625 requires: mumps
626 args: -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
627 output_file: output/ex125_mumps_seq.out
628
629 test:
630 requires: !__float128
631 suffix: mumps_3_nest
632 test:
633 suffix: mumps_3_nest_fp128
634 requires: __float128
635 args: -tol 1e-8
636
637 test:
638 suffix: mumps_4
639 nsize: 3
640 requires: mumps
641 args: -mat_solver_type mumps
642 output_file: output/ex125_mumps_par.out
643
644 testset:
645 nsize: 3
646 requires: mumps
647 args: -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
648 output_file: output/ex125_mumps_par.out
649
650 test:
651 requires: !__float128
652 suffix: mumps_4_nest
653 test:
654 suffix: mumps_4_nest_fp128
655 requires: __float128
656 args: -tol 1e-8
657
658 test:
659 suffix: mumps_5
660 nsize: 3
661 requires: mumps
662 args: -mat_solver_type mumps -cholesky
663 output_file: output/ex125_mumps_par_cholesky.out
664
665 testset:
666 nsize: 3
667 requires: mumps
668 args: -mat_solver_type mumps -cholesky -test_nest -test_nest_bordered {{0 1}}
669 output_file: output/ex125_mumps_par_cholesky.out
670
671 test:
672 requires: !__float128
673 suffix: mumps_5_nest
674 test:
675 suffix: mumps_5_nest_fp128
676 requires: __float128
677 args: -tol 1e-8
678
679 test:
680 nsize: 2
681 requires: mumps
682 args: -mat_solver_type mumps -test_nest -test_nest_bordered -m 13 -n 13
683 output_file: output/ex125_mumps_par.out
684
685 test:
686 requires: !__float128
687 suffix: mumps_6
688 test:
689 suffix: mumps_6_fp128
690 requires: __float128
691 args: -tol 1e-8
692
693 test:
694 suffix: superlu
695 requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu
696 args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu
697 output_file: output/ex125_superlu.out
698
699 test:
700 suffix: superlu_dist
701 nsize: {{1 3}}
702 requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
703 args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
704 output_file: output/ex125_superlu_dist.out
705
706 test:
707 suffix: superlu_dist_2
708 nsize: {{1 3}}
709 requires: superlu_dist !complex
710 args: -n 36 -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
711 output_file: output/ex125_superlu_dist.out
712
713 test:
714 suffix: superlu_dist_3
715 nsize: {{1 3}}
716 requires: superlu_dist !complex
717 requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
718 args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
719 output_file: output/ex125_superlu_dist_nonsymmetric.out
720
721 test:
722 suffix: superlu_dist_complex
723 nsize: 3
724 requires: datafilespath double superlu_dist complex !defined(PETSC_USE_64BIT_INDICES)
725 args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type superlu_dist
726 output_file: output/ex125_superlu_dist_complex.out
727
728 test:
729 suffix: superlu_dist_complex_2
730 nsize: 3
731 requires: superlu_dist complex
732 args: -mat_solver_type superlu_dist
733 output_file: output/ex125_superlu_dist_complex_2.out
734
735 test:
736 suffix: cusparse
737 requires: cuda datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
738 #TODO: fix the bug with cholesky
739 #args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0 1}separate output}
740 args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0}separate output}
741
742 test:
743 suffix: cusparse_2
744 requires: cuda
745 args: -mat_type aijcusparse -mat_solver_type cusparse -cholesky {{0 1}separate output}
746
747 testset:
748 nsize: {{1 2}separate output}
749 requires: double !defined(PETSC_USE_64BIT_INDICES) datafilespath !complex
750 args: -f ${DATAFILESPATH}/matrices/mixed_poisson
751 test:
752 requires: superlu_dist TODO # superlu_dist is broken
753 suffix: saddle_point_superlu_dist
754 args: -mat_solver_type superlu_dist -mat_superlu_dist_rowperm {{norowperm largediag_mc64}} -test_inertia 0
755 test:
756 requires: mumps
757 suffix: saddle_point_mumps_lu
758 args: -mat_solver_type mumps -mat_mumps_icntl_14 100
759 test:
760 requires: mumps
761 suffix: saddle_point_mumps_cholesky
762 args: -cholesky -mat_solver_type mumps
763
764 TEST*/
765