1 static char help[] = "Tests MATHTOOL\n\n";
2
3 #include <petscmat.h>
4
GenEntries(PetscInt sdim,PetscInt M,PetscInt N,const PetscInt * J,const PetscInt * K,PetscScalar * ptr,PetscCtx ctx)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);
9
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 }
20
GenEntriesRectangular(PetscInt sdim,PetscInt M,PetscInt N,const PetscInt * J,const PetscInt * K,PetscScalar * ptr,PetscCtx ctx)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);
25
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 }
36
main(int argc,char ** argv)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 size;
42 PetscScalar *ptr;
43 PetscReal *coords, *gcoords, *scoords, *gscoords, *ctx[2], norm, epsilon = PetscSqrtReal(PETSC_SMALL);
44 MatHtoolKernelFn *kernel = GenEntries;
45 PetscBool flg, sym = PETSC_FALSE, recompression = PETSC_FALSE;
46 PetscRandom rdm;
47 IS iss, ist, is[2];
48 Vec right, left, perm;
49
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 PetscCall(PetscOptionsGetReal(NULL, NULL, "-mat_htool_epsilon", &epsilon, NULL));
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(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
76 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
77 PetscCall(MatViewFromOptions(A, NULL, "-A_view"));
78 PetscCall(MatGetOwnershipIS(A, is, NULL));
79 PetscCall(ISDuplicate(is[0], is + 1));
80 PetscCall(MatIncreaseOverlap(A, 1, is, 2));
81 PetscCall(MatSetBlockSize(A, 2));
82 PetscCall(MatIncreaseOverlap(A, 1, is + 1, 1));
83 PetscCall(ISGetBlockSize(is[1], &bs));
84 PetscCheck(bs == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect block size %" PetscInt_FMT " != 2", bs);
85 PetscCall(MatSetBlockSize(A, 1));
86 PetscCall(ISEqual(is[0], is[1], &flg));
87 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unequal index sets");
88 PetscCall(ISDestroy(is));
89 PetscCall(ISDestroy(is + 1));
90 PetscCall(MatCreateVecs(A, &right, &left));
91 PetscCall(VecSetRandom(right, rdm));
92 PetscCall(MatMult(A, right, left));
93 PetscCall(MatHtoolGetPermutationSource(A, &iss));
94 PetscCall(MatHtoolGetPermutationTarget(A, &ist));
95 PetscCall(VecDuplicate(left, &perm));
96 PetscCall(VecCopy(left, perm));
97 PetscCall(VecPermute(perm, ist, PETSC_FALSE));
98 PetscCall(VecPermute(right, iss, PETSC_FALSE));
99 PetscCall(MatHtoolUsePermutation(A, PETSC_FALSE));
100 PetscCall(MatMult(A, right, left));
101 PetscCall(VecAXPY(left, -1.0, perm));
102 PetscCall(VecNorm(left, NORM_INFINITY, &norm));
103 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);
104 PetscCall(MatHtoolUsePermutation(A, PETSC_TRUE));
105 PetscCall(VecDestroy(&perm));
106 PetscCall(VecDestroy(&left));
107 PetscCall(VecDestroy(&right));
108 PetscCall(ISDestroy(&ist));
109 PetscCall(ISDestroy(&iss));
110 if (PetscAbsReal(epsilon) >= PETSC_SMALL) { /* when there is compression, it is more difficult to check against MATDENSE, so just compare symmetric and nonsymmetric assemblies */
111 PetscReal relative;
112 PetscCall(MatDestroy(&B));
113 PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &B));
114 PetscCall(MatSetOption(B, MAT_SYMMETRIC, (PetscBool)!sym));
115 PetscCall(MatHtoolUseRecompression(B, recompression));
116 PetscCall(MatSetFromOptions(B));
117 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
118 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
119 PetscCall(MatViewFromOptions(B, NULL, "-B_view"));
120 PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &P));
121 PetscCall(MatNorm(P, NORM_FROBENIUS, &relative));
122 PetscCall(MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &R));
123 PetscCall(MatAXPY(R, -1.0, P, SAME_NONZERO_PATTERN));
124 PetscCall(MatNorm(R, NORM_INFINITY, &norm));
125 PetscCheck(PetscAbsReal(norm / relative) <= epsilon, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "||A(!symmetric)-A(symmetric)|| = %g (> %g)", (double)PetscAbsReal(norm / relative), (double)epsilon);
126 PetscCall(MatDestroy(&B));
127 PetscCall(MatDestroy(&R));
128 PetscCall(MatDestroy(&P));
129 } else {
130 PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &D));
131 PetscCall(MatViewFromOptions(D, NULL, "-D_view"));
132 PetscCall(MatMultEqual(A, D, 10, &flg));
133 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Ax != Dx");
134 PetscCall(MatMultTransposeEqual(A, D, 10, &flg));
135 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A^Tx != D^Tx");
136 PetscCall(MatMultAddEqual(A, D, 10, &flg));
137 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "y+Ax != y+Dx");
138 PetscCall(MatGetOwnershipRange(B, &begin, NULL));
139 PetscCall(MatDenseGetArrayWrite(D, &ptr));
140 for (PetscInt i = begin; i < m + begin; ++i)
141 for (PetscInt j = 0; j < M; ++j) PetscCall(GenEntries(dim, 1, 1, &i, &j, ptr + i - begin + j * m, gcoords));
142 PetscCall(MatDenseRestoreArrayWrite(D, &ptr));
143 PetscCall(MatMultEqual(A, D, 10, &flg));
144 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Ax != Dx");
145 PetscCall(MatTranspose(D, MAT_INPLACE_MATRIX, &D));
146 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
147 PetscCall(MatMultEqual(AT, D, 10, &flg));
148 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A^Tx != D^Tx");
149 PetscCall(MatTranspose(A, MAT_REUSE_MATRIX, &AT));
150 PetscCall(MatMultEqual(AT, D, 10, &flg));
151 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A^Tx != D^Tx");
152 PetscCall(MatAXPY(D, -1.0, AT, SAME_NONZERO_PATTERN));
153 PetscCall(MatNorm(D, NORM_INFINITY, &norm));
154 PetscCheck(PetscAbsReal(norm) <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "||A-D|| = %g (> %g)", (double)norm, (double)PETSC_SMALL);
155 PetscCall(MatDestroy(&AT));
156 PetscCall(MatDestroy(&D));
157 PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &P));
158 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
159 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
160 PetscCall(MatMatMultEqual(A, B, P, 10, &flg));
161 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ABx != Px");
162 PetscCall(MatTransposeMatMultEqual(A, B, P, 10, &flg));
163 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A^TBx != P^Tx");
164 PetscCall(MatDestroy(&B));
165 PetscCall(MatDestroy(&P));
166 if (n) {
167 PetscCall(PetscMalloc1(n * dim, &scoords));
168 PetscCall(PetscRandomGetValuesReal(rdm, n * dim, scoords));
169 N = n;
170 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &N, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
171 PetscCall(PetscCalloc1(N * dim, &gscoords));
172 PetscCallMPI(MPI_Exscan(&n, &begin, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
173 PetscCall(PetscArraycpy(gscoords + begin * dim, scoords, n * dim));
174 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, gscoords, N * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
175 kernel = GenEntriesRectangular;
176 ctx[0] = gcoords;
177 ctx[1] = gscoords;
178 PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, n, M, N, dim, coords, scoords, kernel, ctx, &R));
179 PetscCall(MatSetFromOptions(R));
180 PetscCall(MatHtoolUseRecompression(R, recompression));
181 PetscCall(MatAssemblyBegin(R, MAT_FINAL_ASSEMBLY));
182 PetscCall(MatAssemblyEnd(R, MAT_FINAL_ASSEMBLY));
183 PetscCall(MatViewFromOptions(R, NULL, "-R_view"));
184 PetscCall(MatConvert(R, MATDENSE, MAT_INITIAL_MATRIX, &D));
185 PetscCall(MatViewFromOptions(D, NULL, "-D_view"));
186 PetscCall(MatMultEqual(R, D, 10, &flg));
187 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Rx != Dx");
188 PetscCall(MatTranspose(D, MAT_INPLACE_MATRIX, &D));
189 PetscCall(MatTranspose(R, MAT_INITIAL_MATRIX, &RT));
190 PetscCall(MatMultEqual(RT, D, 10, &flg));
191 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "R^Tx != D^Tx");
192 PetscCall(MatTranspose(R, MAT_REUSE_MATRIX, &RT));
193 PetscCall(MatMultEqual(RT, D, 10, &flg));
194 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "R^Tx != D^Tx");
195 PetscCall(MatDestroy(&RT));
196 PetscCall(MatDestroy(&D));
197 PetscCall(MatCreateDense(PETSC_COMM_WORLD, n, PETSC_DECIDE, PETSC_DETERMINE, K, NULL, &B));
198 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
199 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
200 PetscCall(MatSetRandom(B, rdm));
201 PetscCall(MatMatMult(R, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &P));
202 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
203 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
204 PetscCall(MatMatMultEqual(R, B, P, 10, &flg));
205 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "RBx != Px");
206 PetscCall(MatDestroy(&B));
207 PetscCall(MatDestroy(&P));
208 PetscCall(MatCreateVecs(R, &right, &left));
209 PetscCall(VecSetRandom(right, rdm));
210 PetscCall(MatMult(R, right, left));
211 PetscCall(MatHtoolGetPermutationSource(R, &iss));
212 PetscCall(MatHtoolGetPermutationTarget(R, &ist));
213 PetscCall(VecDuplicate(left, &perm));
214 PetscCall(VecCopy(left, perm));
215 PetscCall(VecPermute(perm, ist, PETSC_FALSE));
216 PetscCall(VecPermute(right, iss, PETSC_FALSE));
217 PetscCall(MatHtoolUsePermutation(R, PETSC_FALSE));
218 PetscCall(MatMult(R, right, left));
219 PetscCall(VecAXPY(left, -1.0, perm));
220 PetscCall(VecNorm(left, NORM_INFINITY, &norm));
221 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);
222 PetscCall(MatHtoolUsePermutation(R, PETSC_TRUE));
223 PetscCall(VecDestroy(&perm));
224 PetscCall(VecDestroy(&left));
225 PetscCall(VecDestroy(&right));
226 PetscCall(ISDestroy(&ist));
227 PetscCall(ISDestroy(&iss));
228 PetscCall(MatDestroy(&R));
229 PetscCall(PetscFree(gscoords));
230 PetscCall(PetscFree(scoords));
231 }
232 }
233 PetscCall(PetscRandomDestroy(&rdm));
234 PetscCall(MatDestroy(&A));
235 PetscCall(PetscFree(gcoords));
236 PetscCall(PetscFree(coords));
237 PetscCall(PetscFinalize());
238 return 0;
239 }
240
241 /*TEST
242
243 build:
244 requires: htool
245
246 test:
247 requires: htool
248 suffix: 1
249 nsize: 4
250 args: -m_local 80 -n_local 25 -mat_htool_epsilon 1.0e-11 -symmetric {{false true}shared output} -recompression {{false true}shared output}
251 output_file: output/empty.out
252
253 test:
254 requires: htool
255 suffix: 2
256 nsize: 4
257 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}
258 output_file: output/empty.out
259
260 TEST*/
261