1 static char help[] = "Tests MatForwardSolve and MatBackwardSolve for LU and Cholesky decompositions using C/Pardiso.\n\n";
2
3 #include <petscdmda.h>
4 #include <petscmat.h>
5
main(int argc,char ** args)6 int main(int argc, char **args)
7 {
8 DM da;
9 DMDALocalInfo info;
10 Mat A, F;
11 Vec x, y, b, ytmp;
12 IS rowis, colis;
13 PetscInt i, j, k, n = 5;
14 PetscBool CHOL = PETSC_FALSE;
15 MatStencil row, cols[5];
16 PetscScalar vals[5];
17 PetscReal norm2, tol = 100. * PETSC_MACHINE_EPSILON;
18 PetscRandom rdm;
19 MatFactorInfo finfo;
20
21 PetscFunctionBeginUser;
22 PetscCall(PetscInitialize(&argc, &args, NULL, help));
23 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
24
25 PetscCall(DMDACreate2d(MPI_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, n, n, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
26 PetscCall(DMSetUp(da));
27
28 PetscCall(DMCreateMatrix(da, &A));
29 PetscCall(PetscOptionsHasName(NULL, NULL, "-chol", &CHOL));
30 if (CHOL) PetscCall(MatSetType(A, MATSBAIJ));
31 else PetscCall(MatSetType(A, MATBAIJ));
32 PetscCall(MatSetUp(A));
33
34 PetscCall(DMDAGetLocalInfo(da, &info));
35 for (j = info.ys; j < info.ys + info.ym; j++) {
36 for (i = info.xs; i < info.xs + info.xm; i++) {
37 row.j = j;
38 row.i = i;
39
40 k = 0;
41 if (j != 0) {
42 cols[k].j = j - 1;
43 cols[k].i = i;
44 vals[k] = -1;
45 ++k;
46 }
47 if (i != 0) {
48 cols[k].j = j;
49 cols[k].i = i - 1;
50 vals[k] = -1;
51 ++k;
52 }
53 cols[k].j = j;
54 cols[k].i = i;
55 vals[k] = 4;
56 ++k;
57 if (j != info.my - 1) {
58 cols[k].j = j + 1;
59 cols[k].i = i;
60 vals[k] = -1;
61 ++k;
62 }
63 if (i != info.mx - 1) {
64 cols[k].j = j;
65 cols[k].i = i + 1;
66 vals[k] = -1;
67 ++k;
68 }
69 PetscCall(MatSetValuesStencil(A, 1, &row, k, cols, vals, INSERT_VALUES));
70 }
71 }
72 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
73 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
74
75 if (CHOL) PetscCall(MatSetOption(A, MAT_SPD, PETSC_TRUE));
76
77 /* Create vectors for error checking */
78 PetscCall(MatCreateVecs(A, &x, &b));
79 PetscCall(VecDuplicate(x, &y));
80 PetscCall(VecDuplicate(x, &ytmp));
81 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
82 PetscCall(PetscRandomSetFromOptions(rdm));
83 PetscCall(VecSetRandom(x, rdm));
84 PetscCall(MatMult(A, x, b));
85
86 PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &rowis, &colis));
87 if (CHOL) {
88 PetscCall(PetscPrintf(MPI_COMM_WORLD, "Test Cholesky...\n"));
89 PetscCall(MatGetFactor(A, MATSOLVERMKL_CPARDISO, MAT_FACTOR_CHOLESKY, &F));
90 PetscCall(MatCholeskyFactorSymbolic(F, A, rowis, &finfo));
91 PetscCall(MatCholeskyFactorNumeric(F, A, &finfo));
92 } else {
93 PetscCall(PetscPrintf(MPI_COMM_WORLD, "Test LU...\n"));
94 PetscCall(MatGetFactor(A, MATSOLVERMKL_CPARDISO, MAT_FACTOR_LU, &F));
95 PetscCall(MatLUFactorSymbolic(F, A, rowis, colis, &finfo));
96 PetscCall(MatLUFactorNumeric(F, A, &finfo));
97 }
98
99 PetscCall(PetscPrintf(MPI_COMM_WORLD, "Test MatForwardSolve...\n"));
100 PetscCall(MatForwardSolve(F, b, ytmp));
101 PetscCall(PetscPrintf(MPI_COMM_WORLD, "Test MatBackwardSolve...\n"));
102 PetscCall(MatBackwardSolve(F, ytmp, y));
103 PetscCall(VecAXPY(y, -1.0, x));
104 PetscCall(VecNorm(y, NORM_2, &norm2));
105 if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2));
106
107 /* Free data structures */
108 PetscCall(ISDestroy(&rowis));
109 PetscCall(ISDestroy(&colis));
110 PetscCall(MatDestroy(&F));
111 PetscCall(MatDestroy(&A));
112 PetscCall(DMDestroy(&da));
113 PetscCall(PetscRandomDestroy(&rdm));
114 PetscCall(VecDestroy(&x));
115 PetscCall(VecDestroy(&y));
116 PetscCall(VecDestroy(&ytmp));
117 PetscCall(VecDestroy(&b));
118 PetscCall(PetscFinalize());
119 return 0;
120 }
121
122 /*TEST
123
124 test:
125 requires: mkl_cpardiso
126 nsize: 4
127
128 test:
129 suffix: 2
130 requires: !complex mkl_cpardiso
131 nsize: 4
132 args: -chol
133
134 TEST*/
135