Actual source code: ex82.c

      1: #include <petscksp.h>

      3: static char help[] = "Solves a linear system using PCHPDDM and MATHTOOL.\n\n";

      5: static PetscErrorCode GenEntries(PetscInt sdim, PetscInt M, PetscInt N, const PetscInt *J, const PetscInt *K, PetscScalar *ptr, PetscCtx ctx)
      6: {
      7:   PetscInt  d, j, k;
      8:   PetscReal diff = 0.0, *coords = (PetscReal *)(ctx);

     10: #if !PetscDefined(HAVE_OPENMP)
     11:   PetscFunctionBeginUser;
     12: #endif
     13:   for (j = 0; j < M; j++) {
     14:     for (k = 0; k < N; k++) {
     15:       diff = 0.0;
     16:       for (d = 0; d < sdim; d++) diff += (coords[J[j] * sdim + d] - coords[K[k] * sdim + d]) * (coords[J[j] * sdim + d] - coords[K[k] * sdim + d]);
     17:       ptr[j + M * k] = 1.0 / (1.0e-2 + PetscSqrtReal(diff));
     18:     }
     19:   }
     20: #if !PetscDefined(HAVE_OPENMP)
     21:   PetscFunctionReturn(PETSC_SUCCESS);
     22: #else
     23:   return PETSC_SUCCESS;
     24: #endif
     25: }

     27: int main(int argc, char **argv)
     28: {
     29:   KSP               ksp;
     30:   PC                pc;
     31:   Vec               b, x;
     32:   Mat               A;
     33:   PetscInt          m = 100, dim = 3, M, begin = 0, n = 0, overlap = 1;
     34:   PetscMPIInt       size;
     35:   PetscReal        *coords, *gcoords;
     36:   MatHtoolKernelFn *kernel = GenEntries;
     37:   PetscBool         flg, sym = PETSC_FALSE;
     38:   PetscRandom       rdm;

     40:   PetscFunctionBeginUser;
     41:   PetscCall(PetscInitialize(&argc, &argv, (char *)NULL, help));
     42:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m_local", &m, NULL));
     43:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-symmetric", &sym, NULL));
     44:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL));
     45:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-overlap", &overlap, NULL));
     46:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     47:   M = size * m;
     48:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
     49:   PetscCall(PetscMalloc1(m * dim, &coords));
     50:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
     51:   PetscCall(PetscRandomGetValuesReal(rdm, m * dim, coords));
     52:   PetscCall(PetscCalloc1(M * dim, &gcoords));
     53:   PetscCallMPI(MPI_Exscan(&m, &begin, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
     54:   PetscCall(PetscArraycpy(gcoords + begin * dim, coords, m * dim));
     55:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, gcoords, M * dim, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
     56:   PetscCall(MatCreateHtoolFromKernel(PETSC_COMM_WORLD, m, m, M, M, dim, coords, coords, kernel, gcoords, &A));
     57:   PetscCall(MatSetOption(A, MAT_SYMMETRIC, sym));
     58:   PetscCall(MatSetFromOptions(A));
     59:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     60:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
     61:   PetscCall(MatCreateVecs(A, &b, &x));
     62:   PetscCall(VecSetRandom(b, rdm));
     63:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
     64:   PetscCall(KSPSetOperators(ksp, A, A));
     65:   PetscCall(KSPSetFromOptions(ksp));
     66:   PetscCall(KSPGetPC(ksp, &pc));
     67:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCHPDDM, &flg));
     68:   if (flg) {
     69: #if defined(PETSC_HAVE_HPDDM) && defined(PETSC_HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
     70:     Mat aux;
     71:     IS  is;
     72:     PetscCall(MatGetOwnershipRange(A, &begin, &n));
     73:     n -= begin;
     74:     PetscCall(ISCreateStride(PETSC_COMM_SELF, n, begin, 1, &is));
     75:     PetscCall(MatIncreaseOverlap(A, 1, &is, overlap));
     76:     PetscCall(ISGetLocalSize(is, &n));
     77:     PetscCall(MatCreateDense(PETSC_COMM_SELF, n, n, n, n, NULL, &aux));
     78:     PetscCall(MatSetOption(aux, MAT_SYMMETRIC, sym));
     79:     PetscCall(MatAssemblyBegin(aux, MAT_FINAL_ASSEMBLY));
     80:     PetscCall(MatAssemblyEnd(aux, MAT_FINAL_ASSEMBLY));
     81:     PetscCall(MatShift(aux, 1.0)); /* just the local identity matrix, not very meaningful numerically, but just testing that the necessary plumbing is there */
     82:     PetscCall(PCHPDDMSetAuxiliaryMat(pc, is, aux, NULL, NULL));
     83:     PetscCall(ISDestroy(&is));
     84:     PetscCall(MatDestroy(&aux));
     85: #endif
     86:   }
     87:   PetscCall(KSPSolve(ksp, b, x));
     88:   PetscCall(KSPDestroy(&ksp));
     89:   PetscCall(PetscRandomDestroy(&rdm));
     90:   PetscCall(VecDestroy(&b));
     91:   PetscCall(VecDestroy(&x));
     92:   PetscCall(MatDestroy(&A));
     93:   PetscCall(PetscFree(gcoords));
     94:   PetscCall(PetscFree(coords));
     95:   PetscCall(PetscFinalize());
     96:   return 0;
     97: }

     99: /*TEST

    101:    build:
    102:       requires: htool hpddm

    104:    test:
    105:       requires: htool hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
    106:       nsize: 4
    107:       # different numbers of iterations depending on PetscScalar type
    108:       filter: sed -e "s/symmetry: S/symmetry: N/g" -e "/number of dense/d" -e "s/Linear solve converged due to CONVERGED_RTOL iterations 13/Linear solve converged due to CONVERGED_RTOL iterations 18/g"
    109:       args: -ksp_view -ksp_converged_reason -mat_htool_epsilon 1e-2 -m_local 200 -pc_type hpddm -pc_hpddm_define_subdomains -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 1 -pc_hpddm_coarse_pc_type lu -pc_hpddm_levels_1_eps_gen_non_hermitian -symmetric {{false true}shared output} -overlap 2
    110:       output_file: output/ex82_1.out

    112: TEST*/