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