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;

 38:   PetscInitialize(&argc, &args, (char *)0, help);
 39:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 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++) MatSetValue(Ad, i, j, RBF(sdim, coords + i * sdim, coords + j * sdim, &fctx), INSERT_VALUES);
 75:     }
 76:     MatAssemblyBegin(Ad, MAT_FINAL_ASSEMBLY);
 77:     MatAssemblyEnd(Ad, MAT_FINAL_ASSEMBLY);
 78:   }

 80:   /* Create and assemble the matrix */
 81:   MatCreateH2OpusFromKernel(PETSC_COMM_WORLD, n, n, PETSC_DECIDE, PETSC_DECIDE, sdim, coords, PETSC_FALSE, RBF, &fctx, eta, leafsize, basisord, &A);
 82:   MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
 83:   MatSetOption(A, MAT_SYMMETRY_ETERNAL, PETSC_TRUE);
 84:   PetscObjectSetName((PetscObject)A, "RBF");
 85:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 86:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 87:   MatViewFromOptions(A, NULL, "-rbf_view");

 89:   MatCreateVecs(A, &x, &b);
 90:   VecDuplicate(x, &u);
 91:   VecDuplicate(x, &d);

 93:   {
 94:     PetscReal norm;
 95:     MatComputeOperator(A, MATDENSE, &Ae);
 96:     MatAXPY(Ae, -1.0, Ad, SAME_NONZERO_PATTERN);
 97:     MatGetDiagonal(Ae, d);
 98:     MatViewFromOptions(Ae, NULL, "-A_view");
 99:     MatViewFromOptions(Ae, NULL, "-D_view");
100:     MatNorm(Ae, NORM_FROBENIUS, &norm);
101:     PetscPrintf(PETSC_COMM_WORLD, "Approx err %g\n", norm);
102:     VecNorm(d, NORM_2, &norm);
103:     PetscPrintf(PETSC_COMM_WORLD, "Approx err (diag) %g\n", norm);
104:     MatDestroy(&Ae);
105:   }

107:   VecSet(u, 1.0);
108:   MatMult(Ad, u, b);
109:   MatViewFromOptions(Ad, NULL, "-Ad_view");
110:   KSPCreate(PETSC_COMM_WORLD, &ksp);
111:   KSPSetOperators(ksp, Ad, A);
112:   KSPGetPC(ksp, &pc);
113:   PCSetType(pc, PCH2OPUS);
114:   KSPSetFromOptions(ksp);
115:   /* we can also pass the points coordinates
116:      In this case it is not needed, since the preconditioning
117:      matrix is of type H2OPUS */
118:   PCSetCoordinates(pc, sdim, n, coords);

120:   KSPSolve(ksp, b, x);
121:   VecAXPY(x, -1.0, u);
122:   VecNorm(x, NORM_2, &norm);
123:   KSPGetIterationNumber(ksp, &its);
124:   PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);

126:   /* change lambda and reassemble */
127:   VecSet(x, (scale - 1.) * fctx.lambda);
128:   MatDiagonalSet(Ad, x, ADD_VALUES);
129:   fctx.lambda *= scale;
130:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
131:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
132:   {
133:     PetscReal norm;
134:     MatComputeOperator(A, MATDENSE, &Ae);
135:     MatAXPY(Ae, -1.0, Ad, SAME_NONZERO_PATTERN);
136:     MatGetDiagonal(Ae, d);
137:     MatViewFromOptions(Ae, NULL, "-A_view");
138:     MatViewFromOptions(Ae, NULL, "-D_view");
139:     MatNorm(Ae, NORM_FROBENIUS, &norm);
140:     PetscPrintf(PETSC_COMM_WORLD, "Approx err %g\n", norm);
141:     VecNorm(d, NORM_2, &norm);
142:     PetscPrintf(PETSC_COMM_WORLD, "Approx err (diag) %g\n", norm);
143:     MatDestroy(&Ae);
144:   }
145:   KSPSetOperators(ksp, Ad, A);
146:   MatMult(Ad, u, b);
147:   KSPSolve(ksp, b, x);
148:   MatMult(Ad, x, u);
149:   VecAXPY(u, -1.0, b);
150:   VecNorm(u, NORM_2, &norm);
151:   KSPGetIterationNumber(ksp, &its);
152:   PetscPrintf(PETSC_COMM_WORLD, "Residual norm error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);

154:   PetscFree(coords);
155:   PetscFree(fctx.l);
156:   PetscRandomDestroy(&r);
157:   VecDestroy(&x);
158:   VecDestroy(&u);
159:   VecDestroy(&d);
160:   VecDestroy(&b);
161:   MatDestroy(&A);
162:   MatDestroy(&Ad);
163:   KSPDestroy(&ksp);
164:   PetscFinalize();
165:   return 0;
166: }

168: /*TEST

170:   build:
171:     requires: h2opus

173:   test:
174:     requires: h2opus !single
175:     suffix: 1
176:     args: -ksp_error_if_not_converged -pc_h2opus_monitor

178:   test:
179:     requires: h2opus !single
180:     suffix: 1_ns
181:     output_file: output/ex21_1.out
182:     args: -ksp_error_if_not_converged -pc_h2opus_monitor -pc_h2opus_hyperorder 2

184:   test:
185:     requires: h2opus !single
186:     suffix: 2
187:     args: -ksp_error_if_not_converged -pc_h2opus_monitor -pc_h2opus_hyperorder 4

189: TEST*/