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, (char *)0, 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*/