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*/