Actual source code: ex94.c

  1: static char help[] = "Tests sequential and parallel MatMatMult() and MatPtAP(), MatTransposeMatMult(), sequential MatMatTransposeMult(), MatRARt()\n\
  2: Input arguments are:\n\
  3:   -f0 <input_file> -f1 <input_file> -f2 <input_file> -f3 <input_file> : file to load\n\n";
  4: /* Example of usage:
  5:    ./ex94 -f0 <A_binary> -f1 <B_binary> -matmatmult_mat_view ascii::ascii_info -matmatmulttr_mat_view
  6:    mpiexec -n 3 ./ex94 -f0 medium -f1 medium -f2 arco1 -f3 arco1 -matmatmult_mat_view
  7: */

  9: #include <petscmat.h>

 11: /*
 12:      B = A - B
 13:      norm = norm(B)
 14: */
 15: PetscErrorCode MatNormDifference(Mat A, Mat B, PetscReal *norm)
 16: {
 17:   PetscFunctionBegin;
 18:   PetscCall(MatAXPY(B, -1.0, A, DIFFERENT_NONZERO_PATTERN));
 19:   PetscCall(MatNorm(B, NORM_FROBENIUS, norm));
 20:   PetscFunctionReturn(PETSC_SUCCESS);
 21: }

 23: int main(int argc, char **args)
 24: {
 25:   Mat          A, A_save, B, AT, ATT, BT, BTT, P, R, C, C1;
 26:   Vec          x, v1, v2, v3, v4;
 27:   PetscViewer  viewer;
 28:   PetscMPIInt  size, rank;
 29:   PetscInt     i, m, n, j, *idxn, M, N, nzp, rstart, rend;
 30:   PetscReal    norm, norm_abs, norm_tmp, fill = 4.0;
 31:   PetscRandom  rdm;
 32:   char         file[4][128];
 33:   PetscBool    flg, preload = PETSC_TRUE;
 34:   PetscScalar *a, rval, alpha, none = -1.0;
 35:   PetscBool    Test_MatMatMult = PETSC_TRUE, Test_MatMatTr = PETSC_TRUE, Test_MatPtAP = PETSC_TRUE, Test_MatRARt = PETSC_TRUE, Test_MatMatMatMult = PETSC_TRUE;
 36:   PetscBool    Test_MatAXPY = PETSC_FALSE, view = PETSC_FALSE;
 37:   PetscInt     pm, pn, pM, pN;
 38:   MatInfo      info;
 39:   PetscBool    seqaij;
 40:   MatType      mattype;
 41:   Mat          Cdensetest, Pdense, Cdense, Adense;

 43:   PetscFunctionBeginUser;
 44:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 45:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 46:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

 48:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-fill", &fill, NULL));
 49:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-matops_view", &view, NULL));
 50:   if (view) PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO));

 52:   /*  Load the matrices A_save and B */
 53:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f0", file[0], sizeof(file[0]), &flg));
 54:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate a file name for small matrix A with the -f0 option.");
 55:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f1", file[1], sizeof(file[1]), &flg));
 56:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate a file name for small matrix B with the -f1 option.");
 57:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f2", file[2], sizeof(file[2]), &flg));
 58:   if (!flg) {
 59:     preload = PETSC_FALSE;
 60:   } else {
 61:     PetscCall(PetscOptionsGetString(NULL, NULL, "-f3", file[3], sizeof(file[3]), &flg));
 62:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate a file name for test matrix B with the -f3 option.");
 63:   }

 65:   PetscPreLoadBegin(preload, "Load system");
 66:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[2 * PetscPreLoadIt], FILE_MODE_READ, &viewer));
 67:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A_save));
 68:   PetscCall(MatSetFromOptions(A_save));
 69:   PetscCall(MatLoad(A_save, viewer));
 70:   PetscCall(PetscViewerDestroy(&viewer));

 72:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[2 * PetscPreLoadIt + 1], FILE_MODE_READ, &viewer));
 73:   PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
 74:   PetscCall(MatSetFromOptions(B));
 75:   PetscCall(MatLoad(B, viewer));
 76:   PetscCall(PetscViewerDestroy(&viewer));

 78:   PetscCall(MatGetType(B, &mattype));

 80:   PetscCall(MatGetSize(B, &M, &N));
 81:   nzp = PetscMax((PetscInt)(0.1 * M), 5);
 82:   PetscCall(PetscMalloc2(nzp + 1, &idxn, nzp + 1, &a));

 84:   /* Create vectors v1 and v2 that are compatible with A_save */
 85:   PetscCall(VecCreate(PETSC_COMM_WORLD, &v1));
 86:   PetscCall(MatGetLocalSize(A_save, &m, NULL));
 87:   PetscCall(VecSetSizes(v1, m, PETSC_DECIDE));
 88:   PetscCall(VecSetFromOptions(v1));
 89:   PetscCall(VecDuplicate(v1, &v2));

 91:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
 92:   PetscCall(PetscRandomSetFromOptions(rdm));
 93:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-fill", &fill, NULL));

 95:   /* Test MatAXPY()    */
 96:   /*-------------------*/
 97:   PetscCall(PetscOptionsHasName(NULL, NULL, "-test_MatAXPY", &Test_MatAXPY));
 98:   if (Test_MatAXPY) {
 99:     Mat Btmp;
100:     PetscCall(MatDuplicate(A_save, MAT_COPY_VALUES, &A));
101:     PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &Btmp));
102:     PetscCall(MatAXPY(A, -1.0, B, DIFFERENT_NONZERO_PATTERN)); /* A = -B + A_save */

