1!
2!   Description: Solves a linear system with a block of right-hand sides using KSPHPDDM.
3!
4
5      program main
6#include <petsc/finclude/petscksp.h>
7      use petscksp
8      implicit none
9      Mat                            X,B
10      Mat                            A
11      KSP                            ksp
12      PC                             pc
13      Mat                            F
14      PetscScalar                    alpha
15      PetscReal                      norm
16      PetscInt                       m,K
17      PetscViewer                    viewer
18      character*(PETSC_MAX_PATH_LEN) name
19      PetscBool                      flg
20      PetscErrorCode                 ierr
21
22      call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
23      if (ierr .ne. 0) then
24        print *,'Unable to initialize PETSc'
25        stop
26      endif
27      call PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-f',name,flg,ierr);CHKERRA(ierr)
28      if (flg .eqv. PETSC_FALSE) then
29        SETERRA(PETSC_COMM_WORLD,PETSC_ERR_SUP,'Must provide a binary file for the matrix')
30      endif
31      K = 5
32      call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',K,flg,ierr);CHKERRA(ierr)
33      call MatCreate(PETSC_COMM_WORLD,A,ierr);CHKERRA(ierr)
34      call KSPCreate(PETSC_COMM_WORLD,ksp,ierr);CHKERRA(ierr)
35      call KSPSetOperators(ksp,A,A,ierr);CHKERRA(ierr)
36      call PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,viewer,ierr);CHKERRA(ierr)
37      call MatLoad(A,viewer,ierr);CHKERRA(ierr)
38      call PetscViewerDestroy(viewer,ierr);CHKERRA(ierr)
39      call MatGetLocalSize(A,m,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
40      call MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,K,PETSC_NULL_SCALAR,B,ierr);CHKERRA(ierr)
41      call MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,K,PETSC_NULL_SCALAR,X,ierr);CHKERRA(ierr)
42      call MatSetRandom(B,PETSC_NULL_RANDOM,ierr);CHKERRA(ierr)
43      call KSPSetFromOptions(ksp,ierr);CHKERRA(ierr)
44      call KSPSetUp(ksp,ierr);CHKERRA(ierr)
45      call KSPMatSolve(ksp,B,X,ierr);CHKERRA(ierr)
46      call KSPGetMatSolveBlockSize(ksp,M,ierr);CHKERRA(ierr)
47      if (M .ne. PETSC_DECIDE) then
48        call KSPSetMatSolveBlockSize(ksp,PETSC_DECIDE,ierr);CHKERRA(ierr)
49        call MatZeroEntries(X,ierr);CHKERRA(ierr)
50        call KSPMatSolve(ksp,B,X,ierr);CHKERRA(ierr)
51      endif
52      call KSPGetPC(ksp,pc,ierr);CHKERRA(ierr)
53      call PetscObjectTypeCompare(pc,PCLU,flg,ierr);CHKERRA(ierr)
54      if (flg) then
55        call PCFactorGetMatrix(pc,F,ierr);CHKERRA(ierr)
56        call MatMatSolve(F,B,B,ierr);CHKERRA(ierr)
57        alpha = -1.0
58        call MatAYPX(B,alpha,X,SAME_NONZERO_PATTERN,ierr);CHKERRA(ierr)
59        call MatNorm(B,NORM_INFINITY,norm,ierr);CHKERRA(ierr)
60        if (norm > 100*PETSC_MACHINE_EPSILON) then
61          SETERRA(PETSC_COMM_WORLD,PETSC_ERR_PLIB,'KSPMatSolve() and MatMatSolve() difference has nonzero norm')
62        endif
63      endif
64      call MatDestroy(X,ierr);CHKERRA(ierr)
65      call MatDestroy(B,ierr);CHKERRA(ierr)
66      call MatDestroy(A,ierr);CHKERRA(ierr)
67      call KSPDestroy(ksp,ierr);CHKERRA(ierr)
68      call PetscFinalize(ierr)
69      end
70
71!/*TEST
72!
73!   testset:
74!      nsize: 2
75!      requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
76!      args: -ksp_converged_reason -ksp_max_it 1000 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat
77!      test:
78!         suffix: 1
79!         output_file: output/ex77_1.out
80!         args:
81!      test:
82!         suffix: 2a
83!         requires: hpddm
84!         output_file: output/ex77_2_ksp_hpddm_type-gmres.out
85!         args: -ksp_type hpddm -pc_type asm -ksp_hpddm_type gmres
86!      test:
87!         suffix: 2b
88!         requires: hpddm
89!         output_file: output/ex77_2_ksp_hpddm_type-bgmres.out
90!         args: -ksp_type hpddm -pc_type asm -ksp_hpddm_type bgmres
91!      test:
92!         suffix: 3a
93!         requires: hpddm
94!         output_file: output/ex77_3_ksp_hpddm_type-gcrodr.out
95!         args: -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type gcrodr
96!      test:
97!         suffix: 3b
98!         requires: hpddm
99!         output_file: output/ex77_3_ksp_hpddm_type-bgcrodr.out
100!         args: -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr
101!      test:
102!         nsize: 4
103!         suffix: 4
104!         requires: hpddm
105!         output_file: output/ex77_4.out
106!         args: -ksp_rtol 1e-4 -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr -ksp_view_final_residual -N 12 -ksp_matsolve_block_size 5
107!   test:
108!      nsize: 1
109!      suffix: preonly
110!      requires: hpddm datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
111!      output_file: output/ex77_preonly.out
112!      args: -N 6 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -pc_type lu -ksp_type hpddm -ksp_hpddm_type preonly
113!   test:
114!      nsize: 4
115!      suffix: 4_slepc
116!      requires: hpddm datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) slepc
117!      output_file: output/ex77_4.out
118!      filter: sed "/^ksp_hpddm_recycle_ Linear eigensolve converged/d"
119!      args: -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -ksp_rtol 1e-4 -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr -ksp_view_final_residual -N 12 -ksp_matsolve_block_size 5 -ksp_hpddm_recycle_redistribute 2 -ksp_hpddm_recycle_mat_type dense -ksp_hpddm_recycle_eps_converged_reason -ksp_hpddm_recycle_st_pc_type redundant
120!
121!TEST*/
122