Actual source code: ex221.c

  1: static char help[] = "Tests various routines for MATSHELL\n\n";

  3: #include <petscmat.h>

  5: typedef struct _n_User *User;
  6: struct _n_User {
  7:   Mat B;
  8: };

 10: static PetscErrorCode MatGetDiagonal_User(Mat A, Vec X)
 11: {
 12:   User user;

 14:   MatShellGetContext(A, &user);
 15:   MatGetDiagonal(user->B, X);
 16:   return 0;
 17: }

 19: static PetscErrorCode MatMult_User(Mat A, Vec X, Vec Y)
 20: {
 21:   User user;

 23:   MatShellGetContext(A, &user);
 24:   MatMult(user->B, X, Y);
 25:   return 0;
 26: }

 28: static PetscErrorCode MatMultTranspose_User(Mat A, Vec X, Vec Y)
 29: {
 30:   User user;

 32:   MatShellGetContext(A, &user);
 33:   MatMultTranspose(user->B, X, Y);
 34:   return 0;
 35: }

 37: static PetscErrorCode MatCopy_User(Mat A, Mat X, MatStructure str)
 38: {
 39:   User user, userX;

 41:   MatShellGetContext(A, &user);
 42:   MatShellGetContext(X, &userX);
 44:   PetscObjectReference((PetscObject)user->B);
 45:   return 0;
 46: }

 48: static PetscErrorCode MatDestroy_User(Mat A)
 49: {
 50:   User user;

 52:   MatShellGetContext(A, &user);
 53:   PetscObjectDereference((PetscObject)user->B);
 54:   return 0;
 55: }

 57: int main(int argc, char **args)
 58: {
 59:   User         user;
 60:   Mat          A, S;
 61:   PetscScalar *data, diag = 1.3;
 62:   PetscReal    tol = PETSC_SMALL;
 63:   PetscInt     i, j, m = PETSC_DECIDE, n = PETSC_DECIDE, M = 17, N = 15, s1, s2;
 64:   PetscInt     test, ntest = 2;
 65:   PetscMPIInt  rank, size;
 66:   PetscBool    nc        = PETSC_FALSE, cong, flg;
 67:   PetscBool    ronl      = PETSC_TRUE;
 68:   PetscBool    randomize = PETSC_FALSE, submat = PETSC_FALSE;
 69:   PetscBool    keep         = PETSC_FALSE;
 70:   PetscBool    testzerorows = PETSC_TRUE, testdiagscale = PETSC_TRUE, testgetdiag = PETSC_TRUE, testsubmat = PETSC_TRUE;
 71:   PetscBool    testshift = PETSC_TRUE, testscale = PETSC_TRUE, testdup = PETSC_TRUE, testreset = PETSC_TRUE;
 72:   PetscBool    testaxpy = PETSC_TRUE, testaxpyd = PETSC_TRUE, testaxpyerr = PETSC_FALSE;

 75:   PetscInitialize(&argc, &args, (char *)0, help);
 76:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 77:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 78:   PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
 79:   PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
 80:   PetscOptionsGetInt(NULL, NULL, "-ml", &m, NULL);
 81:   PetscOptionsGetInt(NULL, NULL, "-nl", &n, NULL);
 82:   PetscOptionsGetBool(NULL, NULL, "-square_nc", &nc, NULL);
 83:   PetscOptionsGetBool(NULL, NULL, "-rows_only", &ronl, NULL);
 84:   PetscOptionsGetBool(NULL, NULL, "-randomize", &randomize, NULL);
 85:   PetscOptionsGetBool(NULL, NULL, "-submat", &submat, NULL);
 86:   PetscOptionsGetBool(NULL, NULL, "-test_zerorows", &testzerorows, NULL);
 87:   PetscOptionsGetBool(NULL, NULL, "-test_diagscale", &testdiagscale, NULL);
 88:   PetscOptionsGetBool(NULL, NULL, "-test_getdiag", &testgetdiag, NULL);
 89:   PetscOptionsGetBool(NULL, NULL, "-test_shift", &testshift, NULL);
 90:   PetscOptionsGetBool(NULL, NULL, "-test_scale", &testscale, NULL);
 91:   PetscOptionsGetBool(NULL, NULL, "-test_dup", &testdup, NULL);
 92:   PetscOptionsGetBool(NULL, NULL, "-test_reset", &testreset, NULL);
 93:   PetscOptionsGetBool(NULL, NULL, "-test_submat", &testsubmat, NULL);
 94:   PetscOptionsGetBool(NULL, NULL, "-test_axpy", &testaxpy, NULL);
 95:   PetscOptionsGetBool(NULL, NULL, "-test_axpy_different", &testaxpyd, NULL);
 96:   PetscOptionsGetBool(NULL, NULL, "-test_axpy_error", &testaxpyerr, NULL);
 97:   PetscOptionsGetInt(NULL, NULL, "-loop", &ntest, NULL);
 98:   PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL);
 99:   PetscOptionsGetScalar(NULL, NULL, "-diag", &diag, NULL);
