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