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*/