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