104:     PetscCall(MatScale(A, -1.0));                                 /* A = -A = B - A_save */
105:     PetscCall(MatAXPY(Btmp, -1.0, A, DIFFERENT_NONZERO_PATTERN)); /* Btmp = -A + B = A_save */
106:     PetscCall(MatMultEqual(A_save, Btmp, 10, &flg));
107:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatAXPY() is incorrect");
108:     PetscCall(MatDestroy(&A));
109:     PetscCall(MatDestroy(&Btmp));

111:     Test_MatMatMult    = PETSC_FALSE;
112:     Test_MatMatTr      = PETSC_FALSE;
113:     Test_MatPtAP       = PETSC_FALSE;
114:     Test_MatRARt       = PETSC_FALSE;
115:     Test_MatMatMatMult = PETSC_FALSE;
116:   }

118:   /* 1) Test MatMatMult() */
119:   /* ---------------------*/
120:   if (Test_MatMatMult) {
121:     PetscCall(MatDuplicate(A_save, MAT_COPY_VALUES, &A));
122:     PetscCall(MatCreateTranspose(A, &AT));
123:     PetscCall(MatCreateTranspose(AT, &ATT));
124:     PetscCall(MatCreateTranspose(B, &BT));
125:     PetscCall(MatCreateTranspose(BT, &BTT));

127:     PetscCall(MatMatMult(AT, B, MAT_INITIAL_MATRIX, fill, &C));
128:     PetscCall(MatMatMultEqual(AT, B, C, 10, &flg));
129:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult() for C=AT*B");
130:     PetscCall(MatDestroy(&C));

132:     PetscCall(MatMatMult(ATT, B, MAT_INITIAL_MATRIX, fill, &C));
133:     PetscCall(MatMatMultEqual(ATT, B, C, 10, &flg));
134:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult() for C=ATT*B");
135:     PetscCall(MatDestroy(&C));

137:     PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, fill, &C));
138:     PetscCall(MatMatMultEqual(A, B, C, 10, &flg));
139:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult() for reuse C=A*B");
140:     /* ATT has different matrix type as A (although they have same internal data structure),
141:        we cannot call MatProductReplaceMats(ATT,NULL,NULL,C) and MatMatMult(ATT,B,MAT_REUSE_MATRIX,fill,&C) */
142:     PetscCall(MatDestroy(&C));

144:     PetscCall(MatMatMult(A, BTT, MAT_INITIAL_MATRIX, fill, &C));
145:     PetscCall(MatMatMultEqual(A, BTT, C, 10, &flg));
146:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult() for C=A*BTT");
147:     PetscCall(MatDestroy(&C));

149:     PetscCall(MatMatMult(ATT, BTT, MAT_INITIAL_MATRIX, fill, &C));
150:     PetscCall(MatMatMultEqual(A, B, C, 10, &flg));
151:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult()");
152:     PetscCall(MatDestroy(&C));

154:     PetscCall(MatDestroy(&BTT));
155:     PetscCall(MatDestroy(&BT));
156:     PetscCall(MatDestroy(&ATT));
157:     PetscCall(MatDestroy(&AT));

