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:   PetscFunctionBegin;
 14:   if (a) {
 15:     const PetscScalar *Aa;
 16:     PetscCall(MatDenseGetArrayRead(A, &Aa));
 17:     wA = (PetscBool)(a != Aa);
 18:     PetscCall(MatDenseGetLDA(A, &lda));
 19:     PetscCall(MatGetLocalSize(A, &m, &n));
 20:     for (j = 0; j < n; j++) {
 21:       for (i = m; i < lda; i++) {
 22:         if (Aa[j * lda + i] != MAGIC_NUMBER) wAv = PETSC_TRUE;
 23:       }
 24:     }
 25:     PetscCall(MatDenseRestoreArrayRead(A, &Aa));
 26:   }
 27:   if (b) {
 28:     const PetscScalar *Bb;
 29:     PetscCall(MatDenseGetArrayRead(B, &Bb));
 30:     wB = (PetscBool)(b != Bb);
 31:     PetscCall(MatDenseGetLDA(B, &lda));
 32:     PetscCall(MatGetLocalSize(B, &m, &n));
 33:     for (j = 0; j < n; j++) {
 34:       for (i = m; i < lda; i++) {
 35:         if (Bb[j * lda + i] != MAGIC_NUMBER) wBv = PETSC_TRUE;
 36:       }
 37:     }
 38:     PetscCall(MatDenseRestoreArrayRead(B, &Bb));
 39:   }
 40:   PetscCheck(!wA && !wB, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong array in first Mat? %d, Wrong array in second Mat? %d", wA, wB);
 41:   PetscCheck(!wAv && !wBv, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Wrong data in first Mat? %d, Wrong data in second Mat? %d", wAv, wBv);
 42:   PetscFunctionReturn(PETSC_SUCCESS);
 43: }

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

 51: PetscErrorCode proj_destroy(PetscCtxRt ctx)
 52: {
 53:   proj_data *userdata = *(proj_data **)ctx;

 55:   PetscFunctionBegin;
 56:   PetscCheck(userdata, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing userdata");
 57:   PetscCall(MatDestroy(&userdata->A));
 58:   PetscCall(MatDestroy(&userdata->P));
 59:   PetscCall(MatDestroy(&userdata->R));
 60:   PetscCall(PetscFree(userdata));
 61:   PetscFunctionReturn(PETSC_SUCCESS);
 62: }

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

 71:   PetscFunctionBegin;
 72:   PetscCall(MatShellGetContext(S, &userdata));
 73:   PetscCheck(userdata, PetscObjectComm((PetscObject)S), PETSC_ERR_PLIB, "Missing userdata");
 74:   A = userdata->A;
 75:   R = userdata->R;
 76:   P = userdata->P;
 77:   PetscCheck(A, PetscObjectComm((PetscObject)S), PETSC_ERR_PLIB, "Missing matrix");
 78:   PetscCheck(R || P, PetscObjectComm((PetscObject)S), PETSC_ERR_PLIB, "Missing projectors");
 79:   PetscCheck(!R || !P, PetscObjectComm((PetscObject)S), PETSC_ERR_PLIB, "Both projectors");
 80:   PetscCall(MatCreateVecs(A, &Ax, &Ay));
 81:   if (R) PetscCall(MatCreateVecs(R, &Py, &Px));
 82:   else PetscCall(MatCreateVecs(P, &Px, &Py));
 83:   PetscCall(VecCopy(X, Px));
 84:   if (P) PetscCall(MatMult(P, Px, Py));
 85:   else PetscCall(MatMultTranspose(R, Px, Py));
 86:   PetscCall(VecCopy(Py, Ax));
 87:   PetscCall(MatMult(A, Ax, Ay));
 88:   PetscCall(VecCopy(Ay, Py));
 89:   if (P) PetscCall(MatMultTranspose(P, Py, Px));
 90:   else PetscCall(MatMult(R, Py, Px));
 91:   PetscCall(VecCopy(Px, Y));
 92:   PetscCall(VecDestroy(&Px));
 93:   PetscCall(VecDestroy(&Py));
 94:   PetscCall(VecDestroy(&Ax));
 95:   PetscCall(VecDestroy(&Ay));
 96:   PetscFunctionReturn(PETSC_SUCCESS);
 97: }

 99: PetscErrorCode MyPtShellPMultSymbolic(Mat S, Mat P, Mat PtAP, void **ctx)
100: {
101:   proj_data *userdata;

103:   PetscFunctionBegin;
104:   PetscCall(PetscNew(&userdata));
105:   PetscCall(MatShellSetContext(PtAP, userdata));
106:   *ctx = (void *)userdata;
107:   PetscFunctionReturn(PETSC_SUCCESS);
108: }

110: PetscErrorCode MyPtShellPMultNumeric(Mat S, Mat P, Mat PtAP, PetscCtx ctx)
111: {
112:   Mat        A;
113:   proj_data *userdata = (proj_data *)ctx;

115:   PetscFunctionBegin;
116:   PetscCall(MatShellGetContext(S, &A));
117:   PetscCall(PetscObjectReference((PetscObject)A));
118:   PetscCall(PetscObjectReference((PetscObject)P));
119:   PetscCall(MatDestroy(&userdata->A));
120:   PetscCall(MatDestroy(&userdata->P));
121:   PetscCall(MatDestroy(&userdata->R));
122:   userdata->A = A;
123:   userdata->P = P;
124:   PetscCall(MatShellSetOperation(PtAP, MATOP_MULT, (PetscErrorCodeFn *)proj_mult));
125:   PetscCall(MatSetUp(PtAP));
126:   PetscCall(MatAssemblyBegin(PtAP, MAT_FINAL_ASSEMBLY));
127:   PetscCall(MatAssemblyEnd(PtAP, MAT_FINAL_ASSEMBLY));
128:   PetscFunctionReturn(PETSC_SUCCESS);
129: }

131: PetscErrorCode MyRShellRtMultSymbolic(Mat S, Mat R, Mat RARt, void **ctx)
132: {
133:   proj_data *userdata;

135:   PetscFunctionBegin;
136:   PetscCall(PetscNew(&userdata));
137:   PetscCall(MatShellSetContext(RARt, userdata));
138:   *ctx = (void *)userdata;
139:   PetscFunctionReturn(PETSC_SUCCESS);
140: }

142: PetscErrorCode MyRShellRtMultNumeric(Mat S, Mat R, Mat RARt, PetscCtx ctx)
143: {
144:   Mat        A;
145:   proj_data *userdata = (proj_data *)ctx;

147:   PetscFunctionBegin;
148:   PetscCall(MatShellGetContext(S, &A));
149:   PetscCall(PetscObjectReference((PetscObject)A));
150:   PetscCall(PetscObjectReference((PetscObject)R));
151:   PetscCall(MatDestroy(&userdata->A));
152:   PetscCall(MatDestroy(&userdata->P));
153:   PetscCall(MatDestroy(&userdata->R));
154:   userdata->A = A;
155:   userdata->R = R;
156:   PetscCall(MatShellSetOperation(RARt, MATOP_MULT, (PetscErrorCodeFn *)proj_mult));
157:   PetscCall(MatSetUp(RARt));
158:   PetscCall(MatAssemblyBegin(RARt, MAT_FINAL_ASSEMBLY));
159:   PetscCall(MatAssemblyEnd(RARt, MAT_FINAL_ASSEMBLY));
160:   PetscFunctionReturn(PETSC_SUCCESS);
161: }

163: PetscErrorCode MyMatShellMatMultNumeric(Mat S, Mat B, Mat C, PetscCtx ctx)
164: {
165:   Mat A;

167:   PetscFunctionBegin;
168:   PetscCall(MatShellGetContext(S, &A));
169:   PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &C));
170:   PetscFunctionReturn(PETSC_SUCCESS);
171: }

