1: !
  2: !   Description: Solves a series of linear systems using KSPHPDDM.
  3: !

  5:       program main
  6: #include <petsc/finclude/petscksp.h>
  7:       use petscksp
  8:       implicit none
  9:       Vec                            x,b
 10:       Mat                            A
 11: #if defined(PETSC_HAVE_HPDDM)
 12:       Mat                            U
 13: #endif
 14:       KSP                            ksp
 15:       PetscInt                       i,j,nmat
 16:       PetscViewer                    viewer
 17:       character*(PETSC_MAX_PATH_LEN) dir,name
 18:       character*(8)                  fmt
 19:       character(3)                   cmat
 20:       PetscBool                      flg,reset
 21:       PetscErrorCode                 ierr

 23:       PetscCallA(PetscInitialize(ierr))
 24:       dir = '.'
 25:       PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-load_dir',dir,flg,ierr))
 26:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-nmat',nmat,flg,ierr))
 27:       reset = PETSC_FALSE
 28:       PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-reset',reset,flg,ierr))
 29:       PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
 30:       PetscCallA(KSPCreate(PETSC_COMM_WORLD,ksp,ierr))
 31:       PetscCallA(KSPSetOperators(ksp,A,A,ierr))
 32:       do 50 i=0,nmat-1
 33:         j = i+400
 34:         fmt = '(I3)'
 35:         write (cmat,fmt) j
 36:         write (name,'(a)')trim(dir)//'/A_'//cmat//'.dat'
 37:         PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,viewer,ierr))
 38:         PetscCallA(MatLoad(A,viewer,ierr))
 39:         PetscCallA(PetscViewerDestroy(viewer,ierr))
 40:         if (i .eq. 0) then
 41:           PetscCallA(MatCreateVecs(A,x,b,ierr))
 42:         endif
 43:         write (name,'(a)')trim(dir)//'/rhs_'//cmat//'.dat'
 44:         PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,viewer,ierr))
 45:         PetscCallA(VecLoad(b,viewer,ierr))
 46:         PetscCallA(PetscViewerDestroy(viewer,ierr))
 47:         PetscCallA(KSPSetFromOptions(ksp,ierr))
 48:         PetscCallA(KSPSolve(ksp,b,x,ierr))
 49:         PetscCallA(PetscObjectTypeCompare(ksp,KSPHPDDM,flg,ierr))
 50: #if defined(PETSC_HAVE_HPDDM)
 51:         if (flg .and. reset) then
 52:           PetscCallA(KSPHPDDMGetDeflationMat(ksp,U,ierr))
 53:           PetscCallA(KSPReset(ksp,ierr))
 54:           PetscCallA(KSPSetOperators(ksp,A,A,ierr))
 55:           PetscCallA(KSPSetFromOptions(ksp,ierr))
 56:           PetscCallA(KSPSetUp(ksp,ierr))
 57:           if (.not. PetscObjectIsNull(U)) then
 58:             PetscCallA(KSPHPDDMSetDeflationMat(ksp,U,ierr))
 59:             PetscCallA(MatDestroy(U,ierr))
 60:           endif
 61:         endif
 62: #endif
 63:   50  continue
 64:       PetscCallA(VecDestroy(x,ierr))
 65:       PetscCallA(VecDestroy(b,ierr))
 66:       PetscCallA(MatDestroy(A,ierr))
 67:       PetscCallA(KSPDestroy(ksp,ierr))
 68:       PetscCallA(PetscFinalize(ierr))
 69:       end

 71: !/*TEST
 72: !
 73: !   test:
 74: !      suffix: 1
 75: !      output_file: output/ex75_1.out
 76: !      nsize: 1
 77: !      requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
 78: !      args: -nmat 1 -pc_type none -ksp_converged_reason -ksp_type {{gmres hpddm}shared output} -ksp_max_it 1000 -ksp_gmres_restart 1000 -ksp_rtol 1e-10 -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR
 79: !
 80: !   test:
 81: !      requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
 82: !      suffix: 1_icc
 83: !      output_file: output/ex75_1_icc.out
 84: !      nsize: 1
 85: !      args: -nmat 1 -pc_type icc -ksp_converged_reason -ksp_type {{gmres hpddm}shared output} -ksp_max_it 1000 -ksp_gmres_restart 1000 -ksp_rtol 1e-10 -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR
 86: !
 87: !   testset:
 88: !      requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
 89: !      args: -nmat 3 -pc_type none -ksp_converged_reason -ksp_type hpddm -ksp_max_it 1000 -ksp_gmres_restart 40 -ksp_rtol 1e-10 -ksp_hpddm_type gcrodr -ksp_hpddm_recycle 20 -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR
 90: !      test:
 91: !        nsize: 1
 92: !        suffix: 2_seq
 93: !        output_file: output/ex75_2.out
 94: !      test:
 95: !        nsize: 2
 96: !        suffix: 2_par
 97: !        output_file: output/ex75_2.out
 98: !
 99: !   testset:
100: !      requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
101: !      output_file: output/ex75_2_icc.out
102: !      nsize: 1
103: !      args: -nmat 3 -pc_type icc -ksp_converged_reason -ksp_type hpddm -ksp_max_it 1000 -ksp_gmres_restart 40 -ksp_rtol 1e-10 -ksp_hpddm_type gcrodr -ksp_hpddm_recycle 20 -reset {{false true}shared output} -load_dir ${DATAFILESPATH}/matrices/hpddm/GCRODR
104: !      test:
105: !        suffix: 2_icc
106: !        args:
107: !      test:
108: !        suffix: 2_icc_atol
109: !        output_file: output/ex75_2_icc_atol.out
110: !        args: -ksp_atol 1e-12
111: !
112: !TEST*/