Actual source code: ex76f.F90
petsc-3.14.6 2021-03-30
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: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
27: if (ierr .ne. 0) then
28: print *,'Unable to initialize PETSc'
29: stop
30: endif
31: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
32: if (size .ne. 4) then
33: print *,'This example requires 4 processes'
34: stop
35: endif
36: N = 1
37: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-rhs',N,flg,ierr);CHKERRA(ierr)
38: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
39: call MatCreate(PETSC_COMM_WORLD,A,ierr);CHKERRA(ierr)
40: call MatCreate(PETSC_COMM_SELF,aux,ierr);CHKERRA(ierr)
41: call ISCreate(PETSC_COMM_SELF,is,ierr);CHKERRA(ierr)
42: dir = '.'
43: call PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-load_dir',dir,flg,ierr);CHKERRA(ierr)
44: fmt = '(I1)'
45: write (crank,fmt) rank
46: write (csize,fmt) size
47: write (name,'(a)')trim(dir)//'/sizes_'//crank//'_'//csize//'.dat'
48: call PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ, viewer,ierr);CHKERRA(ierr)
49: call ISCreate(PETSC_COMM_SELF,sizes,ierr);CHKERRA(ierr)
50: call ISLoad(sizes,viewer,ierr);CHKERRA(ierr)
51: call ISGetIndicesF90(sizes,idx,ierr);CHKERRA(ierr)
52: call MatSetSizes(A,idx(1),idx(2),idx(3),idx(4),ierr);CHKERRA(ierr)
53: call ISRestoreIndicesF90(sizes,idx,ierr);CHKERRA(ierr)
54: call ISDestroy(sizes,ierr);CHKERRA(ierr)
55: call PetscViewerDestroy(viewer,ierr);CHKERRA(ierr)
56: call MatSetUp(A,ierr);CHKERRA(ierr)
57: write (name,'(a)')trim(dir)//'/A.dat'
58: call PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,viewer,ierr);CHKERRA(ierr)
59: call MatLoad(A,viewer,ierr);CHKERRA(ierr)
60: call PetscViewerDestroy(viewer,ierr);CHKERRA(ierr)
61: write (name,'(a)')trim(dir)//'/is_'//crank//'_'//csize//'.dat'
62: call PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ,viewer,ierr);CHKERRA(ierr)
63: call ISLoad(is,viewer,ierr);CHKERRA(ierr)
64: call PetscViewerDestroy(viewer,ierr);CHKERRA(ierr)
65: write (name,'(a)')trim(dir)//'/Neumann_'//crank//'_'//csize//'.dat'
66: call PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ,viewer,ierr);CHKERRA(ierr)
67: call MatSetBlockSizesFromMats(aux,A,A,ierr);CHKERRA(ierr)
68: call MatLoad(aux,viewer,ierr);CHKERRA(ierr)
69: call PetscViewerDestroy(viewer,ierr);CHKERRA(ierr)
70: call MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE,ierr);CHKERRA(ierr)
71: call MatSetOption(aux,MAT_SYMMETRIC,PETSC_TRUE,ierr);CHKERRA(ierr)
72: call KSPCreate(PETSC_COMM_WORLD,ksp,ierr);CHKERRA(ierr)
73: call KSPSetOperators(ksp,A,A,ierr);CHKERRA(ierr)
74: call KSPGetPC(ksp,pc,ierr);CHKERRA(ierr)
75: call PCSetType(pc,PCHPDDM,ierr);CHKERRA(ierr)
76: #if defined(PETSC_HAVE_HPDDM)
77: call PCHPDDMSetAuxiliaryMat(pc,is,aux,PETSC_NULL_FUNCTION,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
78: call PCHPDDMHasNeumannMat(pc,PETSC_FALSE,ierr);CHKERRA(ierr)
79: #endif
80: call ISDestroy(is,ierr);CHKERRA(ierr)
81: call MatDestroy(aux,ierr);CHKERRA(ierr)
82: call KSPSetFromOptions(ksp,ierr);CHKERRA(ierr)
83: call MatCreateVecs(A,x,b,ierr);CHKERRA(ierr)
84: one = 1.0
85: call VecSet(b,one,ierr);CHKERRA(ierr)
86: call KSPSolve(ksp,b,x,ierr);CHKERRA(ierr)
87: call VecGetLocalSize(x,m,ierr);CHKERRA(ierr)
88: call VecDestroy(x,ierr);CHKERRA(ierr)
89: call VecDestroy(b,ierr);CHKERRA(ierr)
90: if (N .gt. 1) then
91: call PetscOptionsClearValue(PETSC_NULL_OPTIONS,'-ksp_converged_reason',ierr);CHKERRA(ierr)
92: call KSPSetFromOptions(ksp,ierr);CHKERRA(ierr)
93: call MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,PETSC_NULL_SCALAR,C,ierr);CHKERRA(ierr)
94: call MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,PETSC_NULL_SCALAR,Y,ierr);CHKERRA(ierr)
95: call MatSetRandom(C,PETSC_NULL_RANDOM,ierr);CHKERRA(ierr)
96: call KSPMatSolve(ksp,C,Y,ierr);CHKERRA(ierr)
97: call MatDestroy(Y,ierr);CHKERRA(ierr)
98: call MatDestroy(C,ierr);CHKERRA(ierr)
99: endif
100: call KSPDestroy(ksp,ierr);CHKERRA(ierr)
101: call MatDestroy(A,ierr);CHKERRA(ierr)
102: call PetscFinalize(ierr)
103: end
105: !/*TEST
106: !
107: ! test:
108: ! requires: hpddm slepc datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
109: ! output_file: output/ex76_1.out
110: ! nsize: 4
111: ! 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
112: !
113: ! test:
114: ! requires: hpddm slepc datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
115: ! suffix: geneo
116: ! output_file: output/ex76_geneo_pc_hpddm_levels_1_eps_nev-5.out
117: ! nsize: 4
118: ! 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
119: !
120: ! test:
121: ! requires: hpddm slepc datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
122: ! suffix: fgmres_geneo_20_p_2
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 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
126: !
127: ! test:
128: ! requires: hpddm slepc datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
129: ! suffix: fgmres_geneo_20_p_2_geneo
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 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
133: ! # PCHPDDM + KSPHPDDM test to exercise multilevel + multiple RHS in one go
134: ! test:
135: ! requires: hpddm slepc datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
136: ! suffix: fgmres_geneo_20_p_2_geneo_rhs
137: ! output_file: output/ex76_fgmres_geneo_20_p_2.out
138: ! nsize: 4
139: ! 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
140: !
141: !TEST*/