159:     PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, fill, &C));
160:     PetscCall(MatSetOptionsPrefix(C, "matmatmult_")); /* enable option '-matmatmult_' for matrix C */
161:     PetscCall(MatGetInfo(C, MAT_GLOBAL_SUM, &info));

163:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
164:     alpha = 1.0;
165:     for (i = 0; i < 2; i++) {
166:       alpha -= 0.1;
167:       PetscCall(MatScale(A, alpha));
168:       PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, fill, &C));
169:     }
170:     PetscCall(MatMatMultEqual(A, B, C, 10, &flg));
171:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult()");
172:     PetscCall(MatDestroy(&A));

174:     /* Test MatDuplicate() of C=A*B */
175:     PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &C1));
176:     PetscCall(MatDestroy(&C1));
177:     PetscCall(MatDestroy(&C));
178:   } /* if (Test_MatMatMult) */

180:   /* 2) Test MatTransposeMatMult() and MatMatTransposeMult() */
181:   /* ------------------------------------------------------- */
182:   if (Test_MatMatTr) {
183:     /* Create P */
184:     PetscInt PN, rstart, rend;
185:     PN  = M / 2;
186:     nzp = 5; /* num of nonzeros in each row of P */
187:     PetscCall(MatCreate(PETSC_COMM_WORLD, &P));
188:     PetscCall(MatSetSizes(P, PETSC_DECIDE, PETSC_DECIDE, M, PN));
189:     PetscCall(MatSetType(P, mattype));
190:     PetscCall(MatSeqAIJSetPreallocation(P, nzp, NULL));
191:     PetscCall(MatMPIAIJSetPreallocation(P, nzp, NULL, nzp, NULL));
192:     PetscCall(MatGetOwnershipRange(P, &rstart, &rend));
193:     for (i = 0; i < nzp; i++) PetscCall(PetscRandomGetValue(rdm, &a[i]));
194:     for (i = rstart; i < rend; i++) {
195:       for (j = 0; j < nzp; j++) {
196:         PetscCall(PetscRandomGetValue(rdm, &rval));
197:         idxn[j] = (PetscInt)(PetscRealPart(rval) * PN);
198:       }
199:       PetscCall(MatSetValues(P, 1, &i, nzp, idxn, a, ADD_VALUES));
200:     }
201:     PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
202:     PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));

204:     /* Create R = P^T */
205:     PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &R));

207:     { /* Test R = P^T, C1 = R*B */
208:       PetscCall(MatMatMult(R, B, MAT_INITIAL_MATRIX, fill, &C1));
209:       PetscCall(MatTranspose(P, MAT_REUSE_MATRIX, &R));
210:       PetscCall(MatMatMult(R, B, MAT_REUSE_MATRIX, fill, &C1));
211:       PetscCall(MatDestroy(&C1));
212:     }

214:     /* C = P^T*B */
215:     PetscCall(MatTransposeMatMult(P, B, MAT_INITIAL_MATRIX, fill, &C));
216:     PetscCall(MatGetInfo(C, MAT_GLOBAL_SUM, &info));

218:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
219:     PetscCall(MatTransposeMatMult(P, B, MAT_REUSE_MATRIX, fill, &C));
220:     if (view) {
221:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C = P^T * B:\n"));
222:       PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
223:     }
224:     PetscCall(MatProductClear(C));
225:     if (view) {
226:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nC = P^T * B after MatProductClear():\n"));
227:       PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
228:     }

230:     /* Compare P^T*B and R*B */
231:     PetscCall(MatMatMult(R, B, MAT_INITIAL_MATRIX, fill, &C1));
232:     PetscCall(MatNormDifference(C, C1, &norm));
233:     PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatTransposeMatMult(): %g", (double)norm);
234:     PetscCall(MatDestroy(&C1));

236:     /* Test MatDuplicate() of C=P^T*B */
237:     PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &C1));
238:     PetscCall(MatDestroy(&C1));
239:     PetscCall(MatDestroy(&C));

