1 static char help[] = "Tests MATFACTORHTOOL\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, void *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-1 + PetscSqrtReal(diff)); 16 } 17 } 18 PetscFunctionReturn(PETSC_SUCCESS); 19 } 20 21 int main(int argc, char **argv) 22 { 23 Mat A, Ad, F, Fd, X, Xd, B; 24 Vec x, xd, b; 25 PetscInt m = 100, dim = 3, M, K = 10, begin, n = 0; 26 PetscMPIInt size; 27 PetscReal *coords, *gcoords, norm, epsilon; 28 MatHtoolKernelFn *kernel = GenEntries; 29 PetscBool flg, sym = PETSC_FALSE; 30 PetscRandom rdm; 31 MatSolverType type; 32 33 PetscFunctionBeginUser; 34 PetscCall(PetscInitialize(&argc, &argv, (char *)NULL, help)); 35 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m_local", &m, NULL)); 36 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n_local", &n, NULL)); 37 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL)); 38 PetscCall(PetscOptionsGetInt(NULL, NULL, "-K", &K, NULL)); 39 PetscCall(PetscOptionsGetBool(NULL, NULL, "-symmetric", &sym, NULL)); 40 PetscCall(PetscOptionsGetReal(NULL, NULL, "-mat_htool_epsilon", &epsilon, NULL)); 41 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 42 M = size * m; 43 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 44 PetscCall(PetscMalloc1(m * dim, &coords)); 45 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm)); 46 PetscCall(PetscRandomGetValuesReal(rdm, m * dim, coords)); 47 PetscCall(PetscCalloc1(M * dim, &gcoords)); 48 PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, NULL, &B)); 49 PetscCall(MatSetRandom(B, rdm)); 50 PetscCall(MatGetOwnershipRange(B, &begin, NULL)); 51 PetscCall(PetscArraycpy(gcoords + begin * dim, coords, m * dim)); 52 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, gcoords, M * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD)); 53 PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &A)); 54 PetscCall(MatSetOption(A, MAT_SYMMETRIC, sym)); 55 PetscCall(MatSetFromOptions(A)); 56 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 57 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 58 PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &Ad)); 59 PetscCall(MatMultEqual(A, Ad, 10, &flg)); 60 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Ax != Adx"); 61 PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, NULL, &X)); 62 PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, NULL, &Xd)); 63 PetscCall(MatViewFromOptions(A, NULL, "-A")); 64 PetscCall(MatViewFromOptions(Ad, NULL, "-Ad")); 65 PetscCall(MatViewFromOptions(B, NULL, "-B")); 66 for (PetscInt i = 0; i < 2; ++i) { 67 PetscCall(MatGetFactor(A, MATSOLVERHTOOL, i == 0 ? MAT_FACTOR_LU : MAT_FACTOR_CHOLESKY, &F)); 68 PetscCall(MatGetFactor(Ad, MATSOLVERPETSC, i == 0 ? MAT_FACTOR_LU : MAT_FACTOR_CHOLESKY, &Fd)); 69 PetscCall(MatFactorGetSolverType(F, &type)); 70 PetscCall(PetscStrncmp(type, MATSOLVERHTOOL, 5, &flg)); 71 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MATSOLVERHTOOL != htool"); 72 if (i == 0) { 73 PetscCall(MatLUFactorSymbolic(F, A, NULL, NULL, NULL)); 74 PetscCall(MatLUFactorNumeric(F, A, NULL)); 75 PetscCall(MatLUFactorSymbolic(Fd, Ad, NULL, NULL, NULL)); 76 PetscCall(MatLUFactorNumeric(Fd, Ad, NULL)); 77 } else { 78 PetscCall(MatCholeskyFactorSymbolic(F, A, NULL, NULL)); 79 PetscCall(MatCholeskyFactorNumeric(F, A, NULL)); 80 PetscCall(MatCholeskyFactorSymbolic(Fd, Ad, NULL, NULL)); 81 PetscCall(MatCholeskyFactorNumeric(Fd, Ad, NULL)); 82 } 83 PetscCall(MatMatSolve(F, B, X)); 84 PetscCall(MatMatSolve(Fd, B, Xd)); 85 PetscCall(MatViewFromOptions(X, NULL, "-X")); 86 PetscCall(MatViewFromOptions(Xd, NULL, "-Xd")); 87 PetscCall(MatAXPY(Xd, -1.0, X, SAME_NONZERO_PATTERN)); 88 PetscCall(MatNorm(Xd, NORM_INFINITY, &norm)); 89 PetscCall(MatViewFromOptions(Xd, NULL, "-MatMatSolve")); 90 if (norm > 0.01) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: norm of residual for MatMatSolve %g\n", (double)norm)); 91 if (!PetscDefined(USE_COMPLEX) || i == 0) { 92 PetscCall(MatMatSolveTranspose(F, B, X)); 93 PetscCall(MatMatSolveTranspose(Fd, B, Xd)); 94 PetscCall(MatAXPY(Xd, -1.0, X, SAME_NONZERO_PATTERN)); 95 PetscCall(MatNorm(Xd, NORM_INFINITY, &norm)); 96 PetscCall(MatViewFromOptions(Xd, NULL, "-MatMatSolveTranspose")); 97 if (norm > 0.01) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: norm of residual for MatMatSolveTranspose %g\n", (double)norm)); 98 } 99 PetscCall(MatDenseGetColumnVecRead(B, 0, &b)); 100 PetscCall(MatDenseGetColumnVecWrite(X, 0, &x)); 101 PetscCall(MatDenseGetColumnVecWrite(Xd, 0, &xd)); 102 PetscCall(MatSolve(F, b, x)); 103 PetscCall(MatSolve(Fd, b, xd)); 104 PetscCall(VecAXPY(xd, -1.0, x)); 105 PetscCall(VecNorm(xd, NORM_INFINITY, &norm)); 106 PetscCall(MatViewFromOptions(Xd, NULL, "-MatSolve")); 107 if (norm > 0.01) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: norm of residual for MatSolve %g\n", (double)norm)); 108 if (!PetscDefined(USE_COMPLEX) || i == 0) { 109 PetscCall(MatSolveTranspose(F, b, x)); 110 PetscCall(MatSolveTranspose(Fd, b, xd)); 111 PetscCall(VecAXPY(xd, -1.0, x)); 112 PetscCall(VecNorm(xd, NORM_INFINITY, &norm)); 113 PetscCall(MatViewFromOptions(Xd, NULL, "-MatSolveTranspose")); 114 if (norm > 0.01) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: norm of residual for MatSolveTranspose %g\n", (double)norm)); 115 } 116 PetscCall(MatDenseRestoreColumnVecWrite(Xd, 0, &xd)); 117 PetscCall(MatDenseRestoreColumnVecWrite(X, 0, &x)); 118 PetscCall(MatDenseRestoreColumnVecRead(B, 0, &b)); 119 PetscCall(MatDestroy(&Fd)); 120 PetscCall(MatDestroy(&F)); 121 } 122 PetscCall(MatDestroy(&Xd)); 123 PetscCall(MatDestroy(&X)); 124 PetscCall(PetscRandomDestroy(&rdm)); 125 PetscCall(MatDestroy(&Ad)); 126 PetscCall(MatDestroy(&A)); 127 PetscCall(MatDestroy(&B)); 128 PetscCall(PetscFree(gcoords)); 129 PetscCall(PetscFree(coords)); 130 PetscCall(PetscFinalize()); 131 return 0; 132 } 133 134 /*TEST 135 136 build: 137 requires: htool 138 139 test: 140 requires: htool 141 suffix: 1 142 nsize: 1 143 args: -mat_htool_epsilon 1.0e-11 144 output_file: output/empty.out 145 146 TEST*/ 147