Actual source code: ex70.c

  1: #include <petscmat.h>

  3: static char help[] = "Tests MatMat operations with MAT_REUSE_MATRIX and already allocated dense result.\n\n";

  5: static PetscScalar MAGIC_NUMBER = 12345;

  7: static PetscErrorCode CheckLocal(Mat A, Mat B, PetscScalar *a, PetscScalar *b)
  8: {
  9:   PetscBool wA = PETSC_FALSE, wB = PETSC_FALSE;
 10:   PetscBool wAv = PETSC_FALSE, wBv = PETSC_FALSE;
 11:   PetscInt  lda, i, j, m, n;

 13:   if (a) {
 14:     const PetscScalar *Aa;
 15:     MatDenseGetArrayRead(A, &Aa);
 16:     wA = (PetscBool)(a != Aa);
 17:     MatDenseGetLDA(A, &lda);
 18:     MatGetLocalSize(A, &m, &n);
 19:     for (j = 0; j < n; j++) {
 20:       for (i = m; i < lda; i++) {
 21:         if (Aa[j * lda + i] != MAGIC_NUMBER) wAv = PETSC_TRUE;
 22:       }
 23:     }
 24:     MatDenseRestoreArrayRead(A, &Aa);
 25:   }
 26:   if (b) {
 27:     const PetscScalar *Bb;
 28:     MatDenseGetArrayRead(B, &Bb);
 29:     wB = (PetscBool)(b != Bb);
 30:     MatDenseGetLDA(B, &lda);
 31:     MatGetLocalSize(B, &m, &n);
 32:     for (j = 0; j < n; j++) {
 33:       for (i = m; i < lda; i++) {
 34:         if (Bb[j * lda + i] != MAGIC_NUMBER) wBv = PETSC_TRUE;
 35:       }
 36:     }
 37:     MatDenseRestoreArrayRead(B, &Bb);
 38:   }
 41:   return 0;
 42: }

 44: typedef struct {
 45:   Mat A;
 46:   Mat P;
 47:   Mat R;
 48: } proj_data;

 50: PetscErrorCode proj_destroy(void *ctx)
 51: {
 52:   proj_data *userdata = (proj_data *)ctx;

 55:   MatDestroy(&userdata->A);
 56:   MatDestroy(&userdata->P);
 57:   MatDestroy(&userdata->R);
 58:   PetscFree(userdata);
 59:   return 0;
 60: }

 62: PetscErrorCode proj_mult(Mat S, Vec X, Vec Y)
 63: {
 64:   Mat        A, R, P;
 65:   Vec        Ax, Ay;
 66:   Vec        Px, Py;
 67:   proj_data *userdata;

 69:   MatShellGetContext(S, &userdata);
 71:   A = userdata->A;
 72:   R = userdata->R;
 73:   P = userdata->P;
 77:   MatCreateVecs(A, &Ax, &Ay);
 78:   if (R) {
 79:     MatCreateVecs(R, &Py, &Px);
 80:   } else {
 81:     MatCreateVecs(P, &Px, &Py);
 82:   }
 83:   VecCopy(X, Px);
 84:   if (P) {
 85:     MatMult(P, Px, Py);
 86:   } else {
 87:     MatMultTranspose(R, Px, Py);
 88:   }
 89:   VecCopy(Py, Ax);
 90:   MatMult(A, Ax, Ay);
 91:   VecCopy(Ay, Py);
 92:   if (P) {
 93:     MatMultTranspose(P, Py, Px);
 94:   } else {
 95:     MatMult(R, Py, Px);
 96:   }
 97:   VecCopy(Px, Y);
 98:   VecDestroy(&Px);
 99:   VecDestroy(&Py);
100:   VecDestroy(&Ax);
101:   VecDestroy(&Ay);
102:   return 0;
103: }

105: PetscErrorCode MyPtShellPMultSymbolic(Mat S, Mat P, Mat PtAP, void **ctx)
106: {
107:   proj_data *userdata;

109:   PetscNew(&userdata);
110:   MatShellSetContext(PtAP, userdata);
111:   *ctx = (void *)userdata;
112:   return 0;
113: }

115: PetscErrorCode MyPtShellPMultNumeric(Mat S, Mat P, Mat PtAP, void *ctx)
116: {
117:   Mat        A;
118:   proj_data *userdata = (proj_data *)ctx;

120:   MatShellGetContext(S, &A);
121:   PetscObjectReference((PetscObject)A);
122:   PetscObjectReference((PetscObject)P);
123:   MatDestroy(&userdata->A);
124:   MatDestroy(&userdata->P);
125:   MatDestroy(&userdata->R);
126:   userdata->A = A;
127:   userdata->P = P;
128:   MatShellSetOperation(PtAP, MATOP_MULT, (void (*)(void))proj_mult);
129:   MatSetUp(PtAP);
130:   MatAssemblyBegin(PtAP, MAT_FINAL_ASSEMBLY);
131:   MatAssemblyEnd(PtAP, MAT_FINAL_ASSEMBLY);
132:   return 0;
133: }

