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