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