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