135: PetscErrorCode MyRShellRtMultSymbolic(Mat S, Mat R, Mat RARt, void **ctx)
136: {
137:   proj_data *userdata;

139:   PetscNew(&userdata);
140:   MatShellSetContext(RARt, userdata);
141:   *ctx = (void *)userdata;
142:   return 0;
143: }

145: PetscErrorCode MyRShellRtMultNumeric(Mat S, Mat R, Mat RARt, void *ctx)
146: {
147:   Mat        A;
148:   proj_data *userdata = (proj_data *)ctx;

150:   MatShellGetContext(S, &A);
151:   PetscObjectReference((PetscObject)A);
152:   PetscObjectReference((PetscObject)R);
153:   MatDestroy(&userdata->A);
154:   MatDestroy(&userdata->P);
155:   MatDestroy(&userdata->R);
156:   userdata->A = A;
157:   userdata->R = R;
158:   MatShellSetOperation(RARt, MATOP_MULT, (void (*)(void))proj_mult);
159:   MatSetUp(RARt);
160:   MatAssemblyBegin(RARt, MAT_FINAL_ASSEMBLY);
161:   MatAssemblyEnd(RARt, MAT_FINAL_ASSEMBLY);
162:   return 0;
163: }

165: PetscErrorCode MyMatShellMatMultNumeric(Mat S, Mat B, Mat C, void *ctx)
166: {
167:   Mat A;

169:   MatShellGetContext(S, &A);
170:   MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C);
171:   return 0;
172: }

174: PetscErrorCode MyMatTransposeShellMatMultNumeric(Mat S, Mat B, Mat C, void *ctx)
175: {
176:   Mat A;

178:   MatShellGetContext(S, &A);
179:   MatTransposeMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C);
180:   return 0;
181: }

183: PetscErrorCode MyMatShellMatTransposeMultNumeric(Mat S, Mat B, Mat C, void *ctx)
184: {
185:   Mat A;

187:   MatShellGetContext(S, &A);
188:   MatMatTransposeMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C);
189:   return 0;
190: }

192: int main(int argc, char **args)
193: {
194:   Mat          X, B, A, Bt, T, T2, PtAP = NULL, RARt = NULL, R = NULL;
195:   Vec          r, l, rs, ls;
196:   PetscInt     m, n, k, M = 10, N = 10, K = 5, ldx = 3, ldb = 5, ldr = 4;
197:   const char  *deft = MATAIJ;
198:   char         mattype[256];
199:   PetscBool    flg, symm = PETSC_FALSE, testtt = PETSC_TRUE, testnest = PETSC_TRUE, testtranspose = PETSC_TRUE, testcircular = PETSC_FALSE, local = PETSC_TRUE;
200:   PetscBool    testhtranspose = PETSC_FALSE; /* Hermitian transpose is not handled correctly and generates an error */
201:   PetscBool    xgpu = PETSC_FALSE, bgpu = PETSC_FALSE, testshellops = PETSC_FALSE, testproj = PETSC_TRUE, testrart = PETSC_TRUE, testmatmatt = PETSC_TRUE, testmattmat = PETSC_TRUE;
202:   PetscScalar *dataX = NULL, *dataB = NULL, *dataR = NULL, *dataBt = NULL;
203:   PetscScalar *aX, *aB, *aBt;
204:   PetscReal    err;

207:   PetscInitialize(&argc, &args, NULL, help);
208:   PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
209:   PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
210:   PetscOptionsGetInt(NULL, NULL, "-K", &K, NULL);
211:   PetscOptionsGetBool(NULL, NULL, "-symm", &symm, NULL);
212:   PetscOptionsGetBool(NULL, NULL, "-local", &local, NULL);
213:   PetscOptionsGetInt(NULL, NULL, "-ldx", &ldx, NULL);
214:   PetscOptionsGetInt(NULL, NULL, "-ldb", &ldb, NULL);
215:   PetscOptionsGetInt(NULL, NULL, "-ldr", &ldr, NULL);
216:   PetscOptionsGetBool(NULL, NULL, "-testtranspose", &testtranspose, NULL);
217:   PetscOptionsGetBool(NULL, NULL, "-testnest", &testnest, NULL);
218:   PetscOptionsGetBool(NULL, NULL, "-testtt", &testtt, NULL);
219:   PetscOptionsGetBool(NULL, NULL, "-testcircular", &testcircular, NULL);
220:   PetscOptionsGetBool(NULL, NULL, "-testshellops", &testshellops, NULL);
221:   PetscOptionsGetBool(NULL, NULL, "-testproj", &testproj, NULL);
222:   PetscOptionsGetBool(NULL, NULL, "-testrart", &testrart, NULL);
223:   PetscOptionsGetBool(NULL, NULL, "-testmatmatt", &testmatmatt, NULL);
224:   PetscOptionsGetBool(NULL, NULL, "-testmattmat", &testmattmat, NULL);
225:   PetscOptionsGetBool(NULL, NULL, "-xgpu", &xgpu, NULL);
226:   PetscOptionsGetBool(NULL, NULL, "-bgpu", &bgpu, NULL);
227:   PetscOptionsGetScalar(NULL, NULL, "-magic_number", &MAGIC_NUMBER, NULL);
228:   if (M != N) testproj = PETSC_FALSE;

230:   MatCreate(PETSC_COMM_WORLD, &A);
231:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N);
232:   MatSetType(A, MATAIJ);
233:   MatSetUp(A);
234:   MatSetRandom(A, NULL);
235:   if (M == N && symm) {
236:     Mat AT;

238:     MatTranspose(A, MAT_INITIAL_MATRIX, &AT);
239:     MatAXPY(A, 1.0, AT, DIFFERENT_NONZERO_PATTERN);
240:     MatDestroy(&AT);
241:     MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
242:   }
243:   MatViewFromOptions(A, NULL, "-A_init_view");
244:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "");
245:   PetscOptionsFList("-A_mat_type", "Matrix type", "MatSetType", MatList, deft, mattype, 256, &flg);
246:   PetscOptionsEnd();
247:   if (flg) {
248:     Mat A2;

250:     MatDuplicate(A, MAT_COPY_VALUES, &A2);
251:     MatConvert(A, mattype, MAT_INPLACE_MATRIX, &A);
252:     MatMultEqual(A, A2, 10, &flg);
253:     if (!flg) {
254:       Mat AE, A2E;

256:       PetscPrintf(PETSC_COMM_WORLD, "Error with convert\n");
257:       MatComputeOperator(A, MATDENSE, &AE);
258:       MatComputeOperator(A2, MATDENSE, &A2E);
259:       MatView(AE, NULL);
260:       MatView(A2E, NULL);
261:       MatAXPY(A2E, -1.0, A, SAME_NONZERO_PATTERN);
262:       MatView(A2E, NULL);
263:       MatDestroy(&A2E);
264:       MatDestroy(&AE);
265:     }
266:     MatDestroy(&A2);
267:   }
268:   MatViewFromOptions(A, NULL, "-A_view");

