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);

 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, void *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       size;
 42:   PetscScalar      *ptr;
 43:   PetscReal        *coords, *gcoords, *scoords, *gscoords, *ctx[2], norm, epsilon;
 44:   MatHtoolKernelFn *kernel = GenEntries;
 45:   PetscBool         flg, sym = PETSC_FALSE;
 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(PetscOptionsGetReal(NULL, NULL, "-mat_htool_epsilon", &epsilon, NULL));
 58:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 59:   M = size * m;
 60:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
 61:   PetscCall(PetscMalloc1(m * dim, &coords));
 62:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
 63:   PetscCall(PetscRandomGetValuesReal(rdm, m * dim, coords));
 64:   PetscCall(PetscCalloc1(M * dim, &gcoords));
 65:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, NULL, &B));
 66:   PetscCall(MatSetRandom(B, rdm));
 67:   PetscCall(MatGetOwnershipRange(B, &begin, NULL));
 68:   PetscCall(PetscArraycpy(gcoords + begin * dim, coords, m * dim));
 69:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, gcoords, M * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
 70:   PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &A));
 71:   PetscCall(MatSetOption(A, MAT_SYMMETRIC, sym));
 72:   PetscCall(MatSetFromOptions(A));
 73:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 74:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 75:   PetscCall(MatViewFromOptions(A, NULL, "-A_view"));
 76:   PetscCall(MatGetOwnershipIS(A, is, NULL));
 77:   PetscCall(ISDuplicate(is[0], is + 1));
 78:   PetscCall(MatIncreaseOverlap(A, 1, is, 2));
 79:   PetscCall(MatSetBlockSize(A, 2));
 80:   PetscCall(MatIncreaseOverlap(A, 1, is + 1, 1));
 81:   PetscCall(ISGetBlockSize(is[1], &bs));
 82:   PetscCheck(bs == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect block size %" PetscInt_FMT " != 2", bs);
 83:   PetscCall(MatSetBlockSize(A, 1));
 84:   PetscCall(ISEqual(is[0], is[1], &flg));
 85:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unequal index sets");
 86:   PetscCall(ISDestroy(is));
 87:   PetscCall(ISDestroy(is + 1));
 88:   PetscCall(MatCreateVecs(A, &right, &left));
 89:   PetscCall(VecSetRandom(right, rdm));
 90:   PetscCall(MatMult(A, right, left));
 91:   PetscCall(MatHtoolGetPermutationSource(A, &iss));
 92:   PetscCall(MatHtoolGetPermutationTarget(A, &ist));
 93:   PetscCall(VecDuplicate(left, &perm));
 94:   PetscCall(VecCopy(left, perm));
 95:   PetscCall(VecPermute(perm, ist, PETSC_FALSE));
 96:   PetscCall(VecPermute(right, iss, PETSC_FALSE));
 97:   PetscCall(MatHtoolUsePermutation(A, PETSC_FALSE));
 98:   PetscCall(MatMult(A, right, left));
 99:   PetscCall(VecAXPY(left, -1.0, perm));
100:   PetscCall(VecNorm(left, NORM_INFINITY, &norm));
101:   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);
102:   PetscCall(MatHtoolUsePermutation(A, PETSC_TRUE));
103:   PetscCall(VecDestroy(&perm));
104:   PetscCall(VecDestroy(&left));
105:   PetscCall(VecDestroy(&right));
106:   PetscCall(ISDestroy(&ist));
107:   PetscCall(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:     PetscCall(MatDestroy(&B));
111:     PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &B));
112:     PetscCall(MatSetOption(B, MAT_SYMMETRIC, (PetscBool)!sym));
113:     PetscCall(MatSetFromOptions(B));
114:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
115:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
116:     PetscCall(MatViewFromOptions(B, NULL, "-B_view"));
117:     PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &P));
118:     PetscCall(MatNorm(P, NORM_FROBENIUS, &relative));
119:     PetscCall(MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &R));
120:     PetscCall(MatAXPY(R, -1.0, P, SAME_NONZERO_PATTERN));
121:     PetscCall(MatNorm(R, NORM_INFINITY, &norm));
122:     PetscCheck(PetscAbsReal(norm / relative) <= epsilon, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "||A(!symmetric)-A(symmetric)|| = %g (> %g)", (double)PetscAbsReal(norm / relative), (double)epsilon);
123:     PetscCall(MatDestroy(&B));
124:     PetscCall(MatDestroy(&R));
125:     PetscCall(MatDestroy(&P));
126:   } else {
127:     PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &D));
128:     PetscCall(MatViewFromOptions(D, NULL, "-D_view"));
129:     PetscCall(MatMultEqual(A, D, 10, &flg));
130:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Ax != Dx");
131:     PetscCall(MatMultTransposeEqual(A, D, 10, &flg));
132:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A^Tx != D^Tx");
133:     PetscCall(MatMultAddEqual(A, D, 10, &flg));
134:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "y+Ax != y+Dx");
135:     PetscCall(MatGetOwnershipRange(B, &begin, NULL));
136:     PetscCall(MatDenseGetArrayWrite(D, &ptr));
137:     for (PetscInt i = begin; i < m + begin; ++i)
138:       for (PetscInt j = 0; j < M; ++j) PetscCall(GenEntries(dim, 1, 1, &i, &j, ptr + i - begin + j * m, gcoords));
139:     PetscCall(MatDenseRestoreArrayWrite(D, &ptr));
140:     PetscCall(MatMultEqual(A, D, 10, &flg));
141:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Ax != Dx");
142:     PetscCall(MatTranspose(D, MAT_INPLACE_MATRIX, &D));
143:     PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
144:     PetscCall(MatMultEqual(AT, D, 10, &flg));
145:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A^Tx != D^Tx");
146:     PetscCall(MatTranspose(A, MAT_REUSE_MATRIX, &AT));
147:     PetscCall(MatMultEqual(AT, D, 10, &flg));
148:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A^Tx != D^Tx");
149:     PetscCall(MatAXPY(D, -1.0, AT, SAME_NONZERO_PATTERN));
150:     PetscCall(MatNorm(D, NORM_INFINITY, &norm));
151:     PetscCheck(PetscAbsReal(norm) <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "||A-D|| = %g (> %g)", (double)norm, (double)PETSC_SMALL);
152:     PetscCall(MatDestroy(&AT));
153:     PetscCall(MatDestroy(&D));
154:     PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &P));
155:     PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
156:     PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
157:     PetscCall(MatMatMultEqual(A, B, P, 10, &flg));
158:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ABx != Px");
159:     PetscCall(MatTransposeMatMultEqual(A, B, P, 10, &flg));
160:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "A^TBx != P^Tx");
161:     PetscCall(MatDestroy(&B));
162:     PetscCall(MatDestroy(&P));
163:     if (n) {
164:       PetscCall(PetscMalloc1(n * dim, &scoords));
165:       PetscCall(PetscRandomGetValuesReal(rdm, n * dim, scoords));
166:       N = n;
167:       PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &N, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
168:       PetscCall(PetscCalloc1(N * dim, &gscoords));
169:       PetscCallMPI(MPI_Exscan(&n, &begin, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
170:       PetscCall(PetscArraycpy(gscoords + begin * dim, scoords, n * dim));
171:       PetscCallMPI(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:       PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, n, M, N, dim, coords, scoords, kernel, ctx, &R));
176:       PetscCall(MatSetFromOptions(R));
177:       PetscCall(MatAssemblyBegin(R, MAT_FINAL_ASSEMBLY));
178:       PetscCall(MatAssemblyEnd(R, MAT_FINAL_ASSEMBLY));
179:       PetscCall(MatViewFromOptions(R, NULL, "-R_view"));
180:       PetscCall(MatConvert(R, MATDENSE, MAT_INITIAL_MATRIX, &D));
181:       PetscCall(MatViewFromOptions(D, NULL, "-D_view"));
182:       PetscCall(MatMultEqual(R, D, 10, &flg));
183:       PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Rx != Dx");
184:       PetscCall(MatTranspose(D, MAT_INPLACE_MATRIX, &D));
185:       PetscCall(MatTranspose(R, MAT_INITIAL_MATRIX, &RT));
186:       PetscCall(MatMultEqual(RT, D, 10, &flg));
187:       PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "R^Tx != D^Tx");
188:       PetscCall(MatTranspose(R, MAT_REUSE_MATRIX, &RT));
189:       PetscCall(MatMultEqual(RT, D, 10, &flg));
190:       PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "R^Tx != D^Tx");
191:       PetscCall(MatDestroy(&RT));
192:       PetscCall(MatDestroy(&D));
193:       PetscCall(MatCreateDense(PETSC_COMM_WORLD, n, PETSC_DECIDE, PETSC_DETERMINE, K, NULL, &B));
194:       PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
195:       PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
196:       PetscCall(MatSetRandom(B, rdm));
197:       PetscCall(MatMatMult(R, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &P));
198:       PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
199:       PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
200:       PetscCall(MatMatMultEqual(R, B, P, 10, &flg));
201:       PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "RBx != Px");
202:       PetscCall(MatDestroy(&B));
203:       PetscCall(MatDestroy(&P));
204:       PetscCall(MatCreateVecs(R, &right, &left));
205:       PetscCall(VecSetRandom(right, rdm));
206:       PetscCall(MatMult(R, right, left));
207:       PetscCall(MatHtoolGetPermutationSource(R, &iss));
208:       PetscCall(MatHtoolGetPermutationTarget(R, &ist));
209:       PetscCall(VecDuplicate(left, &perm));
210:       PetscCall(VecCopy(left, perm));
211:       PetscCall(VecPermute(perm, ist, PETSC_FALSE));
212:       PetscCall(VecPermute(right, iss, PETSC_FALSE));
213:       PetscCall(MatHtoolUsePermutation(R, PETSC_FALSE));
214:       PetscCall(MatMult(R, right, left));
215:       PetscCall(VecAXPY(left, -1.0, perm));
216:       PetscCall(VecNorm(left, NORM_INFINITY, &norm));
217:       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);
218:       PetscCall(MatHtoolUsePermutation(R, PETSC_TRUE));
219:       PetscCall(VecDestroy(&perm));
220:       PetscCall(VecDestroy(&left));
221:       PetscCall(VecDestroy(&right));
222:       PetscCall(ISDestroy(&ist));
223:       PetscCall(ISDestroy(&iss));
224:       PetscCall(MatDestroy(&R));
225:       PetscCall(PetscFree(gscoords));
226:       PetscCall(PetscFree(scoords));
227:     }
228:   }
229:   PetscCall(PetscRandomDestroy(&rdm));
230:   PetscCall(MatDestroy(&A));
231:   PetscCall(PetscFree(gcoords));
232:   PetscCall(PetscFree(coords));
233:   PetscCall(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}

248:    test:
249:       requires: htool
250:       suffix: 2
251:       nsize: 4
252:       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}
253:       output_file: output/ex241_1.out

255: TEST*/