173: PetscErrorCode MyMatTransposeShellMatMultNumeric(Mat S, Mat B, Mat C, PetscCtx ctx)
174: {
175:   Mat A;

177:   PetscFunctionBegin;
178:   PetscCall(MatShellGetContext(S, &A));
179:   PetscCall(MatTransposeMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &C));
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

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

187:   PetscFunctionBegin;
188:   PetscCall(MatShellGetContext(S, &A));
189:   PetscCall(MatMatTransposeMult(A, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &C));
190:   PetscFunctionReturn(PETSC_SUCCESS);
191: }

193: int main(int argc, char **args)
194: {
195:   Mat          X, B, A, Bt, T, T2, PtAP = NULL, RARt = NULL, R = NULL;
196:   Vec          r, l, rs, ls;
197:   PetscInt     m, n, k, M = 10, N = 10, K = 5, ldx = 3, ldb = 5, ldr = 4;
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, formt = PETSC_FALSE, testreusemodified = PETSC_FALSE;
202:   PetscScalar *dataX = NULL, *dataB = NULL, *dataR = NULL, *dataBt = NULL;
203:   PetscScalar *aX, *aB, *aBt;
204:   PetscReal    err;

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

232:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
233:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
234:   PetscCall(MatSetType(A, MATAIJ));
235:   PetscCall(MatSeqAIJSetPreallocation(A, PETSC_DEFAULT, NULL));
236:   PetscCall(MatMPIAIJSetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
237:   PetscCall(MatSetUp(A));
238:   PetscCall(MatSetRandom(A, NULL));
239:   if (M == N && symm) {
240:     Mat AT;

242:     PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
243:     PetscCall(MatAXPY(A, 1.0, AT, DIFFERENT_NONZERO_PATTERN));
244:     PetscCall(MatDestroy(&AT));
245:     PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
246:   }
247:   PetscCall(MatViewFromOptions(A, NULL, "-A_init_view"));
248:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "");
249:   PetscCall(PetscOptionsFList("-A_mat_type", "Matrix type", "MatSetType", MatList, MATAIJ, mattype, sizeof(mattype), &flg));
250:   PetscOptionsEnd();
251:   if (flg) {
252:     Mat A2;

254:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
255:     PetscCall(MatConvert(A, mattype, MAT_INPLACE_MATRIX, &A));
256:     PetscCall(MatMultEqual(A, A2, 10, &flg));
257:     if (!flg) {
258:       Mat AE, A2E;

260:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with convert\n"));
261:       PetscCall(MatComputeOperator(A, MATDENSE, &AE));
262:       PetscCall(MatComputeOperator(A2, MATDENSE, &A2E));
263:       PetscCall(MatView(AE, NULL));
264:       PetscCall(MatView(A2E, NULL));
265:       PetscCall(MatAXPY(A2E, -1.0, A, SAME_NONZERO_PATTERN));
266:       PetscCall(MatView(A2E, NULL));
267:       PetscCall(MatDestroy(&A2E));
268:       PetscCall(MatDestroy(&AE));
269:     }
270:     PetscCall(MatDestroy(&A2));
271:   }
272:   if (formt) PetscCall(MatSetOption(A, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE));
273:   PetscCall(MatViewFromOptions(A, NULL, "-A_view"));

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

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

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

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

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

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

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

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

432:   /* Test that MAT_REUSE_MATRIX computes with the current values of A after they have been changed */
433:   if (testreusemodified) {
434:     Mat E;

436:     PetscCall(MatTransposeMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &E));
437:     PetscCall(MatScale(A, 2.0));
438:     PetscCall(MatTransposeMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DETERMINE, &E));
439:     PetscCall(MatTransposeMatMultEqual(A, X, E, 10, &flg));
440:     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with reusage after MatScale() (MatTransposeMat)\n"));
441:     PetscCall(MatScale(A, 0.5));
442:     PetscCall(MatDestroy(&E));
443:   }

