Actual source code: ex70.c
1: #include <petscmat.h>
3: static char help[] = "Tests MatMat operations with MAT_REUSE_MATRIX and already allocated dense result.\n\n";
5: static PetscScalar MAGIC_NUMBER = 12345;
7: static PetscErrorCode CheckLocal(Mat A, Mat B, PetscScalar *a, PetscScalar *b)
8: {
9: PetscBool wA = PETSC_FALSE, wB = PETSC_FALSE;
10: PetscBool wAv = PETSC_FALSE, wBv = PETSC_FALSE;
11: PetscInt lda, i, j, m, n;
13: if (a) {
14: const PetscScalar *Aa;
15: MatDenseGetArrayRead(A, &Aa);
16: wA = (PetscBool)(a != Aa);
17: MatDenseGetLDA(A, &lda);
18: MatGetLocalSize(A, &m, &n);
19: for (j = 0; j < n; j++) {
20: for (i = m; i < lda; i++) {
21: if (Aa[j * lda + i] != MAGIC_NUMBER) wAv = PETSC_TRUE;
22: }
23: }
24: MatDenseRestoreArrayRead(A, &Aa);
25: }
26: if (b) {
27: const PetscScalar *Bb;
28: MatDenseGetArrayRead(B, &Bb);
29: wB = (PetscBool)(b != Bb);
30: MatDenseGetLDA(B, &lda);
31: MatGetLocalSize(B, &m, &n);
32: for (j = 0; j < n; j++) {
33: for (i = m; i < lda; i++) {
34: if (Bb[j * lda + i] != MAGIC_NUMBER) wBv = PETSC_TRUE;
35: }
36: }
37: MatDenseRestoreArrayRead(B, &Bb);
38: }
41: return 0;
42: }
44: typedef struct {
45: Mat A;
46: Mat P;
47: Mat R;
48: } proj_data;
50: PetscErrorCode proj_destroy(void *ctx)
51: {
52: proj_data *userdata = (proj_data *)ctx;
55: MatDestroy(&userdata->A);
56: MatDestroy(&userdata->P);
57: MatDestroy(&userdata->R);
58: PetscFree(userdata);
59: return 0;
60: }
62: PetscErrorCode proj_mult(Mat S, Vec X, Vec Y)
63: {
64: Mat A, R, P;
65: Vec Ax, Ay;
66: Vec Px, Py;
67: proj_data *userdata;
69: MatShellGetContext(S, &userdata);
71: A = userdata->A;
72: R = userdata->R;
73: P = userdata->P;
77: MatCreateVecs(A, &Ax, &Ay);
78: if (R) {
79: MatCreateVecs(R, &Py, &Px);
80: } else {
81: MatCreateVecs(P, &Px, &Py);
82: }
83: VecCopy(X, Px);
84: if (P) {
85: MatMult(P, Px, Py);
86: } else {
87: MatMultTranspose(R, Px, Py);
88: }
89: VecCopy(Py, Ax);
90: MatMult(A, Ax, Ay);
91: VecCopy(Ay, Py);
92: if (P) {
93: MatMultTranspose(P, Py, Px);
94: } else {
95: MatMult(R, Py, Px);
96: }
97: VecCopy(Px, Y);
98: VecDestroy(&Px);
99: VecDestroy(&Py);
100: VecDestroy(&Ax);
101: VecDestroy(&Ay);
102: return 0;
103: }
105: PetscErrorCode MyPtShellPMultSymbolic(Mat S, Mat P, Mat PtAP, void **ctx)
106: {
107: proj_data *userdata;
109: PetscNew(&userdata);
110: MatShellSetContext(PtAP, userdata);
111: *ctx = (void *)userdata;
112: return 0;
113: }
115: PetscErrorCode MyPtShellPMultNumeric(Mat S, Mat P, Mat PtAP, void *ctx)
116: {
117: Mat A;
118: proj_data *userdata = (proj_data *)ctx;
120: MatShellGetContext(S, &A);
121: PetscObjectReference((PetscObject)A);
122: PetscObjectReference((PetscObject)P);
123: MatDestroy(&userdata->A);
124: MatDestroy(&userdata->P);
125: MatDestroy(&userdata->R);
126: userdata->A = A;
127: userdata->P = P;
128: MatShellSetOperation(PtAP, MATOP_MULT, (void (*)(void))proj_mult);
129: MatSetUp(PtAP);
130: MatAssemblyBegin(PtAP, MAT_FINAL_ASSEMBLY);
131: MatAssemblyEnd(PtAP, MAT_FINAL_ASSEMBLY);
132: return 0;
133: }
135: PetscErrorCode MyRShellRtMultSymbolic(Mat S, Mat R, Mat RARt, void **ctx)
136: {
137: proj_data *userdata;
139: PetscNew(&userdata);
140: MatShellSetContext(RARt, userdata);
141: *ctx = (void *)userdata;
142: return 0;
143: }
145: PetscErrorCode MyRShellRtMultNumeric(Mat S, Mat R, Mat RARt, void *ctx)
146: {
147: Mat A;
148: proj_data *userdata = (proj_data *)ctx;
150: MatShellGetContext(S, &A);
151: PetscObjectReference((PetscObject)A);
152: PetscObjectReference((PetscObject)R);
153: MatDestroy(&userdata->A);
154: MatDestroy(&userdata->P);
155: MatDestroy(&userdata->R);
156: userdata->A = A;
157: userdata->R = R;
158: MatShellSetOperation(RARt, MATOP_MULT, (void (*)(void))proj_mult);
159: MatSetUp(RARt);
160: MatAssemblyBegin(RARt, MAT_FINAL_ASSEMBLY);
161: MatAssemblyEnd(RARt, MAT_FINAL_ASSEMBLY);
162: return 0;
163: }
165: PetscErrorCode MyMatShellMatMultNumeric(Mat S, Mat B, Mat C, void *ctx)
166: {
167: Mat A;
169: MatShellGetContext(S, &A);
170: MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C);
171: return 0;
172: }
174: PetscErrorCode MyMatTransposeShellMatMultNumeric(Mat S, Mat B, Mat C, void *ctx)
175: {
176: Mat A;
178: MatShellGetContext(S, &A);
179: MatTransposeMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C);
180: return 0;
181: }
183: PetscErrorCode MyMatShellMatTransposeMultNumeric(Mat S, Mat B, Mat C, void *ctx)
184: {
185: Mat A;
187: MatShellGetContext(S, &A);
188: MatMatTransposeMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C);
189: return 0;
190: }
192: int main(int argc, char **args)
193: {
194: Mat X, B, A, Bt, T, T2, PtAP = NULL, RARt = NULL, R = NULL;
195: Vec r, l, rs, ls;
196: PetscInt m, n, k, M = 10, N = 10, K = 5, ldx = 3, ldb = 5, ldr = 4;
197: const char *deft = MATAIJ;
198: char mattype[256];
199: PetscBool flg, symm = PETSC_FALSE, testtt = PETSC_TRUE, testnest = PETSC_TRUE, testtranspose = PETSC_TRUE, testcircular = PETSC_FALSE, local = PETSC_TRUE;
200: PetscBool testhtranspose = PETSC_FALSE; /* Hermitian transpose is not handled correctly and generates an error */
201: PetscBool xgpu = PETSC_FALSE, bgpu = PETSC_FALSE, testshellops = PETSC_FALSE, testproj = PETSC_TRUE, testrart = PETSC_TRUE, testmatmatt = PETSC_TRUE, testmattmat = PETSC_TRUE;
202: PetscScalar *dataX = NULL, *dataB = NULL, *dataR = NULL, *dataBt = NULL;
203: PetscScalar *aX, *aB, *aBt;
204: PetscReal err;
207: PetscInitialize(&argc, &args, NULL, help);
208: PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
209: PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
210: PetscOptionsGetInt(NULL, NULL, "-K", &K, NULL);
211: PetscOptionsGetBool(NULL, NULL, "-symm", &symm, NULL);
212: PetscOptionsGetBool(NULL, NULL, "-local", &local, NULL);
213: PetscOptionsGetInt(NULL, NULL, "-ldx", &ldx, NULL);
214: PetscOptionsGetInt(NULL, NULL, "-ldb", &ldb, NULL);
215: PetscOptionsGetInt(NULL, NULL, "-ldr", &ldr, NULL);
216: PetscOptionsGetBool(NULL, NULL, "-testtranspose", &testtranspose, NULL);
217: PetscOptionsGetBool(NULL, NULL, "-testnest", &testnest, NULL);
218: PetscOptionsGetBool(NULL, NULL, "-testtt", &testtt, NULL);
219: PetscOptionsGetBool(NULL, NULL, "-testcircular", &testcircular, NULL);
220: PetscOptionsGetBool(NULL, NULL, "-testshellops", &testshellops, NULL);
221: PetscOptionsGetBool(NULL, NULL, "-testproj", &testproj, NULL);
222: PetscOptionsGetBool(NULL, NULL, "-testrart", &testrart, NULL);
223: PetscOptionsGetBool(NULL, NULL, "-testmatmatt", &testmatmatt, NULL);
224: PetscOptionsGetBool(NULL, NULL, "-testmattmat", &testmattmat, NULL);
225: PetscOptionsGetBool(NULL, NULL, "-xgpu", &xgpu, NULL);
226: PetscOptionsGetBool(NULL, NULL, "-bgpu", &bgpu, NULL);
227: PetscOptionsGetScalar(NULL, NULL, "-magic_number", &MAGIC_NUMBER, NULL);
228: if (M != N) testproj = PETSC_FALSE;
230: MatCreate(PETSC_COMM_WORLD, &A);
231: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N);
232: MatSetType(A, MATAIJ);
233: MatSetUp(A);
234: MatSetRandom(A, NULL);
235: if (M == N && symm) {
236: Mat AT;
238: MatTranspose(A, MAT_INITIAL_MATRIX, &AT);
239: MatAXPY(A, 1.0, AT, DIFFERENT_NONZERO_PATTERN);
240: MatDestroy(&AT);
241: MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
242: }
243: MatViewFromOptions(A, NULL, "-A_init_view");
244: PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "");
245: PetscOptionsFList("-A_mat_type", "Matrix type", "MatSetType", MatList, deft, mattype, 256, &flg);
246: PetscOptionsEnd();
247: if (flg) {
248: Mat A2;
250: MatDuplicate(A, MAT_COPY_VALUES, &A2);
251: MatConvert(A, mattype, MAT_INPLACE_MATRIX, &A);
252: MatMultEqual(A, A2, 10, &flg);
253: if (!flg) {
254: Mat AE, A2E;
256: PetscPrintf(PETSC_COMM_WORLD, "Error with convert\n");
257: MatComputeOperator(A, MATDENSE, &AE);
258: MatComputeOperator(A2, MATDENSE, &A2E);
259: MatView(AE, NULL);
260: MatView(A2E, NULL);
261: MatAXPY(A2E, -1.0, A, SAME_NONZERO_PATTERN);
262: MatView(A2E, NULL);
263: MatDestroy(&A2E);
264: MatDestroy(&AE);
265: }
266: MatDestroy(&A2);
267: }
268: MatViewFromOptions(A, NULL, "-A_view");
270: MatGetLocalSize(A, &m, &n);
271: if (local) {
272: PetscInt i;
274: PetscMalloc1((m + ldx) * K, &dataX);
275: PetscMalloc1((n + ldb) * K, &dataB);
276: for (i = 0; i < (m + ldx) * K; i++) dataX[i] = MAGIC_NUMBER;
277: for (i = 0; i < (n + ldb) * K; i++) dataB[i] = MAGIC_NUMBER;
278: }
279: MatCreateDense(PETSC_COMM_WORLD, n, PETSC_DECIDE, N, K, dataB, &B);
280: MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, dataX, &X);
281: if (local) {
282: MatDenseSetLDA(X, m + ldx);
283: MatDenseSetLDA(B, n + ldb);
284: }
285: MatGetLocalSize(B, NULL, &k);
286: if (local) {
287: PetscInt i;
289: PetscMalloc1((k + ldr) * N, &dataBt);
290: for (i = 0; i < (k + ldr) * N; i++) dataBt[i] = MAGIC_NUMBER;
291: }
292: MatCreateDense(PETSC_COMM_WORLD, k, n, K, N, dataBt, &Bt);
293: if (local) MatDenseSetLDA(Bt, k + ldr);
295: /* store pointer to dense data for testing */
296: MatDenseGetArrayRead(B, (const PetscScalar **)&dataB);
297: MatDenseGetArrayRead(X, (const PetscScalar **)&dataX);
298: MatDenseGetArrayRead(Bt, (const PetscScalar **)&dataBt);
299: aX = dataX;
300: aB = dataB;
301: aBt = dataBt;
302: MatDenseRestoreArrayRead(Bt, (const PetscScalar **)&dataBt);
303: MatDenseRestoreArrayRead(B, (const PetscScalar **)&dataB);
304: MatDenseRestoreArrayRead(X, (const PetscScalar **)&dataX);
305: if (local) {
306: dataX = aX;
307: dataB = aB;
308: dataBt = aBt;
309: }
311: MatSetRandom(X, NULL);
312: MatSetRandom(B, NULL);
313: MatSetRandom(Bt, NULL);
314: CheckLocal(X, NULL, aX, NULL);
315: CheckLocal(Bt, B, aBt, aB);
317: /* convert to CUDA if needed */
318: if (bgpu) {
319: MatConvert(B, MATDENSECUDA, MAT_INPLACE_MATRIX, &B);
320: MatConvert(Bt, MATDENSECUDA, MAT_INPLACE_MATRIX, &Bt);
321: }
322: if (xgpu) MatConvert(X, MATDENSECUDA, MAT_INPLACE_MATRIX, &X);
323: CheckLocal(B, X, aB, aX);
325: /* Test MatDenseGetSubMatrix */
326: {
327: Mat B2, T3, T4;
329: MatDuplicate(B, MAT_COPY_VALUES, &B2);
330: MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &T4);
331: MatSetRandom(T4, NULL);
332: MatAXPY(B2, 1.0, T4, SAME_NONZERO_PATTERN);
333: MatDenseGetSubMatrix(B, PETSC_DECIDE, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T);
334: MatDenseGetSubMatrix(T4, PETSC_DECIDE, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T2);
335: MatDenseGetSubMatrix(B2, PETSC_DECIDE, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T3);
336: MatAXPY(T, 1.0, T2, SAME_NONZERO_PATTERN);
337: MatAXPY(T3, -1.0, T, SAME_NONZERO_PATTERN);
338: MatNorm(T3, NORM_FROBENIUS, &err);
339: if (err > PETSC_SMALL) {
340: PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetSubMatrix\n");
341: MatView(T3, NULL);
342: }
343: MatDenseRestoreSubMatrix(B, &T);
344: MatDenseRestoreSubMatrix(T4, &T2);
345: MatDenseRestoreSubMatrix(B2, &T3);
346: CheckLocal(B, NULL, aB, NULL);
347: MatDestroy(&B2);
348: MatDestroy(&T4);
349: if (N >= 2) {
350: MatDuplicate(B, MAT_COPY_VALUES, &B2);
351: MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &T4);
352: MatSetRandom(T4, NULL);
353: MatAXPY(B2, 1.0, T4, SAME_NONZERO_PATTERN);
354: MatDenseGetSubMatrix(B, N - 2, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T);
355: MatDenseGetSubMatrix(T4, N - 2, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T2);
356: MatDenseGetSubMatrix(B2, N - 2, PETSC_DECIDE, PetscMin(1, K - 1), PetscMin(2, K), &T3);
357: MatAXPY(T, 1.0, T2, SAME_NONZERO_PATTERN);
358: MatAXPY(T3, -1.0, T, SAME_NONZERO_PATTERN);
359: MatNorm(T3, NORM_FROBENIUS, &err);
360: if (err > PETSC_SMALL) {
361: PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetSubMatrix\n");
362: MatView(T3, NULL);
363: }
364: MatDenseRestoreSubMatrix(B, &T);
365: MatDenseRestoreSubMatrix(T4, &T2);
366: MatDenseRestoreSubMatrix(B2, &T3);
367: CheckLocal(B, NULL, aB, NULL);
368: MatDestroy(&B2);
369: MatDestroy(&T4);
370: MatDuplicate(B, MAT_COPY_VALUES, &B2);
371: MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &T4);
372: MatSetRandom(T4, NULL);
373: MatAXPY(B2, 1.0, T4, SAME_NONZERO_PATTERN);
374: MatDenseGetSubMatrix(B, PETSC_DECIDE, 2, PetscMin(1, K - 1), PetscMin(2, K), &T);
375: MatDenseGetSubMatrix(T4, PETSC_DECIDE, 2, PetscMin(1, K - 1), PetscMin(2, K), &T2);
376: MatDenseGetSubMatrix(B2, PETSC_DECIDE, 2, PetscMin(1, K - 1), PetscMin(2, K), &T3);
377: MatAXPY(T, 1.0, T2, SAME_NONZERO_PATTERN);
378: MatAXPY(T3, -1.0, T, SAME_NONZERO_PATTERN);
379: MatNorm(T3, NORM_FROBENIUS, &err);
380: if (err > PETSC_SMALL) {
381: PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetSubMatrix\n");
382: MatView(T3, NULL);
383: }
384: MatDenseRestoreSubMatrix(B, &T);
385: MatDenseRestoreSubMatrix(T4, &T2);
386: MatDenseRestoreSubMatrix(B2, &T3);
387: CheckLocal(B, NULL, aB, NULL);
388: MatDestroy(&B2);
389: MatDestroy(&T4);
390: }
391: }
393: /* Test reusing a previously allocated dense buffer */
394: MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
395: CheckLocal(B, X, aB, aX);
396: MatMatMultEqual(A, B, X, 10, &flg);
397: if (!flg) {
398: PetscPrintf(PETSC_COMM_WORLD, "Error with reusage\n");
399: MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
400: MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN);
401: MatView(T, NULL);
402: MatDestroy(&T);
403: }
405: /* Test MatTransposeMat and MatMatTranspose */
406: if (testmattmat) {
407: MatTransposeMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
408: CheckLocal(B, X, aB, aX);
409: MatTransposeMatMultEqual(A, X, B, 10, &flg);
410: if (!flg) {
411: PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatTransposeMat)\n");
412: MatTransposeMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &B);
413: MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN);
414: MatView(T, NULL);
415: MatDestroy(&T);
416: }
417: }
418: if (testmatmatt) {
419: MatMatTransposeMult(A, Bt, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
420: CheckLocal(Bt, X, aBt, aX);
421: MatMatTransposeMultEqual(A, Bt, X, 10, &flg);
422: if (!flg) {
423: PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatMatTranspose)\n");
424: MatMatTransposeMult(A, Bt, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
425: MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN);
426: MatView(T, NULL);
427: MatDestroy(&T);
428: }
429: }
431: /* Test projection operations (PtAP and RARt) */
432: if (testproj) {
433: MatPtAP(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &PtAP);
434: MatPtAPMultEqual(A, B, PtAP, 10, &flg);
435: if (!flg) {
436: PetscPrintf(PETSC_COMM_WORLD, "Error with PtAP\n");
437: MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
438: MatTransposeMatMult(B, T, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T2);
439: MatAXPY(T2, -1.0, PtAP, SAME_NONZERO_PATTERN);
440: MatView(T2, NULL);
441: MatDestroy(&T2);
442: MatDestroy(&T);
443: }
444: PetscMalloc1((k + ldr) * M, &dataR);
445: MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, m, K, M, dataR, &R);
446: MatDenseSetLDA(R, k + ldr);
447: MatSetRandom(R, NULL);
448: if (testrart) { /* fails for AIJCUSPARSE because RA operation is not defined */
449: MatRARt(A, R, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &RARt);
450: MatRARtMultEqual(A, R, RARt, 10, &flg);
451: if (!flg) {
452: PetscPrintf(PETSC_COMM_WORLD, "Error with RARt\n");
453: MatMatTransposeMult(A, R, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
454: MatMatMult(R, T, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T2);
455: MatAXPY(T2, -1.0, RARt, SAME_NONZERO_PATTERN);
456: MatView(T2, NULL);
457: MatDestroy(&T2);
458: MatDestroy(&T);
459: }
460: }
461: }
463: /* Test MatDenseGetColumnVec and friends */
464: MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
465: MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
466: MatDuplicate(T, MAT_DO_NOT_COPY_VALUES, &T2);
467: for (k = 0; k < K; k++) {
468: Vec Xv, Tv, T2v;
470: MatDenseGetColumnVecRead(X, k, &Xv);
471: MatDenseGetColumnVec(T, k, &Tv);
472: MatDenseGetColumnVecWrite(T2, k, &T2v);
473: VecCopy(Xv, T2v);
474: VecAXPY(Tv, -1., Xv);
475: MatDenseRestoreColumnVecRead(X, k, &Xv);
476: MatDenseRestoreColumnVec(T, k, &Tv);
477: MatDenseRestoreColumnVecWrite(T2, k, &T2v);
478: }
479: MatNorm(T, NORM_FROBENIUS, &err);
480: if (err > PETSC_SMALL) {
481: PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetColumnVec\n");
482: MatView(T, NULL);
483: }
484: MatAXPY(T2, -1., X, SAME_NONZERO_PATTERN);
485: MatNorm(T2, NORM_FROBENIUS, &err);
486: if (err > PETSC_SMALL) {
487: PetscPrintf(PETSC_COMM_WORLD, "Error with MatDenseGetColumnVecWrite\n");
488: MatView(T2, NULL);
489: }
490: MatDestroy(&T);
491: MatDestroy(&T2);
493: /* Test with MatShell */
494: MatDuplicate(A, MAT_COPY_VALUES, &T);
495: MatConvert(T, MATSHELL, MAT_INITIAL_MATRIX, &T2);
496: MatDestroy(&T);
498: /* scale matrix */
499: MatScale(A, 2.0);
500: MatScale(T2, 2.0);
501: MatCreateVecs(A, &r, &l);
502: VecSetRandom(r, NULL);
503: VecSetRandom(l, NULL);
504: MatCreateVecs(T2, &rs, &ls);
505: VecCopy(r, rs);
506: VecCopy(l, ls);
507: if (testproj) {
508: MatDiagonalScale(A, r, r);
509: MatDiagonalScale(T2, rs, rs);
510: } else {
511: MatDiagonalScale(A, l, r);
512: MatDiagonalScale(T2, ls, rs);
513: }
514: MatDuplicate(A, MAT_COPY_VALUES, &T);
515: MatAXPY(A, 4.5, T, SAME_NONZERO_PATTERN);
516: MatAXPY(T2, 4.5, T, DIFFERENT_NONZERO_PATTERN);
517: MatMultEqual(T2, A, 10, &flg);
518: if (!flg) PetscPrintf(PETSC_COMM_WORLD, "Error with MATSHELL (MatMult)\n");
519: MatMultTransposeEqual(T2, A, 10, &flg);
520: if (!flg) PetscPrintf(PETSC_COMM_WORLD, "Error with MATSHELL (MatMultTranspose)\n");
521: MatDestroy(&T);
522: VecDestroy(&ls);
523: VecDestroy(&rs);
524: VecDestroy(&l);
525: VecDestroy(&r);
527: /* recompute projections, test reusage */
528: if (PtAP) MatPtAP(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &PtAP);
529: if (RARt) MatRARt(A, R, MAT_REUSE_MATRIX, PETSC_DEFAULT, &RARt);
530: if (testshellops) { /* test callbacks for user defined MatProducts */
531: MatShellSetMatProductOperation(T2, MATPRODUCT_AB, NULL, MyMatShellMatMultNumeric, NULL, MATDENSE, MATDENSE);
532: MatShellSetMatProductOperation(T2, MATPRODUCT_AB, NULL, MyMatShellMatMultNumeric, NULL, MATDENSECUDA, MATDENSECUDA);
533: MatShellSetMatProductOperation(T2, MATPRODUCT_AtB, NULL, MyMatTransposeShellMatMultNumeric, NULL, MATDENSE, MATDENSE);
534: MatShellSetMatProductOperation(T2, MATPRODUCT_AtB, NULL, MyMatTransposeShellMatMultNumeric, NULL, MATDENSECUDA, MATDENSECUDA);
535: MatShellSetMatProductOperation(T2, MATPRODUCT_ABt, NULL, MyMatShellMatTransposeMultNumeric, NULL, MATDENSE, MATDENSE);
536: MatShellSetMatProductOperation(T2, MATPRODUCT_ABt, NULL, MyMatShellMatTransposeMultNumeric, NULL, MATDENSECUDA, MATDENSECUDA);
537: if (testproj) {
538: MatShellSetMatProductOperation(T2, MATPRODUCT_PtAP, MyPtShellPMultSymbolic, MyPtShellPMultNumeric, proj_destroy, MATDENSE, MATSHELL);
539: MatShellSetMatProductOperation(T2, MATPRODUCT_PtAP, MyPtShellPMultSymbolic, MyPtShellPMultNumeric, proj_destroy, MATDENSECUDA, MATSHELL);
540: MatShellSetMatProductOperation(T2, MATPRODUCT_RARt, MyRShellRtMultSymbolic, MyRShellRtMultNumeric, proj_destroy, MATDENSE, MATSHELL);
541: MatShellSetMatProductOperation(T2, MATPRODUCT_RARt, MyRShellRtMultSymbolic, MyRShellRtMultNumeric, proj_destroy, MATDENSECUDA, MATSHELL);
542: }
543: }
544: CheckLocal(B, X, aB, aX);
545: /* we either use the shell operations or the loop over columns code, applying the operator */
546: MatMatMult(T2, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
547: CheckLocal(B, X, aB, aX);
548: MatMatMultEqual(T2, B, X, 10, &flg);
549: if (!flg) {
550: PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MATSHELL)\n");
551: MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
552: MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN);
553: MatView(T, NULL);
554: MatDestroy(&T);
555: }
556: if (testproj) {
557: MatPtAPMultEqual(T2, B, PtAP, 10, &flg);
558: if (!flg) PetscPrintf(PETSC_COMM_WORLD, "Error with PtAP (MATSHELL)\n");
559: if (testshellops) { /* projections fail if the product operations are not specified */
560: MatPtAP(T2, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
561: MatPtAP(T2, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &T);
562: MatPtAPMultEqual(T2, B, T, 10, &flg);
563: if (!flg) {
564: Mat TE;
566: PetscPrintf(PETSC_COMM_WORLD, "Error with PtAP (MATSHELL user defined)\n");
567: MatComputeOperator(T, MATDENSE, &TE);
568: MatView(TE, NULL);
569: MatView(PtAP, NULL);
570: MatAXPY(TE, -1.0, PtAP, SAME_NONZERO_PATTERN);
571: MatView(TE, NULL);
572: MatDestroy(&TE);
573: }
574: MatDestroy(&T);
575: }
576: if (RARt) {
577: MatRARtMultEqual(T2, R, RARt, 10, &flg);
578: if (!flg) PetscPrintf(PETSC_COMM_WORLD, "Error with RARt (MATSHELL)\n");
579: }
580: if (testshellops) {
581: MatRARt(T2, R, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
582: MatRARt(T2, R, MAT_REUSE_MATRIX, PETSC_DEFAULT, &T);
583: MatRARtMultEqual(T2, R, T, 10, &flg);
584: if (!flg) {
585: Mat TE;
587: PetscPrintf(PETSC_COMM_WORLD, "Error with RARt (MATSHELL user defined)\n");
588: MatComputeOperator(T, MATDENSE, &TE);
589: MatView(TE, NULL);
590: if (RARt) {
591: MatView(RARt, NULL);
592: MatAXPY(TE, -1.0, RARt, SAME_NONZERO_PATTERN);
593: MatView(TE, NULL);
594: }
595: MatDestroy(&TE);
596: }
597: MatDestroy(&T);
598: }
599: }
601: if (testmattmat) { /* we either use the shell operations or the loop over columns code applying the transposed operator */
602: MatTransposeMatMult(T2, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
603: CheckLocal(B, X, aB, aX);
604: MatTransposeMatMultEqual(T2, X, B, 10, &flg);
605: if (!flg) {
606: PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatTranspose, MATSHELL)\n");
607: MatTransposeMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
608: MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN);
609: MatView(T, NULL);
610: MatDestroy(&T);
611: }
612: }
613: if (testmatmatt && testshellops) { /* only when shell operations are set */
614: MatMatTransposeMult(T2, Bt, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
615: CheckLocal(Bt, X, aBt, aX);
616: MatMatTransposeMultEqual(T2, Bt, X, 10, &flg);
617: if (!flg) {
618: PetscPrintf(PETSC_COMM_WORLD, "Error with reusage (MatMatTranspose, MATSHELL)\n");
619: MatMatTransposeMult(A, Bt, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
620: MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN);
621: MatView(T, NULL);
622: MatDestroy(&T);
623: }
624: }
625: MatDestroy(&T2);
627: if (testnest) { /* test with MatNest */
628: Mat NA;
630: MatCreateNest(PETSC_COMM_WORLD, 1, NULL, 1, NULL, &A, &NA);
631: MatViewFromOptions(NA, NULL, "-NA_view");
632: MatMatMult(NA, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
633: CheckLocal(B, X, aB, aX);
634: MatMatMultEqual(NA, B, X, 10, &flg);
635: if (!flg) {
636: PetscPrintf(PETSC_COMM_WORLD, "Error with Nest\n");
637: MatMatMult(NA, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
638: MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN);
639: MatView(T, NULL);
640: MatDestroy(&T);
641: }
642: MatDestroy(&NA);
643: }
645: if (testtranspose) { /* test with Transpose */
646: Mat TA;
648: MatCreateTranspose(A, &TA);
649: MatMatMult(TA, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
650: CheckLocal(B, X, aB, aX);
651: MatMatMultEqual(TA, X, B, 10, &flg);
652: if (!flg) {
653: PetscPrintf(PETSC_COMM_WORLD, "Error with Transpose\n");
654: MatMatMult(TA, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
655: MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN);
656: MatView(T, NULL);
657: MatDestroy(&T);
658: }
659: MatDestroy(&TA);
660: }
662: if (testhtranspose) { /* test with Hermitian Transpose */
663: Mat TA;
665: MatCreateHermitianTranspose(A, &TA);
666: MatMatMult(TA, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
667: CheckLocal(B, X, aB, aX);
668: MatMatMultEqual(TA, X, B, 10, &flg);
669: if (!flg) {
670: PetscPrintf(PETSC_COMM_WORLD, "Error with Transpose\n");
671: MatMatMult(TA, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
672: MatAXPY(T, -1.0, B, SAME_NONZERO_PATTERN);
673: MatView(T, NULL);
674: MatDestroy(&T);
675: }
676: MatDestroy(&TA);
677: }
679: if (testtt) { /* test with Transpose(Transpose) */
680: Mat TA, TTA;
682: MatCreateTranspose(A, &TA);
683: MatCreateTranspose(TA, &TTA);
684: MatMatMult(TTA, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
685: CheckLocal(B, X, aB, aX);
686: MatMatMultEqual(TTA, B, X, 10, &flg);
687: if (!flg) {
688: PetscPrintf(PETSC_COMM_WORLD, "Error with Transpose(Transpose)\n");
689: MatMatMult(TTA, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
690: MatAXPY(T, -1.0, X, SAME_NONZERO_PATTERN);
691: MatView(T, NULL);
692: MatDestroy(&T);
693: }
694: MatDestroy(&TA);
695: MatDestroy(&TTA);
696: }
698: if (testcircular) { /* test circular */
699: Mat AB;
701: MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AB);
702: MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &X);
703: CheckLocal(B, X, aB, aX);
704: if (M == N && N == K) {
705: MatMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
706: } else {
707: MatTransposeMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
708: }
709: CheckLocal(B, X, aB, aX);
710: MatDestroy(&AB);
711: }
713: /* Test by Pierre Jolivet */
714: {
715: Mat C, D, D2, AtA;
716: MatCreateNormal(A, &AtA);
717: MatDuplicate(X, MAT_DO_NOT_COPY_VALUES, &C);
718: MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &D);
719: MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &D2);
720: MatSetRandom(B, NULL);
721: MatSetRandom(C, NULL);
722: MatSetRandom(D, NULL);
723: MatSetRandom(D2, NULL);
724: MatProductCreateWithMat(A, B, NULL, C);
725: MatProductSetType(C, MATPRODUCT_AB);
726: MatProductSetFromOptions(C);
727: MatProductSymbolic(C);
728: MatProductCreateWithMat(A, C, NULL, D);
729: MatProductSetType(D, MATPRODUCT_AtB);
730: MatProductSetFromOptions(D);
731: MatProductSymbolic(D);
732: MatProductNumeric(C);
733: MatMatMultEqual(A, B, C, 10, &flg);
734: if (!flg) {
735: PetscPrintf(PETSC_COMM_WORLD, "Error with Normal (AB != C)\n");
736: MatView(A, NULL);
737: MatView(B, NULL);
738: MatView(C, NULL);
739: }
740: MatProductNumeric(D);
741: MatMatMultEqual(AtA, B, D, 10, &flg);
742: if (!flg) {
743: PetscPrintf(PETSC_COMM_WORLD, "Error with Normal (2)\n");
744: MatMatMult(AtA, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &T);
745: MatView(D, NULL);
746: MatView(T, NULL);
747: MatAXPY(T, -1.0, D, SAME_NONZERO_PATTERN);
748: MatView(T, NULL);
749: MatDestroy(&T);
750: }
751: MatDestroy(&C);
752: MatDestroy(&D);
753: MatDestroy(&D2);
754: MatDestroy(&AtA);
755: }
757: MatDestroy(&X);
758: MatDestroy(&Bt);
759: MatDestroy(&B);
760: MatDestroy(&A);
761: MatDestroy(&R);
762: MatDestroy(&PtAP);
763: MatDestroy(&RARt);
764: PetscFree(dataX);
765: PetscFree(dataB);
766: PetscFree(dataR);
767: PetscFree(dataBt);
768: PetscFinalize();
769: return 0;
770: }
772: /*TEST
774: test:
775: suffix: 1
776: args: -local {{0 1}} -testshellops
778: test:
779: output_file: output/ex70_1.out
780: requires: cuda
781: suffix: 1_cuda
782: args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{seqaijcusparse seqaij}} -testshellops {{0 1}}
784: test:
785: output_file: output/ex70_1.out
786: nsize: 2
787: suffix: 1_par
788: args: -local {{0 1}} -testmatmatt 0
790: test:
791: output_file: output/ex70_1.out
792: requires: cuda
793: nsize: 2
794: suffix: 1_par_cuda
795: args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{mpiaijcusparse mpiaij}} -testnest 0 -testmatmatt 0 -matmatmult_Bbn 3
797: test:
798: output_file: output/ex70_1.out
799: suffix: 2
800: nsize: 1
801: args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}}
803: testset:
804: requires: cuda
805: output_file: output/ex70_1.out
806: nsize: 1
807: args: -M 7 -N 9 -K 2 -local {{0 1}} -testnest 0 -A_mat_type {{seqdensecuda seqdense}} -xgpu {{0 1}} -bgpu {{0 1}}
808: test:
809: requires: !complex
810: suffix: 2_cuda_real
811: test:
812: # complex+single gives a little bigger error in the MatDenseGetColumnVec test
813: requires: complex !single
814: suffix: 2_cuda_complex
816: test:
817: output_file: output/ex70_1.out
818: suffix: 2_par
819: nsize: 2
820: args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} -testcircular -testmatmatt 0
822: test:
823: requires: cuda
824: output_file: output/ex70_1.out
825: suffix: 2_par_cuda
826: nsize: 2
827: args: -M 11 -N 9 -K 1 -local {{0 1}} -testcircular 0 -A_mat_type mpiaijcusparse -xgpu -bgpu -testnest 0 -testmatmatt 0
829: test:
830: output_file: output/ex70_1.out
831: suffix: 3
832: nsize: {{1 3}}
833: args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type sbaij -symm -testproj 0 -testmatmatt 0
835: test:
836: output_file: output/ex70_1.out
837: suffix: 4
838: nsize: 1
839: args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular
841: test:
842: output_file: output/ex70_1.out
843: suffix: 5
844: nsize: {{2 4}}
845: args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular -testmatmatt 0
847: test:
848: output_file: output/ex70_1.out
849: suffix: 6
850: nsize: 1
851: args: -M {{1 3}} -N {{2 5}} -K {{1 2}} -local {{0 1}} -testcircular
853: test:
854: output_file: output/ex70_1.out
855: suffix: 7
856: nsize: 1
857: args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type dense -testnest -testcircular
859: TEST*/