100:   PetscOptionsGetBool(NULL, NULL, "-keep", &keep, NULL);
101:   /* This tests square matrices with different row/col layout */
102:   if (nc && size > 1) {
103:     M = PetscMax(PetscMax(N, M), 1);
104:     N = M;
105:     m = n = 0;
106:     if (rank == 0) {
107:       m = M - 1;
108:       n = 1;
109:     } else if (rank == 1) {
110:       m = 1;
111:       n = N - 1;
112:     }
113:   }
114:   MatCreateDense(PETSC_COMM_WORLD, m, n, M, N, NULL, &A);
115:   MatGetLocalSize(A, &m, &n);
116:   MatGetSize(A, &M, &N);
117:   MatGetOwnershipRange(A, &s1, NULL);
118:   s2 = 1;
119:   while (s2 < M) s2 *= 10;
120:   MatDenseGetArray(A, &data);
121:   for (j = 0; j < N; j++) {
122:     for (i = 0; i < m; i++) data[j * m + i] = s2 * j + i + s1 + 1;
123:   }
124:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
125:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

127:   if (submat) {
128:     Mat      A2;
129:     IS       r, c;
130:     PetscInt rst, ren, cst, cen;

132:     MatGetOwnershipRange(A, &rst, &ren);
133:     MatGetOwnershipRangeColumn(A, &cst, &cen);
134:     ISCreateStride(PetscObjectComm((PetscObject)A), (ren - rst) / 2, rst, 1, &r);
135:     ISCreateStride(PetscObjectComm((PetscObject)A), (cen - cst) / 2, cst, 1, &c);
136:     MatCreateSubMatrix(A, r, c, MAT_INITIAL_MATRIX, &A2);
137:     ISDestroy(&r);
138:     ISDestroy(&c);
139:     MatDestroy(&A);
140:     A = A2;
141:   }

143:   MatGetSize(A, &M, &N);
144:   MatGetLocalSize(A, &m, &n);
145:   MatHasCongruentLayouts(A, &cong);

147:   MatConvert(A, MATAIJ, MAT_INPLACE_MATRIX, &A);
148:   MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, keep);
149:   PetscObjectSetName((PetscObject)A, "initial");
150:   MatViewFromOptions(A, NULL, "-view_mat");

152:   PetscNew(&user);
153:   MatCreateShell(PETSC_COMM_WORLD, m, n, M, N, user, &S);
154:   MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_User);
155:   MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_User);
156:   if (cong) MatShellSetOperation(S, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_User);
157:   MatShellSetOperation(S, MATOP_COPY, (void (*)(void))MatCopy_User);
158:   MatShellSetOperation(S, MATOP_DESTROY, (void (*)(void))MatDestroy_User);
159:   MatDuplicate(A, MAT_COPY_VALUES, &user->B);

161:   /* Square and rows only scaling */
162:   ronl = cong ? ronl : PETSC_TRUE;