445:   /* Test projection operations (PtAP and RARt) */
446:   if (testproj) {
447:     PetscCall(MatPtAP(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &PtAP));
448:     PetscCall(MatPtAPMultEqual(A, B, PtAP, 10, &flg));
449:     if (!flg) {
450:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with PtAP\n"));
451:       PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T));
452:       PetscCall(MatTransposeMatMult(B, T, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T2));
453:       PetscCall(MatAXPY(T2, -1.0, PtAP, SAME_NONZERO_PATTERN));
454:       PetscCall(MatView(T2, NULL));
455:       PetscCall(MatDestroy(&T2));
456:       PetscCall(MatDestroy(&T));
457:     }
458:     PetscCall(PetscMalloc1((k + ldr) * M, &dataR));
459:     PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, m, K, M, dataR, &R));
460:     PetscCall(MatDenseSetLDA(R, k + ldr));
461:     PetscCall(MatSetRandom(R, NULL));
462:     if (testrart) { /* fails for AIJCUSPARSE because RA operation is not defined */
463:       PetscCall(MatRARt(A, R, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &RARt));
464:       PetscCall(MatRARtMultEqual(A, R, RARt, 10, &flg));
465:       if (!flg) {
466:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with RARt\n"));
467:         PetscCall(MatMatTransposeMult(A, R, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T));
468:         PetscCall(MatMatMult(R, T, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T2));
469:         PetscCall(MatAXPY(T2, -1.0, RARt, SAME_NONZERO_PATTERN));
470:         PetscCall(MatView(T2, NULL));
471:         PetscCall(MatDestroy(&T2));
472:         PetscCall(MatDestroy(&T));
473:       }
474:     }
475:   }

