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: PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
21: /* 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 */
22: PetscRandomSetInterval(rctx, zero, one);
23: PetscRandomSetFromOptions(rctx);
24: for (i = 0; i < (nr * nc); i++) MatCreateSeqDense(PETSC_COMM_WORLD, arow[i], acol[i], NULL, &subs[i]);
25: MatCreateNest(PETSC_COMM_WORLD, nr, NULL, nc, NULL, subs, &A);
26: MatCreateVecs(A, &x, NULL);
27: MatCreateVecs(A, NULL, &y);
28: VecDuplicate(x, &ATy);
29: VecDuplicate(y, &Ax);
30: MatSetRandom(A, rctx);
31: MatTranspose(A, MAT_INITIAL_MATRIX, &Atranspose);
33: MatView(A, PETSC_VIEWER_STDOUT_WORLD);
34: MatNestGetSubMats(A, NULL, NULL, &block);
35: for (i = 0; i < nr; i++) {
36: for (j = 0; j < nc; j++) MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
37: }
39: MatView(Atranspose, PETSC_VIEWER_STDOUT_WORLD);
40: MatNestGetSubMats(Atranspose, NULL, NULL, &block);
41: for (i = 0; i < nc; i++) {
42: for (j = 0; j < nr; j++) MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
43: }
45: /* Check <Ax, y> = <x, A^Ty> */
46: for (i = 0; i < 10; i++) {
47: VecSetRandom(x, rctx);
48: VecSetRandom(y, rctx);
50: MatMult(A, x, Ax);
51: VecDot(Ax, y, &dot1);
52: MatMult(Atranspose, y, ATy);
53: VecDot(ATy, x, &dot2);
55: PetscPrintf(PETSC_COMM_WORLD, "<Ax, y> = %g\n", (double)PetscRealPart(dot1));
56: PetscPrintf(PETSC_COMM_WORLD, "<x, A^Ty> = %g\n", (double)PetscRealPart(dot2));
57: }
58: VecDestroy(&x);
59: VecDestroy(&y);
60: VecDestroy(&Ax);
62: MatCreateSeqDense(PETSC_COMM_WORLD, acol[0] + acol[nr] + acol[2 * nr], nk, NULL, &B);
63: MatSetRandom(B, rctx);
64: MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C);
65: MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C);
66: MatMatMultEqual(A, B, C, 10, &equal);
69: MatGetSize(A, &M, &N);
70: MatGetLocalSize(A, &m, &n);
71: for (i = 0; i < nk; i++) {
72: MatDenseGetColumn(B, i, &valsB);
73: VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, valsB, &x);
74: MatCreateVecs(A, NULL, &Ax);
75: MatMult(A, x, Ax);
76: VecDestroy(&x);
77: MatDenseGetColumn(C, i, &valsC);
78: VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, m, M, valsC, &y);
79: VecAXPY(y, -1.0, Ax);
80: VecDestroy(&Ax);
81: VecDestroy(&y);
82: MatDenseRestoreColumn(C, &valsC);
83: MatDenseRestoreColumn(B, &valsB);
84: }
85: MatNorm(C, NORM_INFINITY, &norm);
87: MatDestroy(&C);
88: MatDestroy(&B);
90: for (i = 0; i < (nr * nc); i++) MatDestroy(&subs[i]);
91: MatDestroy(&A);
92: MatDestroy(&Atranspose);
93: VecDestroy(&ATy);
94: PetscRandomDestroy(&rctx);
95: return 0;
96: }
98: PetscErrorCode TestReuseMatrix(void)
99: {
100: const PetscInt n = 2;
101: Mat A;
102: Mat subs[2 * 2], **block;
103: PetscInt i, j;
104: PetscRandom rctx;
105: PetscScalar zero = 0.0, one = 1.0;
107: PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
108: PetscRandomSetInterval(rctx, zero, one);
109: PetscRandomSetFromOptions(rctx);
110: for (i = 0; i < (n * n); i++) MatCreateSeqDense(PETSC_COMM_WORLD, n, n, NULL, &subs[i]);
111: MatCreateNest(PETSC_COMM_WORLD, n, NULL, n, NULL, subs, &A);
112: MatSetRandom(A, rctx);
114: MatView(A, PETSC_VIEWER_STDOUT_WORLD);
115: MatNestGetSubMats(A, NULL, NULL, &block);
116: for (i = 0; i < n; i++) {
117: for (j = 0; j < n; j++) MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
118: }
119: MatTranspose(A, MAT_INPLACE_MATRIX, &A);
120: MatView(A, PETSC_VIEWER_STDOUT_WORLD);
121: MatNestGetSubMats(A, NULL, NULL, &block);
122: for (i = 0; i < n; i++) {
123: for (j = 0; j < n; j++) MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
124: }
126: for (i = 0; i < (n * n); i++) MatDestroy(&subs[i]);
127: MatDestroy(&A);
128: PetscRandomDestroy(&rctx);
129: return 0;
130: }
132: int main(int argc, char **args)
133: {
135: PetscInitialize(&argc, &args, (char *)0, help);
136: TestInitialMatrix();
137: TestReuseMatrix();
138: PetscFinalize();
139: return 0;
140: }
142: /*TEST
144: test:
145: args: -malloc_dump
147: TEST*/