1 static char help[] = "Tests MATHTOOL\n\n"; 2 3 #include <petscmat.h> 4 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 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 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