477:   /* Test MatDenseGetColumnVec and friends */
478:   PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &X));
479:   PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T));
480:   PetscCall(MatDuplicate(T, MAT_DO_NOT_COPY_VALUES, &T2));
481:   for (k = 0; k < K; k++) {
482:     Vec Xv, Tv, T2v;

484:     PetscCall(MatDenseGetColumnVecRead(X, k, &Xv));
485:     PetscCall(MatDenseGetColumnVec(T, k, &Tv));
486:     PetscCall(MatDenseGetColumnVecWrite(T2, k, &T2v));
487:     PetscCall(VecCopy(Xv, T2v));
488:     PetscCall(VecAXPY(Tv, -1., Xv));
489:     PetscCall(MatDenseRestoreColumnVecRead(X, k, &Xv));
490:     PetscCall(MatDenseRestoreColumnVec(T, k, &Tv));
491:     PetscCall(MatDenseRestoreColumnVecWrite(T2, k, &T2v));
492:   }
493:   PetscCall(MatNorm(T, NORM_FROBENIUS, &err));
494:   if (err > PETSC_SMALL) {
495:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetColumnVec\n"));
496:     PetscCall(MatView(T, NULL));
497:   }
498:   PetscCall(MatAXPY(T2, -1., X, SAME_NONZERO_PATTERN));
499:   PetscCall(MatNorm(T2, NORM_FROBENIUS, &err));
500:   if (err > PETSC_SMALL) {
501:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetColumnVecWrite\n"));
502:     PetscCall(MatView(T2, NULL));
503:   }
504:   PetscCall(MatDestroy(&T));
505:   PetscCall(MatDestroy(&T2));

507:   /* Test with MatShell */
508:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &T));
509:   PetscCall(MatConvert(T, MATSHELL, MAT_INITIAL_MATRIX, &T2));
510:   PetscCall(MatDestroy(&T));

512:   /* scale matrix */
513:   PetscCall(MatScale(A, 2.0));
514:   PetscCall(MatScale(T2, 2.0));
515:   PetscCall(MatCreateVecs(A, &r, &l));
516:   PetscCall(VecSetRandom(r, NULL));
517:   PetscCall(VecSetRandom(l, NULL));
518:   PetscCall(MatCreateVecs(T2, &rs, &ls));
519:   PetscCall(VecCopy(r, rs));
520:   PetscCall(VecCopy(l, ls));
521:   if (testproj) {
522:     PetscCall(MatDiagonalScale(A, r, r));
523:     PetscCall(MatDiagonalScale(T2, rs, rs));
524:   } else {
525:     PetscCall(MatDiagonalScale(A, l, r));
526:     PetscCall(MatDiagonalScale(T2, ls, rs));
527:   }
528:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &T));
529:   PetscCall(MatAXPY(A, 4.5, T, SAME_NONZERO_PATTERN));
530:   PetscCall(MatAXPY(T2, 4.5, T, DIFFERENT_NONZERO_PATTERN));
531:   PetscCall(MatMultEqual(T2, A, 10, &flg));
532:   if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with MATSHELL (MatMult)\n"));
533:   PetscCall(MatMultTransposeEqual(T2, A, 10, &flg));
534:   if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with MATSHELL (MatMultTranspose)\n"));
535:   PetscCall(MatDestroy(&T));
536:   PetscCall(VecDestroy(&ls));
537:   PetscCall(VecDestroy(&rs));
538:   PetscCall(VecDestroy(&l));
539:   PetscCall(VecDestroy(&r));