241:     /* C = B*R^T */
242:     PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &seqaij));
243:     if (size == 1 && seqaij) {
244:       PetscCall(MatMatTransposeMult(B, R, MAT_INITIAL_MATRIX, fill, &C));
245:       PetscCall(MatSetOptionsPrefix(C, "matmatmulttr_")); /* enable '-matmatmulttr_' for matrix C */
246:       PetscCall(MatGetInfo(C, MAT_GLOBAL_SUM, &info));

248:       /* Test MAT_REUSE_MATRIX - reuse symbolic C */
249:       PetscCall(MatMatTransposeMult(B, R, MAT_REUSE_MATRIX, fill, &C));

251:       /* Check */
252:       PetscCall(MatMatMult(B, P, MAT_INITIAL_MATRIX, fill, &C1));
253:       PetscCall(MatNormDifference(C, C1, &norm));
254:       PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatMatTransposeMult() %g", (double)norm);
255:       PetscCall(MatDestroy(&C1));
256:       PetscCall(MatDestroy(&C));
257:     }
258:     PetscCall(MatDestroy(&P));
259:     PetscCall(MatDestroy(&R));
260:   }

262:   /* 3) Test MatPtAP() */
263:   /*-------------------*/
264:   if (Test_MatPtAP) {
265:     PetscInt PN;
266:     Mat      Cdup;

268:     PetscCall(MatDuplicate(A_save, MAT_COPY_VALUES, &A));
269:     PetscCall(MatGetSize(A, &M, &N));
270:     PetscCall(MatGetLocalSize(A, &m, &n));

272:     PN  = M / 2;
273:     nzp = (PetscInt)(0.1 * PN + 1); /* num of nonzeros in each row of P */
274:     PetscCall(MatCreate(PETSC_COMM_WORLD, &P));
275:     PetscCall(MatSetSizes(P, PETSC_DECIDE, PETSC_DECIDE, N, PN));
276:     PetscCall(MatSetType(P, mattype));
277:     PetscCall(MatSeqAIJSetPreallocation(P, nzp, NULL));
278:     PetscCall(MatMPIAIJSetPreallocation(P, nzp, NULL, nzp, NULL));
279:     for (i = 0; i < nzp; i++) PetscCall(PetscRandomGetValue(rdm, &a[i]));
280:     PetscCall(MatGetOwnershipRange(P, &rstart, &rend));
281:     for (i = rstart; i < rend; i++) {
282:       for (j = 0; j < nzp; j++) {
283:         PetscCall(PetscRandomGetValue(rdm, &rval));
284:         idxn[j] = (PetscInt)(PetscRealPart(rval) * PN);
285:       }
286:       PetscCall(MatSetValues(P, 1, &i, nzp, idxn, a, ADD_VALUES));
287:     }
288:     PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
289:     PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));

291:     /* PetscCall(MatView(P,PETSC_VIEWER_STDOUT_WORLD)); */
292:     PetscCall(MatGetSize(P, &pM, &pN));
293:     PetscCall(MatGetLocalSize(P, &pm, &pn));
294:     PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, fill, &C));

296:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
297:     alpha = 1.0;
298:     for (i = 0; i < 2; i++) {
299:       alpha -= 0.1;
300:       PetscCall(MatScale(A, alpha));
301:       PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, fill, &C));
302:     }

304:     /* Test PtAP ops with P Dense and A either AIJ or SeqDense (it assumes MatPtAP_XAIJ_XAIJ is fine) */
305:     PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &seqaij));
306:     if (seqaij) {
307:       PetscCall(MatConvert(C, MATSEQDENSE, MAT_INITIAL_MATRIX, &Cdensetest));
308:       PetscCall(MatConvert(P, MATSEQDENSE, MAT_INITIAL_MATRIX, &Pdense));
309:     } else {
310:       PetscCall(MatConvert(C, MATMPIDENSE, MAT_INITIAL_MATRIX, &Cdensetest));
311:       PetscCall(MatConvert(P, MATMPIDENSE, MAT_INITIAL_MATRIX, &Pdense));
312:     }

314:     /* test with A(AIJ), Pdense -- call MatPtAP_Basic() when np>1 */
315:     PetscCall(MatPtAP(A, Pdense, MAT_INITIAL_MATRIX, fill, &Cdense));
316:     PetscCall(MatPtAP(A, Pdense, MAT_REUSE_MATRIX, fill, &Cdense));
317:     PetscCall(MatPtAPMultEqual(A, Pdense, Cdense, 10, &flg));
318:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatPtAP with A AIJ and P Dense");
319:     PetscCall(MatDestroy(&Cdense));

