Actual source code: ex129.c
1: /*
2: Laplacian in 3D. Use for testing MatSolve routines.
3: Modeled by the partial differential equation
5: - Laplacian u = 1,0 < x,y,z < 1,
7: with boundary conditions
8: u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
9: */
11: static char help[] = "This example is for testing different MatSolve routines :MatSolve(), MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd(), MatMatSolve(), and MatMatSolveTranspose(), including how they flag a solution they could not compute.\n\
12: Example usage: ./ex129 -mat_type aij -dof 2\n\n";
14: #include <petscdm.h>
15: #include <petscdmda.h>
17: extern PetscErrorCode ComputeMatrix(DM, Mat);
18: extern PetscErrorCode ComputeRHS(DM, Vec);
19: extern PetscErrorCode ComputeRHSMatrix(PetscInt, PetscInt, Mat *);
21: /*
22: Every entry of a solution that could not be computed must have been flagged with positive infinity, see MatFlag() and VecFlag()
23: */
24: static PetscErrorCode CheckInf(Mat X, const char name[])
25: {
26: const PetscScalar *x;
27: PetscInt i, j, m, N, lda;
29: PetscFunctionBeginUser;
30: PetscCall(MatGetLocalSize(X, &m, NULL));
31: PetscCall(MatGetSize(X, NULL, &N));
32: PetscCall(MatDenseGetLDA(X, &lda));
33: PetscCall(MatDenseGetArrayRead(X, &x));
34: for (j = 0; j < N; j++)
35: for (i = 0; i < m; i++)
36: PetscCheck(PetscIsInfReal(PetscRealPart(x[i + j * lda])) && PetscRealPart(x[i + j * lda]) > 0.0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s left entry (%" PetscInt_FMT ",%" PetscInt_FMT ") unflagged after a failed factorization", name, i, j);
37: PetscCall(MatDenseRestoreArrayRead(X, &x));
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: /*
42: A factorization that hit a zero pivot must make MatMatSolve() and MatMatSolveTranspose() flag every column of X, exactly as MatSolve() flags x
43: */
44: static PetscErrorCode TestZeroPivot(PetscInt n, PetscInt nrhs)
45: {
46: Mat A, F, B, X;
47: Vec b, x;
48: IS perm, iperm;
49: MatFactorInfo info;
50: const PetscScalar *xx;
51: PetscInt i;
53: PetscFunctionBeginUser;
54: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 1, NULL, &A));
55: /* an explicit zero on the diagonal gives a numeric, rather than structural, zero pivot */
56: for (i = 0; i < n; i++) PetscCall(MatSetValue(A, i, i, i == n / 2 ? 0.0 : 1.0, INSERT_VALUES));
57: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
58: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
60: PetscCall(MatFactorInfoInitialize(&info));
61: PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &perm, &iperm));
62: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F));
63: PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
64: /* the zero pivot is divided by during the factorization, so do not trap floating point exceptions until the solves are done */
65: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
66: PetscCall(MatLUFactorNumeric(F, A, &info));
68: PetscCall(MatCreateVecs(A, &x, &b));
69: PetscCall(VecSet(b, 1.0));
70: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, n, nrhs, NULL, &B));
71: PetscCall(MatZeroEntries(B));
72: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &X)); /* start from zeros so that an untouched X cannot pass the checks below */
74: /* the reference behaviour that the block solves must reproduce */
75: PetscCall(MatSolve(F, b, x));
76: PetscCall(VecGetArrayRead(x, &xx));
77: for (i = 0; i < n; i++) PetscCheck(PetscIsInfReal(PetscRealPart(xx[i])) && PetscRealPart(xx[i]) > 0.0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatSolve() left entry %" PetscInt_FMT " unflagged after a failed factorization", i);
78: PetscCall(VecRestoreArrayRead(x, &xx));
80: PetscCall(MatMatSolve(F, B, X));
81: PetscCall(CheckInf(X, "MatMatSolve()"));
83: PetscCall(MatZeroEntries(X));
84: PetscCall(MatMatSolveTranspose(F, B, X));
85: PetscCall(CheckInf(X, "MatMatSolveTranspose()"));
86: PetscCall(PetscFPTrapPop());
88: PetscCall(MatDestroy(&A));
89: PetscCall(MatDestroy(&F));
90: PetscCall(MatDestroy(&B));
91: PetscCall(MatDestroy(&X));
92: PetscCall(VecDestroy(&x));
93: PetscCall(VecDestroy(&b));
94: PetscCall(ISDestroy(&perm));
95: PetscCall(ISDestroy(&iperm));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: int main(int argc, char **args)
100: {
101: PetscMPIInt size;
102: Vec x, b, y, b1;
103: DM da;
104: Mat A, F, RHS, X, C1;
105: MatFactorInfo info;
106: IS perm, iperm;
107: PetscInt dof = 1, M = 8, m, n, nrhs;
108: PetscScalar one = 1.0;
109: PetscReal norm, tol = 1000 * PETSC_MACHINE_EPSILON;
110: PetscBool InplaceLU = PETSC_FALSE;
112: PetscFunctionBeginUser;
113: PetscCall(PetscInitialize(&argc, &args, NULL, help));
114: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
115: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only");
116: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL));
117: PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
119: PetscCall(TestZeroPivot(4, 2));
121: PetscCall(DMDACreate(PETSC_COMM_WORLD, &da));
122: PetscCall(DMSetDimension(da, 3));
123: PetscCall(DMDASetBoundaryType(da, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE));
124: PetscCall(DMDASetStencilType(da, DMDA_STENCIL_STAR));
125: PetscCall(DMDASetSizes(da, M, M, M));
126: PetscCall(DMDASetNumProcs(da, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
127: PetscCall(DMDASetDof(da, dof));
128: PetscCall(DMDASetStencilWidth(da, 1));
129: PetscCall(DMDASetOwnershipRanges(da, NULL, NULL, NULL));
130: PetscCall(DMSetMatType(da, MATBAIJ));
131: PetscCall(DMSetFromOptions(da));
132: PetscCall(DMSetUp(da));
134: PetscCall(DMCreateGlobalVector(da, &x));
135: PetscCall(DMCreateGlobalVector(da, &b));
136: PetscCall(VecDuplicate(b, &y));
137: PetscCall(ComputeRHS(da, b));
138: PetscCall(VecSet(y, one));
139: PetscCall(DMCreateMatrix(da, &A));
140: PetscCall(ComputeMatrix(da, A));
141: PetscCall(MatGetSize(A, &m, &n));
142: nrhs = 2;
143: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
144: PetscCall(ComputeRHSMatrix(m, nrhs, &RHS));
145: PetscCall(MatDuplicate(RHS, MAT_DO_NOT_COPY_VALUES, &X));
147: PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm));
149: PetscCall(PetscOptionsGetBool(NULL, NULL, "-inplacelu", &InplaceLU, NULL));
150: PetscCall(MatFactorInfoInitialize(&info));
151: if (!InplaceLU) {
152: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F));
153: info.fill = 5.0;
154: PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
155: PetscCall(MatLUFactorNumeric(F, A, &info));
156: } else { /* Test inplace factorization */
157: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &F));
158: PetscCall(MatLUFactor(F, perm, iperm, &info));
159: }
161: PetscCall(VecDuplicate(y, &b1));
163: /* MatSolve */
164: PetscCall(MatSolve(F, b, x));
165: PetscCall(MatMult(A, x, b1));
166: PetscCall(VecAXPY(b1, -1.0, b));
167: PetscCall(VecNorm(b1, NORM_2, &norm));
168: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolve : Error of norm %g\n", (double)norm));
170: /* MatSolveTranspose */
171: PetscCall(MatSolveTranspose(F, b, x));
172: PetscCall(MatMultTranspose(A, x, b1));
173: PetscCall(VecAXPY(b1, -1.0, b));
174: PetscCall(VecNorm(b1, NORM_2, &norm));
175: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveTranspose : Error of norm %g\n", (double)norm));
177: /* MatSolveAdd */
178: PetscCall(MatSolveAdd(F, b, y, x));
179: PetscCall(MatMult(A, y, b1));
180: PetscCall(VecScale(b1, -1.0));
181: PetscCall(MatMultAdd(A, x, b1, b1));
182: PetscCall(VecAXPY(b1, -1.0, b));
183: PetscCall(VecNorm(b1, NORM_2, &norm));
184: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveAdd : Error of norm %g\n", (double)norm));
186: /* MatSolveTransposeAdd */
187: PetscCall(MatSolveTransposeAdd(F, b, y, x));
188: PetscCall(MatMultTranspose(A, y, b1));
189: PetscCall(VecScale(b1, -1.0));
190: PetscCall(MatMultTransposeAdd(A, x, b1, b1));
191: PetscCall(VecAXPY(b1, -1.0, b));
192: PetscCall(VecNorm(b1, NORM_2, &norm));
193: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveTransposeAdd : Error of norm %g\n", (double)norm));
195: /* MatMatSolve */
196: PetscCall(MatMatSolve(F, RHS, X));
197: PetscCall(MatMatMult(A, X, MAT_INITIAL_MATRIX, 2.0, &C1));
198: PetscCall(MatAXPY(C1, -1.0, RHS, SAME_NONZERO_PATTERN));
199: PetscCall(MatNorm(C1, NORM_FROBENIUS, &norm));
200: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve : Error of norm %g\n", (double)norm));
202: PetscCall(VecDestroy(&x));
203: PetscCall(VecDestroy(&b));
204: PetscCall(VecDestroy(&b1));
205: PetscCall(VecDestroy(&y));
206: PetscCall(MatDestroy(&A));
207: PetscCall(MatDestroy(&F));
208: PetscCall(MatDestroy(&RHS));
209: PetscCall(MatDestroy(&C1));
210: PetscCall(MatDestroy(&X));
211: PetscCall(ISDestroy(&perm));
212: PetscCall(ISDestroy(&iperm));
213: PetscCall(DMDestroy(&da));
214: PetscCall(PetscFinalize());
215: return 0;
216: }
218: PetscErrorCode ComputeRHS(DM da, Vec b)
219: {
220: PetscInt mx, my, mz;
221: PetscScalar h;
223: PetscFunctionBegin;
224: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
225: h = 1.0 / ((mx - 1) * (my - 1) * (mz - 1));
226: PetscCall(VecSet(b, h));
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: PetscErrorCode ComputeRHSMatrix(PetscInt m, PetscInt nrhs, Mat *C)
231: {
232: PetscRandom rand;
233: Mat RHS;
234: PetscScalar *array, rval;
235: PetscInt i, k;
237: PetscFunctionBegin;
238: PetscCall(MatCreate(PETSC_COMM_WORLD, &RHS));
239: PetscCall(MatSetSizes(RHS, m, PETSC_DECIDE, PETSC_DECIDE, nrhs));
240: PetscCall(MatSetType(RHS, MATSEQDENSE));
241: PetscCall(MatSetUp(RHS));
243: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
244: PetscCall(PetscRandomSetFromOptions(rand));
245: PetscCall(MatDenseGetArray(RHS, &array));
246: for (i = 0; i < m; i++) {
247: PetscCall(PetscRandomGetValue(rand, &rval));
248: array[i] = rval;
249: }
250: if (nrhs > 1) {
251: for (k = 1; k < nrhs; k++) {
252: for (i = 0; i < m; i++) array[m * k + i] = array[i];
253: }
254: }
255: PetscCall(MatDenseRestoreArray(RHS, &array));
256: PetscCall(MatAssemblyBegin(RHS, MAT_FINAL_ASSEMBLY));
257: PetscCall(MatAssemblyEnd(RHS, MAT_FINAL_ASSEMBLY));
258: *C = RHS;
259: PetscCall(PetscRandomDestroy(&rand));
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: PetscErrorCode ComputeMatrix(DM da, Mat B)
264: {
265: PetscInt i, j, k, mx, my, mz, xm, ym, zm, xs, ys, zs, dof, k1, k2, k3;
266: PetscScalar *v, *v_neighbor, Hx, Hy, Hz, HxHydHz, HyHzdHx, HxHzdHy, r1, r2;
267: MatStencil row, col;
268: PetscRandom rand;
270: PetscFunctionBegin;
271: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
272: PetscCall(PetscRandomSetSeed(rand, 1));
273: PetscCall(PetscRandomSetInterval(rand, -.001, .001));
274: PetscCall(PetscRandomSetFromOptions(rand));
276: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
277: /* For simplicity, this example only works on mx=my=mz */
278: PetscCheck(mx == my && mx == mz, PETSC_COMM_SELF, PETSC_ERR_SUP, "This example only works with mx %" PetscInt_FMT " = my %" PetscInt_FMT " = mz %" PetscInt_FMT, mx, my, mz);
280: Hx = 1.0 / (PetscReal)(mx - 1);
281: Hy = 1.0 / (PetscReal)(my - 1);
282: Hz = 1.0 / (PetscReal)(mz - 1);
283: HxHydHz = Hx * Hy / Hz;
284: HxHzdHy = Hx * Hz / Hy;
285: HyHzdHx = Hy * Hz / Hx;
287: PetscCall(PetscMalloc1(2 * dof * dof + 1, &v));
288: v_neighbor = v + dof * dof;
289: PetscCall(PetscArrayzero(v, 2 * dof * dof + 1));
290: k3 = 0;
291: for (k1 = 0; k1 < dof; k1++) {
292: for (k2 = 0; k2 < dof; k2++) {
293: if (k1 == k2) {
294: v[k3] = 2.0 * (HxHydHz + HxHzdHy + HyHzdHx);
295: v_neighbor[k3] = -HxHydHz;
296: } else {
297: PetscCall(PetscRandomGetValue(rand, &r1));
298: PetscCall(PetscRandomGetValue(rand, &r2));
300: v[k3] = r1;
301: v_neighbor[k3] = r2;
302: }
303: k3++;
304: }
305: }
306: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
308: for (k = zs; k < zs + zm; k++) {
309: for (j = ys; j < ys + ym; j++) {
310: for (i = xs; i < xs + xm; i++) {
311: row.i = i;
312: row.j = j;
313: row.k = k;
314: if (i == 0 || j == 0 || k == 0 || i == mx - 1 || j == my - 1 || k == mz - 1) { /* boundary points */
315: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &row, v, INSERT_VALUES));
316: } else { /* interior points */
317: /* center */
318: col.i = i;
319: col.j = j;
320: col.k = k;
321: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v, INSERT_VALUES));
323: /* x neighbors */
324: col.i = i - 1;
325: col.j = j;
326: col.k = k;
327: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
328: col.i = i + 1;
329: col.j = j;
330: col.k = k;
331: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
333: /* y neighbors */
334: col.i = i;
335: col.j = j - 1;
336: col.k = k;
337: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
338: col.i = i;
339: col.j = j + 1;
340: col.k = k;
341: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
343: /* z neighbors */
344: col.i = i;
345: col.j = j;
346: col.k = k - 1;
347: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
348: col.i = i;
349: col.j = j;
350: col.k = k + 1;
351: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
352: }
353: }
354: }
355: }
356: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
357: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
358: PetscCall(PetscFree(v));
359: PetscCall(PetscRandomDestroy(&rand));
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: /*TEST
365: test:
366: args: -dm_mat_type aij -dof 1
367: output_file: output/empty.out
369: test:
370: suffix: 2
371: args: -dm_mat_type aij -dof 1 -inplacelu
372: output_file: output/empty.out
374: TEST*/