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,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*/