270:   MatGetLocalSize(A, &m, &n);
271:   if (local) {
272:     PetscInt i;

274:     PetscMalloc1((m + ldx) * K, &dataX);
275:     PetscMalloc1((n + ldb) * K, &dataB);
276:     for (i = 0; i < (m + ldx) * K; i++) dataX[i] = MAGIC_NUMBER;
277:     for (i = 0; i < (n + ldb) * K; i++) dataB[i] = MAGIC_NUMBER;
278:   }
279:   MatCreateDense(PETSC_COMM_WORLD, n, PETSC_DECIDE, N, K, dataB, &B);
280:   MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, dataX, &X);
281:   if (local) {
282:     MatDenseSetLDA(X, m + ldx);
283:     MatDenseSetLDA(B, n + ldb);
284:   }
285:   MatGetLocalSize(B, NULL, &k);
286:   if (local) {
287:     PetscInt i;

289:     PetscMalloc1((k + ldr) * N, &dataBt);
290:     for (i = 0; i < (k + ldr) * N; i++) dataBt[i] = MAGIC_NUMBER;
291:   }
292:   MatCreateDense(PETSC_COMM_WORLD, k, n, K, N, dataBt, &Bt);
293:   if (local) MatDenseSetLDA(Bt, k + ldr);

295:   /* store pointer to dense data for testing */
296:   MatDenseGetArrayRead(B, (const PetscScalar **)&dataB);
297:   MatDenseGetArrayRead(X, (const PetscScalar **)&dataX);
298:   MatDenseGetArrayRead(Bt, (const PetscScalar **)&dataBt);
299:   aX  = dataX;
300:   aB  = dataB;
301:   aBt = dataBt;
302:   MatDenseRestoreArrayRead(Bt, (const PetscScalar **)&dataBt);
303:   MatDenseRestoreArrayRead(B, (const PetscScalar **)&dataB);
304:   MatDenseRestoreArrayRead(X, (const PetscScalar **)&dataX);
305:   if (local) {
306:     dataX  = aX;
307:     dataB  = aB;
308:     dataBt = aBt;
309:   }

311:   MatSetRandom(X, NULL);
312:   MatSetRandom(B, NULL);
313:   MatSetRandom(Bt, NULL);
314:   CheckLocal(X, NULL, aX, NULL);
315:   CheckLocal(Bt, B, aBt, aB);

317:   /* convert to CUDA if needed */
318:   if (bgpu) {
319:     MatConvert(B, MATDENSECUDA, MAT_INPLACE_MATRIX, &B);
320:     MatConvert(Bt, MATDENSECUDA, MAT_INPLACE_MATRIX, &Bt);
321:   }
322:   if (xgpu) MatConvert(X, MATDENSECUDA, MAT_INPLACE_MATRIX, &X);
323:   CheckLocal(B, X, aB, aX);

