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*/