Actual source code: ex263.c
1: static char help[] = "Tests MatForwardSolve and MatBackwardSolve for LU and Cholesky decompositions using C/Pardiso.\n\n";
3: #include <petscdmda.h>
4: #include <petscmat.h>
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;
21: PetscFunctionBeginUser;
22: PetscCall(PetscInitialize(&argc, &args, NULL, help));
23: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
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));
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));
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;
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));
75: if (CHOL) PetscCall(MatSetOption(A, MAT_SPD, PETSC_TRUE));
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));
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: }
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));
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: }
122: /*TEST
124: test:
125: requires: mkl_cpardiso
126: nsize: 4
128: test:
129: suffix: 2
130: requires: !complex mkl_cpardiso
131: nsize: 4
132: args: -chol
134: TEST*/