1 static char help[] = "Tests MATFACTORHTOOL\n\n";
2
3 #include <petscmat.h>
4
GenEntries(PetscInt sdim,PetscInt M,PetscInt N,const PetscInt * J,const PetscInt * K,PetscScalar * ptr,PetscCtx ctx)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-1 + PetscSqrtReal(diff));
16 }
17 }
18 PetscFunctionReturn(PETSC_SUCCESS);
19 }
20
main(int argc,char ** argv)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