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