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