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, (char *)0, 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(PetscMalloc((nzp + 1) * (sizeof(PetscInt) + sizeof(PetscScalar)), &idxn));
 83:   a = (PetscScalar *)(idxn + nzp);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

425:   /* Destroy objects */
426:   PetscCall(VecDestroy(&v1));
427:   PetscCall(VecDestroy(&v2));
428:   PetscCall(PetscRandomDestroy(&rdm));
429:   PetscCall(PetscFree(idxn));

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

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

439: /*TEST

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

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

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

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

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

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

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

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

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

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

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

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

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

527: TEST*/