541:   /* recompute projections, test reusage */
542:   if (PtAP) PetscCall(MatPtAP(A, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &PtAP));
543:   if (RARt) PetscCall(MatRARt(A, R, MAT_REUSE_MATRIX, PETSC_DETERMINE, &RARt));
544:   if (testshellops) { /* test callbacks for user defined MatProducts */
545:     PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_AB, NULL, MyMatShellMatMultNumeric, NULL, MATDENSE, MATDENSE));
546:     PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_AB, NULL, MyMatShellMatMultNumeric, NULL, MATDENSECUDA, MATDENSECUDA));
547:     PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_AtB, NULL, MyMatTransposeShellMatMultNumeric, NULL, MATDENSE, MATDENSE));
548:     PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_AtB, NULL, MyMatTransposeShellMatMultNumeric, NULL, MATDENSECUDA, MATDENSECUDA));
549:     PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_ABt, NULL, MyMatShellMatTransposeMultNumeric, NULL, MATDENSE, MATDENSE));
550:     PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_ABt, NULL, MyMatShellMatTransposeMultNumeric, NULL, MATDENSECUDA, MATDENSECUDA));
551:     if (testproj) {
552:       PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_PtAP, MyPtShellPMultSymbolic, MyPtShellPMultNumeric, proj_destroy, MATDENSE, MATSHELL));
553:       PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_PtAP, MyPtShellPMultSymbolic, MyPtShellPMultNumeric, proj_destroy, MATDENSECUDA, MATSHELL));
554:       PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_RARt, MyRShellRtMultSymbolic, MyRShellRtMultNumeric, proj_destroy, MATDENSE, MATSHELL));
555:       PetscCall(MatShellSetMatProductOperation(T2, MATPRODUCT_RARt, MyRShellRtMultSymbolic, MyRShellRtMultNumeric, proj_destroy, MATDENSECUDA, MATSHELL));
556:     }
557:   }
558:   PetscCall(CheckLocal(B, X, aB, aX));
559:   /* we either use the shell operations or the loop over columns code, applying the operator */
560:   PetscCall(MatMatMult(T2, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &X));
561:   PetscCall(CheckLocal(B, X, aB, aX));
562:   PetscCall(MatMatMultEqual(T2, B, X, 10, &flg));
563:   if (!flg) {
564:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MATSHELL)\n"));
565:     PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T));
566:     PetscCall(MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN));
567:     PetscCall(MatView(T, NULL));
568:     PetscCall(MatDestroy(&T));
569:   }
570:   if (testproj) {
571:     PetscCall(MatPtAPMultEqual(T2, B, PtAP, 10, &flg));
572:     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with PtAP (MATSHELL)\n"));
573:     if (testshellops) { /* projections fail if the product operations are not specified */
574:       PetscCall(MatPtAP(T2, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T));
575:       PetscCall(MatPtAP(T2, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &T));
576:       PetscCall(MatPtAPMultEqual(T2, B, T, 10, &flg));
577:       if (!flg) {
578:         Mat TE;

580:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with PtAP (MATSHELL user defined)\n"));
581:         PetscCall(MatComputeOperator(T, MATDENSE, &TE));
582:         PetscCall(MatView(TE, NULL));
583:         PetscCall(MatView(PtAP, NULL));
584:         PetscCall(MatAXPY(TE, -1.0, PtAP, SAME_NONZERO_PATTERN));
585:         PetscCall(MatView(TE, NULL));
586:         PetscCall(MatDestroy(&TE));
587:       }
588:       PetscCall(MatDestroy(&T));
589:     }
590:     if (RARt) {
591:       PetscCall(MatRARtMultEqual(T2, R, RARt, 10, &flg));
592:       if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with RARt (MATSHELL)\n"));
593:     }
594:     if (testshellops) {
595:       PetscCall(MatRARt(T2, R, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T));
596:       PetscCall(MatRARt(T2, R, MAT_REUSE_MATRIX, PETSC_DETERMINE, &T));
597:       PetscCall(MatRARtMultEqual(T2, R, T, 10, &flg));
598:       if (!flg) {
599:         Mat TE;

601:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with RARt (MATSHELL user defined)\n"));
602:         PetscCall(MatComputeOperator(T, MATDENSE, &TE));
603:         PetscCall(MatView(TE, NULL));
604:         if (RARt) {
605:           PetscCall(MatView(RARt, NULL));
606:           PetscCall(MatAXPY(TE, -1.0, RARt, SAME_NONZERO_PATTERN));
607:           PetscCall(MatView(TE, NULL));
608:         }
609:         PetscCall(MatDestroy(&TE));
610:       }
611:       PetscCall(MatDestroy(&T));
612:     }
613:   }

