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;
34: PetscInt basisord,leafsize,sdim,n,its,i;
35: PetscMPIInt size;
36: RBFCtx fctx;
38: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
39: MPI_Comm_size(PETSC_COMM_WORLD,&size);
40: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");
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++) {
75: MatSetValue(Ad,i,j,RBF(sdim,coords + i*sdim,coords + j*sdim,&fctx),INSERT_VALUES);
76: }
77: }
78: MatAssemblyBegin(Ad,MAT_FINAL_ASSEMBLY);
79: MatAssemblyEnd(Ad,MAT_FINAL_ASSEMBLY);
80: }
82: /* Create and assemble the matrix */
83: MatCreateH2OpusFromKernel(PETSC_COMM_WORLD,n,n,PETSC_DECIDE,PETSC_DECIDE,sdim,coords,PETSC_FALSE,RBF,&fctx,eta,leafsize,basisord,&A);
84: MatSetOption(A,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
85: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
86: PetscObjectSetName((PetscObject)A,"RBF");
87: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
88: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
89: MatViewFromOptions(A,NULL,"-rbf_view");
91: MatCreateVecs(A,&x,&b);
92: VecDuplicate(x,&u);
93: VecDuplicate(x,&d);
95: {
96: PetscReal norm;
97: MatComputeOperator(A,MATDENSE,&Ae);
98: MatAXPY(Ae,-1.0,Ad,SAME_NONZERO_PATTERN);
99: MatGetDiagonal(Ae,d);
100: MatViewFromOptions(Ae,NULL,"-A_view");
101: MatViewFromOptions(Ae,NULL,"-D_view");
102: MatNorm(Ae,NORM_FROBENIUS,&norm);
103: PetscPrintf(PETSC_COMM_WORLD,"Approx err %g\n",norm);
104: VecNorm(d,NORM_2,&norm);
105: PetscPrintf(PETSC_COMM_WORLD,"Approx err (diag) %g\n",norm);
106: MatDestroy(&Ae);
107: }
109: VecSet(u,1.0);
110: MatMult(Ad,u,b);
111: MatViewFromOptions(Ad,NULL,"-Ad_view");
112: KSPCreate(PETSC_COMM_WORLD,&ksp);
113: KSPSetOperators(ksp,Ad,A);
114: KSPGetPC(ksp,&pc);
115: PCSetType(pc,PCH2OPUS);
116: KSPSetFromOptions(ksp);
117: /* we can also pass the points coordinates
118: In this case it is not needed, since the preconditioning
119: matrix is of type H2OPUS */
120: PCSetCoordinates(pc,sdim,n,coords);
122: KSPSolve(ksp,b,x);
123: VecAXPY(x,-1.0,u);
124: VecNorm(x,NORM_2,&norm);
125: KSPGetIterationNumber(ksp,&its);
126: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
128: /* change lambda and reassemble */
129: VecSet(x,(scale-1.)*fctx.lambda);
130: MatDiagonalSet(Ad,x,ADD_VALUES);
131: fctx.lambda *= scale;
132: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
133: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
134: {
135: PetscReal norm;
136: MatComputeOperator(A,MATDENSE,&Ae);
137: MatAXPY(Ae,-1.0,Ad,SAME_NONZERO_PATTERN);
138: MatGetDiagonal(Ae,d);
139: MatViewFromOptions(Ae,NULL,"-A_view");
140: MatViewFromOptions(Ae,NULL,"-D_view");
141: MatNorm(Ae,NORM_FROBENIUS,&norm);
142: PetscPrintf(PETSC_COMM_WORLD,"Approx err %g\n",norm);
143: VecNorm(d,NORM_2,&norm);
144: PetscPrintf(PETSC_COMM_WORLD,"Approx err (diag) %g\n",norm);
145: MatDestroy(&Ae);
146: }
147: KSPSetOperators(ksp,Ad,A);
148: MatMult(Ad,u,b);
149: KSPSolve(ksp,b,x);
150: MatMult(Ad,x,u);
151: VecAXPY(u,-1.0,b);
152: VecNorm(u,NORM_2,&norm);
153: KSPGetIterationNumber(ksp,&its);
154: PetscPrintf(PETSC_COMM_WORLD,"Residual norm error %g, Iterations %D\n",(double)norm,its);
156: PetscFree(coords);
157: PetscFree(fctx.l);
158: PetscRandomDestroy(&r);
159: VecDestroy(&x);
160: VecDestroy(&u);
161: VecDestroy(&d);
162: VecDestroy(&b);
163: MatDestroy(&A);
164: MatDestroy(&Ad);
165: KSPDestroy(&ksp);
166: PetscFinalize();
167: return ierr;
168: }
170: /*TEST
172: build:
173: requires: h2opus
175: test:
176: requires: h2opus !single
177: suffix: 1
178: args: -ksp_error_if_not_converged -pc_h2opus_monitor
180: test:
181: requires: h2opus !single
182: suffix: 1_ns
183: output_file: output/ex21_1.out
184: args: -ksp_error_if_not_converged -pc_h2opus_monitor -pc_h2opus_hyperorder 2
186: test:
187: requires: h2opus !single
188: suffix: 2
189: args: -ksp_error_if_not_converged -pc_h2opus_monitor -pc_h2opus_hyperorder 4
191: TEST*/