Actual source code: ex202.c
1: static char help[] = "Tests the use of MatTranspose_Nest and MatMatMult_Nest_Dense\n";
3: #include <petscmat.h>
5: PetscErrorCode TestInitialMatrix(void)
6: {
7: const PetscInt nr = 2, nc = 3, nk = 10;
8: PetscInt n, N, m, M;
9: const PetscInt arow[2 * 3] = {2, 2, 2, 3, 3, 3};
10: const PetscInt acol[2 * 3] = {3, 2, 4, 3, 2, 4};
11: Mat A, Atranspose, B, C;
12: Mat subs[2 * 3], **block;
13: Vec x, y, Ax, ATy;
14: PetscInt i, j;
15: PetscScalar dot1, dot2, zero = 0.0, one = 1.0, *valsB, *valsC;
16: PetscReal norm;
17: PetscRandom rctx;
18: PetscBool equal;
20: PetscFunctionBegin;
21: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
22: /* Force the random numbers to have imaginary part 0 so printed results are the same for --with-scalar-type=real or --with-scalar-type=complex */
23: PetscCall(PetscRandomSetInterval(rctx, zero, one));
24: PetscCall(PetscRandomSetFromOptions(rctx));
25: for (i = 0; i < (nr * nc); i++) PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD, arow[i], acol[i], NULL, &subs[i]));
26: PetscCall(MatCreateNest(PETSC_COMM_WORLD, nr, NULL, nc, NULL, subs, &A));
27: PetscCall(MatCreateVecs(A, &x, NULL));
28: PetscCall(MatCreateVecs(A, NULL, &y));
29: PetscCall(VecDuplicate(x, &ATy));
30: PetscCall(VecDuplicate(y, &Ax));
31: PetscCall(MatSetRandom(A, rctx));
32: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atranspose));
34: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
35: PetscCall(MatNestGetSubMats(A, NULL, NULL, &block));
36: for (i = 0; i < nr; i++) {
37: for (j = 0; j < nc; j++) PetscCall(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD));
38: }
40: PetscCall(MatView(Atranspose, PETSC_VIEWER_STDOUT_WORLD));
41: PetscCall(MatNestGetSubMats(Atranspose, NULL, NULL, &block));
42: for (i = 0; i < nc; i++) {
43: for (j = 0; j < nr; j++) PetscCall(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD));
44: }
46: /* Check <Ax, y> = <x, A^Ty> */
47: for (i = 0; i < 10; i++) {
48: PetscCall(VecSetRandom(x, rctx));
49: PetscCall(VecSetRandom(y, rctx));
51: PetscCall(MatMult(A, x, Ax));
52: PetscCall(VecDot(Ax, y, &dot1));
53: PetscCall(MatMult(Atranspose, y, ATy));
54: PetscCall(VecDot(ATy, x, &dot2));
56: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "<Ax, y> = %g\n", (double)PetscRealPart(dot1)));
57: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "<x, A^Ty> = %g\n", (double)PetscRealPart(dot2)));
58: }
59: PetscCall(VecDestroy(&x));
60: PetscCall(VecDestroy(&y));
61: PetscCall(VecDestroy(&Ax));
63: PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD, acol[0] + acol[nr] + acol[2 * nr], nk, NULL, &B));
64: PetscCall(MatSetRandom(B, rctx));
65: PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
66: PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &C));
67: PetscCall(MatMatMultEqual(A, B, C, 10, &equal));
68: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in C != A*B");
70: PetscCall(MatGetSize(A, &M, &N));
71: PetscCall(MatGetLocalSize(A, &m, &n));
72: for (i = 0; i < nk; i++) {
73: PetscCall(MatDenseGetColumn(B, i, &valsB));
74: PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, valsB, &x));
75: PetscCall(MatCreateVecs(A, NULL, &Ax));
76: PetscCall(MatMult(A, x, Ax));
77: PetscCall(VecDestroy(&x));
78: PetscCall(MatDenseGetColumn(C, i, &valsC));
79: PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, m, M, valsC, &y));
80: PetscCall(VecAXPY(y, -1.0, Ax));
81: PetscCall(VecDestroy(&Ax));
82: PetscCall(VecDestroy(&y));
83: PetscCall(MatDenseRestoreColumn(C, &valsC));
84: PetscCall(MatDenseRestoreColumn(B, &valsB));
85: }
86: PetscCall(MatNorm(C, NORM_INFINITY, &norm));
87: PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatMatMult(): %g", (double)norm);
88: PetscCall(MatDestroy(&C));
89: PetscCall(MatDestroy(&B));
91: for (i = 0; i < (nr * nc); i++) PetscCall(MatDestroy(&subs[i]));
92: PetscCall(MatDestroy(&A));
93: PetscCall(MatDestroy(&Atranspose));
94: PetscCall(VecDestroy(&ATy));
95: PetscCall(PetscRandomDestroy(&rctx));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: PetscErrorCode TestReuseMatrix(void)
100: {
101: const PetscInt n = 2;
102: Mat A;
103: Mat subs[2 * 2], **block;
104: PetscInt i, j;
105: PetscRandom rctx;
106: PetscScalar zero = 0.0, one = 1.0;
108: PetscFunctionBegin;
109: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
110: PetscCall(PetscRandomSetInterval(rctx, zero, one));
111: PetscCall(PetscRandomSetFromOptions(rctx));
112: for (i = 0; i < (n * n); i++) PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD, n, n, NULL, &subs[i]));
113: PetscCall(MatCreateNest(PETSC_COMM_WORLD, n, NULL, n, NULL, subs, &A));
114: PetscCall(MatSetRandom(A, rctx));
116: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
117: PetscCall(MatNestGetSubMats(A, NULL, NULL, &block));
118: for (i = 0; i < n; i++) {
119: for (j = 0; j < n; j++) PetscCall(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD));
120: }
121: PetscCall(MatTranspose(A, MAT_INPLACE_MATRIX, &A));
122: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
123: PetscCall(MatNestGetSubMats(A, NULL, NULL, &block));
124: for (i = 0; i < n; i++) {
125: for (j = 0; j < n; j++) PetscCall(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD));
126: }
128: for (i = 0; i < (n * n); i++) PetscCall(MatDestroy(&subs[i]));
129: PetscCall(MatDestroy(&A));
130: PetscCall(PetscRandomDestroy(&rctx));
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: int main(int argc, char **args)
135: {
136: PetscFunctionBeginUser;
137: PetscCall(PetscInitialize(&argc, &args, NULL, help));
138: PetscCall(TestInitialMatrix());
139: PetscCall(TestReuseMatrix());
140: PetscCall(PetscFinalize());
141: return 0;
142: }
144: /*TEST
146: test:
147: args: -malloc_dump
149: TEST*/