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,void *ctx)
  6: {
  7:   PetscInt  d,j,k;
  8:   PetscReal diff = 0.0,*coords = (PetscReal*)(ctx);

 11:   for (j = 0; j < M; j++) {
 12:     for (k = 0; k < N; k++) {
 13:       diff = 0.0;
 14:       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]);
 15:       ptr[j+M*k] = 1.0/(1.0e-2 + PetscSqrtReal(diff));
 16:     }
 17:   }
 18:   return 0;
 19: }

 21: int main(int argc,char **argv)
 22: {
 23:   KSP            ksp;
 24:   PC             pc;
 25:   Vec            b,x;
 26:   Mat            A;
 27:   PetscInt       m = 100,dim = 3,M,begin = 0,n = 0,overlap = 1;
 28:   PetscMPIInt    size;
 29:   PetscReal      *coords,*gcoords;
 30:   MatHtoolKernel kernel = GenEntries;
 31:   PetscBool      flg,sym = PETSC_FALSE;
 32:   PetscRandom    rdm;

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

 92: /*TEST

 94:    build:
 95:       requires: htool hpddm

 97:    test:
 98:       requires: htool hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
 99:       nsize: 4
100:       # different numbers of iterations depending on PetscScalar type
101:       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"
102:       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
103:       output_file: output/ex82_1.out

105: TEST*/