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, void *ctx)
6: {
7: PetscInt d, j, k;
8: PetscReal diff = 0.0, *coords = (PetscReal *)(ctx);
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: return 0;
19: }
21: static PetscErrorCode GenEntriesRectangular(PetscInt sdim, PetscInt M, PetscInt N, const PetscInt *J, const PetscInt *K, PetscScalar *ptr, void *ctx)
22: {
23: PetscInt d, j, k;
24: PetscReal diff = 0.0, **coords = (PetscReal **)(ctx);
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: return 0;
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 size;
42: PetscScalar *ptr;
43: PetscReal *coords, *gcoords, *scoords, *gscoords, *(ctx[2]), norm, epsilon;
44: MatHtoolKernel kernel = GenEntries;
45: PetscBool flg, sym = PETSC_FALSE;
46: PetscRandom rdm;
47: IS iss, ist, is[2];
48: Vec right, left, perm;
51: PetscInitialize(&argc, &argv, (char *)NULL, help);
52: PetscOptionsGetInt(NULL, NULL, "-m_local", &m, NULL);
53: PetscOptionsGetInt(NULL, NULL, "-n_local", &n, NULL);
54: PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL);
55: PetscOptionsGetInt(NULL, NULL, "-K", &K, NULL);
56: PetscOptionsGetBool(NULL, NULL, "-symmetric", &sym, NULL);
57: PetscOptionsGetReal(NULL, NULL, "-mat_htool_epsilon", &epsilon, NULL);
58: MPI_Comm_size(PETSC_COMM_WORLD, &size);
59: M = size * m;
60: PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
61: PetscMalloc1(m * dim, &coords);
62: PetscRandomCreate(PETSC_COMM_WORLD, &rdm);
63: PetscRandomGetValuesReal(rdm, m * dim, coords);
64: PetscCalloc1(M * dim, &gcoords);
65: MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, NULL, &B);
66: MatSetRandom(B, rdm);
67: MatGetOwnershipRange(B, &begin, NULL);
68: PetscArraycpy(gcoords + begin * dim, coords, m * dim);
69: MPIU_Allreduce(MPI_IN_PLACE, gcoords, M * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD);
70: MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &A);
71: MatSetOption(A, MAT_SYMMETRIC, sym);
72: MatSetFromOptions(A);
73: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
74: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
75: MatViewFromOptions(A, NULL, "-A_view");
76: MatGetOwnershipIS(A, is, NULL);
77: ISDuplicate(is[0], is + 1);
78: MatIncreaseOverlap(A, 1, is, 2);
79: MatSetBlockSize(A, 2);
80: MatIncreaseOverlap(A, 1, is + 1, 1);
81: ISGetBlockSize(is[1], &bs);
83: MatSetBlockSize(A, 1);
84: ISEqual(is[0], is[1], &flg);
86: ISDestroy(is);
87: ISDestroy(is + 1);
88: MatCreateVecs(A, &right, &left);
89: VecSetRandom(right, rdm);
90: MatMult(A, right, left);
91: MatHtoolGetPermutationSource(A, &iss);
92: MatHtoolGetPermutationTarget(A, &ist);
93: VecDuplicate(left, &perm);
94: VecCopy(left, perm);
95: VecPermute(perm, ist, PETSC_FALSE);
96: VecPermute(right, iss, PETSC_FALSE);
97: MatHtoolUsePermutation(A, PETSC_FALSE);
98: MatMult(A, right, left);
99: VecAXPY(left, -1.0, perm);
100: VecNorm(left, NORM_INFINITY, &norm);
102: MatHtoolUsePermutation(A, PETSC_TRUE);
103: VecDestroy(&perm);
104: VecDestroy(&left);
105: VecDestroy(&right);
106: ISDestroy(&ist);
107: ISDestroy(&iss);
108: if (PetscAbsReal(epsilon) >= PETSC_SMALL) { /* when there is compression, it is more difficult to check against MATDENSE, so just compare symmetric and nonsymmetric assemblies */
109: PetscReal relative;
110: MatDestroy(&B);
111: MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &B);
112: MatSetOption(B, MAT_SYMMETRIC, (PetscBool)!sym);
113: MatSetFromOptions(B);
114: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
115: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
116: MatViewFromOptions(B, NULL, "-B_view");
117: MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &P);
118: MatNorm(P, NORM_FROBENIUS, &relative);
119: MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &R);
120: MatAXPY(R, -1.0, P, SAME_NONZERO_PATTERN);
121: MatNorm(R, NORM_INFINITY, &norm);
123: MatDestroy(&B);
124: MatDestroy(&R);
125: MatDestroy(&P);
126: } else {
127: MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &D);
128: MatViewFromOptions(D, NULL, "-D_view");
129: MatMultEqual(A, D, 10, &flg);
131: MatMultTransposeEqual(A, D, 10, &flg);
133: MatMultAddEqual(A, D, 10, &flg);
135: MatGetOwnershipRange(B, &begin, NULL);
136: MatDenseGetArrayWrite(D, &ptr);
137: for (PetscInt i = begin; i < m + begin; ++i)
138: for (PetscInt j = 0; j < M; ++j) GenEntries(dim, 1, 1, &i, &j, ptr + i - begin + j * m, gcoords);
139: MatDenseRestoreArrayWrite(D, &ptr);
140: MatMultEqual(A, D, 10, &flg);
142: MatTranspose(D, MAT_INPLACE_MATRIX, &D);
143: MatTranspose(A, MAT_INITIAL_MATRIX, &AT);
144: MatMultEqual(AT, D, 10, &flg);
146: MatTranspose(A, MAT_REUSE_MATRIX, &AT);
147: MatMultEqual(AT, D, 10, &flg);
149: MatAXPY(D, -1.0, AT, SAME_NONZERO_PATTERN);
150: MatNorm(D, NORM_INFINITY, &norm);
152: MatDestroy(&AT);
153: MatDestroy(&D);
154: MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &P);
155: MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
156: MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
157: MatMatMultEqual(A, B, P, 10, &flg);
159: MatTransposeMatMultEqual(A, B, P, 10, &flg);
161: MatDestroy(&B);
162: MatDestroy(&P);
163: if (n) {
164: PetscMalloc1(n * dim, &scoords);
165: PetscRandomGetValuesReal(rdm, n * dim, scoords);
166: N = n;
167: MPIU_Allreduce(MPI_IN_PLACE, &N, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD);
168: PetscCalloc1(N * dim, &gscoords);
169: MPI_Exscan(&n, &begin, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD);
170: PetscArraycpy(gscoords + begin * dim, scoords, n * dim);
171: MPIU_Allreduce(MPI_IN_PLACE, gscoords, N * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD);
172: kernel = GenEntriesRectangular;
173: ctx[0] = gcoords;
174: ctx[1] = gscoords;
175: MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, n, M, N, dim, coords, scoords, kernel, ctx, &R);
176: MatSetFromOptions(R);
177: MatAssemblyBegin(R, MAT_FINAL_ASSEMBLY);
178: MatAssemblyEnd(R, MAT_FINAL_ASSEMBLY);
179: MatViewFromOptions(R, NULL, "-R_view");
180: MatConvert(R, MATDENSE, MAT_INITIAL_MATRIX, &D);
181: MatViewFromOptions(D, NULL, "-D_view");
182: MatMultEqual(R, D, 10, &flg);
184: MatTranspose(D, MAT_INPLACE_MATRIX, &D);
185: MatTranspose(R, MAT_INITIAL_MATRIX, &RT);
186: MatMultEqual(RT, D, 10, &flg);
188: MatTranspose(R, MAT_REUSE_MATRIX, &RT);
189: MatMultEqual(RT, D, 10, &flg);
191: MatDestroy(&RT);
192: MatDestroy(&D);
193: MatCreateDense(PETSC_COMM_WORLD, n, PETSC_DECIDE, PETSC_DETERMINE, K, NULL, &B);
194: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
195: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
196: MatSetRandom(B, rdm);
197: MatMatMult(R, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &P);
198: MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
199: MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
200: MatMatMultEqual(R, B, P, 10, &flg);
202: MatDestroy(&B);
203: MatDestroy(&P);
204: MatCreateVecs(R, &right, &left);
205: VecSetRandom(right, rdm);
206: MatMult(R, right, left);
207: MatHtoolGetPermutationSource(R, &iss);
208: MatHtoolGetPermutationTarget(R, &ist);
209: VecDuplicate(left, &perm);
210: VecCopy(left, perm);
211: VecPermute(perm, ist, PETSC_FALSE);
212: VecPermute(right, iss, PETSC_FALSE);
213: MatHtoolUsePermutation(R, PETSC_FALSE);
214: MatMult(R, right, left);
215: VecAXPY(left, -1.0, perm);
216: VecNorm(left, NORM_INFINITY, &norm);
218: MatHtoolUsePermutation(R, PETSC_TRUE);
219: VecDestroy(&perm);
220: VecDestroy(&left);
221: VecDestroy(&right);
222: ISDestroy(&ist);
223: ISDestroy(&iss);
224: MatDestroy(&R);
225: PetscFree(gscoords);
226: PetscFree(scoords);
227: }
228: }
229: PetscRandomDestroy(&rdm);
230: MatDestroy(&A);
231: PetscFree(gcoords);
232: PetscFree(coords);
233: PetscFinalize();
234: return 0;
235: }
237: /*TEST
239: build:
240: requires: htool
242: test:
243: requires: htool
244: suffix: 1
245: nsize: 4
246: args: -m_local 80 -n_local 25 -mat_htool_epsilon 1.0e-11 -symmetric {{false true}shared output}
247: output_file: output/ex101.out
249: test:
250: requires: htool
251: suffix: 2
252: nsize: 4
253: 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}
254: output_file: output/ex101.out
256: TEST*/