325:   /* Test MatDenseGetSubMatrix */
326:   {
327:     Mat B2, T3, T4;

329:     MatDuplicate(B, MAT_COPY_VALUES, &B2);
330:     MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &T4);
331:     MatSetRandom(T4, NULL);
332:     MatAXPY(B2, 1.0, T4, SAME_NONZERO_PATTERN);
333:     MatDenseGetSubMatrix(B, PETSC_DECIDE, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T);
334:     MatDenseGetSubMatrix(T4, PETSC_DECIDE, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T2);
335:     MatDenseGetSubMatrix(B2, PETSC_DECIDE, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T3);
336:     MatAXPY(T, 1.0, T2, SAME_NONZERO_PATTERN);
337:     MatAXPY(T3, -1.0, T, SAME_NONZERO_PATTERN);
338:     MatNorm(T3, NORM_FROBENIUS, &err);
339:     if (err > PETSC_SMALL) {
340:       PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetSubMatrix\n");
341:       MatView(T3, NULL);
342:     }
343:     MatDenseRestoreSubMatrix(B, &T);
344:     MatDenseRestoreSubMatrix(T4, &T2);
345:     MatDenseRestoreSubMatrix(B2, &T3);
346:     CheckLocal(B, NULL, aB, NULL);
347:     MatDestroy(&B2);
348:     MatDestroy(&T4);
349:     if (N >= 2) {
350:       MatDuplicate(B, MAT_COPY_VALUES, &B2);
351:       MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &T4);
352:       MatSetRandom(T4, NULL);
353:       MatAXPY(B2, 1.0, T4, SAME_NONZERO_PATTERN);
354:       MatDenseGetSubMatrix(B, N - 2, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T);
355:       MatDenseGetSubMatrix(T4, N - 2, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T2);
356:       MatDenseGetSubMatrix(B2, N - 2, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T3);
357:       MatAXPY(T, 1.0, T2, SAME_NONZERO_PATTERN);
358:       MatAXPY(T3, -1.0, T, SAME_NONZERO_PATTERN);
359:       MatNorm(T3, NORM_FROBENIUS, &err);
360:       if (err > PETSC_SMALL) {
361:         PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetSubMatrix\n");
362:         MatView(T3, NULL);
363:       }
364:       MatDenseRestoreSubMatrix(B, &T);
365:       MatDenseRestoreSubMatrix(T4, &T2);
366:       MatDenseRestoreSubMatrix(B2, &T3);
367:       CheckLocal(B, NULL, aB, NULL);
368:       MatDestroy(&B2);
369:       MatDestroy(&T4);
370:       MatDuplicate(B, MAT_COPY_VALUES, &B2);
371:       MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &T4);
372:       MatSetRandom(T4, NULL);
373:       MatAXPY(B2, 1.0, T4, SAME_NONZERO_PATTERN);
374:       MatDenseGetSubMatrix(B, PETSC_DECIDE, 2, PetscMin(1, K - 1), PetscMin(2, K), &T);
375:       MatDenseGetSubMatrix(T4, PETSC_DECIDE, 2, PetscMin(1, K - 1), PetscMin(2, K), &T2);
376:       MatDenseGetSubMatrix(B2, PETSC_DECIDE, 2, PetscMin(1, K - 1), PetscMin(2, K), &T3);
377:       MatAXPY(T, 1.0, T2, SAME_NONZERO_PATTERN);
378:       MatAXPY(T3, -1.0, T, SAME_NONZERO_PATTERN);
379:       MatNorm(T3, NORM_FROBENIUS, &err);
380:       if (err > PETSC_SMALL) {
381:         PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetSubMatrix\n");
382:         MatView(T3, NULL);
383:       }
384:       MatDenseRestoreSubMatrix(B, &T);
385:       MatDenseRestoreSubMatrix(T4, &T2);
386:       MatDenseRestoreSubMatrix(B2, &T3);
387:       CheckLocal(B, NULL, aB, NULL);
388:       MatDestroy(&B2);
389:       MatDestroy(&T4);
390:     }
391:   }

393:   /* Test reusing a previously allocated dense buffer */
394:   MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
395:   CheckLocal(B, X, aB, aX);
396:   MatMatMultEqual(A, B, X, 10, &flg);
397:   if (!flg) {
398:     PetscPrintf(PETSC_COMM_WORLD, "Error with reusage\n");
399:     MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
400:     MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN);
401:     MatView(T, NULL);
402:     MatDestroy(&T);
403:   }