164:   for (test = 0; test < ntest; test++) {
165:     PetscReal err;

167:     MatMultAddEqual(A, S, 10, &flg);
168:     if (!flg) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mult add\n", test);
169:     MatMultTransposeAddEqual(A, S, 10, &flg);
170:     if (!flg) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mult add (T)\n", test);
171:     if (testzerorows) {
172:       Mat       ST, B, C, BT, BTT;
173:       IS        zr;
174:       Vec       x = NULL, b1 = NULL, b2 = NULL;
175:       PetscInt *idxs = NULL, nr = 0;

177:       if (rank == (test % size)) {
178:         nr = 1;
179:         PetscMalloc1(nr, &idxs);
180:         if (test % 2) {
181:           idxs[0] = (2 * M - 1 - test / 2) % M;
182:         } else {
183:           idxs[0] = (test / 2) % M;
184:         }
185:         idxs[0] = PetscMax(idxs[0], 0);
186:       }
187:       ISCreateGeneral(PETSC_COMM_WORLD, nr, idxs, PETSC_OWN_POINTER, &zr);
188:       PetscObjectSetName((PetscObject)zr, "ZR");
189:       ISViewFromOptions(zr, NULL, "-view_is");
190:       MatCreateVecs(A, &x, &b1);
191:       if (randomize) {
192:         VecSetRandom(x, NULL);
193:         VecSetRandom(b1, NULL);
194:       } else {
195:         VecSet(x, 11.4);
196:         VecSet(b1, -14.2);
197:       }
198:       VecDuplicate(b1, &b2);
199:       VecCopy(b1, b2);
200:       PetscObjectSetName((PetscObject)b1, "A_B1");
201:       PetscObjectSetName((PetscObject)b2, "A_B2");
202:       if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */
203:         VecDestroy(&b1);
204:       }
205:       if (ronl) {
206:         MatZeroRowsIS(A, zr, diag, x, b1);
207:         MatZeroRowsIS(S, zr, diag, x, b2);
208:       } else {
209:         MatZeroRowsColumnsIS(A, zr, diag, x, b1);
210:         MatZeroRowsColumnsIS(S, zr, diag, x, b2);
211:         ISDestroy(&zr);
212:         /* Mix zerorows and zerorowscols */
213:         nr   = 0;
214:         idxs = NULL;
215:         if (rank == 0) {
216:           nr = 1;
217:           PetscMalloc1(nr, &idxs);
218:           if (test % 2) {
219:             idxs[0] = (3 * M - 2 - test / 2) % M;
220:           } else {
221:             idxs[0] = (test / 2 + 1) % M;
222:           }
223:           idxs[0] = PetscMax(idxs[0], 0);
224:         }
225:         ISCreateGeneral(PETSC_COMM_WORLD, nr, idxs, PETSC_OWN_POINTER, &zr);
226:         PetscObjectSetName((PetscObject)zr, "ZR2");
227:         ISViewFromOptions(zr, NULL, "-view_is");
228:         MatZeroRowsIS(A, zr, diag * 2.0 + PETSC_SMALL, NULL, NULL);
229:         MatZeroRowsIS(S, zr, diag * 2.0 + PETSC_SMALL, NULL, NULL);
230:       }
231:       ISDestroy(&zr);

233:       if (b1) {
234:         Vec b;

236:         VecViewFromOptions(b1, NULL, "-view_b");
237:         VecViewFromOptions(b2, NULL, "-view_b");
238:         VecDuplicate(b1, &b);
239:         VecCopy(b1, b);
240:         VecAXPY(b, -1.0, b2);
241:         VecNorm(b, NORM_INFINITY, &err);
242:         if (err >= tol) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error b %g\n", test, (double)err);
243:         VecDestroy(&b);
244:       }
245:       VecDestroy(&b1);
246:       VecDestroy(&b2);
247:       VecDestroy(&x);
248:       MatConvert(S, MATDENSE, MAT_INITIAL_MATRIX, &B);

250:       MatCreateTranspose(S, &ST);
251:       MatComputeOperator(ST, MATDENSE, &BT);
252:       MatTranspose(BT, MAT_INITIAL_MATRIX, &BTT);
253:       PetscObjectSetName((PetscObject)B, "S");
254:       PetscObjectSetName((PetscObject)BTT, "STT");
255:       MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &C);
256:       PetscObjectSetName((PetscObject)C, "A");

258:       MatViewFromOptions(C, NULL, "-view_mat");
259:       MatViewFromOptions(B, NULL, "-view_mat");
260:       MatViewFromOptions(BTT, NULL, "-view_mat");

262:       MatAXPY(C, -1.0, B, SAME_NONZERO_PATTERN);
263:       MatNorm(C, NORM_FROBENIUS, &err);
264:       if (err >= tol) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat mult after %s %g\n", test, ronl ? "MatZeroRows" : "MatZeroRowsColumns", (double)err);

