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