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 */
8
9 #include <petscmat.h>
10
11 /*
12 B = A - B
13 norm = norm(B)
14 */
MatNormDifference(Mat A,Mat B,PetscReal * norm)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 }
22
main(int argc,char ** args)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;
42
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));
47
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));
51
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 }
64
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));
71
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));
77
78 PetscCall(MatGetType(B, &mattype));
79
80 PetscCall(MatGetSize(B, &M, &N));
81 nzp = PetscMax((PetscInt)(0.1 * M), 5);
82 PetscCall(PetscMalloc2(nzp + 1, &idxn, nzp + 1, &a));
83
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));
90
91 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
92 PetscCall(PetscRandomSetFromOptions(rdm));
93 PetscCall(PetscOptionsGetReal(NULL, NULL, "-fill", &fill, NULL));
94
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 */
103
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));
110
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 }
117
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));
126
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));
131
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));
136
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));
143
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));
148
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));
153
154 PetscCall(MatDestroy(&BTT));
155 PetscCall(MatDestroy(&BT));
156 PetscCall(MatDestroy(&ATT));
157 PetscCall(MatDestroy(&AT));
158
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));
162
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));
173
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) */
179
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));
203
204 /* Create R = P^T */
205 PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &R));
206
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 }
213
214 /* C = P^T*B */
215 PetscCall(MatTransposeMatMult(P, B, MAT_INITIAL_MATRIX, fill, &C));
216 PetscCall(MatGetInfo(C, MAT_GLOBAL_SUM, &info));
217
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 }
229
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));
235
236 /* Test MatDuplicate() of C=P^T*B */
237 PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &C1));
238 PetscCall(MatDestroy(&C1));
239 PetscCall(MatDestroy(&C));
240
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));
247
248 /* Test MAT_REUSE_MATRIX - reuse symbolic C */
249 PetscCall(MatMatTransposeMult(B, R, MAT_REUSE_MATRIX, fill, &C));
250
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 }
261
262 /* 3) Test MatPtAP() */
263 /*-------------------*/
264 if (Test_MatPtAP) {
265 PetscInt PN;
266 Mat Cdup;
267
268 PetscCall(MatDuplicate(A_save, MAT_COPY_VALUES, &A));
269 PetscCall(MatGetSize(A, &M, &N));
270 PetscCall(MatGetLocalSize(A, &m, &n));
271
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));
290
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));
295
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 }
303
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 }
313
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));
320
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));
333
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));
339
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));
343
344 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nCdup:\n"));
345 PetscCall(MatView(Cdup, PETSC_VIEWER_STDOUT_WORLD));
346 }
347 PetscCall(MatDestroy(&Cdup));
348
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));
355
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));
360
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));
367
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 }
376
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 }
387
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));
393
394 PetscCall(VecCreate(PETSC_COMM_WORLD, &v3));
395 PetscCall(VecSetSizes(v3, n, PETSC_DECIDE));
396 PetscCall(VecSetFromOptions(v3));
397 PetscCall(VecDuplicate(v3, &v4));
398
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 */
404
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));
410
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);
415
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 }
423
424 /* Destroy objects */
425 PetscCall(VecDestroy(&v1));
426 PetscCall(VecDestroy(&v2));
427 PetscCall(PetscRandomDestroy(&rdm));
428 PetscCall(PetscFree2(idxn, a));
429
430 PetscCall(MatDestroy(&A_save));
431 PetscCall(MatDestroy(&B));
432
433 PetscPreLoadEnd();
434 PetscCall(PetscFinalize());
435 return 0;
436 }
437
438 /*TEST
439
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
445 output_file: output/empty.out
446
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/empty.out
453
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/empty.out
460
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/empty.out
467
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/empty.out
474
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/empty.out
480
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/empty.out
486
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/empty.out
493
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/empty.out
500
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/empty.out
507
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/empty.out
513
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/empty.out
519
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
526
527 TEST*/
528