Actual source code: ex221.c
petsc-3.14.6 2021-03-30
1: static char help[] = "Tests various routines for MATSHELL\n\n";
3: #include <petscmat.h>
5: typedef struct _n_User *User;
6: struct _n_User {
7: Mat B;
8: };
10: static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X)
11: {
12: User user;
16: MatShellGetContext(A,&user);
17: MatGetDiagonal(user->B,X);
18: return(0);
19: }
21: static PetscErrorCode MatMult_User(Mat A,Vec X,Vec Y)
22: {
23: User user;
27: MatShellGetContext(A,&user);
28: MatMult(user->B,X,Y);
29: return(0);
30: }
32: static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y)
33: {
34: User user;
38: MatShellGetContext(A,&user);
39: MatMultTranspose(user->B,X,Y);
40: return(0);
41: }
43: static PetscErrorCode MatCopy_User(Mat A,Mat X,MatStructure str)
44: {
45: User user,userX;
49: MatShellGetContext(A,&user);
50: MatShellGetContext(X,&userX);
51: if (user != userX) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen");
52: PetscObjectReference((PetscObject)user->B);
53: return(0);
54: }
56: static PetscErrorCode MatDestroy_User(Mat A)
57: {
58: User user;
62: MatShellGetContext(A,&user);
63: PetscObjectDereference((PetscObject)user->B);
64: return(0);
65: }
67: int main(int argc,char **args)
68: {
69: User user;
70: Mat A,S;
71: PetscScalar *data,diag = 1.3;
72: PetscReal tol = PETSC_SMALL;
73: PetscInt i,j,m = PETSC_DECIDE,n = PETSC_DECIDE,M = 17,N = 15,s1,s2;
74: PetscInt test, ntest = 2;
75: PetscMPIInt rank,size;
76: PetscBool nc = PETSC_FALSE, cong, flg;
77: PetscBool ronl = PETSC_TRUE;
78: PetscBool randomize = PETSC_FALSE, submat = PETSC_FALSE;
79: PetscBool keep = PETSC_FALSE;
80: PetscBool testzerorows = PETSC_TRUE, testdiagscale = PETSC_TRUE, testgetdiag = PETSC_TRUE, testsubmat = PETSC_TRUE;
81: PetscBool testshift = PETSC_TRUE, testscale = PETSC_TRUE, testdup = PETSC_TRUE, testreset = PETSC_TRUE;
82: PetscBool testaxpy = PETSC_TRUE, testaxpyd = PETSC_TRUE, testaxpyerr = PETSC_FALSE;
85: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
86: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
87: MPI_Comm_size(PETSC_COMM_WORLD,&size);
88: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
89: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
90: PetscOptionsGetInt(NULL,NULL,"-ml",&m,NULL);
91: PetscOptionsGetInt(NULL,NULL,"-nl",&n,NULL);
92: PetscOptionsGetBool(NULL,NULL,"-square_nc",&nc,NULL);
93: PetscOptionsGetBool(NULL,NULL,"-rows_only",&ronl,NULL);
94: PetscOptionsGetBool(NULL,NULL,"-randomize",&randomize,NULL);
95: PetscOptionsGetBool(NULL,NULL,"-submat",&submat,NULL);
96: PetscOptionsGetBool(NULL,NULL,"-test_zerorows",&testzerorows,NULL);
97: PetscOptionsGetBool(NULL,NULL,"-test_diagscale",&testdiagscale,NULL);
98: PetscOptionsGetBool(NULL,NULL,"-test_getdiag",&testgetdiag,NULL);
99: PetscOptionsGetBool(NULL,NULL,"-test_shift",&testshift,NULL);
100: PetscOptionsGetBool(NULL,NULL,"-test_scale",&testscale,NULL);
101: PetscOptionsGetBool(NULL,NULL,"-test_dup",&testdup,NULL);
102: PetscOptionsGetBool(NULL,NULL,"-test_reset",&testreset,NULL);
103: PetscOptionsGetBool(NULL,NULL,"-test_submat",&testsubmat,NULL);
104: PetscOptionsGetBool(NULL,NULL,"-test_axpy",&testaxpy,NULL);
105: PetscOptionsGetBool(NULL,NULL,"-test_axpy_different",&testaxpyd,NULL);
106: PetscOptionsGetBool(NULL,NULL,"-test_axpy_error",&testaxpyerr,NULL);
107: PetscOptionsGetInt(NULL,NULL,"-loop",&ntest,NULL);
108: PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL);
109: PetscOptionsGetScalar(NULL,NULL,"-diag",&diag,NULL);
110: PetscOptionsGetBool(NULL,NULL,"-keep",&keep,NULL);
111: /* This tests square matrices with different row/col layout */
112: if (nc && size > 1) {
113: M = PetscMax(PetscMax(N,M),1);
114: N = M;
115: m = n = 0;
116: if (rank == 0) { m = M-1; n = 1; }
117: else if (rank == 1) { m = 1; n = N-1; }
118: }
119: MatCreateDense(PETSC_COMM_WORLD,m,n,M,N,NULL,&A);
120: MatGetLocalSize(A,&m,&n);
121: MatGetSize(A,&M,&N);
122: MatGetOwnershipRange(A,&s1,NULL);
123: s2 = 1;
124: while (s2 < M) s2 *= 10;
125: MatDenseGetArray(A,&data);
126: for (j = 0; j < N; j++) {
127: for (i = 0; i < m; i++) {
128: data[j*m + i] = s2*j + i + s1 + 1;
129: }
130: }
131: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
132: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
134: if (submat) {
135: Mat A2;
136: IS r,c;
137: PetscInt rst,ren,cst,cen;
139: MatGetOwnershipRange(A,&rst,&ren);
140: MatGetOwnershipRangeColumn(A,&cst,&cen);
141: ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r);
142: ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c);
143: MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&A2);
144: ISDestroy(&r);
145: ISDestroy(&c);
146: MatDestroy(&A);
147: A = A2;
148: }
150: MatGetSize(A,&M,&N);
151: MatGetLocalSize(A,&m,&n);
152: MatHasCongruentLayouts(A,&cong);
154: MatConvert(A,MATAIJ,MAT_INPLACE_MATRIX,&A);
155: MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,keep);
156: PetscObjectSetName((PetscObject)A,"initial");
157: MatViewFromOptions(A,NULL,"-view_mat");
159: PetscNew(&user);
160: MatCreateShell(PETSC_COMM_WORLD,m,n,M,N,user,&S);
161: MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User);
162: MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User);
163: if (cong) {
164: MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User);
165: }
166: MatShellSetOperation(S,MATOP_COPY,(void (*)(void))MatCopy_User);
167: MatShellSetOperation(S,MATOP_DESTROY,(void (*)(void))MatDestroy_User);
168: MatDuplicate(A,MAT_COPY_VALUES,&user->B);
170: /* Square and rows only scaling */
171: ronl = cong ? ronl : PETSC_TRUE;
173: for (test = 0; test < ntest; test++) {
174: PetscReal err;
176: MatMultAddEqual(A,S,10,&flg);
177: if (!flg) {
178: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mult add\n",test);
179: }
180: MatMultTransposeAddEqual(A,S,10,&flg);
181: if (!flg) {
182: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mult add (T)\n",test);
183: }
184: if (testzerorows) {
185: Mat ST,B,C,BT,BTT;
186: IS zr;
187: Vec x = NULL, b1 = NULL, b2 = NULL;
188: PetscInt *idxs = NULL, nr = 0;
190: if (rank == (test%size)) {
191: nr = 1;
192: PetscMalloc1(nr,&idxs);
193: if (test%2) {
194: idxs[0] = (2*M - 1 - test/2)%M;
195: } else {
196: idxs[0] = (test/2)%M;
197: }
198: idxs[0] = PetscMax(idxs[0],0);
199: }
200: ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr);
201: PetscObjectSetName((PetscObject)zr,"ZR");
202: ISViewFromOptions(zr,NULL,"-view_is");
203: MatCreateVecs(A,&x,&b1);
204: if (randomize) {
205: VecSetRandom(x,NULL);
206: VecSetRandom(b1,NULL);
207: } else {
208: VecSet(x,11.4);
209: VecSet(b1,-14.2);
210: }
211: VecDuplicate(b1,&b2);
212: VecCopy(b1,b2);
213: PetscObjectSetName((PetscObject)b1,"A_B1");
214: PetscObjectSetName((PetscObject)b2,"A_B2");
215: if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */
216: VecDestroy(&b1);
217: }
218: if (ronl) {
219: MatZeroRowsIS(A,zr,diag,x,b1);
220: MatZeroRowsIS(S,zr,diag,x,b2);
221: } else {
222: MatZeroRowsColumnsIS(A,zr,diag,x,b1);
223: MatZeroRowsColumnsIS(S,zr,diag,x,b2);
224: ISDestroy(&zr);
225: /* Mix zerorows and zerorowscols */
226: nr = 0;
227: idxs = NULL;
228: if (!rank) {
229: nr = 1;
230: PetscMalloc1(nr,&idxs);
231: if (test%2) {
232: idxs[0] = (3*M - 2 - test/2)%M;
233: } else {
234: idxs[0] = (test/2+1)%M;
235: }
236: idxs[0] = PetscMax(idxs[0],0);
237: }
238: ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr);
239: PetscObjectSetName((PetscObject)zr,"ZR2");
240: ISViewFromOptions(zr,NULL,"-view_is");
241: MatZeroRowsIS(A,zr,diag*2.0+PETSC_SMALL,NULL,NULL);
242: MatZeroRowsIS(S,zr,diag*2.0+PETSC_SMALL,NULL,NULL);
243: }
244: ISDestroy(&zr);
246: if (b1) {
247: Vec b;
249: VecViewFromOptions(b1,NULL,"-view_b");
250: VecViewFromOptions(b2,NULL,"-view_b");
251: VecDuplicate(b1,&b);
252: VecCopy(b1,b);
253: VecAXPY(b,-1.0,b2);
254: VecNorm(b,NORM_INFINITY,&err);
255: if (err >= tol) {
256: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error b %g\n",test,(double)err);
257: }
258: VecDestroy(&b);
259: }
260: VecDestroy(&b1);
261: VecDestroy(&b2);
262: VecDestroy(&x);
263: MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&B);
265: MatCreateTranspose(S,&ST);
266: MatComputeOperator(ST,MATDENSE,&BT);
267: MatTranspose(BT,MAT_INITIAL_MATRIX,&BTT);
268: PetscObjectSetName((PetscObject)B,"S");
269: PetscObjectSetName((PetscObject)BTT,"STT");
270: MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&C);
271: PetscObjectSetName((PetscObject)C,"A");
273: MatViewFromOptions(C,NULL,"-view_mat");
274: MatViewFromOptions(B,NULL,"-view_mat");
275: MatViewFromOptions(BTT,NULL,"-view_mat");
277: MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN);
278: MatNorm(C,NORM_FROBENIUS,&err);
279: if (err >= tol) {
280: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat mult after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err);
281: }
283: MatConvert(A,MATDENSE,MAT_REUSE_MATRIX,&C);
284: MatAXPY(C,-1.0,BTT,SAME_NONZERO_PATTERN);
285: MatNorm(C,NORM_FROBENIUS,&err);
286: if (err >= tol) {
287: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat mult transpose after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err);
288: }
290: MatDestroy(&ST);
291: MatDestroy(&BTT);
292: MatDestroy(&BT);
293: MatDestroy(&B);
294: MatDestroy(&C);
295: }
296: if (testdiagscale) { /* MatDiagonalScale() */
297: Vec vr,vl;
299: MatCreateVecs(A,&vr,&vl);
300: if (randomize) {
301: VecSetRandom(vr,NULL);
302: VecSetRandom(vl,NULL);
303: } else {
304: VecSet(vr,test%2 ? 0.15 : 1.0/0.15);
305: VecSet(vl,test%2 ? -1.2 : 1.0/-1.2);
306: }
307: MatDiagonalScale(A,vl,vr);
308: MatDiagonalScale(S,vl,vr);
309: VecDestroy(&vr);
310: VecDestroy(&vl);
311: }
313: if (testscale) { /* MatScale() */
314: MatScale(A,test%2 ? 1.4 : 1.0/1.4);
315: MatScale(S,test%2 ? 1.4 : 1.0/1.4);
316: }
318: if (testshift && cong) { /* MatShift() : MATSHELL shift is broken when row/cols layout are not congruent and left/right scaling have been applied */
319: MatShift(A,test%2 ? -77.5 : 77.5);
320: MatShift(S,test%2 ? -77.5 : 77.5);
321: }
323: if (testgetdiag && cong) { /* MatGetDiagonal() */
324: Vec dA,dS;
326: MatCreateVecs(A,&dA,NULL);
327: MatCreateVecs(S,&dS,NULL);
328: MatGetDiagonal(A,dA);
329: MatGetDiagonal(S,dS);
330: VecAXPY(dA,-1.0,dS);
331: VecNorm(dA,NORM_INFINITY,&err);
332: if (err >= tol) {
333: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error diag %g\n",test,(double)err);
334: }
335: VecDestroy(&dA);
336: VecDestroy(&dS);
337: }
339: if (testdup && !test) {
340: Mat A2, S2;
342: MatDuplicate(A,MAT_COPY_VALUES,&A2);
343: MatDuplicate(S,MAT_COPY_VALUES,&S2);
344: MatDestroy(&A);
345: MatDestroy(&S);
346: A = A2;
347: S = S2;
348: }
350: if (testsubmat) {
351: Mat sA,sS,dA,dS,At,St;
352: IS r,c;
353: PetscInt rst,ren,cst,cen;
355: MatGetOwnershipRange(A,&rst,&ren);
356: MatGetOwnershipRangeColumn(A,&cst,&cen);
357: ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r);
358: ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c);
359: MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&sA);
360: MatCreateSubMatrix(S,r,c,MAT_INITIAL_MATRIX,&sS);
361: MatMultAddEqual(sA,sS,10,&flg);
362: if (!flg) {
363: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error submatrix mult add\n",test);
364: }
365: MatMultTransposeAddEqual(sA,sS,10,&flg);
366: if (!flg) {
367: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error submatrix mult add (T)\n",test);
368: }
369: MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA);
370: MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS);
371: MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN);
372: MatNorm(dA,NORM_FROBENIUS,&err);
373: if (err >= tol) {
374: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat submatrix %g\n",test,(double)err);
375: }
376: MatDestroy(&sA);
377: MatDestroy(&sS);
378: MatDestroy(&dA);
379: MatDestroy(&dS);
380: MatCreateTranspose(A,&At);
381: MatCreateTranspose(S,&St);
382: MatCreateSubMatrix(At,c,r,MAT_INITIAL_MATRIX,&sA);
383: MatCreateSubMatrix(St,c,r,MAT_INITIAL_MATRIX,&sS);
384: MatMultAddEqual(sA,sS,10,&flg);
385: if (!flg) {
386: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error submatrix (T) mult add\n",test);
387: }
388: MatMultTransposeAddEqual(sA,sS,10,&flg);
389: if (!flg) {
390: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error submatrix (T) mult add (T)\n",test);
391: }
392: MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA);
393: MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS);
394: MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN);
395: MatNorm(dA,NORM_FROBENIUS,&err);
396: if (err >= tol) {
397: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat submatrix (T) %g\n",test,(double)err);
398: }
399: MatDestroy(&sA);
400: MatDestroy(&sS);
401: MatDestroy(&dA);
402: MatDestroy(&dS);
403: MatDestroy(&At);
404: MatDestroy(&St);
405: ISDestroy(&r);
406: ISDestroy(&c);
407: }
409: if (testaxpy) {
410: Mat tA,tS,dA,dS;
411: MatStructure str[3] = { SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN };
413: MatDuplicate(A,MAT_COPY_VALUES,&tA);
414: if (testaxpyd && !(test%2)) {
415: PetscObjectReference((PetscObject)tA);
416: tS = tA;
417: } else {
418: PetscObjectReference((PetscObject)S);
419: tS = S;
420: }
421: MatAXPY(A,0.5,tA,str[test%3]);
422: MatAXPY(S,0.5,tS,str[test%3]);
423: /* this will trigger an error the next MatMult or MatMultTranspose call for S */
424: if (testaxpyerr) { MatScale(tA,0); }
425: MatDestroy(&tA);
426: MatDestroy(&tS);
427: MatMultAddEqual(A,S,10,&flg);
428: if (!flg) {
429: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error axpy mult add\n",test);
430: }
431: MatMultTransposeAddEqual(A,S,10,&flg);
432: if (!flg) {
433: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error axpy mult add (T)\n",test);
434: }
435: MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&dA);
436: MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&dS);
437: MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN);
438: MatNorm(dA,NORM_FROBENIUS,&err);
439: if (err >= tol) {
440: PetscPrintf(PETSC_COMM_WORLD,"[test %D] Error mat submatrix %g\n",test,(double)err);
441: }
442: MatDestroy(&dA);
443: MatDestroy(&dS);
444: }
446: if (testreset && (ntest == 1 || test == ntest-2)) {
447: /* reset MATSHELL */
448: MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY);
449: MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY);
450: /* reset A */
451: MatCopy(user->B,A,DIFFERENT_NONZERO_PATTERN);
452: }
453: }
455: MatDestroy(&A);
456: MatDestroy(&S);
457: PetscFree(user);
458: PetscFinalize();
459: return ierr;
460: }
463: /*TEST
465: testset:
466: suffix: rect
467: requires: !single
468: output_file: output/ex221_1.out
469: nsize: {{1 3}}
470: args: -loop 3 -keep {{0 1}} -M {{12 19}} -N {{19 12}} -submat {{0 1}} -test_axpy_different {{0 1}}
472: testset:
473: suffix: square
474: requires: !single
475: output_file: output/ex221_1.out
476: nsize: {{1 3}}
477: args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}} -submat {{0 1}} -test_axpy_different {{0 1}}
478: TEST*/