Actual source code: ex70.c
petsc-3.13.6 2020-09-29
1: #include <petscmat.h>
3: static char help[] = "Tests MatMatMult with MAT_REUSE_MATRIX and already allocated dense result.\n\n";
5: static PetscErrorCode CheckLocal(Mat A, Mat B, PetscScalar *a, PetscScalar *b)
6: {
8: PetscBool wA = PETSC_FALSE, wB = PETSC_FALSE;
11: if (a) {
12: const PetscScalar *Aa;
13: MatDenseGetArrayRead(A,&Aa);
14: wA = (PetscBool)(a != Aa);
15: MatDenseRestoreArrayRead(A,&Aa);
16: }
17: if (b) {
18: const PetscScalar *Bb;
19: MatDenseGetArrayRead(B,&Bb);
20: wB = (PetscBool)(b != Bb);
21: MatDenseRestoreArrayRead(B,&Bb);
22: }
23: if (wA || wB) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong array in A Mat %d, Wrong array in B Mat %d",wA,wB);
24: return(0);
25: }
27: int main(int argc,char **args)
28: {
29: Mat X,B,A,T;
30: PetscInt m,n,M = 10,N = 10,K = 5, ldx = 0, ldb = 0;
31: const char *deft = MATAIJ;
32: char mattype[256];
33: PetscBool flg,symm = PETSC_FALSE,testtt = PETSC_TRUE, testnest = PETSC_TRUE, testtranspose = PETSC_TRUE, testcircular = PETSC_FALSE, local = PETSC_TRUE;
34: PetscScalar *dataX = NULL,*dataB = NULL;
35: PetscScalar *aX,*aB;
38: PetscInitialize(&argc,&args,NULL,help);if (ierr) return ierr;
39: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
40: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
41: PetscOptionsGetInt(NULL,NULL,"-K",&K,NULL);
42: PetscOptionsGetBool(NULL,NULL,"-symm",&symm,NULL);
43: PetscOptionsGetBool(NULL,NULL,"-local",&local,NULL);
44: PetscOptionsGetInt(NULL,NULL,"-ldx",&ldx,NULL);
45: PetscOptionsGetInt(NULL,NULL,"-ldb",&ldb,NULL);
46: PetscOptionsGetBool(NULL,NULL,"-testtranspose",&testtranspose,NULL);
47: PetscOptionsGetBool(NULL,NULL,"-testnest",&testnest,NULL);
48: PetscOptionsGetBool(NULL,NULL,"-testtt",&testtt,NULL);
49: PetscOptionsGetBool(NULL,NULL,"-testcircular",&testcircular,NULL);
50: #if defined(PETSC_USE_COMPLEX)
51: testtranspose = PETSC_FALSE;
52: testtt = PETSC_FALSE;
53: #endif
54: MatCreate(PETSC_COMM_WORLD,&A);
55: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
56: MatSetType(A,MATAIJ);
57: MatSetUp(A);
58: MatSetRandom(A,NULL);
59: if (M==N && symm) {
60: Mat AT;
62: MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&AT);
63: MatAXPY(A,1.0,AT,DIFFERENT_NONZERO_PATTERN);
64: MatDestroy(&AT);
65: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
66: }
67: MatViewFromOptions(A,NULL,"-A_init_view");
68: PetscOptionsBegin(PETSC_COMM_WORLD,"","","");
69: PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,mattype,256,&flg);
70: PetscOptionsEnd();
71: if (flg) {
72: Mat A2;
74: MatConvert(A,mattype,MAT_INITIAL_MATRIX,&A2);
75: MatMultEqual(A,A2,10,&flg);
76: if (!flg) {
77: Mat AE,A2E;
79: PetscPrintf(PETSC_COMM_WORLD,"Error with convert\n");
80: MatComputeOperator(A,MATDENSE,&AE);
81: MatComputeOperator(A2,MATDENSE,&A2E);
82: MatView(AE,NULL);
83: MatView(A2E,NULL);
84: MatAXPY(A2E,-1.0,A,SAME_NONZERO_PATTERN);
85: MatView(A2E,NULL);
86: MatDestroy(&A2E);
87: MatDestroy(&AE);
88: }
89: MatDestroy(&A);
90: A = A2;
91: }
92: MatViewFromOptions(A,NULL,"-A_view");
94: MatGetLocalSize(A,&m,&n);
95: if (local) {
96: PetscMalloc1(PetscMax((m+ldx)*K,1),&dataX);
97: PetscMalloc1(PetscMax((n+ldb)*K,1),&dataB);
98: }
99: MatCreateDense(PETSC_COMM_WORLD,n,PETSC_DECIDE,N,K,dataB,&B);
100: MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,M,K,dataX,&X);
102: /* store pointer to dense data for testing */
103: MatDenseGetArrayRead(B,(const PetscScalar**)&dataB);
104: MatDenseGetArrayRead(X,(const PetscScalar**)&dataX);
105: aX = dataX;
106: aB = dataB;
107: MatDenseRestoreArrayRead(B,(const PetscScalar**)&dataB);
108: MatDenseRestoreArrayRead(X,(const PetscScalar**)&dataX);
109: if (local) {
110: dataX = aX;
111: dataB = aB;
112: }
113: MatSetRandom(B,NULL);
115: /* Test reusing a previously allocated dense buffer */
116: MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);
117: CheckLocal(B,X,aB,aX);
118: MatMatMultEqual(A,B,X,10,&flg);
119: if (!flg) {
120: PetscPrintf(PETSC_COMM_WORLD,"Error with reusage\n");
121: MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);
122: MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);
123: MatView(T,NULL);
124: MatDestroy(&T);
125: }
127: if (testnest) { /* test with MatNest */
128: Mat NA;
130: MatCreateNest(PETSC_COMM_WORLD,1,NULL,1,NULL,&A,&NA);
131: MatViewFromOptions(NA,NULL,"-NA_view");
132: MatMatMult(NA,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);
133: CheckLocal(B,X,aB,aX);
134: MatMatMultEqual(NA,B,X,10,&flg);
135: if (!flg) {
136: PetscPrintf(PETSC_COMM_WORLD,"Error with Nest\n");
137: MatMatMult(NA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);
138: MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);
139: MatView(T,NULL);
140: MatDestroy(&T);
141: }
142: MatDestroy(&NA);
143: }
145: if (testtranspose) { /* test with Transpose */
146: Mat TA;
148: MatCreateHermitianTranspose(A,&TA);
149: MatMatMult(TA,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);
150: CheckLocal(B,X,aB,aX);
151: MatMatMultEqual(TA,X,B,10,&flg);
152: if (!flg) {
153: PetscPrintf(PETSC_COMM_WORLD,"Error with Transpose\n");
154: MatMatMult(TA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);
155: MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);
156: MatView(T,NULL);
157: MatDestroy(&T);
158: }
159: MatDestroy(&TA);
160: }
162: if (testtt) { /* test with Transpose(Transpose) */
163: Mat TA, TTA;
165: MatCreateHermitianTranspose(A,&TA);
166: MatCreateHermitianTranspose(TA,&TTA);
167: MatMatMult(TTA,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);
168: CheckLocal(B,X,aB,aX);
169: MatMatMultEqual(TTA,B,X,10,&flg);
170: if (!flg) {
171: PetscPrintf(PETSC_COMM_WORLD,"Error with Transpose(Transpose)\n");
172: MatMatMult(TTA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);
173: MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);
174: MatView(T,NULL);
175: MatDestroy(&T);
176: }
177: MatDestroy(&TA);
178: MatDestroy(&TTA);
179: }
181: if (testcircular) { /* test circular */
182: Mat AB;
184: MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AB);
185: MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);
186: CheckLocal(B,X,aB,aX);
187: if (M == N && N == K) {
188: MatMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);
189: } else {
190: MatTransposeMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);
191: }
192: CheckLocal(B,X,aB,aX);
193: MatDestroy(&AB);
194: }
195: MatDestroy(&X);
196: MatDestroy(&B);
197: MatDestroy(&A);
198: PetscFree(dataX);
199: PetscFree(dataB);
200: PetscFinalize();
201: return ierr;
202: }
204: /*TEST
206: test:
207: suffix: 1
208: args: -local {{0 1}}
210: test:
211: output_file: output/ex70_1.out
212: nsize: 2
213: suffix: 1_par
214: args: -testtranspose 0 -local {{0 1}}
216: test:
217: TODO: MPIAIJ x MPIDENSE broken for MatTransposeMatMult
218: output_file: output/ex70_1.out
219: nsize: 2
220: suffix: 1_par_broken
221: args: -testtranspose -local {{0 1}}
223: test:
224: output_file: output/ex70_1.out
225: suffix: 2
226: nsize: 1
227: args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}}
229: test:
230: output_file: output/ex70_1.out
231: suffix: 2_par
232: nsize: 2
233: args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} -testcircular 0
235: test:
236: TODO: MatTransposeMatMultSymbolic_SeqAIJ_SeqDense plays with the destroy routine
237: output_file: output/ex70_1.out
238: suffix: 2_broken
239: nsize: 2
240: args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} -testcircular 1
242: test:
243: output_file: output/ex70_1.out
244: suffix: 3
245: nsize: 1
246: args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -mat_type sbaij -symm -testtranspose 0
248: test:
249: output_file: output/ex70_1.out
250: suffix: 4
251: nsize: 1
252: args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular
254: test:
255: output_file: output/ex70_1.out
256: suffix: 5
257: nsize: {{2 4}}
258: args: -M 3 -N 3 -K 3 -local 1 -testcircular -testtranspose 0
260: test:
261: TODO: MatCreateDense broken with some processors not having local rows
262: output_file: output/ex70_1.out
263: suffix: 5_broken
264: nsize: {{2 4}}
265: args: -M 3 -N 3 -K 3 -local 0 -testcircular -testtranspose 0
267: test:
268: output_file: output/ex70_1.out
269: suffix: 6
270: nsize: 1
271: args: -M {{1 3}} -N {{2 5}} -K {{1 2}} -local {{0 1}} -testcircular -testtranspose 0
273: test:
274: TODO: MatTransposeMatMultSymbolic_SeqAIJ_SeqDense plays with the destroy routine
275: output_file: output/ex70_1.out
276: suffix: 6_broken
277: nsize: 1
278: args: -M {{1 3}} -N {{2 5}} -K {{1 2}} -local {{0 1}} -testcircular -testtranspose
280: test:
281: output_file: output/ex70_1.out
282: suffix: 7
283: nsize: 1
284: args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -mat_type dense -testnest 0 -testcircular
286: test:
287: TODO: NEST x DENSE with dense nested matrices seems broken in this case
288: output_file: output/ex70_1.out
289: suffix: 7_broken
290: nsize: 1
291: args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -mat_type dense -testnest -testcircular
292: TEST*/