321:     /* test with A SeqDense */
322:     if (seqaij) {
323:       PetscCall(MatConvert(A, MATSEQDENSE, MAT_INITIAL_MATRIX, &Adense));
324:       PetscCall(MatPtAP(Adense, Pdense, MAT_INITIAL_MATRIX, fill, &Cdense));
325:       PetscCall(MatPtAP(Adense, Pdense, MAT_REUSE_MATRIX, fill, &Cdense));
326:       PetscCall(MatPtAPMultEqual(Adense, Pdense, Cdense, 10, &flg));
327:       PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in MatPtAP with A SeqDense and P SeqDense");
328:       PetscCall(MatDestroy(&Cdense));
329:       PetscCall(MatDestroy(&Adense));
330:     }
331:     PetscCall(MatDestroy(&Cdensetest));
332:     PetscCall(MatDestroy(&Pdense));

334:     /* Test MatDuplicate() of C=PtAP and MatView(Cdup,...) */
335:     PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &Cdup));
336:     if (view) {
337:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nC = P^T * A * P:\n"));
338:       PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));

340:       PetscCall(MatProductClear(C));
341:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nC = P^T * A * P after MatProductClear():\n"));
342:       PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));

344:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nCdup:\n"));
345:       PetscCall(MatView(Cdup, PETSC_VIEWER_STDOUT_WORLD));
346:     }
347:     PetscCall(MatDestroy(&Cdup));

349:     if (size > 1 || !seqaij) Test_MatRARt = PETSC_FALSE;
350:     /* 4) Test MatRARt() */
351:     /* ----------------- */
352:     if (Test_MatRARt) {
353:       Mat R, RARt, Rdense, RARtdense;
354:       PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &R));

356:       /* Test MatRARt_Basic(), MatMatMatMult_Basic() */
357:       PetscCall(MatConvert(R, MATDENSE, MAT_INITIAL_MATRIX, &Rdense));
358:       PetscCall(MatRARt(A, Rdense, MAT_INITIAL_MATRIX, 2.0, &RARtdense));
359:       PetscCall(MatRARt(A, Rdense, MAT_REUSE_MATRIX, 2.0, &RARtdense));

361:       PetscCall(MatConvert(RARtdense, MATAIJ, MAT_INITIAL_MATRIX, &RARt));
362:       PetscCall(MatNormDifference(C, RARt, &norm));
363:       PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "|PtAP - RARtdense| = %g", (double)norm);
364:       PetscCall(MatDestroy(&Rdense));
365:       PetscCall(MatDestroy(&RARtdense));
366:       PetscCall(MatDestroy(&RARt));

368:       /* Test MatRARt() for aij matrices */
369:       PetscCall(MatRARt(A, R, MAT_INITIAL_MATRIX, 2.0, &RARt));
370:       PetscCall(MatRARt(A, R, MAT_REUSE_MATRIX, 2.0, &RARt));
371:       PetscCall(MatNormDifference(C, RARt, &norm));
372:       PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "|PtAP - RARt| = %g", (double)norm);
373:       PetscCall(MatDestroy(&R));
374:       PetscCall(MatDestroy(&RARt));
375:     }

377:     if (Test_MatMatMatMult && size == 1) {
378:       Mat R, RAP;
379:       PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &R));
380:       PetscCall(MatMatMatMult(R, A, P, MAT_INITIAL_MATRIX, 2.0, &RAP));
381:       PetscCall(MatMatMatMult(R, A, P, MAT_REUSE_MATRIX, 2.0, &RAP));
382:       PetscCall(MatNormDifference(C, RAP, &norm));
383:       PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PtAP != RAP %g", (double)norm);
384:       PetscCall(MatDestroy(&R));
385:       PetscCall(MatDestroy(&RAP));
386:     }

388:     /* Create vector x that is compatible with P */
389:     PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
390:     PetscCall(MatGetLocalSize(P, &m, &n));
391:     PetscCall(VecSetSizes(x, n, PETSC_DECIDE));
392:     PetscCall(VecSetFromOptions(x));

394:     PetscCall(VecCreate(PETSC_COMM_WORLD, &v3));
395:     PetscCall(VecSetSizes(v3, n, PETSC_DECIDE));
396:     PetscCall(VecSetFromOptions(v3));
397:     PetscCall(VecDuplicate(v3, &v4));

