Actual source code: ex94.c


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

 10: #include <petscmat.h>

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

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

 44:   PetscFunctionBeginUser;
 45:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 46:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 47:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

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

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

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

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

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

 81:   PetscCall(MatGetSize(B, &M, &N));
 82:   nzp = PetscMax((PetscInt)(0.1 * M), 5);
 83:   PetscCall(PetscMalloc((nzp + 1) * (sizeof(PetscInt) + sizeof(PetscScalar)), &idxn));
 84:   a = (PetscScalar *)(idxn + nzp);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

440: /*TEST

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

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

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

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

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

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

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

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

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

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

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

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

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

528: TEST*/