Actual source code: ex268.c
1: static char help[] = "Tests MATFACTORHTOOL\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-1 + PetscSqrtReal(diff));
16: }
17: }
18: PetscFunctionReturn(PETSC_SUCCESS);
19: }
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;
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: }
134: /*TEST
136: build:
137: requires: htool
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
146: TEST*/