1 static char help[] = "Tests the use of MatTranspose_Nest and MatMatMult_Nest_Dense\n";
2
3 #include <petscmat.h>
4
TestInitialMatrix(void)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;
19
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));
33
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 }
39
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 }
45
46 /* Check <Ax, y> = <x, A^Ty> */
47 for (i = 0; i < 10; i++) {
48 PetscCall(VecSetRandom(x, rctx));
49 PetscCall(VecSetRandom(y, rctx));
50
51 PetscCall(MatMult(A, x, Ax));
52 PetscCall(VecDot(Ax, y, &dot1));
53 PetscCall(MatMult(Atranspose, y, ATy));
54 PetscCall(VecDot(ATy, x, &dot2));
55
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));
62
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");
69
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));
90
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 }
98
TestReuseMatrix(void)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;
107
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));
115
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 }
127
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 }
133
main(int argc,char ** args)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 }
143
144 /*TEST
145
146 test:
147 args: -malloc_dump
148
149 TEST*/
150