266:       MatConvert(A, MATDENSE, MAT_REUSE_MATRIX, &C);
267:       MatAXPY(C, -1.0, BTT, SAME_NONZERO_PATTERN);
268:       MatNorm(C, NORM_FROBENIUS, &err);
269:       if (err >= tol) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat mult transpose after %s %g\n", test, ronl ? "MatZeroRows" : "MatZeroRowsColumns", (double)err);

271:       MatDestroy(&ST);
272:       MatDestroy(&BTT);
273:       MatDestroy(&BT);
274:       MatDestroy(&B);
275:       MatDestroy(&C);
276:     }
277:     if (testdiagscale) { /* MatDiagonalScale() */
278:       Vec vr, vl;

280:       MatCreateVecs(A, &vr, &vl);
281:       if (randomize) {
282:         VecSetRandom(vr, NULL);
283:         VecSetRandom(vl, NULL);
284:       } else {
285:         VecSet(vr, test % 2 ? 0.15 : 1.0 / 0.15);
286:         VecSet(vl, test % 2 ? -1.2 : 1.0 / -1.2);
287:       }
288:       MatDiagonalScale(A, vl, vr);
289:       MatDiagonalScale(S, vl, vr);
290:       VecDestroy(&vr);
291:       VecDestroy(&vl);
292:     }

294:     if (testscale) { /* MatScale() */
295:       MatScale(A, test % 2 ? 1.4 : 1.0 / 1.4);
296:       MatScale(S, test % 2 ? 1.4 : 1.0 / 1.4);
297:     }

299:     if (testshift && cong) { /* MatShift() : MATSHELL shift is broken when row/cols layout are not congruent and left/right scaling have been applied */
300:       MatShift(A, test % 2 ? -77.5 : 77.5);
301:       MatShift(S, test % 2 ? -77.5 : 77.5);
302:     }

304:     if (testgetdiag && cong) { /* MatGetDiagonal() */
305:       Vec dA, dS;

307:       MatCreateVecs(A, &dA, NULL);
308:       MatCreateVecs(S, &dS, NULL);
309:       MatGetDiagonal(A, dA);
310:       MatGetDiagonal(S, dS);
311:       VecAXPY(dA, -1.0, dS);
312:       VecNorm(dA, NORM_INFINITY, &err);
313:       if (err >= tol) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error diag %g\n", test, (double)err);
314:       VecDestroy(&dA);
315:       VecDestroy(&dS);
316:     }

318:     if (testdup && !test) {
319:       Mat A2, S2;

321:       MatDuplicate(A, MAT_COPY_VALUES, &A2);
322:       MatDuplicate(S, MAT_COPY_VALUES, &S2);
323:       MatDestroy(&A);
324:       MatDestroy(&S);
325:       A = A2;
326:       S = S2;
327:     }

