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: PetscInitialize(&argc,&args,(char*)0,help);
38: MPI_Comm_size(PETSC_COMM_WORLD,&size);
40: PetscRandomCreate(PETSC_COMM_WORLD,&r);
41: PetscRandomSetFromOptions(r);
43: sdim = 2;
44: PetscOptionsGetInt(NULL,NULL,"-sdim",&sdim,NULL);
45: n = 32;
46: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
47: eta = 0.6;
48: PetscOptionsGetReal(NULL,NULL,"-eta",&eta,NULL);
49: leafsize = 32;
50: PetscOptionsGetInt(NULL,NULL,"-leafsize",&leafsize,NULL);
51: basisord = 8;
52: PetscOptionsGetInt(NULL,NULL,"-basisord",&basisord,NULL);
54: /* Create random points */
55: PetscMalloc1(sdim*n,&coords);
56: PetscRandomGetValuesReal(r,sdim*n,coords);
58: fctx.lambda = 0.01;
59: PetscOptionsGetReal(NULL,NULL,"-lambda",&fctx.lambda,NULL);
60: PetscRandomGetValueReal(r,&fctx.sigma);
61: PetscOptionsGetReal(NULL,NULL,"-sigma",&fctx.sigma,NULL);
62: PetscMalloc1(sdim,&fctx.l);
63: PetscRandomGetValuesReal(r,sdim,fctx.l);
64: PetscOptionsGetRealArray(NULL,NULL,"-l",fctx.l,(i=sdim,&i),NULL);
65: PetscOptionsGetReal(NULL,NULL,"-scale",&scale,NULL);
67: /* Populate dense matrix for comparisons */
68: {
69: PetscInt i,j;
71: MatCreateDense(PETSC_COMM_WORLD,n,n,PETSC_DECIDE,PETSC_DECIDE,NULL,&Ad);
72: for (i = 0; i < n; i++) {
73: for (j = 0; j < n; j++) {
74: MatSetValue(Ad,i,j,RBF(sdim,coords + i*sdim,coords + j*sdim,&fctx),INSERT_VALUES);
75: }
76: }
77: MatAssemblyBegin(Ad,MAT_FINAL_ASSEMBLY);
78: MatAssemblyEnd(Ad,MAT_FINAL_ASSEMBLY);
79: }
81: /* Create and assemble the matrix */
82: MatCreateH2OpusFromKernel(PETSC_COMM_WORLD,n,n,PETSC_DECIDE,PETSC_DECIDE,sdim,coords,PETSC_FALSE,RBF,&fctx,eta,leafsize,basisord,&A);
83: MatSetOption(A,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
84: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
85: PetscObjectSetName((PetscObject)A,"RBF");
86: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
87: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
88: MatViewFromOptions(A,NULL,"-rbf_view");
90: MatCreateVecs(A,&x,&b);
91: VecDuplicate(x,&u);
92: VecDuplicate(x,&d);
94: {
95: PetscReal norm;
96: MatComputeOperator(A,MATDENSE,&Ae);
97: MatAXPY(Ae,-1.0,Ad,SAME_NONZERO_PATTERN);
98: MatGetDiagonal(Ae,d);
99: MatViewFromOptions(Ae,NULL,"-A_view");
100: MatViewFromOptions(Ae,NULL,"-D_view");
101: MatNorm(Ae,NORM_FROBENIUS,&norm);
102: PetscPrintf(PETSC_COMM_WORLD,"Approx err %g\n",norm);
103: VecNorm(d,NORM_2,&norm);
104: PetscPrintf(PETSC_COMM_WORLD,"Approx err (diag) %g\n",norm);
105: MatDestroy(&Ae);
106: }
108: VecSet(u,1.0);
109: MatMult(Ad,u,b);
110: MatViewFromOptions(Ad,NULL,"-Ad_view");
111: KSPCreate(PETSC_COMM_WORLD,&ksp);
112: KSPSetOperators(ksp,Ad,A);
113: KSPGetPC(ksp,&pc);
114: PCSetType(pc,PCH2OPUS);
115: KSPSetFromOptions(ksp);
116: /* we can also pass the points coordinates
117: In this case it is not needed, since the preconditioning
118: matrix is of type H2OPUS */
119: PCSetCoordinates(pc,sdim,n,coords);
121: KSPSolve(ksp,b,x);
122: VecAXPY(x,-1.0,u);
123: VecNorm(x,NORM_2,&norm);
124: KSPGetIterationNumber(ksp,&its);
125: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
127: /* change lambda and reassemble */
128: VecSet(x,(scale-1.)*fctx.lambda);
129: MatDiagonalSet(Ad,x,ADD_VALUES);
130: fctx.lambda *= scale;
131: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
132: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
133: {
134: PetscReal norm;
135: MatComputeOperator(A,MATDENSE,&Ae);
136: MatAXPY(Ae,-1.0,Ad,SAME_NONZERO_PATTERN);
137: MatGetDiagonal(Ae,d);
138: MatViewFromOptions(Ae,NULL,"-A_view");
139: MatViewFromOptions(Ae,NULL,"-D_view");
140: MatNorm(Ae,NORM_FROBENIUS,&norm);
141: PetscPrintf(PETSC_COMM_WORLD,"Approx err %g\n",norm);
142: VecNorm(d,NORM_2,&norm);
143: PetscPrintf(PETSC_COMM_WORLD,"Approx err (diag) %g\n",norm);
144: MatDestroy(&Ae);
145: }
146: KSPSetOperators(ksp,Ad,A);
147: MatMult(Ad,u,b);
148: KSPSolve(ksp,b,x);
149: MatMult(Ad,x,u);
150: VecAXPY(u,-1.0,b);
151: VecNorm(u,NORM_2,&norm);
152: KSPGetIterationNumber(ksp,&its);
153: PetscPrintf(PETSC_COMM_WORLD,"Residual norm error %g, Iterations %D\n",(double)norm,its);
155: PetscFree(coords);
156: PetscFree(fctx.l);
157: PetscRandomDestroy(&r);
158: VecDestroy(&x);
159: VecDestroy(&u);
160: VecDestroy(&d);
161: VecDestroy(&b);
162: MatDestroy(&A);
163: MatDestroy(&Ad);
164: KSPDestroy(&ksp);
165: PetscFinalize();
166: return 0;
167: }
169: /*TEST
171: build:
172: requires: h2opus
174: test:
175: requires: h2opus !single
176: suffix: 1
177: args: -ksp_error_if_not_converged -pc_h2opus_monitor
179: test:
180: requires: h2opus !single
181: suffix: 1_ns
182: output_file: output/ex21_1.out
183: args: -ksp_error_if_not_converged -pc_h2opus_monitor -pc_h2opus_hyperorder 2
185: test:
186: requires: h2opus !single
187: suffix: 2
188: args: -ksp_error_if_not_converged -pc_h2opus_monitor -pc_h2opus_hyperorder 4
190: TEST*/