615:   if (testmattmat) { /* we either use the shell operations or the loop over columns code applying the transposed operator */
616:     PetscCall(MatTransposeMatMult(T2, X, MAT_REUSE_MATRIX, PETSC_DETERMINE, &B));
617:     PetscCall(CheckLocal(B, X, aB, aX));
618:     PetscCall(MatTransposeMatMultEqual(T2, X, B, 10, &flg));
619:     if (!flg) {
620:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatTranspose, MATSHELL)\n"));
621:       PetscCall(MatTransposeMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T));
622:       PetscCall(MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN));
623:       PetscCall(MatView(T, NULL));
624:       PetscCall(MatDestroy(&T));
625:     }
626:   }
627:   if (testmatmatt && testshellops) { /* only when shell operations are set */
628:     PetscCall(MatMatTransposeMult(T2, Bt, MAT_REUSE_MATRIX, PETSC_DETERMINE, &X));
629:     PetscCall(CheckLocal(Bt, X, aBt, aX));
630:     PetscCall(MatMatTransposeMultEqual(T2, Bt, X, 10, &flg));
631:     if (!flg) {
632:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatMatTranspose, MATSHELL)\n"));
633:       PetscCall(MatMatTransposeMult(A, Bt, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T));
634:       PetscCall(MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN));
635:       PetscCall(MatView(T, NULL));
636:       PetscCall(MatDestroy(&T));
637:     }
638:   }
639:   PetscCall(MatDestroy(&T2));

641:   if (testnest) { /* test with MatNest */
642:     Mat NA;

644:     PetscCall(MatCreateNest(PETSC_COMM_WORLD, 1, NULL, 1, NULL, &A, &NA));
645:     PetscCall(MatViewFromOptions(NA, NULL, "-NA_view"));
646:     PetscCall(MatMatMult(NA, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &X));
647:     PetscCall(CheckLocal(B, X, aB, aX));
648:     PetscCall(MatMatMultEqual(NA, B, X, 10, &flg));
649:     if (!flg) {
650:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with Nest\n"));
651:       PetscCall(MatMatMult(NA, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T));
652:       PetscCall(MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN));
653:       PetscCall(MatView(T, NULL));
654:       PetscCall(MatDestroy(&T));
655:     }
656:     PetscCall(MatDestroy(&NA));
657:   }

659:   if (testtranspose) { /* test with Transpose */
660:     Mat TA;

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

676:   if (testhtranspose) { /* test with Hermitian Transpose */
677:     Mat TA;

679:     PetscCall(MatCreateHermitianTranspose(A, &TA));
680:     PetscCall(MatMatMult(TA, X, MAT_REUSE_MATRIX, PETSC_DETERMINE, &B));
681:     PetscCall(CheckLocal(B, X, aB, aX));
682:     PetscCall(MatMatMultEqual(TA, X, B, 10, &flg));
683:     if (!flg) {
684:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with Transpose\n"));
685:       PetscCall(MatMatMult(TA, X, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T));
686:       PetscCall(MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN));
687:       PetscCall(MatView(T, NULL));
688:       PetscCall(MatDestroy(&T));
689:     }
690:     PetscCall(MatDestroy(&TA));
691:   }

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

696:     PetscCall(MatCreateTranspose(A, &TA));
697:     PetscCall(MatCreateTranspose(TA, &TTA));
698:     PetscCall(MatMatMult(TTA, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &X));
699:     PetscCall(CheckLocal(B, X, aB, aX));
700:     PetscCall(MatMatMultEqual(TTA, B, X, 10, &flg));
701:     if (!flg) {
702:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with Transpose(Transpose)\n"));
703:       PetscCall(MatMatMult(TTA, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T));
704:       PetscCall(MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN));
705:       PetscCall(MatView(T, NULL));
706:       PetscCall(MatDestroy(&T));
707:     }
708:     PetscCall(MatDestroy(&TA));
709:     PetscCall(MatDestroy(&TTA));
710:   }