329:     if (testsubmat) {
330:       Mat      sA, sS, dA, dS, At, St;
331:       IS       r, c;
332:       PetscInt rst, ren, cst, cen;

334:       MatGetOwnershipRange(A, &rst, &ren);
335:       MatGetOwnershipRangeColumn(A, &cst, &cen);
336:       ISCreateStride(PetscObjectComm((PetscObject)A), (ren - rst) / 2, rst, 1, &r);
337:       ISCreateStride(PetscObjectComm((PetscObject)A), (cen - cst) / 2, cst, 1, &c);
338:       MatCreateSubMatrix(A, r, c, MAT_INITIAL_MATRIX, &sA);
339:       MatCreateSubMatrix(S, r, c, MAT_INITIAL_MATRIX, &sS);
340:       MatMultAddEqual(sA, sS, 10, &flg);
341:       if (!flg) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix mult add\n", test);
342:       MatMultTransposeAddEqual(sA, sS, 10, &flg);
343:       if (!flg) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix mult add (T)\n", test);
344:       MatConvert(sA, MATDENSE, MAT_INITIAL_MATRIX, &dA);
345:       MatConvert(sS, MATDENSE, MAT_INITIAL_MATRIX, &dS);
346:       MatAXPY(dA, -1.0, dS, SAME_NONZERO_PATTERN);
347:       MatNorm(dA, NORM_FROBENIUS, &err);
348:       if (err >= tol) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat submatrix %g\n", test, (double)err);
349:       MatDestroy(&sA);
350:       MatDestroy(&sS);
351:       MatDestroy(&dA);
352:       MatDestroy(&dS);
353:       MatCreateTranspose(A, &At);
354:       MatCreateTranspose(S, &St);
355:       MatCreateSubMatrix(At, c, r, MAT_INITIAL_MATRIX, &sA);
356:       MatCreateSubMatrix(St, c, r, MAT_INITIAL_MATRIX, &sS);
357:       MatMultAddEqual(sA, sS, 10, &flg);
358:       if (!flg) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix (T) mult add\n", test);
359:       MatMultTransposeAddEqual(sA, sS, 10, &flg);
360:       if (!flg) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error submatrix (T) mult add (T)\n", test);
361:       MatConvert(sA, MATDENSE, MAT_INITIAL_MATRIX, &dA);
362:       MatConvert(sS, MATDENSE, MAT_INITIAL_MATRIX, &dS);
363:       MatAXPY(dA, -1.0, dS, SAME_NONZERO_PATTERN);
364:       MatNorm(dA, NORM_FROBENIUS, &err);
365:       if (err >= tol) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat submatrix (T) %g\n", test, (double)err);
366:       MatDestroy(&sA);
367:       MatDestroy(&sS);
368:       MatDestroy(&dA);
369:       MatDestroy(&dS);
370:       MatDestroy(&At);
371:       MatDestroy(&St);
372:       ISDestroy(&r);
373:       ISDestroy(&c);
374:     }

376:     if (testaxpy) {
377:       Mat          tA, tS, dA, dS;
378:       MatStructure str[3] = {SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN};

380:       MatDuplicate(A, MAT_COPY_VALUES, &tA);
381:       if (testaxpyd && !(test % 2)) {
382:         PetscObjectReference((PetscObject)tA);
383:         tS = tA;
384:       } else {
385:         PetscObjectReference((PetscObject)S);
386:         tS = S;
387:       }
388:       MatAXPY(A, 0.5, tA, str[test % 3]);
389:       MatAXPY(S, 0.5, tS, str[test % 3]);
390:       /* this will trigger an error the next MatMult or MatMultTranspose call for S */
391:       if (testaxpyerr) MatScale(tA, 0);
392:       MatDestroy(&tA);
393:       MatDestroy(&tS);
394:       MatMultAddEqual(A, S, 10, &flg);
395:       if (!flg) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error axpy mult add\n", test);
396:       MatMultTransposeAddEqual(A, S, 10, &flg);
397:       if (!flg) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error axpy mult add (T)\n", test);
398:       MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &dA);
399:       MatConvert(S, MATDENSE, MAT_INITIAL_MATRIX, &dS);
400:       MatAXPY(dA, -1.0, dS, SAME_NONZERO_PATTERN);
401:       MatNorm(dA, NORM_FROBENIUS, &err);
402:       if (err >= tol) PetscPrintf(PETSC_COMM_WORLD, "[test %" PetscInt_FMT "] Error mat submatrix %g\n", test, (double)err);
403:       MatDestroy(&dA);
404:       MatDestroy(&dS);
405:     }

407:     if (testreset && (ntest == 1 || test == ntest - 2)) {
408:       /* reset MATSHELL */
409:       MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
410:       MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
411:       /* reset A */
412:       MatCopy(user->B, A, DIFFERENT_NONZERO_PATTERN);
413:     }
414:   }

416:   MatDestroy(&A);
417:   MatDestroy(&S);
418:   PetscFree(user);
419:   PetscFinalize();
420:   return 0;
421: }

423: /*TEST

425:    testset:
426:      suffix: rect
427:      requires: !single
428:      output_file: output/ex221_1.out
429:      nsize: {{1 3}}
430:      args: -loop 3 -keep {{0 1}} -M {{12 19}} -N {{19 12}} -submat {{0 1}} -test_axpy_different {{0 1}}

432:    testset:
433:      suffix: square
434:      requires: !single
435:      output_file: output/ex221_1.out
436:      nsize: {{1 3}}
437:      args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}} -submat {{0 1}} -test_axpy_different {{0 1}}
438: TEST*/