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;

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

 93: /*TEST

 95:    build:
 96:       requires: htool hpddm

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

106: TEST*/