712:   if (testcircular) { /* test circular */
713:     Mat AB;

715:     PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &AB));
716:     PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &X));
717:     PetscCall(CheckLocal(B, X, aB, aX));
718:     if (M == N && N == K) PetscCall(MatMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DETERMINE, &B));
719:     else PetscCall(MatTransposeMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DETERMINE, &B));
720:     PetscCall(CheckLocal(B, X, aB, aX));
721:     PetscCall(MatDestroy(&AB));
722:   }

724:   /* Test by Pierre Jolivet */
725:   {
726:     Mat C, D, D2, AtA;
727:     PetscCall(MatCreateNormal(A, &AtA));
728:     PetscCall(MatDuplicate(X, MAT_DO_NOT_COPY_VALUES, &C));
729:     PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &D));
730:     PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &D2));
731:     PetscCall(MatSetRandom(B, NULL));
732:     PetscCall(MatSetRandom(C, NULL));
733:     PetscCall(MatSetRandom(D, NULL));
734:     PetscCall(MatSetRandom(D2, NULL));
735:     PetscCall(MatProductCreateWithMat(A, B, NULL, C));
736:     PetscCall(MatProductSetType(C, MATPRODUCT_AB));
737:     PetscCall(MatProductSetFromOptions(C));
738:     PetscCall(MatProductSymbolic(C));
739:     PetscCall(MatProductCreateWithMat(A, C, NULL, D));
740:     PetscCall(MatProductSetType(D, MATPRODUCT_AtB));
741:     PetscCall(MatProductSetFromOptions(D));
742:     PetscCall(MatProductSymbolic(D));
743:     PetscCall(MatProductNumeric(C));
744:     PetscCall(MatMatMultEqual(A, B, C, 10, &flg));
745:     if (!flg) {
746:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with Normal (AB != C)\n"));
747:       PetscCall(MatView(A, NULL));
748:       PetscCall(MatView(B, NULL));
749:       PetscCall(MatView(C, NULL));
750:     }
751:     PetscCall(MatProductNumeric(D));
752:     PetscCall(MatMatMultEqual(AtA, B, D, 10, &flg));
753:     if (!flg) {
754:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with Normal (2)\n"));
755:       PetscCall(MatMatMult(AtA, C, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &T));
756:       PetscCall(MatView(D, NULL));
757:       PetscCall(MatView(T, NULL));
758:       PetscCall(MatAXPY(T, -1.0, D, SAME_NONZERO_PATTERN));
759:       PetscCall(MatView(T, NULL));
760:       PetscCall(MatDestroy(&T));
761:     }
762:     PetscCall(MatDestroy(&C));
763:     PetscCall(MatDestroy(&D));
764:     PetscCall(MatDestroy(&D2));
765:     PetscCall(MatDestroy(&AtA));
766:   }

768:   PetscCall(MatDestroy(&X));
769:   PetscCall(MatDestroy(&Bt));
770:   PetscCall(MatDestroy(&B));
771:   PetscCall(MatDestroy(&A));
772:   PetscCall(MatDestroy(&R));
773:   PetscCall(MatDestroy(&PtAP));
774:   PetscCall(MatDestroy(&RARt));
775:   PetscCall(PetscFree(dataX));
776:   PetscCall(PetscFree(dataB));
777:   PetscCall(PetscFree(dataR));
778:   PetscCall(PetscFree(dataBt));
779:   PetscCall(PetscFinalize());
780:   return 0;
781: }

