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(), and MatMatSolve().\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: int main(int argc, char **args)
22: {
23: PetscMPIInt size;
24: Vec x, b, y, b1;
25: DM da;
26: Mat A, F, RHS, X, C1;
27: MatFactorInfo info;
28: IS perm, iperm;
29: PetscInt dof = 1, M = 8, m, n, nrhs;
30: PetscScalar one = 1.0;
31: PetscReal norm, tol = 1000 * PETSC_MACHINE_EPSILON;
32: PetscBool InplaceLU = PETSC_FALSE;
34: PetscFunctionBeginUser;
35: PetscCall(PetscInitialize(&argc, &args, NULL, help));
36: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
37: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only");
38: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL));
39: PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
41: PetscCall(DMDACreate(PETSC_COMM_WORLD, &da));
42: PetscCall(DMSetDimension(da, 3));
43: PetscCall(DMDASetBoundaryType(da, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE));
44: PetscCall(DMDASetStencilType(da, DMDA_STENCIL_STAR));
45: PetscCall(DMDASetSizes(da, M, M, M));
46: PetscCall(DMDASetNumProcs(da, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
47: PetscCall(DMDASetDof(da, dof));
48: PetscCall(DMDASetStencilWidth(da, 1));
49: PetscCall(DMDASetOwnershipRanges(da, NULL, NULL, NULL));
50: PetscCall(DMSetMatType(da, MATBAIJ));
51: PetscCall(DMSetFromOptions(da));
52: PetscCall(DMSetUp(da));
54: PetscCall(DMCreateGlobalVector(da, &x));
55: PetscCall(DMCreateGlobalVector(da, &b));
56: PetscCall(VecDuplicate(b, &y));
57: PetscCall(ComputeRHS(da, b));
58: PetscCall(VecSet(y, one));
59: PetscCall(DMCreateMatrix(da, &A));
60: PetscCall(ComputeMatrix(da, A));
61: PetscCall(MatGetSize(A, &m, &n));
62: nrhs = 2;
63: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
64: PetscCall(ComputeRHSMatrix(m, nrhs, &RHS));
65: PetscCall(MatDuplicate(RHS, MAT_DO_NOT_COPY_VALUES, &X));
67: PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm));
69: PetscCall(PetscOptionsGetBool(NULL, NULL, "-inplacelu", &InplaceLU, NULL));
70: PetscCall(MatFactorInfoInitialize(&info));
71: if (!InplaceLU) {
72: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F));
73: info.fill = 5.0;
74: PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
75: PetscCall(MatLUFactorNumeric(F, A, &info));
76: } else { /* Test inplace factorization */
77: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &F));
78: PetscCall(MatLUFactor(F, perm, iperm, &info));
79: }
81: PetscCall(VecDuplicate(y, &b1));
83: /* MatSolve */
84: PetscCall(MatSolve(F, b, x));
85: PetscCall(MatMult(A, x, b1));
86: PetscCall(VecAXPY(b1, -1.0, b));
87: PetscCall(VecNorm(b1, NORM_2, &norm));
88: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolve : Error of norm %g\n", (double)norm));
90: /* MatSolveTranspose */
91: PetscCall(MatSolveTranspose(F, b, x));
92: PetscCall(MatMultTranspose(A, x, b1));
93: PetscCall(VecAXPY(b1, -1.0, b));
94: PetscCall(VecNorm(b1, NORM_2, &norm));
95: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveTranspose : Error of norm %g\n", (double)norm));
97: /* MatSolveAdd */
98: PetscCall(MatSolveAdd(F, b, y, x));
99: PetscCall(MatMult(A, y, b1));
100: PetscCall(VecScale(b1, -1.0));
101: PetscCall(MatMultAdd(A, x, b1, b1));
102: PetscCall(VecAXPY(b1, -1.0, b));
103: PetscCall(VecNorm(b1, NORM_2, &norm));
104: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveAdd : Error of norm %g\n", (double)norm));
106: /* MatSolveTransposeAdd */
107: PetscCall(MatSolveTransposeAdd(F, b, y, x));
108: PetscCall(MatMultTranspose(A, y, b1));
109: PetscCall(VecScale(b1, -1.0));
110: PetscCall(MatMultTransposeAdd(A, x, b1, b1));
111: PetscCall(VecAXPY(b1, -1.0, b));
112: PetscCall(VecNorm(b1, NORM_2, &norm));
113: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveTransposeAdd : Error of norm %g\n", (double)norm));
115: /* MatMatSolve */
116: PetscCall(MatMatSolve(F, RHS, X));
117: PetscCall(MatMatMult(A, X, MAT_INITIAL_MATRIX, 2.0, &C1));
118: PetscCall(MatAXPY(C1, -1.0, RHS, SAME_NONZERO_PATTERN));
119: PetscCall(MatNorm(C1, NORM_FROBENIUS, &norm));
120: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve : Error of norm %g\n", (double)norm));
122: PetscCall(VecDestroy(&x));
123: PetscCall(VecDestroy(&b));
124: PetscCall(VecDestroy(&b1));
125: PetscCall(VecDestroy(&y));
126: PetscCall(MatDestroy(&A));
127: PetscCall(MatDestroy(&F));
128: PetscCall(MatDestroy(&RHS));
129: PetscCall(MatDestroy(&C1));
130: PetscCall(MatDestroy(&X));
131: PetscCall(ISDestroy(&perm));
132: PetscCall(ISDestroy(&iperm));
133: PetscCall(DMDestroy(&da));
134: PetscCall(PetscFinalize());
135: return 0;
136: }
138: PetscErrorCode ComputeRHS(DM da, Vec b)
139: {
140: PetscInt mx, my, mz;
141: PetscScalar h;
143: PetscFunctionBegin;
144: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
145: h = 1.0 / ((mx - 1) * (my - 1) * (mz - 1));
146: PetscCall(VecSet(b, h));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: PetscErrorCode ComputeRHSMatrix(PetscInt m, PetscInt nrhs, Mat *C)
151: {
152: PetscRandom rand;
153: Mat RHS;
154: PetscScalar *array, rval;
155: PetscInt i, k;
157: PetscFunctionBegin;
158: PetscCall(MatCreate(PETSC_COMM_WORLD, &RHS));
159: PetscCall(MatSetSizes(RHS, m, PETSC_DECIDE, PETSC_DECIDE, nrhs));
160: PetscCall(MatSetType(RHS, MATSEQDENSE));
161: PetscCall(MatSetUp(RHS));
163: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
164: PetscCall(PetscRandomSetFromOptions(rand));
165: PetscCall(MatDenseGetArray(RHS, &array));
166: for (i = 0; i < m; i++) {
167: PetscCall(PetscRandomGetValue(rand, &rval));
168: array[i] = rval;
169: }
170: if (nrhs > 1) {
171: for (k = 1; k < nrhs; k++) {
172: for (i = 0; i < m; i++) array[m * k + i] = array[i];
173: }
174: }
175: PetscCall(MatDenseRestoreArray(RHS, &array));
176: PetscCall(MatAssemblyBegin(RHS, MAT_FINAL_ASSEMBLY));
177: PetscCall(MatAssemblyEnd(RHS, MAT_FINAL_ASSEMBLY));
178: *C = RHS;
179: PetscCall(PetscRandomDestroy(&rand));
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: PetscErrorCode ComputeMatrix(DM da, Mat B)
184: {
185: PetscInt i, j, k, mx, my, mz, xm, ym, zm, xs, ys, zs, dof, k1, k2, k3;
186: PetscScalar *v, *v_neighbor, Hx, Hy, Hz, HxHydHz, HyHzdHx, HxHzdHy, r1, r2;
187: MatStencil row, col;
188: PetscRandom rand;
190: PetscFunctionBegin;
191: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
192: PetscCall(PetscRandomSetSeed(rand, 1));
193: PetscCall(PetscRandomSetInterval(rand, -.001, .001));
194: PetscCall(PetscRandomSetFromOptions(rand));
196: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
197: /* For simplicity, this example only works on mx=my=mz */
198: 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);
200: Hx = 1.0 / (PetscReal)(mx - 1);
201: Hy = 1.0 / (PetscReal)(my - 1);
202: Hz = 1.0 / (PetscReal)(mz - 1);
203: HxHydHz = Hx * Hy / Hz;
204: HxHzdHy = Hx * Hz / Hy;
205: HyHzdHx = Hy * Hz / Hx;
207: PetscCall(PetscMalloc1(2 * dof * dof + 1, &v));
208: v_neighbor = v + dof * dof;
209: PetscCall(PetscArrayzero(v, 2 * dof * dof + 1));
210: k3 = 0;
211: for (k1 = 0; k1 < dof; k1++) {
212: for (k2 = 0; k2 < dof; k2++) {
213: if (k1 == k2) {
214: v[k3] = 2.0 * (HxHydHz + HxHzdHy + HyHzdHx);
215: v_neighbor[k3] = -HxHydHz;
216: } else {
217: PetscCall(PetscRandomGetValue(rand, &r1));
218: PetscCall(PetscRandomGetValue(rand, &r2));
220: v[k3] = r1;
221: v_neighbor[k3] = r2;
222: }
223: k3++;
224: }
225: }
226: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
228: for (k = zs; k < zs + zm; k++) {
229: for (j = ys; j < ys + ym; j++) {
230: for (i = xs; i < xs + xm; i++) {
231: row.i = i;
232: row.j = j;
233: row.k = k;
234: if (i == 0 || j == 0 || k == 0 || i == mx - 1 || j == my - 1 || k == mz - 1) { /* boundary points */
235: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &row, v, INSERT_VALUES));
236: } else { /* interior points */
237: /* center */
238: col.i = i;
239: col.j = j;
240: col.k = k;
241: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v, INSERT_VALUES));
243: /* x neighbors */
244: col.i = i - 1;
245: col.j = j;
246: col.k = k;
247: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
248: col.i = i + 1;
249: col.j = j;
250: col.k = k;
251: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
253: /* y neighbors */
254: col.i = i;
255: col.j = j - 1;
256: col.k = k;
257: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
258: col.i = i;
259: col.j = j + 1;
260: col.k = k;
261: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
263: /* z neighbors */
264: col.i = i;
265: col.j = j;
266: col.k = k - 1;
267: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
268: col.i = i;
269: col.j = j;
270: col.k = k + 1;
271: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
272: }
273: }
274: }
275: }
276: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
277: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
278: PetscCall(PetscFree(v));
279: PetscCall(PetscRandomDestroy(&rand));
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: /*TEST
285: test:
286: args: -dm_mat_type aij -dof 1
287: output_file: output/ex129.out
289: test:
290: suffix: 2
291: args: -dm_mat_type aij -dof 1 -inplacelu
292: output_file: output/ex129.out
294: TEST*/