399:     norm = 0.0;
400:     for (i = 0; i < 10; i++) {
401:       PetscCall(VecSetRandom(x, rdm));
402:       PetscCall(MatMult(P, x, v1));
403:       PetscCall(MatMult(A, v1, v2)); /* v2 = A*P*x */

405:       PetscCall(MatMultTranspose(P, v2, v3)); /* v3 = Pt*A*P*x */
406:       PetscCall(MatMult(C, x, v4));           /* v3 = C*x   */
407:       PetscCall(VecNorm(v4, NORM_2, &norm_abs));
408:       PetscCall(VecAXPY(v4, none, v3));
409:       PetscCall(VecNorm(v4, NORM_2, &norm_tmp));

411:       norm_tmp /= norm_abs;
412:       if (norm_tmp > norm) norm = norm_tmp;
413:     }
414:     PetscCheck(norm < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatPtAP(), |v1 - v2|: %g", (double)norm);

416:     PetscCall(MatDestroy(&A));
417:     PetscCall(MatDestroy(&P));
418:     PetscCall(MatDestroy(&C));
419:     PetscCall(VecDestroy(&v3));
420:     PetscCall(VecDestroy(&v4));
421:     PetscCall(VecDestroy(&x));
422:   }

424:   /* Destroy objects */
425:   PetscCall(VecDestroy(&v1));
426:   PetscCall(VecDestroy(&v2));
427:   PetscCall(PetscRandomDestroy(&rdm));
428:   PetscCall(PetscFree2(idxn, a));

430:   PetscCall(MatDestroy(&A_save));
431:   PetscCall(MatDestroy(&B));

433:   PetscPreLoadEnd();
434:   PetscCall(PetscFinalize());
435:   return 0;
436: }

438: /*TEST

440:    test:
441:       suffix: 2_mattransposematmult_matmatmult
442:       nsize: 3
443:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
444:       args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium -mattransposematmult_via at*b> ex94_2.tmp 2>&1

446:    test:
447:       suffix: 2_mattransposematmult_scalable
448:       nsize: 3
449:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
450:       args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium -mattransposematmult_via scalable> ex94_2.tmp 2>&1
451:       output_file: output/ex94_1.out

453:    test:
454:       suffix: axpy_mpiaij
455:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
456:       nsize: 8
457:       args: -f0 ${DATAFILESPATH}/matrices/poisson_2d5p -f1 ${DATAFILESPATH}/matrices/poisson_2d13p -test_MatAXPY
458:       output_file: output/ex94_1.out

460:    test:
461:       suffix: axpy_mpibaij
462:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
463:       nsize: 8
464:       args: -f0 ${DATAFILESPATH}/matrices/poisson_2d5p -f1 ${DATAFILESPATH}/matrices/poisson_2d13p -test_MatAXPY -mat_type baij
465:       output_file: output/ex94_1.out

467:    test:
468:       suffix: axpy_mpisbaij
469:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
470:       nsize: 8
471:       args: -f0 ${DATAFILESPATH}/matrices/poisson_2d5p -f1 ${DATAFILESPATH}/matrices/poisson_2d13p -test_MatAXPY -mat_type sbaij
472:       output_file: output/ex94_1.out

474:    test:
475:       suffix: matmatmult
476:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
477:       args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -viewer_binary_skip_info
478:       output_file: output/ex94_1.out

480:    test:
481:       suffix: matmatmult_2
482:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
483:       args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -mat_type mpiaij -viewer_binary_skip_info
484:       output_file: output/ex94_1.out

486:    test:
487:       suffix: matmatmult_scalable
488:       nsize: 4
489:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
490:       args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -matmatmult_via scalable
491:       output_file: output/ex94_1.out

493:    test:
494:       suffix: ptap
495:       nsize: 3
496:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
497:       args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium -matptap_via scalable
498:       output_file: output/ex94_1.out

500:    test:
501:       suffix: rap
502:       nsize: 3
503:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
504:       args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium
505:       output_file: output/ex94_1.out

507:    test:
508:       suffix: scalable0
509:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
510:       args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -viewer_binary_skip_info
511:       output_file: output/ex94_1.out

513:    test:
514:       suffix: scalable1
515:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
516:       args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -viewer_binary_skip_info -matptap_via scalable
517:       output_file: output/ex94_1.out

519:    test:
520:       suffix: view
521:       nsize: 2
522:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
523:       args: -f0 ${DATAFILESPATH}/matrices/tiny -f1 ${DATAFILESPATH}/matrices/tiny -viewer_binary_skip_info -matops_view
524:       output_file: output/ex94_2.out

526: TEST*/