783: /*TEST

785:   test:
786:     output_file: output/empty.out
787:     suffix: 1
788:     args: -local {{0 1}} -testshellops

790:   test:
791:     output_file: output/empty.out
792:     requires: cuda
793:     suffix: 1_cuda
794:     args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{seqaijcusparse seqaij}} -testshellops {{0 1}}

796:   # RARt is skipped for aijkokkos because the R*A products with a dense R it needs are not defined
797:   # -xgpu and -bgpu are not used with aijkokkos because the checks in MatMatMultEqual() and friends cannot
798:   # copy between the VECCUDA of a MATDENSECUDA operand and the VECKOKKOS of an aijkokkos matrix
799:   test:
800:     output_file: output/empty.out
801:     requires: kokkos_kernels
802:     suffix: 1_kokkos
803:     args: -local {{0 1}} -A_mat_type aijkokkos -testrart 0 -A_form_explicit_transpose {{0 1}}

805:   test:
806:     output_file: output/empty.out
807:     nsize: 2
808:     suffix: 1_par
809:     args: -local {{0 1}} -testmatmatt 0

811:   test:
812:     output_file: output/empty.out
813:     nsize: 2
814:     requires: kokkos_kernels
815:     suffix: 1_par_kokkos
816:     args: -local {{0 1}} -testmatmatt 0 -A_mat_type aijkokkos -testrart 0

818:   # MatCreateMAIJ() converts its result to MATAIJKOKKOS or MATAIJCUSPARSE, so the MAIJ matrix
819:   # cached by MatTransposeMatMultSymbolic_MPIAIJ_MPIDense() holds a copy of A that goes stale
820:   # when A's values change before a MAT_REUSE_MATRIX product
821:   test:
822:     output_file: output/empty.out
823:     nsize: 2
824:     requires: kokkos_kernels
825:     suffix: 1_par_kokkos_reuse
826:     args: -testmatmatt 0 -A_mat_type aijkokkos -testrart 0 -testreusemodified
827:     TODO: MATPRODUCT_AtB with MATMPIAIJKOKKOS computes with stale values of A on MAT_REUSE_MATRIX

829:   test:
830:     output_file: output/empty.out
831:     nsize: 2
832:     requires: cuda
833:     suffix: 1_par_cuda_reuse
834:     args: -testnest 0 -testmatmatt 0 -matproduct_batch_size 3 -A_mat_type mpiaijcusparse -testreusemodified
835:     TODO: MATPRODUCT_AtB with MATMPIAIJCUSPARSE computes with stale values of A on MAT_REUSE_MATRIX

837:   test:
838:     output_file: output/empty.out
839:     requires: cuda
840:     nsize: 2
841:     suffix: 1_par_cuda
842:     args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{mpiaijcusparse mpiaij}} -testnest 0 -testmatmatt 0 -matproduct_batch_size 3

844:   test:
845:     output_file: output/empty.out
846:     suffix: 2
847:     nsize: 1
848:     args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}}

850:   testset:
851:     requires: cuda
852:     output_file: output/empty.out
853:     nsize: 1
854:     args: -M 7 -N 9 -K 2 -local {{0 1}} -testnest 0 -A_mat_type {{seqdensecuda seqdense}} -xgpu {{0 1}} -bgpu {{0 1}}
855:     test:
856:       requires: !complex
857:       suffix: 2_cuda_real
858:     test:
859:       # complex+single gives a little bigger error in the MatDenseGetColumnVec test
860:       requires: complex !single
861:       suffix: 2_cuda_complex

863:   test:
864:     output_file: output/empty.out
865:     suffix: 2_par
866:     nsize: 2
867:     args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} -testcircular -testmatmatt 0

869:   test:
870:     requires: cuda
871:     output_file: output/empty.out
872:     suffix: 2_par_cuda
873:     nsize: 2
874:     args: -M 11 -N 9 -K 1 -local {{0 1}} -testcircular 0 -A_mat_type mpiaijcusparse -xgpu -bgpu -testnest 0 -testmatmatt 0

876:   test:
877:     output_file: output/empty.out
878:     suffix: 3
879:     nsize: {{1 3}}
880:     args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type sbaij -symm -testproj 0 -testmatmatt 0

882:   test:
883:     output_file: output/empty.out
884:     suffix: 4
885:     nsize: 1
886:     args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular

888:   test:
889:     output_file: output/empty.out
890:     suffix: 5
891:     nsize: {{2 4}}
892:     args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular -testmatmatt 0

894:   test:
895:     output_file: output/empty.out
896:     suffix: 6
897:     nsize: 1
898:     args: -M {{1 3}} -N {{2 5}} -K {{1 2}} -local {{0 1}} -testcircular

900:   test:
901:     output_file: output/empty.out
902:     suffix: 7
903:     nsize: 1
904:     args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type dense -testnest -testcircular

906: TEST*/