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;
38: PetscInitialize(&argc, &args, (char *)0, help);
39: MPI_Comm_size(PETSC_COMM_WORLD, &size);
41: PetscRandomCreate(PETSC_COMM_WORLD, &r);
42: PetscRandomSetFromOptions(r);
44: sdim = 2;
45: PetscOptionsGetInt(NULL, NULL, "-sdim", &sdim, NULL);
46: n = 32;
47: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
48: eta = 0.6;
49: PetscOptionsGetReal(NULL, NULL, "-eta", &eta, NULL);
50: leafsize = 32;
51: PetscOptionsGetInt(NULL, NULL, "-leafsize", &leafsize, NULL);
52: basisord = 8;
53: PetscOptionsGetInt(NULL, NULL, "-basisord", &basisord, NULL);
55: /* Create random points */
56: PetscMalloc1(sdim * n, &coords);
57: PetscRandomGetValuesReal(r, sdim * n, coords);
59: fctx.lambda = 0.01;
60: PetscOptionsGetReal(NULL, NULL, "-lambda", &fctx.lambda, NULL);
61: PetscRandomGetValueReal(r, &fctx.sigma);
62: PetscOptionsGetReal(NULL, NULL, "-sigma", &fctx.sigma, NULL);
63: PetscMalloc1(sdim, &fctx.l);
64: PetscRandomGetValuesReal(r, sdim, fctx.l);
65: PetscOptionsGetRealArray(NULL, NULL, "-l", fctx.l, (i = sdim, &i), NULL);
66: PetscOptionsGetReal(NULL, NULL, "-scale", &scale, NULL);
68: /* Populate dense matrix for comparisons */
69: {
70: PetscInt i, j;
72: 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++) MatSetValue(Ad, i, j, RBF(sdim, coords + i * sdim, coords + j * sdim, &fctx), INSERT_VALUES);
75: }
76: MatAssemblyBegin(Ad, MAT_FINAL_ASSEMBLY);
77: MatAssemblyEnd(Ad, MAT_FINAL_ASSEMBLY);
78: }
80: /* Create and assemble the matrix */
81: MatCreateH2OpusFromKernel(PETSC_COMM_WORLD, n, n, PETSC_DECIDE, PETSC_DECIDE, sdim, coords, PETSC_FALSE, RBF, &fctx, eta, leafsize, basisord, &A);
82: MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
83: MatSetOption(A, MAT_SYMMETRY_ETERNAL, PETSC_TRUE);
84: PetscObjectSetName((PetscObject)A, "RBF");
85: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
86: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
87: MatViewFromOptions(A, NULL, "-rbf_view");
89: MatCreateVecs(A, &x, &b);
90: VecDuplicate(x, &u);
91: VecDuplicate(x, &d);
93: {
94: PetscReal norm;
95: MatComputeOperator(A, MATDENSE, &Ae);
96: MatAXPY(Ae, -1.0, Ad, SAME_NONZERO_PATTERN);
97: MatGetDiagonal(Ae, d);
98: MatViewFromOptions(Ae, NULL, "-A_view");
99: MatViewFromOptions(Ae, NULL, "-D_view");
100: MatNorm(Ae, NORM_FROBENIUS, &norm);
101: PetscPrintf(PETSC_COMM_WORLD, "Approx err %g\n", norm);
102: VecNorm(d, NORM_2, &norm);
103: PetscPrintf(PETSC_COMM_WORLD, "Approx err (diag) %g\n", norm);
104: MatDestroy(&Ae);
105: }
107: VecSet(u, 1.0);
108: MatMult(Ad, u, b);
109: MatViewFromOptions(Ad, NULL, "-Ad_view");
110: KSPCreate(PETSC_COMM_WORLD, &ksp);
111: KSPSetOperators(ksp, Ad, A);
112: KSPGetPC(ksp, &pc);
113: PCSetType(pc, PCH2OPUS);
114: 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: PCSetCoordinates(pc, sdim, n, coords);
120: KSPSolve(ksp, b, x);
121: VecAXPY(x, -1.0, u);
122: VecNorm(x, NORM_2, &norm);
123: KSPGetIterationNumber(ksp, &its);
124: PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);
126: /* change lambda and reassemble */
127: VecSet(x, (scale - 1.) * fctx.lambda);
128: MatDiagonalSet(Ad, x, ADD_VALUES);
129: fctx.lambda *= scale;
130: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
131: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
132: {
133: PetscReal norm;
134: MatComputeOperator(A, MATDENSE, &Ae);
135: MatAXPY(Ae, -1.0, Ad, SAME_NONZERO_PATTERN);
136: MatGetDiagonal(Ae, d);
137: MatViewFromOptions(Ae, NULL, "-A_view");
138: MatViewFromOptions(Ae, NULL, "-D_view");
139: MatNorm(Ae, NORM_FROBENIUS, &norm);
140: PetscPrintf(PETSC_COMM_WORLD, "Approx err %g\n", norm);
141: VecNorm(d, NORM_2, &norm);
142: PetscPrintf(PETSC_COMM_WORLD, "Approx err (diag) %g\n", norm);
143: MatDestroy(&Ae);
144: }
145: KSPSetOperators(ksp, Ad, A);
146: MatMult(Ad, u, b);
147: KSPSolve(ksp, b, x);
148: MatMult(Ad, x, u);
149: VecAXPY(u, -1.0, b);
150: VecNorm(u, NORM_2, &norm);
151: KSPGetIterationNumber(ksp, &its);
152: PetscPrintf(PETSC_COMM_WORLD, "Residual norm error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);
154: PetscFree(coords);
155: PetscFree(fctx.l);
156: PetscRandomDestroy(&r);
157: VecDestroy(&x);
158: VecDestroy(&u);
159: VecDestroy(&d);
160: VecDestroy(&b);
161: MatDestroy(&A);
162: MatDestroy(&Ad);
163: KSPDestroy(&ksp);
164: 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*/