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