Actual source code: ex21.c

  1: static char help[] = "Solves a RBF kernel matrix with KSP and PCH2OPUS.\n\n";

  3: #include <petscksp.h>

  5: typedef struct {
  6:   PetscReal  sigma;
  7:   PetscReal *l;
  8:   PetscReal  lambda;
  9: } RBFCtx;

 11: static PetscScalar RBF(PetscInt sdim, PetscReal x[], PetscReal y[], void *ctx)
 12: {
 13:   RBFCtx    *rbfctx = (RBFCtx *)ctx;
 14:   PetscInt   d;
 15:   PetscReal  diff   = 0.0;
 16:   PetscReal  s      = rbfctx->sigma;
 17:   PetscReal *l      = rbfctx->l;
 18:   PetscReal  lambda = rbfctx->lambda;

 20:   for (d = 0; d < sdim; d++) diff += (x[d] - y[d]) * (x[d] - y[d]) / (l[d] * l[d]);
 21:   return s * s * PetscExpReal(-0.5 * diff) + (diff != 0.0 ? 0.0 : lambda);
 22: }

 24: int main(int argc, char **args)
 25: {
 26:   Vec         x, b, u, d;
 27:   Mat         A, Ae = NULL, Ad = NULL;
 28:   KSP         ksp;
 29:   PetscRandom r;
 30:   PC          pc;
 31:   PetscReal   norm, *coords, eta, scale = 0.5;
 32:   PetscInt    basisord, leafsize, sdim, n, its, i;
 33:   PetscMPIInt size;
 34:   RBFCtx      fctx;

 36:   PetscFunctionBeginUser;
 37:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 38:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 39:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
 40:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &r));
 41:   PetscCall(PetscRandomSetFromOptions(r));

 43:   sdim = 2;
 44:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-sdim", &sdim, NULL));
 45:   n = 32;
 46:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 47:   eta = 0.6;
 48:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-eta", &eta, NULL));
 49:   leafsize = 32;
 50:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-leafsize", &leafsize, NULL));
 51:   basisord = 8;
 52:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-basisord", &basisord, NULL));

 54:   /* Create random points */
 55:   PetscCall(PetscMalloc1(sdim * n, &coords));
 56:   PetscCall(PetscRandomGetValuesReal(r, sdim * n, coords));

 58:   fctx.lambda = 0.01;
 59:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-lambda", &fctx.lambda, NULL));
 60:   PetscCall(PetscRandomGetValueReal(r, &fctx.sigma));
 61:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma", &fctx.sigma, NULL));
 62:   PetscCall(PetscMalloc1(sdim, &fctx.l));
 63:   PetscCall(PetscRandomGetValuesReal(r, sdim, fctx.l));
 64:   PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-l", fctx.l, (i = sdim, &i), NULL));
 65:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-scale", &scale, NULL));

 67:   /* Populate dense matrix for comparisons */
 68:   {
 69:     PetscInt i, j;

 71:     PetscCall(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++) PetscCall(MatSetValue(Ad, i, j, RBF(sdim, coords + i * sdim, coords + j * sdim, &fctx), INSERT_VALUES));
 74:     }
 75:     PetscCall(MatAssemblyBegin(Ad, MAT_FINAL_ASSEMBLY));
 76:     PetscCall(MatAssemblyEnd(Ad, MAT_FINAL_ASSEMBLY));
 77:   }

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

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

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

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

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

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

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

167: /*TEST

169:   build:
170:     requires: h2opus

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

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

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

188: TEST*/