Actual source code: ex241.c
1: static char help[] = "Tests MATHTOOL\n\n";
3: #include <petscmat.h>
5: static PetscErrorCode GenEntries(PetscInt sdim, PetscInt M, PetscInt N, const PetscInt *J, const PetscInt *K, PetscScalar *ptr, PetscCtx ctx)
6: {
7: PetscInt d, j, k;
8: PetscReal diff = 0.0, *coords = (PetscReal *)(ctx);
10: PetscFunctionBeginUser;
11: for (j = 0; j < M; j++) {
12: for (k = 0; k < N; k++) {
13: diff = 0.0;
14: for (d = 0; d < sdim; d++) diff += (coords[J[j] * sdim + d] - coords[K[k] * sdim + d]) * (coords[J[j] * sdim + d] - coords[K[k] * sdim + d]);
15: ptr[j + M * k] = 1.0 / (1.0e-2 + PetscSqrtReal(diff));
16: }
17: }
18: PetscFunctionReturn(PETSC_SUCCESS);
19: }
21: static PetscErrorCode GenEntriesRectangular(PetscInt sdim, PetscInt M, PetscInt N, const PetscInt *J, const PetscInt *K, PetscScalar *ptr, PetscCtx ctx)
22: {
23: PetscInt d, j, k;
24: PetscReal diff = 0.0, **coords = (PetscReal **)(ctx);
26: PetscFunctionBeginUser;
27: for (j = 0; j < M; j++) {
28: for (k = 0; k < N; k++) {
29: diff = 0.0;
30: for (d = 0; d < sdim; d++) diff += (coords[0][J[j] * sdim + d] - coords[1][K[k] * sdim + d]) * (coords[0][J[j] * sdim + d] - coords[1][K[k] * sdim + d]);
31: ptr[j + M * k] = 1.0 / (1.0e-2 + PetscSqrtReal(diff));
32: }
33: }
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: int main(int argc, char **argv)
38: {
39: Mat A, AT, D, B, P, R, RT;
40: PetscInt m = 100, dim = 3, M, K = 10, begin, n = 0, N, bs;
41: PetscMPIInt rank, size;
42: PetscScalar *ptr;
43: PetscReal *coords, *gcoords, *scoords, *gscoords, *ctx[2], norm, epsilon;
44: MatHtoolKernelFn *kernel = GenEntries;
45: PetscBool flg, sym = PETSC_FALSE, recompression = PETSC_FALSE, consistent;
46: PetscRandom rdm;
47: IS iss, ist, is[2];
48: Vec right, left, perm;
50: PetscFunctionBeginUser;
51: PetscCall(PetscInitialize(&argc, &argv, (char *)NULL, help));
52: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m_local", &m, NULL));
53: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n_local", &n, NULL));
54: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL));
55: PetscCall(PetscOptionsGetInt(NULL, NULL, "-K", &K, NULL));
56: PetscCall(PetscOptionsGetBool(NULL, NULL, "-symmetric", &sym, NULL));
57: PetscCall(PetscOptionsGetBool(NULL, NULL, "-recompression", &recompression, NULL));
58: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
59: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
60: M = size * m;
61: PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
62: PetscCall(PetscMalloc1(m * dim, &coords));
63: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
64: PetscCall(PetscRandomGetValuesReal(rdm, m * dim, coords));
65: PetscCall(PetscCalloc1(M * dim, &gcoords));
66: PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, NULL, &B));
67: PetscCall(MatSetRandom(B, rdm));
68: PetscCall(MatGetOwnershipRange(B, &begin, NULL));
69: PetscCall(PetscArraycpy(gcoords + begin * dim, coords, m * dim));
70: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, gcoords, M * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
71: PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &A));
72: PetscCall(MatSetOption(A, MAT_SYMMETRIC, sym));
73: PetscCall(MatHtoolUseRecompression(A, recompression));
74: PetscCall(MatSetFromOptions(A));
75: PetscCall(MatHtoolGetEpsilon(A, &epsilon));
76: PetscCall(MatHtoolGetBlockTreeConsistency(A, &consistent));
77: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
78: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
79: PetscCall(MatViewFromOptions(A, NULL, "-A_view"));
80: PetscCall(MatGetOwnershipIS(A, is, NULL));
81: PetscCall(ISDuplicate(is[0], is + 1));
82: PetscCall(MatIncreaseOverlap(A, 1, is, 2));
83: PetscCall(MatSetBlockSize(A, 2));
84: PetscCall(MatIncreaseOverlap(A, 1, is + 1, 1));
85: PetscCall(ISGetBlockSize(is[1], &bs));
86: PetscCheck(bs == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect block size %" PetscInt_FMT " != 2", bs);
87: PetscCall(MatSetBlockSize(A, 1));
88: PetscCall(ISEqual(is[0], is[1], &flg));
89: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unequal index sets");
90: PetscCall(ISDestroy(is));
91: PetscCall(ISDestroy(is + 1));
92: PetscCall(MatCreateVecs(A, &right, &left));
93: PetscCall(VecSetRandom(right, rdm));
94: PetscCall(MatMult(A, right, left));
95: PetscCall(MatHtoolGetPermutationSource(A, &iss));
96: PetscCall(MatHtoolGetPermutationTarget(A, &ist));
97: PetscCall(VecDuplicate(left, &perm));
98: PetscCall(VecCopy(left, perm));
99: PetscCall(VecPermute(perm, ist, PETSC_FALSE));
100: PetscCall(VecPermute(right, iss, PETSC_FALSE));
101: PetscCall(MatHtoolUsePermutation(A, PETSC_FALSE));
102: PetscCall(MatMult(A, right, left));
103: PetscCall(VecAXPY(left, -1.0, perm));
104: PetscCall(VecNorm(left, NORM_INFINITY, &norm));
105: PetscCheck(PetscAbsReal(norm) <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "||y(with permutation)-y(without permutation)|| = %g (> %g)", (double)PetscAbsReal(norm), (double)PETSC_SMALL);
106: PetscCall(MatHtoolUsePermutation(A, PETSC_TRUE));
107: PetscCall(VecDestroy(&perm));
108: PetscCall(VecDestroy(&left));
109: PetscCall(VecDestroy(&right));
110: PetscCall(ISDestroy(&ist));
111: PetscCall(ISDestroy(&iss));
112: if (PetscAbsReal(epsilon) >= PETSC_SMALL) { /* when there is compression, it is more difficult to check against MATDENSE, so just compare symmetric and nonsymmetric assemblies */
113: PetscReal relative;
114: PetscCall(MatDestroy(&B));
115: PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &B));
116: PetscCall(MatSetOption(B, MAT_SYMMETRIC, (PetscBool)!sym));
117: PetscCall(MatHtoolUseRecompression(B, recompression));
118: PetscCall(MatSetFromOptions(B));
119: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
120: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
121: PetscCall(MatViewFromOptions(B, NULL, "-B_view"));
122: PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &P));
123: PetscCall(MatNorm(P, NORM_FROBENIUS, &relative));
124: PetscCall(MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &R));
125: PetscCall(MatAXPY(R, -1.0, P, SAME_NONZERO_PATTERN));
126: PetscCall(MatNorm(R, NORM_INFINITY, &norm));
127: PetscCheck(PetscAbsReal(norm / relative) <= epsilon, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "||A(!symmetric)-A(symmetric)|| = %g (> %g)", (double)PetscAbsReal(norm / relative), (double)epsilon);
128: PetscCall(MatDestroy(&B));
129: PetscCall(MatDestroy(&R));
130: PetscCall(MatDestroy(&P));
131: } else {
132: PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &D));
133: PetscCall(MatViewFromOptions(D, NULL, "-D_view"));
134: PetscCall(MatMultEqual(A, D, 10, &flg));
135: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Ax != Dx");
136: PetscCall(MatMultTransposeEqual(A, D, 10, &flg));
137: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A^Tx != D^Tx");
138: PetscCall(MatMultAddEqual(A, D, 10, &flg));
139: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "y+Ax != y+Dx");
140: PetscCall(MatGetOwnershipRange(B, &begin, NULL));
141: PetscCall(MatDenseGetArrayWrite(D, &ptr));
142: for (PetscInt i = begin; i < m + begin; ++i)
143: for (PetscInt j = 0; j < M; ++j) PetscCall(GenEntries(dim, 1, 1, &i, &j, ptr + i - begin + j * m, gcoords));
144: PetscCall(MatDenseRestoreArrayWrite(D, &ptr));
145: PetscCall(MatMultEqual(A, D, 10, &flg));
146: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Ax != Dx");
147: PetscCall(MatTranspose(D, MAT_INPLACE_MATRIX, &D));
148: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
149: PetscCall(MatMultEqual(AT, D, 10, &flg));
150: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A^Tx != D^Tx");
151: PetscCall(MatTranspose(A, MAT_REUSE_MATRIX, &AT));
152: PetscCall(MatMultEqual(AT, D, 10, &flg));
153: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A^Tx != D^Tx");
154: PetscCall(MatAXPY(D, -1.0, AT, SAME_NONZERO_PATTERN));
155: PetscCall(MatNorm(D, NORM_INFINITY, &norm));
156: PetscCheck(PetscAbsReal(norm) <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "||A-D|| = %g (> %g)", (double)norm, (double)PETSC_SMALL);
157: PetscCall(MatDestroy(&AT));
158: PetscCall(MatDestroy(&D));
159: PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &P));
160: PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
161: PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
162: PetscCall(MatMatMultEqual(A, B, P, 10, &flg));
163: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ABx != Px");
164: PetscCall(MatTransposeMatMultEqual(A, B, P, 10, &flg));
165: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A^TBx != P^Tx");
166: PetscCall(MatDestroy(&B));
167: PetscCall(MatDestroy(&P));
168: if (n) {
169: PetscCall(PetscMalloc1(n * dim, &scoords));
170: PetscCall(PetscRandomGetValuesReal(rdm, n * dim, scoords));
171: N = n;
172: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &N, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
173: PetscCall(PetscCalloc1(N * dim, &gscoords));
174: PetscCallMPI(MPI_Exscan(&n, &begin, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
175: PetscCall(PetscArraycpy(gscoords + begin * dim, scoords, n * dim));
176: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, gscoords, N * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
177: kernel = GenEntriesRectangular;
178: ctx[0] = gcoords;
179: ctx[1] = gscoords;
180: PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, n, M, N, dim, coords, scoords, kernel, ctx, &R));
181: PetscCall(MatSetFromOptions(R));
182: PetscCall(MatHtoolUseRecompression(R, recompression));
183: PetscCall(MatAssemblyBegin(R, MAT_FINAL_ASSEMBLY));
184: PetscCall(MatAssemblyEnd(R, MAT_FINAL_ASSEMBLY));
185: PetscCall(MatViewFromOptions(R, NULL, "-R_view"));
186: PetscCall(MatConvert(R, MATDENSE, MAT_INITIAL_MATRIX, &D));
187: PetscCall(MatViewFromOptions(D, NULL, "-D_view"));
188: PetscCall(MatMultEqual(R, D, 10, &flg));
189: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Rx != Dx");
190: PetscCall(MatTranspose(D, MAT_INPLACE_MATRIX, &D));
191: PetscCall(MatTranspose(R, MAT_INITIAL_MATRIX, &RT));
192: PetscCall(MatMultEqual(RT, D, 10, &flg));
193: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "R^Tx != D^Tx");
194: PetscCall(MatTranspose(R, MAT_REUSE_MATRIX, &RT));
195: PetscCall(MatMultEqual(RT, D, 10, &flg));
196: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "R^Tx != D^Tx");
197: PetscCall(MatDestroy(&RT));
198: PetscCall(MatDestroy(&D));
199: PetscCall(MatCreateDense(PETSC_COMM_WORLD, n, PETSC_DECIDE, PETSC_DETERMINE, K, NULL, &B));
200: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
201: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
202: PetscCall(MatSetRandom(B, rdm));
203: PetscCall(MatMatMult(R, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &P));
204: PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
205: PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
206: PetscCall(MatMatMultEqual(R, B, P, 10, &flg));
207: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "RBx != Px");
208: PetscCall(MatDestroy(&B));
209: PetscCall(MatDestroy(&P));
210: PetscCall(MatCreateVecs(R, &right, &left));
211: PetscCall(VecSetRandom(right, rdm));
212: PetscCall(MatMult(R, right, left));
213: PetscCall(MatHtoolGetPermutationSource(R, &iss));
214: PetscCall(MatHtoolGetPermutationTarget(R, &ist));
215: PetscCall(VecDuplicate(left, &perm));
216: PetscCall(VecCopy(left, perm));
217: PetscCall(VecPermute(perm, ist, PETSC_FALSE));
218: PetscCall(VecPermute(right, iss, PETSC_FALSE));
219: PetscCall(MatHtoolUsePermutation(R, PETSC_FALSE));
220: PetscCall(MatMult(R, right, left));
221: PetscCall(VecAXPY(left, -1.0, perm));
222: PetscCall(VecNorm(left, NORM_INFINITY, &norm));
223: PetscCheck(PetscAbsReal(norm) <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "||y(with permutation)-y(without permutation)|| = %g (> %g)", (double)PetscAbsReal(norm), (double)PETSC_SMALL);
224: PetscCall(MatHtoolUsePermutation(R, PETSC_TRUE));
225: PetscCall(VecDestroy(&perm));
226: PetscCall(VecDestroy(&left));
227: PetscCall(VecDestroy(&right));
228: PetscCall(ISDestroy(&ist));
229: PetscCall(ISDestroy(&iss));
230: PetscCall(MatDestroy(&R));
231: PetscCall(PetscFree(gscoords));
232: PetscCall(PetscFree(scoords));
233: }
234: }
235: /* verify that MatGetDiagonalBlock() returns the proper matrix */
236: if (consistent) {
237: Mat D, B, C;
239: PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, size, NULL, &B));
240: for (PetscMPIInt i = 0; i < size; ++i) {
241: Vec col, local;
243: PetscCall(MatDenseGetColumnVecWrite(B, i, &col));
244: PetscCall(VecCreateLocalVector(col, &local));
245: PetscCall(VecGetLocalVector(col, local));
246: if (i == rank) PetscCall(VecSetRandom(local, rdm));
247: else PetscCall(VecSet(local, 0.0));
248: PetscCall(VecRestoreLocalVector(col, local));
249: PetscCall(VecDestroy(&local));
250: PetscCall(MatDenseRestoreColumnVecWrite(B, i, &col));
251: }
252: PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
253: PetscCall(MatGetDiagonalBlock(A, &D));
254: PetscCall(PetscObjectTypeCompare((PetscObject)D, MATHTOOL, &flg));
255: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MatGetDiagonalBlock() did not return a MATHTOOL");
256: for (PetscMPIInt i = 0; i < size; ++i) {
257: Vec col[2], local[2];
259: PetscCall(MatDenseGetColumnVecRead(B, i, col));
260: PetscCall(MatDenseGetColumnVecRead(C, i, col + 1));
261: if (i == rank) {
262: Vec y;
264: PetscCall(VecCreateLocalVector(col[0], local));
265: PetscCall(VecGetLocalVectorRead(col[0], local[0]));
266: PetscCall(VecCreateLocalVector(col[1], local + 1));
267: PetscCall(VecGetLocalVectorRead(col[1], local[1]));
268: PetscCall(MatCreateVecs(D, NULL, &y));
269: PetscCall(MatMult(D, local[0], y));
270: PetscCall(VecAXPY(y, -1.0, local[1]));
271: PetscCall(VecNorm(y, NORM_INFINITY, &norm));
272: PetscCheck(norm < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Difference (= %g) greater than PETSC_SMALL (= %g)", (double)norm, (double)PETSC_SMALL);
273: PetscCall(VecDestroy(&y));
274: PetscCall(VecRestoreLocalVectorRead(col[1], local[1]));
275: PetscCall(VecDestroy(local + 1));
276: PetscCall(VecRestoreLocalVectorRead(col[0], local[0]));
277: PetscCall(VecDestroy(local));
278: }
279: PetscCall(MatDenseRestoreColumnVecRead(C, i, col + 1));
280: PetscCall(MatDenseRestoreColumnVecRead(B, i, col));
281: }
282: PetscCall(MatDestroy(&C));
283: PetscCall(MatDestroy(&B));
284: }
285: PetscCall(PetscRandomDestroy(&rdm));
286: PetscCall(MatDestroy(&A));
287: PetscCall(PetscFree(gcoords));
288: PetscCall(PetscFree(coords));
289: PetscCall(PetscFinalize());
290: return 0;
291: }
293: /*TEST
295: build:
296: requires: htool
298: test:
299: suffix: 1
300: nsize: 4
301: args: -m_local 80 -n_local 25 -mat_htool_epsilon 1.0e-11 -symmetric {{false true}shared output} -recompression {{false true}shared output}
302: output_file: output/empty.out
304: test:
305: suffix: 2
306: nsize: 4
307: args: -m_local 120 -mat_htool_epsilon 1.0e-2 -mat_htool_compressor {{sympartialACA fullACA SVD}shared output} -mat_htool_clustering {{PCARegular PCAGeometric BoundingBox1Regular BoundingBox1Geometric}shared output}
308: output_file: output/empty.out
310: test:
311: suffix: 3
312: nsize: 4
313: args: -m_local 80 -n_local 25 -mat_htool_epsilon 1.0e-11 -mat_htool_block_tree_consistency false
314: output_file: output/empty.out
316: TEST*/