405:   /* Test MatTransposeMat and MatMatTranspose */
406:   if (testmattmat) {
407:     MatTransposeMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
408:     CheckLocal(B, X, aB, aX);
409:     MatTransposeMatMultEqual(A, X, B, 10, &flg);
410:     if (!flg) {
411:       PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatTransposeMat)\n");
412:       MatTransposeMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &B);
413:       MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN);
414:       MatView(T, NULL);
415:       MatDestroy(&T);
416:     }
417:   }
418:   if (testmatmatt) {
419:     MatMatTransposeMult(A, Bt, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
420:     CheckLocal(Bt, X, aBt, aX);
421:     MatMatTransposeMultEqual(A, Bt, X, 10, &flg);
422:     if (!flg) {
423:       PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatMatTranspose)\n");
424:       MatMatTransposeMult(A, Bt, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
425:       MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN);
426:       MatView(T, NULL);
427:       MatDestroy(&T);
428:     }
429:   }

431:   /* Test projection operations (PtAP and RARt) */
432:   if (testproj) {
433:     MatPtAP(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &PtAP);
434:     MatPtAPMultEqual(A, B, PtAP, 10, &flg);
435:     if (!flg) {
436:       PetscPrintf(PETSC_COMM_WORLD, "Error with PtAP\n");
437:       MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
438:       MatTransposeMatMult(B, T, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T2);
439:       MatAXPY(T2, -1.0, PtAP, SAME_NONZERO_PATTERN);
440:       MatView(T2, NULL);
441:       MatDestroy(&T2);
442:       MatDestroy(&T);
443:     }
444:     PetscMalloc1((k + ldr) * M, &dataR);
445:     MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, m, K, M, dataR, &R);
446:     MatDenseSetLDA(R, k + ldr);
447:     MatSetRandom(R, NULL);
448:     if (testrart) { /* fails for AIJCUSPARSE because RA operation is not defined */
449:       MatRARt(A, R, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &RARt);
450:       MatRARtMultEqual(A, R, RARt, 10, &flg);
451:       if (!flg) {
452:         PetscPrintf(PETSC_COMM_WORLD, "Error with RARt\n");
453:         MatMatTransposeMult(A, R, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
454:         MatMatMult(R, T, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T2);
455:         MatAXPY(T2, -1.0, RARt, SAME_NONZERO_PATTERN);
456:         MatView(T2, NULL);
457:         MatDestroy(&T2);
458:         MatDestroy(&T);
459:       }
460:     }
461:   }

463:   /* Test MatDenseGetColumnVec and friends */
464:   MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
465:   MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
466:   MatDuplicate(T, MAT_DO_NOT_COPY_VALUES, &T2);
467:   for (k = 0; k < K; k++) {
468:     Vec Xv, Tv, T2v;

470:     MatDenseGetColumnVecRead(X, k, &Xv);
471:     MatDenseGetColumnVec(T, k, &Tv);
472:     MatDenseGetColumnVecWrite(T2, k, &T2v);
473:     VecCopy(Xv, T2v);
474:     VecAXPY(Tv, -1., Xv);
475:     MatDenseRestoreColumnVecRead(X, k, &Xv);
476:     MatDenseRestoreColumnVec(T, k, &Tv);
477:     MatDenseRestoreColumnVecWrite(T2, k, &T2v);
478:   }
479:   MatNorm(T, NORM_FROBENIUS, &err);
480:   if (err > PETSC_SMALL) {
481:     PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetColumnVec\n");
482:     MatView(T, NULL);
483:   }
484:   MatAXPY(T2, -1., X, SAME_NONZERO_PATTERN);
485:   MatNorm(T2, NORM_FROBENIUS, &err);
486:   if (err > PETSC_SMALL) {
487:     PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetColumnVecWrite\n");
488:     MatView(T2, NULL);
489:   }
490:   MatDestroy(&T);
491:   MatDestroy(&T2);

493:   /* Test with MatShell */
494:   MatDuplicate(A, MAT_COPY_VALUES, &T);
495:   MatConvert(T, MATSHELL, MAT_INITIAL_MATRIX, &T2);
496:   MatDestroy(&T);

498:   /* scale matrix */
499:   MatScale(A, 2.0);
500:   MatScale(T2, 2.0);
501:   MatCreateVecs(A, &r, &l);
502:   VecSetRandom(r, NULL);
503:   VecSetRandom(l, NULL);
504:   MatCreateVecs(T2, &rs, &ls);
505:   VecCopy(r, rs);
506:   VecCopy(l, ls);
507:   if (testproj) {
508:     MatDiagonalScale(A, r, r);
509:     MatDiagonalScale(T2, rs, rs);
510:   } else {
511:     MatDiagonalScale(A, l, r);
512:     MatDiagonalScale(T2, ls, rs);
513:   }
514:   MatDuplicate(A, MAT_COPY_VALUES, &T);
515:   MatAXPY(A, 4.5, T, SAME_NONZERO_PATTERN);
516:   MatAXPY(T2, 4.5, T, DIFFERENT_NONZERO_PATTERN);
517:   MatMultEqual(T2, A, 10, &flg);
518:   if (!flg) PetscPrintf(PETSC_COMM_WORLD, "Error with MATSHELL (MatMult)\n");
519:   MatMultTransposeEqual(T2, A, 10, &flg);
520:   if (!flg) PetscPrintf(PETSC_COMM_WORLD, "Error with MATSHELL (MatMultTranspose)\n");
521:   MatDestroy(&T);
522:   VecDestroy(&ls);
523:   VecDestroy(&rs);
524:   VecDestroy(&l);
525:   VecDestroy(&r);

