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 == 0) then
 41:       PetscCallA(MatCreateVecs(A, x, b, ierr))
 42:     end if
 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:       end if
 61:     end if
 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*/