Actual source code: ex21.c
2: static char help[] = "Solves a RBF kernel matrix with KSP and PCH2OPUS.\n\n";
4: #include <petscksp.h>
6: typedef struct {
7: PetscReal sigma;
8: PetscReal *l;
9: PetscReal lambda;
10: } RBFCtx;
12: static PetscScalar RBF(PetscInt sdim, PetscReal x[], PetscReal y[], void *ctx)
13: {
14: RBFCtx *rbfctx = (RBFCtx *)ctx;
15: PetscInt d;
16: PetscReal diff = 0.0;
17: PetscReal s = rbfctx->sigma;
18: PetscReal *l = rbfctx->l;
19: PetscReal lambda = rbfctx->lambda;
21: for (d = 0; d < sdim; d++) diff += (x[d] - y[d]) * (x[d] - y[d]) / (l[d] * l[d]);
22: return s * s * PetscExpReal(-0.5 * diff) + (diff != 0.0 ? 0.0 : lambda);
23: }
25: int main(int argc, char **args)
26: {
27: Vec x, b, u, d;
28: Mat A, Ae = NULL, Ad = NULL;
29: KSP ksp;
30: PetscRandom r;
31: PC pc;
32: PetscReal norm, *coords, eta, scale = 0.5;
33: PetscInt basisord, leafsize, sdim, n, its, i;
34: PetscMPIInt size;
35: RBFCtx fctx;
37: PetscFunctionBeginUser;
38: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
39: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
40: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
41: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &r));
42: PetscCall(PetscRandomSetFromOptions(r));
44: sdim = 2;
45: PetscCall(PetscOptionsGetInt(NULL, NULL, "-sdim", &sdim, NULL));
46: n = 32;
47: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
48: eta = 0.6;
49: PetscCall(PetscOptionsGetReal(NULL, NULL, "-eta", &eta, NULL));
50: leafsize = 32;
51: PetscCall(PetscOptionsGetInt(NULL, NULL, "-leafsize", &leafsize, NULL));
52: basisord = 8;
53: PetscCall(PetscOptionsGetInt(NULL, NULL, "-basisord", &basisord, NULL));
55: /* Create random points */
56: PetscCall(PetscMalloc1(sdim * n, &coords));
57: PetscCall(PetscRandomGetValuesReal(r, sdim * n, coords));
59: fctx.lambda = 0.01;
60: PetscCall(PetscOptionsGetReal(NULL, NULL, "-lambda", &fctx.lambda, NULL));
61: PetscCall(PetscRandomGetValueReal(r, &fctx.sigma));
62: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma", &fctx.sigma, NULL));
63: PetscCall(PetscMalloc1(sdim, &fctx.l));
64: PetscCall(PetscRandomGetValuesReal(r, sdim, fctx.l));
65: PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-l", fctx.l, (i = sdim, &i), NULL));
66: PetscCall(PetscOptionsGetReal(NULL, NULL, "-scale", &scale, NULL));
68: /* Populate dense matrix for comparisons */
69: {
70: PetscInt i, j;
72: PetscCall(MatCreateDense(PETSC_COMM_WORLD, n, n, PETSC_DECIDE, PETSC_DECIDE, NULL, &Ad));
73: for (i = 0; i < n; i++) {
74: for (j = 0; j < n; j++) PetscCall(MatSetValue(Ad, i, j, RBF(sdim, coords + i * sdim, coords + j * sdim, &fctx), INSERT_VALUES));
75: }
76: PetscCall(MatAssemblyBegin(Ad, MAT_FINAL_ASSEMBLY));
77: PetscCall(MatAssemblyEnd(Ad, MAT_FINAL_ASSEMBLY));
78: }
80: /* Create and assemble the matrix */
81: PetscCall(MatCreateH2OpusFromKernel(PETSC_COMM_WORLD, n, n, PETSC_DECIDE, PETSC_DECIDE, sdim, coords, PETSC_FALSE, RBF, &fctx, eta, leafsize, basisord, &A));
82: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
83: PetscCall(MatSetOption(A, MAT_SYMMETRY_ETERNAL, PETSC_TRUE));
84: PetscCall(PetscObjectSetName((PetscObject)A, "RBF"));
85: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
86: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
87: PetscCall(MatViewFromOptions(A, NULL, "-rbf_view"));
89: PetscCall(MatCreateVecs(A, &x, &b));
90: PetscCall(VecDuplicate(x, &u));
91: PetscCall(VecDuplicate(x, &d));
93: {
94: PetscReal norm;
95: PetscCall(MatComputeOperator(A, MATDENSE, &Ae));
96: PetscCall(MatAXPY(Ae, -1.0, Ad, SAME_NONZERO_PATTERN));
97: PetscCall(MatGetDiagonal(Ae, d));
98: PetscCall(MatViewFromOptions(Ae, NULL, "-A_view"));
99: PetscCall(MatViewFromOptions(Ae, NULL, "-D_view"));
100: PetscCall(MatNorm(Ae, NORM_FROBENIUS, &norm));
101: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Approx err %g\n", norm));
102: PetscCall(VecNorm(d, NORM_2, &norm));
103: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Approx err (diag) %g\n", norm));
104: PetscCall(MatDestroy(&Ae));
105: }
107: PetscCall(VecSet(u, 1.0));
108: PetscCall(MatMult(Ad, u, b));
109: PetscCall(MatViewFromOptions(Ad, NULL, "-Ad_view"));
110: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
111: PetscCall(KSPSetOperators(ksp, Ad, A));
112: PetscCall(KSPGetPC(ksp, &pc));
113: PetscCall(PCSetType(pc, PCH2OPUS));
114: PetscCall(KSPSetFromOptions(ksp));
115: /* we can also pass the points coordinates
116: In this case it is not needed, since the preconditioning
117: matrix is of type H2OPUS */
118: PetscCall(PCSetCoordinates(pc, sdim, n, coords));
120: PetscCall(KSPSolve(ksp, b, x));
121: PetscCall(VecAXPY(x, -1.0, u));
122: PetscCall(VecNorm(x, NORM_2, &norm));
123: PetscCall(KSPGetIterationNumber(ksp, &its));
124: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
126: /* change lambda and reassemble */
127: PetscCall(VecSet(x, (scale - 1.) * fctx.lambda));
128: PetscCall(MatDiagonalSet(Ad, x, ADD_VALUES));
129: fctx.lambda *= scale;
130: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
131: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
132: {
133: PetscReal norm;
134: PetscCall(MatComputeOperator(A, MATDENSE, &Ae));
135: PetscCall(MatAXPY(Ae, -1.0, Ad, SAME_NONZERO_PATTERN));
136: PetscCall(MatGetDiagonal(Ae, d));
137: PetscCall(MatViewFromOptions(Ae, NULL, "-A_view"));
138: PetscCall(MatViewFromOptions(Ae, NULL, "-D_view"));
139: PetscCall(MatNorm(Ae, NORM_FROBENIUS, &norm));
140: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Approx err %g\n", norm));
141: PetscCall(VecNorm(d, NORM_2, &norm));
142: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Approx err (diag) %g\n", norm));
143: PetscCall(MatDestroy(&Ae));
144: }
145: PetscCall(KSPSetOperators(ksp, Ad, A));
146: PetscCall(MatMult(Ad, u, b));
147: PetscCall(KSPSolve(ksp, b, x));
148: PetscCall(MatMult(Ad, x, u));
149: PetscCall(VecAXPY(u, -1.0, b));
150: PetscCall(VecNorm(u, NORM_2, &norm));
151: PetscCall(KSPGetIterationNumber(ksp, &its));
152: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
154: PetscCall(PetscFree(coords));
155: PetscCall(PetscFree(fctx.l));
156: PetscCall(PetscRandomDestroy(&r));
157: PetscCall(VecDestroy(&x));
158: PetscCall(VecDestroy(&u));
159: PetscCall(VecDestroy(&d));
160: PetscCall(VecDestroy(&b));
161: PetscCall(MatDestroy(&A));
162: PetscCall(MatDestroy(&Ad));
163: PetscCall(KSPDestroy(&ksp));
164: PetscCall(PetscFinalize());
165: return 0;
166: }
168: /*TEST
170: build:
171: requires: h2opus
173: test:
174: requires: h2opus !single
175: suffix: 1
176: args: -ksp_error_if_not_converged -pc_h2opus_monitor
178: test:
179: requires: h2opus !single
180: suffix: 1_ns
181: output_file: output/ex21_1.out
182: args: -ksp_error_if_not_converged -pc_h2opus_monitor -pc_h2opus_hyperorder 2
184: test:
185: requires: h2opus !single
186: suffix: 2
187: args: -ksp_error_if_not_converged -pc_h2opus_monitor -pc_h2opus_hyperorder 4
189: TEST*/