2: static char help[] = "Solves a RBF kernel matrix with KSP and PCHARA.\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;
 28:   Mat            A,Ae = NULL, Ad = NULL;
 29:   KSP            ksp;
 30:   PetscRandom    r;
 31:   PC             pc;
 32:   PetscReal      norm,*coords,eta;
 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);
 54:   PetscMalloc1(sdim*n,&coords);
 55:   for (i=0;i<sdim*n;i++) {
 56:     PetscRandomGetValueReal(r,coords + i);
 57:   }
 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:   for (i=0;i<sdim;i++) {
 64:     PetscRandomGetValueReal(r,&fctx.l[i]);
 65:   }
 66:   PetscOptionsGetRealArray(NULL,NULL,"-l",fctx.l,(i=sdim,&i),NULL);
 67:   {
 68:     PetscInt i,j;

 70:     MatCreateDense(PETSC_COMM_WORLD,n,n,PETSC_DECIDE,PETSC_DECIDE,NULL,&Ad);
 71:     for (i = 0; i < n; i++) {
 72:       for (j = 0; j < n; j++) {
 73:         MatSetValue(Ad,i,j,RBF(sdim,coords + i*sdim,coords + j*sdim,&fctx),INSERT_VALUES);
 74:       }
 75:     }
 76:     MatAssemblyBegin(Ad,MAT_FINAL_ASSEMBLY);
 77:     MatAssemblyEnd(Ad,MAT_FINAL_ASSEMBLY);
 78:   }
 79:   MatCreateHaraFromKernel(PETSC_COMM_WORLD,n,n,PETSC_DECIDE,PETSC_DECIDE,sdim,coords,RBF,&fctx,eta,leafsize,basisord,&A);
 80:   MatSetOption(A,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
 81:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
 82:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 83:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 84:   {
 85:     PetscReal norm;
 86:     MatComputeOperator(A,MATDENSE,&Ae);
 87:     MatAXPY(Ae,-1.0,Ad,SAME_NONZERO_PATTERN);
 88:     MatViewFromOptions(Ae,NULL,"-A_view");
 89:     MatNorm(Ae,NORM_FROBENIUS,&norm);
 90:     PetscPrintf(PETSC_COMM_WORLD,"Approx err %g\n",norm);
 91:   }
 92:   MatCreateVecs(A,&x,&b);
 93:   VecDuplicate(x,&u);

 95:   VecSet(u,1.0);
 96:   MatMult(A,u,b);
 97:   MatViewFromOptions(Ad,NULL,"-Ad_view");
 98:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 99:   KSPSetOperators(ksp,Ad ? Ad : A,A);
100:   KSPGetPC(ksp,&pc);
101:   PCSetType(pc,PCHARA);
102:   KSPSetFromOptions(ksp);
103:   KSPGetPC(ksp,&pc);
104:   PCSetCoordinates(pc,sdim,n,coords);
105:   KSPSolve(ksp,b,x);
106:   VecAXPY(x,-1.0,u);
107:   VecNorm(x,NORM_2,&norm);
108:   KSPGetIterationNumber(ksp,&its);
109:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);

111:   PetscFree(coords);
112:   PetscFree(fctx.l);
113:   PetscRandomDestroy(&r);
114:   VecDestroy(&x);
115:   VecDestroy(&u);
116:   VecDestroy(&b);
117:   MatDestroy(&A);
118:   MatDestroy(&Ae);
119:   MatDestroy(&Ad);
120:   KSPDestroy(&ksp);
121:   PetscFinalize();
122:   return ierr;
123: }

125: /*TEST

127:   build:
128:     requires: hara

130:   test:
131:     requires: hara
132:     suffix: 1
133:     args: -pc_hara_monitor

135:   test:
136:     requires: hara
137:     suffix: 1_ns
138:     output_file: output/ex21_1.out
139:     args: -pc_hara_monitor -pc_hara_hyperorder 2

141:   test:
142:     requires: hara
143:     suffix: 2
144:     args: -pc_hara_monitor -pc_hara_hyperorder 4

146: TEST*/