Actual source code: ex76f.F90
1: !
2: ! Description: Solves a linear systems using PCHPDDM.
3: !
5: program main
6: #include <petsc/finclude/petscksp.h>
7: use petscksp
8: use petscisdef
9: implicit none
10: Vec x,b
11: Mat A,aux,Y,C
12: KSP ksp
13: PC pc
14: IS is,sizes
15: PetscScalar one
16: PetscInt, pointer :: idx(:)
17: PetscMPIInt rank,size
18: PetscInt m,N
19: PetscViewer viewer
20: character*(PETSC_MAX_PATH_LEN) dir,name
21: character*(8) fmt
22: character(1) crank,csize
23: PetscBool flg
24: PetscErrorCode ierr
26: PetscCallA(PetscInitialize(ierr))
28: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
29: N = 1
30: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-rhs',N,flg,ierr))
31: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
32: PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
33: PetscCallA(MatCreate(PETSC_COMM_SELF,aux,ierr))
34: PetscCallA(ISCreate(PETSC_COMM_SELF,is,ierr))
35: dir = '.'
36: PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-load_dir',dir,flg,ierr))
37: fmt = '(I1)'
38: write (crank,fmt) rank
39: write (csize,fmt) size
40: write (name,'(a)')trim(dir)//'/sizes_'//crank//'_'//csize//'.dat'
41: PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ, viewer,ierr))
42: PetscCallA(ISCreate(PETSC_COMM_SELF,sizes,ierr))
43: PetscCallA(ISLoad(sizes,viewer,ierr))
44: PetscCallA(ISGetIndicesF90(sizes,idx,ierr))
45: PetscCallA(MatSetSizes(A,idx(1),idx(2),idx(3),idx(4),ierr))
46: PetscCallA(ISRestoreIndicesF90(sizes,idx,ierr))
47: PetscCallA(ISDestroy(sizes,ierr))
48: PetscCallA(PetscViewerDestroy(viewer,ierr))
49: PetscCallA(MatSetUp(A,ierr))
50: write (name,'(a)')trim(dir)//'/A.dat'
51: PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,viewer,ierr))
52: PetscCallA(MatLoad(A,viewer,ierr))
53: PetscCallA(PetscViewerDestroy(viewer,ierr))
54: write (name,'(a)')trim(dir)//'/is_'//crank//'_'//csize//'.dat'
55: PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ,viewer,ierr))
56: PetscCallA(ISLoad(is,viewer,ierr))
57: PetscCallA(PetscViewerDestroy(viewer,ierr))
58: write (name,'(a)')trim(dir)//'/Neumann_'//crank//'_'//csize//'.dat'
59: PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ,viewer,ierr))
60: PetscCallA(MatSetBlockSizesFromMats(aux,A,A,ierr))
61: PetscCallA(MatLoad(aux,viewer,ierr))
62: PetscCallA(PetscViewerDestroy(viewer,ierr))
63: PetscCallA(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE,ierr))
64: PetscCallA(MatSetOption(aux,MAT_SYMMETRIC,PETSC_TRUE,ierr))
65: PetscCallA(KSPCreate(PETSC_COMM_WORLD,ksp,ierr))
66: PetscCallA(KSPSetOperators(ksp,A,A,ierr))
67: PetscCallA(KSPGetPC(ksp,pc,ierr))
68: PetscCallA(PCSetType(pc,PCHPDDM,ierr))
69: #if defined(PETSC_HAVE_HPDDM) && defined(PETSC_HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
70: PetscCallA(PCHPDDMSetAuxiliaryMat(pc,is,aux,PETSC_NULL_FUNCTION,PETSC_NULL_INTEGER,ierr))
71: PetscCallA(PCHPDDMHasNeumannMat(pc,PETSC_FALSE,ierr))
72: #endif
73: PetscCallA(ISDestroy(is,ierr))
74: PetscCallA(MatDestroy(aux,ierr))
75: PetscCallA(KSPSetFromOptions(ksp,ierr))
76: PetscCallA(MatCreateVecs(A,x,b,ierr))
77: one = 1.0
78: PetscCallA(VecSet(b,one,ierr))
79: PetscCallA(KSPSolve(ksp,b,x,ierr))
80: PetscCallA(VecGetLocalSize(x,m,ierr))
81: PetscCallA(VecDestroy(x,ierr))
82: PetscCallA(VecDestroy(b,ierr))
83: if (N .gt. 1) then
84: PetscCallA(PetscOptionsClearValue(PETSC_NULL_OPTIONS,'-ksp_converged_reason',ierr))
85: PetscCallA(KSPSetFromOptions(ksp,ierr))
86: PetscCallA(MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,PETSC_NULL_SCALAR,C,ierr))
87: PetscCallA(MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,PETSC_NULL_SCALAR,Y,ierr))
88: PetscCallA(MatSetRandom(C,PETSC_NULL_RANDOM,ierr))
89: PetscCallA(KSPMatSolve(ksp,C,Y,ierr))
90: PetscCallA(MatDestroy(Y,ierr))
91: PetscCallA(MatDestroy(C,ierr))
92: endif
93: PetscCallA(KSPDestroy(ksp,ierr))
94: PetscCallA(MatDestroy(A,ierr))
95: PetscCallA(PetscFinalize(ierr))
96: end
98: !/*TEST
99: !
100: ! test:
101: ! requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
102: ! output_file: output/ex76_1.out
103: ! nsize: 4
104: ! args: -ksp_rtol 1e-3 -ksp_converged_reason -pc_type {{bjacobi hpddm}shared output} -pc_hpddm_coarse_sub_pc_type lu -sub_pc_type lu -options_left no -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO
105: !
106: ! test:
107: ! requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
108: ! suffix: geneo
109: ! output_file: output/ex76_geneo_pc_hpddm_levels_1_eps_nev-5.out
110: ! nsize: 4
111: ! args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev 5 -pc_hpddm_levels_1_st_pc_type cholesky -pc_hpddm_coarse_p {{1 2}shared output} -pc_hpddm_coarse_pc_type redundant -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO
112: !
113: ! test:
114: ! requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
115: ! suffix: fgmres_geneo_20_p_2
116: ! output_file: output/ex76_fgmres_geneo_20_p_2.out
117: ! nsize: 4
118: ! args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_pc_type redundant -ksp_type fgmres -pc_hpddm_coarse_mat_type {{baij sbaij}shared output} -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO
119: !
120: ! test:
121: ! requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
122: ! suffix: fgmres_geneo_20_p_2_geneo
123: ! output_file: output/ex76_fgmres_geneo_20_p_2.out
124: ! nsize: 4
125: ! args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_2_p 2 -pc_hpddm_levels_2_mat_type {{baij sbaij}shared output} -pc_hpddm_levels_2_eps_nev {{5 20}shared output} -pc_hpddm_levels_2_sub_pc_type cholesky -pc_hpddm_levels_2_ksp_type gmres -ksp_type fgmres -pc_hpddm_coarse_mat_type {{baij sbaij}shared output} -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO
126: ! # PCHPDDM + KSPHPDDM test to exercise multilevel + multiple RHS in one go
127: ! test:
128: ! requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
129: ! suffix: fgmres_geneo_20_p_2_geneo_rhs
130: ! output_file: output/ex76_fgmres_geneo_20_p_2.out
131: ! nsize: 4
132: ! args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_2_p 2 -pc_hpddm_levels_2_mat_type baij -pc_hpddm_levels_2_eps_nev 5 -pc_hpddm_levels_2_sub_pc_type cholesky -pc_hpddm_levels_2_ksp_max_it 10 -pc_hpddm_levels_2_ksp_type hpddm -pc_hpddm_levels_2_ksp_hpddm_type gmres -ksp_type hpddm -ksp_hpddm_variant flexible -pc_hpddm_coarse_mat_type baij -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO -rhs 4
133: !
134: !TEST*/