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:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 24:       if (ierr .ne. 0) then
 25:         print *,'Unable to initialize PETSc'
 26:         stop
 27:       endif
 28:       dir = '.'
 29:       call PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-load_dir',dir,flg,ierr);CHKERRA(ierr)
 30:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-nmat',nmat,flg,ierr);CHKERRA(ierr)
 31:       reset = PETSC_FALSE
 32:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-reset',reset,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:       do 50 i=0,nmat-1
 37:         j = i+400
 38:         fmt = '(I3)'
 39:         write (cmat,fmt) j
 40:         write (name,'(a)')trim(dir)//'/A_'//cmat//'.dat'
 41:         call PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,viewer,ierr);CHKERRA(ierr)
 42:         call MatLoad(A,viewer,ierr);CHKERRA(ierr)
 43:         call PetscViewerDestroy(viewer,ierr);CHKERRA(ierr)
 44:         if (i .eq. 0) then
 45:           call MatCreateVecs(A,x,b,ierr);CHKERRA(ierr)
 46:         endif
 47:         write (name,'(a)')trim(dir)//'/rhs_'//cmat//'.dat'
 48:         call PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,viewer,ierr);CHKERRA(ierr)
 49:         call VecLoad(b,viewer,ierr);CHKERRA(ierr)
 50:         call PetscViewerDestroy(viewer,ierr);CHKERRA(ierr)
 51:         call KSPSetFromOptions(ksp,ierr);CHKERRA(ierr)
 52:         call KSPSolve(ksp,b,x,ierr);CHKERRA(ierr)
 53:         call PetscObjectTypeCompare(ksp,KSPHPDDM,flg,ierr);CHKERRA(ierr)
 54: #if defined(PETSC_HAVE_HPDDM)
 55:         if (flg .and. reset) then
 56:           call KSPHPDDMGetDeflationSpace(ksp,U,ierr);CHKERRA(ierr)
 57:           call KSPReset(ksp,ierr);CHKERRA(ierr)
 58:           call KSPSetOperators(ksp,A,A,ierr);CHKERRA(ierr)
 59:           call KSPSetFromOptions(ksp,ierr);CHKERRA(ierr)
 60:           call KSPSetUp(ksp,ierr);CHKERRA(ierr)
 61:           if (U .ne. PETSC_NULL_MAT) then
 62:             call KSPHPDDMSetDeflationSpace(ksp,U,ierr);CHKERRA(ierr)
 63:             call MatDestroy(U,ierr);CHKERRA(ierr)
 64:           endif
 65:         endif
 66: #endif
 67:   50  continue
 68:       call VecDestroy(x,ierr);CHKERRA(ierr)
 69:       call VecDestroy(b,ierr);CHKERRA(ierr)
 70:       call MatDestroy(A,ierr);CHKERRA(ierr)
 71:       call KSPDestroy(ksp,ierr);CHKERRA(ierr)
 72:       call PetscFinalize(ierr)
 73:       end

 75: !/*TEST
 76: !
 77: !   test:
 78: !      suffix: 1
 79: !      output_file: output/ex75_1.out
 80: !      nsize: 1
 81: !      requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
 82: !      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
 83: !
 84: !   test:
 85: !      requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
 86: !      suffix: 1_icc
 87: !      output_file: output/ex75_1_icc.out
 88: !      nsize: 1
 89: !      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
 90: !
 91: !   testset:
 92: !      requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
 93: !      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
 94: !      test:
 95: !        nsize: 1
 96: !        suffix: 2_seq
 97: !        output_file: output/ex75_2.out
 98: !      test:
 99: !        nsize: 2
100: !        suffix: 2_par
101: !        output_file: output/ex75_2.out
102: !
103: !   testset:
104: !      requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
105: !      output_file: output/ex75_2_icc.out
106: !      nsize: 1
107: !      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
108: !      test:
109: !        suffix: 2_icc
110: !        args:
111: !      test:
112: !        suffix: 2_icc_atol
113: !        output_file: output/ex75_2_icc_atol.out
114: !        args: -ksp_atol 1e-12
115: !
116: !TEST*/