Actual source code: ex21.c

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