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*/