527:   /* recompute projections, test reusage */
528:   if (PtAP) MatPtAP(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &PtAP);
529:   if (RARt) MatRARt(A, R, MAT_REUSE_MATRIX, PETSC_DEFAULT, &RARt);
530:   if (testshellops) { /* test callbacks for user defined MatProducts */
531:     MatShellSetMatProductOperation(T2, MATPRODUCT_AB, NULL, MyMatShellMatMultNumeric, NULL, MATDENSE, MATDENSE);
532:     MatShellSetMatProductOperation(T2, MATPRODUCT_AB, NULL, MyMatShellMatMultNumeric, NULL, MATDENSECUDA, MATDENSECUDA);
533:     MatShellSetMatProductOperation(T2, MATPRODUCT_AtB, NULL, MyMatTransposeShellMatMultNumeric, NULL, MATDENSE, MATDENSE);
534:     MatShellSetMatProductOperation(T2, MATPRODUCT_AtB, NULL, MyMatTransposeShellMatMultNumeric, NULL, MATDENSECUDA, MATDENSECUDA);
535:     MatShellSetMatProductOperation(T2, MATPRODUCT_ABt, NULL, MyMatShellMatTransposeMultNumeric, NULL, MATDENSE, MATDENSE);
536:     MatShellSetMatProductOperation(T2, MATPRODUCT_ABt, NULL, MyMatShellMatTransposeMultNumeric, NULL, MATDENSECUDA, MATDENSECUDA);
537:     if (testproj) {
538:       MatShellSetMatProductOperation(T2, MATPRODUCT_PtAP, MyPtShellPMultSymbolic, MyPtShellPMultNumeric, proj_destroy, MATDENSE, MATSHELL);
539:       MatShellSetMatProductOperation(T2, MATPRODUCT_PtAP, MyPtShellPMultSymbolic, MyPtShellPMultNumeric, proj_destroy, MATDENSECUDA, MATSHELL);
540:       MatShellSetMatProductOperation(T2, MATPRODUCT_RARt, MyRShellRtMultSymbolic, MyRShellRtMultNumeric, proj_destroy, MATDENSE, MATSHELL);
541:       MatShellSetMatProductOperation(T2, MATPRODUCT_RARt, MyRShellRtMultSymbolic, MyRShellRtMultNumeric, proj_destroy, MATDENSECUDA, MATSHELL);
542:     }
543:   }
544:   CheckLocal(B, X, aB, aX);
545:   /* we either use the shell operations or the loop over columns code, applying the operator */
546:   MatMatMult(T2, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
547:   CheckLocal(B, X, aB, aX);
548:   MatMatMultEqual(T2, B, X, 10, &flg);
549:   if (!flg) {
550:     PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MATSHELL)\n");
551:     MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
552:     MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN);
553:     MatView(T, NULL);
554:     MatDestroy(&T);
555:   }
556:   if (testproj) {
557:     MatPtAPMultEqual(T2, B, PtAP, 10, &flg);
558:     if (!flg) PetscPrintf(PETSC_COMM_WORLD, "Error with PtAP (MATSHELL)\n");
559:     if (testshellops) { /* projections fail if the product operations are not specified */
560:       MatPtAP(T2, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
561:       MatPtAP(T2, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &T);
562:       MatPtAPMultEqual(T2, B, T, 10, &flg);
563:       if (!flg) {
564:         Mat TE;

566:         PetscPrintf(PETSC_COMM_WORLD, "Error with PtAP (MATSHELL user defined)\n");
567:         MatComputeOperator(T, MATDENSE, &TE);
568:         MatView(TE, NULL);
569:         MatView(PtAP, NULL);
570:         MatAXPY(TE, -1.0, PtAP, SAME_NONZERO_PATTERN);
571:         MatView(TE, NULL);
572:         MatDestroy(&TE);
573:       }
574:       MatDestroy(&T);
575:     }
576:     if (RARt) {
577:       MatRARtMultEqual(T2, R, RARt, 10, &flg);
578:       if (!flg) PetscPrintf(PETSC_COMM_WORLD, "Error with RARt (MATSHELL)\n");
579:     }
580:     if (testshellops) {
581:       MatRARt(T2, R, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
582:       MatRARt(T2, R, MAT_REUSE_MATRIX, PETSC_DEFAULT, &T);
583:       MatRARtMultEqual(T2, R, T, 10, &flg);
584:       if (!flg) {
585:         Mat TE;

587:         PetscPrintf(PETSC_COMM_WORLD, "Error with RARt (MATSHELL user defined)\n");
588:         MatComputeOperator(T, MATDENSE, &TE);
589:         MatView(TE, NULL);
590:         if (RARt) {
591:           MatView(RARt, NULL);
592:           MatAXPY(TE, -1.0, RARt, SAME_NONZERO_PATTERN);
593:           MatView(TE, NULL);
594:         }
595:         MatDestroy(&TE);
596:       }
597:       MatDestroy(&T);
598:     }
599:   }

601:   if (testmattmat) { /* we either use the shell operations or the loop over columns code applying the transposed operator */
602:     MatTransposeMatMult(T2, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
603:     CheckLocal(B, X, aB, aX);
604:     MatTransposeMatMultEqual(T2, X, B, 10, &flg);
605:     if (!flg) {
606:       PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatTranspose, MATSHELL)\n");
607:       MatTransposeMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
608:       MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN);
609:       MatView(T, NULL);
610:       MatDestroy(&T);
611:     }
612:   }
613:   if (testmatmatt && testshellops) { /* only when shell operations are set */
614:     MatMatTransposeMult(T2, Bt, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
615:     CheckLocal(Bt, X, aBt, aX);
616:     MatMatTransposeMultEqual(T2, Bt, X, 10, &flg);
617:     if (!flg) {
618:       PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatMatTranspose, MATSHELL)\n");
619:       MatMatTransposeMult(A, Bt, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
620:       MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN);
621:       MatView(T, NULL);
622:       MatDestroy(&T);
623:     }
624:   }
625:   MatDestroy(&T2);

627:   if (testnest) { /* test with MatNest */
628:     Mat NA;

630:     MatCreateNest(PETSC_COMM_WORLD, 1, NULL, 1, NULL, &A, &NA);
631:     MatViewFromOptions(NA, NULL, "-NA_view");
632:     MatMatMult(NA, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
633:     CheckLocal(B, X, aB, aX);
634:     MatMatMultEqual(NA, B, X, 10, &flg);
635:     if (!flg) {
636:       PetscPrintf(PETSC_COMM_WORLD, "Error with Nest\n");
637:       MatMatMult(NA, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
638:       MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN);
639:       MatView(T, NULL);
640:       MatDestroy(&T);
641:     }
642:     MatDestroy(&NA);
643:   }

645:   if (testtranspose) { /* test with Transpose */
646:     Mat TA;

648:     MatCreateTranspose(A, &TA);
649:     MatMatMult(TA, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
650:     CheckLocal(B, X, aB, aX);
651:     MatMatMultEqual(TA, X, B, 10, &flg);
652:     if (!flg) {
653:       PetscPrintf(PETSC_COMM_WORLD, "Error with Transpose\n");
654:       MatMatMult(TA, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
655:       MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN);
656:       MatView(T, NULL);
657:       MatDestroy(&T);
658:     }
659:     MatDestroy(&TA);
660:   }

662:   if (testhtranspose) { /* test with Hermitian Transpose */
663:     Mat TA;

665:     MatCreateHermitianTranspose(A, &TA);
666:     MatMatMult(TA, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
667:     CheckLocal(B, X, aB, aX);
668:     MatMatMultEqual(TA, X, B, 10, &flg);
669:     if (!flg) {
670:       PetscPrintf(PETSC_COMM_WORLD, "Error with Transpose\n");
671:       MatMatMult(TA, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
672:       MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN);
673:       MatView(T, NULL);
674:       MatDestroy(&T);
675:     }
676:     MatDestroy(&TA);
677:   }

679:   if (testtt) { /* test with Transpose(Transpose) */
680:     Mat TA, TTA;

682:     MatCreateTranspose(A, &TA);
683:     MatCreateTranspose(TA, &TTA);
684:     MatMatMult(TTA, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
685:     CheckLocal(B, X, aB, aX);
686:     MatMatMultEqual(TTA, B, X, 10, &flg);
687:     if (!flg) {
688:       PetscPrintf(PETSC_COMM_WORLD, "Error with Transpose(Transpose)\n");
689:       MatMatMult(TTA, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
690:       MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN);
691:       MatView(T, NULL);
692:       MatDestroy(&T);
693:     }
694:     MatDestroy(&TA);
695:     MatDestroy(&TTA);
696:   }

698:   if (testcircular) { /* test circular */
699:     Mat AB;

701:     MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AB);
702:     MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
703:     CheckLocal(B, X, aB, aX);
704:     if (M == N && N == K) {
705:       MatMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
706:     } else {
707:       MatTransposeMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
708:     }
709:     CheckLocal(B, X, aB, aX);
710:     MatDestroy(&AB);
711:   }

713:   /* Test by Pierre Jolivet */
714:   {
715:     Mat C, D, D2, AtA;
716:     MatCreateNormal(A, &AtA);
717:     MatDuplicate(X, MAT_DO_NOT_COPY_VALUES, &C);
718:     MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &D);
719:     MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &D2);
720:     MatSetRandom(B, NULL);
721:     MatSetRandom(C, NULL);
722:     MatSetRandom(D, NULL);
723:     MatSetRandom(D2, NULL);
724:     MatProductCreateWithMat(A, B, NULL, C);
725:     MatProductSetType(C, MATPRODUCT_AB);
726:     MatProductSetFromOptions(C);
727:     MatProductSymbolic(C);
728:     MatProductCreateWithMat(A, C, NULL, D);
729:     MatProductSetType(D, MATPRODUCT_AtB);
730:     MatProductSetFromOptions(D);
731:     MatProductSymbolic(D);
732:     MatProductNumeric(C);
733:     MatMatMultEqual(A, B, C, 10, &flg);
734:     if (!flg) {
735:       PetscPrintf(PETSC_COMM_WORLD, "Error with Normal (AB != C)\n");
736:       MatView(A, NULL);
737:       MatView(B, NULL);
738:       MatView(C, NULL);
739:     }
740:     MatProductNumeric(D);
741:     MatMatMultEqual(AtA, B, D, 10, &flg);
742:     if (!flg) {
743:       PetscPrintf(PETSC_COMM_WORLD, "Error with Normal (2)\n");
744:       MatMatMult(AtA, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
745:       MatView(D, NULL);
746:       MatView(T, NULL);
747:       MatAXPY(T, -1.0, D, SAME_NONZERO_PATTERN);
748:       MatView(T, NULL);
749:       MatDestroy(&T);
750:     }
751:     MatDestroy(&C);
752:     MatDestroy(&D);
753:     MatDestroy(&D2);
754:     MatDestroy(&AtA);
755:   }

757:   MatDestroy(&X);
758:   MatDestroy(&Bt);
759:   MatDestroy(&B);
760:   MatDestroy(&A);
761:   MatDestroy(&R);
762:   MatDestroy(&PtAP);
763:   MatDestroy(&RARt);
764:   PetscFree(dataX);
765:   PetscFree(dataB);
766:   PetscFree(dataR);
767:   PetscFree(dataBt);
768:   PetscFinalize();
769:   return 0;
770: }

772: /*TEST

774:   test:
775:     suffix: 1
776:     args: -local {{0 1}} -testshellops

778:   test:
779:     output_file: output/ex70_1.out
780:     requires: cuda
781:     suffix: 1_cuda
782:     args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{seqaijcusparse seqaij}} -testshellops {{0 1}}

784:   test:
785:     output_file: output/ex70_1.out
786:     nsize: 2
787:     suffix: 1_par
788:     args: -local {{0 1}} -testmatmatt 0

790:   test:
791:     output_file: output/ex70_1.out
792:     requires: cuda
793:     nsize: 2
794:     suffix: 1_par_cuda
795:     args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{mpiaijcusparse mpiaij}} -testnest 0 -testmatmatt 0 -matmatmult_Bbn 3

797:   test:
798:     output_file: output/ex70_1.out
799:     suffix: 2
800:     nsize: 1
801:     args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}}

803:   testset:
804:     requires: cuda
805:     output_file: output/ex70_1.out
806:     nsize: 1
807:     args: -M 7 -N 9 -K 2 -local {{0 1}} -testnest 0 -A_mat_type {{seqdensecuda seqdense}} -xgpu {{0 1}} -bgpu {{0 1}}
808:     test:
809:       requires: !complex
810:       suffix: 2_cuda_real
811:     test:
812:       # complex+single gives a little bigger error in the MatDenseGetColumnVec test
813:       requires: complex !single
814:       suffix: 2_cuda_complex

816:   test:
817:     output_file: output/ex70_1.out
818:     suffix: 2_par
819:     nsize: 2
820:     args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} -testcircular -testmatmatt 0

822:   test:
823:     requires: cuda
824:     output_file: output/ex70_1.out
825:     suffix: 2_par_cuda
826:     nsize: 2
827:     args: -M 11 -N 9 -K 1 -local {{0 1}} -testcircular 0 -A_mat_type mpiaijcusparse -xgpu -bgpu -testnest 0 -testmatmatt 0

829:   test:
830:     output_file: output/ex70_1.out
831:     suffix: 3
832:     nsize: {{1 3}}
833:     args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type sbaij -symm -testproj 0 -testmatmatt 0

835:   test:
836:     output_file: output/ex70_1.out
837:     suffix: 4
838:     nsize: 1
839:     args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular

841:   test:
842:     output_file: output/ex70_1.out
843:     suffix: 5
844:     nsize: {{2 4}}
845:     args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular -testmatmatt 0

847:   test:
848:     output_file: output/ex70_1.out
849:     suffix: 6
850:     nsize: 1
851:     args: -M {{1 3}} -N {{2 5}} -K {{1 2}} -local {{0 1}} -testcircular

853:   test:
854:     output_file: output/ex70_1.out
855:     suffix: 7
856:     nsize: 1
857:     args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type dense -testnest -testcircular

859: TEST*/