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, PetscCtx 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: #if PetscDefined(USE_COMPLEX)
22: /* GenEntries() scaled by the unitary diag(exp(i x_0)) on both sides: Hermitian, but not symmetric */
23: static PetscErrorCode GenEntriesHermitian(PetscInt sdim, PetscInt M, PetscInt N, const PetscInt *J, const PetscInt *K, PetscScalar *ptr, PetscCtx ctx)
24: {
25: PetscReal *coords = (PetscReal *)(ctx);
27: PetscFunctionBeginUser;
28: PetscCall(GenEntries(sdim, M, N, J, K, ptr, ctx));
29: for (PetscInt j = 0; j < M; j++) {
30: for (PetscInt k = 0; k < N; k++) ptr[j + M * k] *= PetscExpComplex(PETSC_i * (coords[J[j] * sdim] - coords[K[k] * sdim]));
31: }
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
34: #endif
36: int main(int argc, char **argv)
37: {
38: Mat A, Ad, F, Fd, X, Xd, B;
39: Vec x, xd, b;
40: PetscInt m = 100, dim = 3, M, K = 10, begin, n = 0;
41: PetscMPIInt size;
42: PetscReal *coords, *gcoords, norm, epsilon;
43: MatHtoolKernelFn *kernel = GenEntries;
44: PetscBool flg, sym = PETSC_FALSE, set_sym, herm = PETSC_FALSE, set_herm, spd = PETSC_FALSE, set_spd;
45: PetscRandom rdm;
46: MatSolverType type;
48: PetscFunctionBeginUser;
49: PetscCall(PetscInitialize(&argc, &argv, (char *)NULL, help));
50: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m_local", &m, NULL));
51: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n_local", &n, NULL));
52: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL));
53: PetscCall(PetscOptionsGetInt(NULL, NULL, "-K", &K, NULL));
54: PetscCall(PetscOptionsGetBool(NULL, NULL, "-symmetric", &sym, &set_sym));
55: PetscCall(PetscOptionsGetBool(NULL, NULL, "-hermitian", &herm, &set_herm));
56: PetscCall(PetscOptionsGetBool(NULL, NULL, "-spd", &spd, &set_spd));
57: #if PetscDefined(USE_COMPLEX)
58: if (herm) kernel = GenEntriesHermitian;
59: #endif
60: PetscCall(PetscOptionsGetReal(NULL, NULL, "-mat_htool_epsilon", &epsilon, NULL));
61: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
62: M = size * m;
63: PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
64: PetscCall(PetscMalloc1(m * dim, &coords));
65: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
66: PetscCall(PetscRandomGetValuesReal(rdm, m * dim, coords));
67: PetscCall(PetscCalloc1(M * dim, &gcoords));
68: PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, NULL, &B));
69: PetscCall(MatSetRandom(B, rdm));
70: PetscCall(MatGetOwnershipRange(B, &begin, NULL));
71: PetscCall(PetscArraycpy(gcoords + begin * dim, coords, m * dim));
72: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, gcoords, M * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
73: PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &A));
74: if (set_sym) PetscCall(MatSetOption(A, MAT_SYMMETRIC, sym));
75: if (set_herm) PetscCall(MatSetOption(A, MAT_HERMITIAN, herm));
76: if (set_spd) PetscCall(MatSetOption(A, MAT_SPD, spd));
77: PetscCall(MatSetFromOptions(A));
78: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
79: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
80: PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &Ad));
81: PetscCall(MatPropagateSymmetryOptions(A, Ad));
82: PetscCall(MatMultEqual(A, Ad, 10, &flg));
83: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Ax != Adx");
84: PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, NULL, &X));
85: PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, M, K, NULL, &Xd));
86: PetscCall(MatViewFromOptions(A, NULL, "-A"));
87: PetscCall(MatViewFromOptions(Ad, NULL, "-Ad"));
88: PetscCall(MatViewFromOptions(B, NULL, "-B"));
89: for (PetscInt i = sym || herm || spd ? 1 : 0; i < 2; ++i) { // LU requires full storage, i.e., a MATHTOOL not flagged symmetric
90: PetscCall(MatGetFactor(A, MATSOLVERHTOOL, i == 0 ? MAT_FACTOR_LU : MAT_FACTOR_CHOLESKY, &F));
91: PetscCall(MatGetFactor(Ad, MATSOLVERPETSC, i == 0 ? MAT_FACTOR_LU : MAT_FACTOR_CHOLESKY, &Fd));
92: PetscCall(MatFactorGetSolverType(F, &type));
93: PetscCall(PetscStrncmp(type, MATSOLVERHTOOL, 5, &flg));
94: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MATSOLVERHTOOL != htool");
95: if (i == 0) {
96: PetscCall(MatLUFactorSymbolic(F, A, NULL, NULL, NULL));
97: PetscCall(MatLUFactorNumeric(F, A, NULL));
98: PetscCall(MatLUFactorSymbolic(Fd, Ad, NULL, NULL, NULL));
99: PetscCall(MatLUFactorNumeric(Fd, Ad, NULL));
100: } else {
101: PetscCall(MatCholeskyFactorSymbolic(F, A, NULL, NULL));
102: PetscCall(MatCholeskyFactorNumeric(F, A, NULL));
103: PetscCall(MatCholeskyFactorSymbolic(Fd, Ad, NULL, NULL));
104: PetscCall(MatCholeskyFactorNumeric(Fd, Ad, NULL));
105: }
106: PetscCall(MatMatSolve(F, B, X));
107: PetscCall(MatMatSolve(Fd, B, Xd));
108: PetscCall(MatViewFromOptions(X, NULL, "-X"));
109: PetscCall(MatViewFromOptions(Xd, NULL, "-Xd"));
110: PetscCall(MatAXPY(Xd, -1.0, X, SAME_NONZERO_PATTERN));
111: PetscCall(MatNorm(Xd, NORM_INFINITY, &norm));
112: PetscCall(MatViewFromOptions(Xd, NULL, "-MatMatSolve"));
113: if (norm > 0.01) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: norm of residual for MatMatSolve %g\n", (double)norm));
114: if (!PetscDefined(USE_COMPLEX) || i == 0) {
115: PetscCall(MatMatSolveTranspose(F, B, X));
116: PetscCall(MatMatSolveTranspose(Fd, B, Xd));
117: PetscCall(MatAXPY(Xd, -1.0, X, SAME_NONZERO_PATTERN));
118: PetscCall(MatNorm(Xd, NORM_INFINITY, &norm));
119: PetscCall(MatViewFromOptions(Xd, NULL, "-MatMatSolveTranspose"));
120: if (norm > 0.01) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: norm of residual for MatMatSolveTranspose %g\n", (double)norm));
121: }
122: PetscCall(MatDenseGetColumnVecRead(B, 0, &b));
123: PetscCall(MatDenseGetColumnVecWrite(X, 0, &x));
124: PetscCall(MatDenseGetColumnVecWrite(Xd, 0, &xd));
125: PetscCall(MatSolve(F, b, x));
126: PetscCall(MatSolve(Fd, b, xd));
127: PetscCall(VecAXPY(xd, -1.0, x));
128: PetscCall(VecNorm(xd, NORM_INFINITY, &norm));
129: PetscCall(MatViewFromOptions(Xd, NULL, "-MatSolve"));
130: if (norm > 0.01) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: norm of residual for MatSolve %g\n", (double)norm));
131: if (!PetscDefined(USE_COMPLEX) || i == 0) {
132: PetscCall(MatSolveTranspose(F, b, x));
133: PetscCall(MatSolveTranspose(Fd, b, xd));
134: PetscCall(VecAXPY(xd, -1.0, x));
135: PetscCall(VecNorm(xd, NORM_INFINITY, &norm));
136: PetscCall(MatViewFromOptions(Xd, NULL, "-MatSolveTranspose"));
137: if (norm > 0.01) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: norm of residual for MatSolveTranspose %g\n", (double)norm));
138: }
139: PetscCall(MatDenseRestoreColumnVecWrite(Xd, 0, &xd));
140: PetscCall(MatDenseRestoreColumnVecWrite(X, 0, &x));
141: PetscCall(MatDenseRestoreColumnVecRead(B, 0, &b));
142: PetscCall(MatDestroy(&Fd));
143: PetscCall(MatDestroy(&F));
144: }
145: PetscCall(MatDestroy(&Xd));
146: PetscCall(MatDestroy(&X));
147: PetscCall(PetscRandomDestroy(&rdm));
148: PetscCall(MatDestroy(&Ad));
149: PetscCall(MatDestroy(&A));
150: PetscCall(MatDestroy(&B));
151: PetscCall(PetscFree(gcoords));
152: PetscCall(PetscFree(coords));
153: PetscCall(PetscFinalize());
154: return 0;
155: }
157: /*TEST
159: build:
160: requires: htool
162: test:
163: requires: htool
164: suffix: 1
165: nsize: 1
166: args: -mat_htool_epsilon 1.0e-11 -symmetric {{false true}shared output}
167: output_file: output/empty.out
169: test:
170: requires: htool
171: suffix: ldlt
172: nsize: 1
173: args: -mat_htool_epsilon 1.0e-11 -spd false -symmetric {{false true}shared output}
174: output_file: output/empty.out
176: test:
177: requires: htool
178: suffix: spd
179: nsize: 1
180: args: -mat_htool_epsilon 1.0e-11 -spd
181: output_file: output/empty.out
183: test:
184: requires: htool complex
185: suffix: hermitian
186: nsize: 1
187: args: -mat_htool_epsilon 1.0e-11 -hermitian -spd {{false true}shared output}
188: output_file: output/empty.out
190: TEST*/