1 /*
2 Laplacian in 3D. Use for testing MatSolve routines.
3 Modeled by the partial differential equation
4
5 - Laplacian u = 1,0 < x,y,z < 1,
6
7 with boundary conditions
8 u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
9 */
10
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";
13
14 #include <petscdm.h>
15 #include <petscdmda.h>
16
17 extern PetscErrorCode ComputeMatrix(DM, Mat);
18 extern PetscErrorCode ComputeRHS(DM, Vec);
19 extern PetscErrorCode ComputeRHSMatrix(PetscInt, PetscInt, Mat *);
20
main(int argc,char ** args)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;
33
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));
40
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));
53
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));
66
67 PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm));
68
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 }
80
81 PetscCall(VecDuplicate(y, &b1));
82
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));
89
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));
96
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));
105
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));
114
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));
121
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 }
137
ComputeRHS(DM da,Vec b)138 PetscErrorCode ComputeRHS(DM da, Vec b)
139 {
140 PetscInt mx, my, mz;
141 PetscScalar h;
142
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 }
149
ComputeRHSMatrix(PetscInt m,PetscInt nrhs,Mat * C)150 PetscErrorCode ComputeRHSMatrix(PetscInt m, PetscInt nrhs, Mat *C)
151 {
152 PetscRandom rand;
153 Mat RHS;
154 PetscScalar *array, rval;
155 PetscInt i, k;
156
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));
162
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 }
182
ComputeMatrix(DM da,Mat B)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;
189
190 PetscFunctionBegin;
191 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
192 PetscCall(PetscRandomSetSeed(rand, 1));
193 PetscCall(PetscRandomSetInterval(rand, -.001, .001));
194 PetscCall(PetscRandomSetFromOptions(rand));
195
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);
199
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;
206
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));
219
220 v[k3] = r1;
221 v_neighbor[k3] = r2;
222 }
223 k3++;
224 }
225 }
226 PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
227
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));
242
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));
252
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));
262
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 }
282
283 /*TEST
284
285 test:
286 args: -dm_mat_type aij -dof 1
287 output_file: output/empty.out
288
289 test:
290 suffix: 2
291 args: -dm_mat_type aij -dof 1 -inplacelu
292 output_file: output/empty.out
293
294 TEST*/
295