eigen_mod.F90

Go to the documentation of this file.

module eigen_mod

  use choral_constants
  use real_type
  use choral_variables
  use abstract_interfaces, only: rntorn, rntornxl
  use basic_tools

  implicit none
  private



  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%
  !! TYPE
  public :: eigen

  !! GENERIC SUBROUTINES
  public :: clear
  public :: set, solve
  public :: print


  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%
  type eigen

     logical :: created = .false.

     !! parameters provided by the constructor 'eig = eigen(...)'
     !!
     integer  :: size

     integer  :: nev

     logical :: standard

     logical :: symmetric

     integer :: mode

     integer  :: which

     !! parameters provided by 'call set(eig, ...)'
     !!
     integer  :: type

     real(RP) :: tol 

     integer  :: itmax

     real(RP) :: shift

     integer  :: verb

     !! Parameter defined by the code with 'call set(eig, ...)'
     !!
     character(10) :: name

     !! Data computed with 'call solve(eig, ...)'
     !!
     real(RP), dimension(:), allocatable :: lambda

     real(RP), dimension(:,:), allocatable :: v

     !! Parameter returned by the code after solving
     !!
     integer  :: iter
     integer  :: aeval
     integer  :: beval
     integer  :: opeval
     integer  :: ierr

    contains
        final :: eigen_clear
  end type eigen

  !       %----------------------------------------%
  !       |                                        |
  !       |          GENERIC INTERFACES            |
  !       |                                        |
  !       %----------------------------------------%
  interface print
     module procedure eigen_print
  end interface print

  interface clear
     module procedure eigen_clear
  end interface clear

  !! constructor
  interface eigen
     module procedure eigen_create
  end interface eigen

  interface set
     module procedure eigen_set
  end interface set

  interface solve
     module procedure eigen_solve
  end interface solve

contains


  subroutine eigen_clear(eig)
    type(eigen), intent(inout) :: eig

    eig%created   = .false.
    call freemem(eig%lambda)
    call freemem(eig%v)

  end subroutine eigen_clear

  function eigen_create(size, nev, sym, std) result(eig)
    type(eigen) :: eig
    integer, intent(in) :: size, nev
    logical, intent(in) :: sym, std

    if (choral_verb>1) write(*,*) &
         & "eigen_mod       : eigen_create"

    !! checking size
    if (size<=0) call quit(&
         & "eigen_mod: eigen_create: 'size' <= 0")

    !! checking nev
    if ((nev<=0).OR.(nev>size)) call quit(&
         & "eigen_mod: eigen_create: wrong argument 'nev'")

    call clear(eig)

    eig%size      = size
    eig%nev       = nev
    eig%symmetric = sym
    eig%standard  = std

    ! set parameters to default ones
    eig%type    = eig_arpack
    eig%mode    = eig_regular 
    eig%shift   = 0.0_rp
    eig%tol     = 1e-8_rp
    eig%itMax   = 100
    eig%verb    = 0
    eig%name    = 'ARPACK'
    eig%which   = eig_which_sm

    eig%created = .true.

  end function eigen_create


  subroutine eigen_set(eig, type, mode, which, tol, itMax, &
       & shift, verb)
    type(eigen) , intent(inout)    :: eig
    integer , intent(in), optional :: type, mode, itMax, verb, which
    real(RP), intent(in), optional :: tol, shift

    if (.NOT.eig%created) call quit("eigen_mod: eigen_set:&
         & not created, 'eig=eigen()' missing")

    if (choral_verb>1) write(*,*) &
         & "eigen_mod       : eigen_set"

    if(present(type )) eig%type  = type     
    if(present(mode )) eig%mode  = mode
    if(present(tol  )) eig%tol   = tol
    if(present(itmax)) eig%itMax = itmax 
    if(present(shift)) eig%shift = shift
    if(present(verb )) eig%verb  = verb
    if(present(which)) eig%which = which

    !! checking which
    if(present(which)) then
       select case(eig%which)
       case(eig_which_sm, eig_which_lm)
       case default
          call quit("eigen_mod: eigen_set: 'which' = uncorrect")
       end select
    end if

    !! checking type
    if (present(type)) then
       select case(eig%type)
       case(eig_arpack)
          eig%name = 'ARPACK'
       case default
          call quit("eigen_mod: eigen_set:  'type' = uncorrect")
       end select
    end if

    !! checking mode
    if (present(mode)) then
       select case(mode)
       case(eig_regular, eig_inverse, eig_shift_invert)
       case default
          call quit("eigen_mod: eigen_set:  'mode' = uncorrect")
       end select
    end if

    !! checking tol
    if (present(tol)) then
       if (tol<=0.0_rp) call quit(&
         & "eigen_mod: eigen_set: 'tol' <= 0")
    end if

    !! checking itMax
    if (present(itmax)) then
       if (itmax<=0) call quit(&
         & "eigen_mod: eigen_set: 'itMax' <= 0")
    end if


  end subroutine eigen_set


  subroutine eigen_print (eig)
    type(eigen), intent(in) :: eig

    write(*,*)"eigen_mod       : eigen_print"
    write(*,*)"  Created                        =", eig%created
    if (.NOT.eig%created) return
    if (eig%standard) then
       write(*,*)"  Standard eigenproblem"
    else
       write(*,*)"  Generalised eigenproblem"
    end if
    if (eig%symmetric) then
       write(*,*)"  Symmetric case"
    else
       write(*,*)"  Non symmetric case"
    end if
    write(*,*)"  Solver Name                    = ", trim(eig%name)
    select case(eig%mode)
    case(eig_regular)
       write(*,*)"  Mode                           = regular"
    case(eig_inverse)
       write(*,*)"  Mode                           = inverse"
    case(eig_shift_invert)
       write(*,*)"  Mode                           = shift-invert"
       write(*,*)"  Shift factor                   =", eig%shift
    case default
       write(*,*)"  Mode                           = undefined"
    end select
    select case(eig%which)
    case(eig_which_sm)
       write(*,*)"  Which eigenvalues              = smallest magnitude"
    case(eig_which_lm)
       write(*,*)"  Which eigenvalues              = largest magnitude"
    case default
       write(*,*)"  Wich eigenvalues               = undefined"
    end select
    write(*,*)"  Number of eigenvalues/vectors  =", eig%nev    
    write(*,*)"  Problem size                   =", eig%size   
    write(*,*)"  Tolerance                      =", eig%tol    
    write(*,*)"  ItMax                          =", eig%itMax  
    write(*,*)"  Verbosity                      =", eig%verb

  end subroutine eigen_print



  subroutine eigen_solve(eig, A, B, Op)
    type(eigen), intent(inout)    :: eig
    procedure(rntorn)  , optional :: A, B
    procedure(rntornxl), optional :: Op

    real(DP) :: cpu

    if (.NOT.eig%created) call quit(&
         & "eigen: mode: eigen_solve: 'eig' not created yet,&
         & try 'eig = eigen(...)' first")

    if ( (eig%nev <= 0) .OR. (eig%nev > eig%size) ) call quit(&
         & "eigen: mode: eigen_solve: wrong eig%nev =&
         &  number of eigenvalues to be computed")


    if (eig%verb>0) write(*,*)"eigen_mod       :&
         & eigen_solve    = ", trim(eig%name)

    !! Check argument list
    !!
    if (eig%standard) then

       if (eig%verb>1) write(*,*) "eigen_mod       :&
            & eigen_solve    = standard eigen-problem"

       if (present(b)) call quit(&
            & "eigen_mod: eigen_solve:&
            &  standard eigen-problem: unconsistent arg 'B'")

       select case(eig%mode)
       case(eig_regular)

          if (eig%verb>1) write(*,*) "eigen_mod       :&
               & eigen_solve    = regular mode"

          if (.NOT.present(a)) call quit(&
               & "eigen_mod: eigen_solve:&
               & standard + EIG_REGULAR: operator 'A' missing")
          if (present(op)) call quit(&
               & "eigen_mod: eigen_solve:&
               & standard + EIG_REGULAR: unconsistent arg 'Op'")

       case(eig_inverse)
          call quit("eigen_mod: eigen_solve:&
               & standard + EIG_INVERSE not compatible")

       case(eig_shift_invert)
          if (eig%verb>1) write(*,*) "eigen_mod       :&
               & eigen_solve    = shift-invert mode"

          if (.NOT.present(op)) call quit(&
               & "eigen_mod: eigen_solve:&
               & standard + EIG_SHIFT_INVERT: arg 'Op' missing")
          if (present(a)) call quit(&
               & "eigen_mod: eigen_solve:&
               & standard + EIG_SHIFT_INVERT: unconsistent arg  'A'")
          
       end select
       
    else  !! generalised problem

       if (eig%verb>1) write(*,*) "eigen_mod       :&
            & eigen_solve    = generalised eigen-problem"

       select case(eig%mode)
       case(eig_regular)
          call quit("eigen_mod: eigen_solve:&
               & generalised + EIG_REGULAR: not compatible")
          
       case(eig_inverse)
          if (eig%verb>1) write(*,*) "eigen_mod       :&
               & eigen_solve    = inverse mode"

          if (.NOT.present(op)) call quit(&
               & "eigen_mod: eigen_solve:&
               & generalised + EIG_INVERSE: operator 'Op' missing")
          if (.NOT.present(b)) call quit(&
               & "eigen_mod: eigen_solve:&
               & generalised + EIG_INVERSE: operator 'B' missing")
          if (.NOT.present(a)) call quit(&
               & "eigen_mod: eigen_solve:&
               & generalised + EIG_INVERSE: operator 'A' missing")
          
       case(eig_shift_invert)
          if (eig%verb>1) write(*,*) "eigen_mod       :&
               & eigen_solve    = shift-invert mode"
          if (.NOT.present(op)) call quit(&
               & "eigen_mod: eigen_solve:&
               & generalised + EIG_SHIFT_INVERT: operator 'Op' missing")
          if (.NOT.present(b)) call quit(&
               & "eigen_mod: eigen_solve:&
               & generalised + EIG_SHIFT_INVERT: operator 'B' missing")
          if (present(a)) call quit(&
               & "eigen_mod: eigen_solve:&
               & generalised + EIG_SHIFT_INVERT: unconsistent arg  'A'")
          
       end select
    end if


    !! initialise 
    cpu       = clock()
    eig%ierr  = 0
    eig%AEval = 0
    eig%BEval = 0
    eig%OpEval = 0
    call allocmem(eig%lambda, eig%nev)
    call allocmem(eig%v, eig%size, eig%nev)

    select case(eig%type)
    case(eig_arpack)
#ifdef WARPACK
       call arpack_solve(eig, a, b, op)  
#else
       call quit("eigen_mod: solve_standard: set 'WITH_ARPACK=1' in the&
            & file 'CONFIG.in', re-configure (ARPACK must be installed)")
#endif      
    end select

    !! Display
    if (eig%verb>1) then
       if (eig%verb>1) write(*,*)"eigen_mod       :&
            & eigen_solve    = finished"
       write(*,*) "  Iterations                     =", eig%iter
       if (eig%AEval>0) then
          write(*,*) "  Nb of mat-vec products x-->A*x =", eig%AEval
       end if
       if (eig%BEval>0) then
          write(*,*) "  Nb of mat-vec products x-->B*x =", eig%BEval
       end if
       if (eig%OpEval>0) then
          write(*,*) "  Nb of linear system inversions =", eig%OpEval
       end if
       cpu = clock() - cpu
       write(*,*) '  CPU                            =', real(cpu, sp)
       select case(eig%ierr)
       case(0)
          write(*,*) "  Converged"
       case default
          write(*,*) "  Not converged, ierr            =", eig%ierr
       end select
    end if

  end subroutine eigen_solve


  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !!           WRAPPER TO ARPACK
  !!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
#ifdef WARPACK


  subroutine arpack_solve(eig, A, B, Op)
    type(eigen), intent(inout)    :: eig
    procedure(rntorn)  , optional :: A, B
    procedure(rntornxl), optional :: Op

    !! double precision only
    !!
    if (rp/=dp) call quit("eigen_mod:&
         & arpack_solve: double precision only")
    
    if (eig%symmetric) then
       if (eig%verb>1) write(*,*)"eigen_mod       :&
            & arpack_solve   = symmetric"
       call arpack_ds(eig, a, b, op)
    else
       call quit("eigen_mod: arpack_solve:&
            & non-symmetric case not implemented")
    end if

  end subroutine arpack_solve


  subroutine arpack_ds(eig, A, B, Op)
    type(eigen), intent(inout)    :: eig
    procedure(rntorn)  , optional :: A, B
    procedure(rntornxl), optional :: Op

    real(RP), dimension(:)  , allocatable :: resid, workd, workl
    real(RP), dimension(:,:), allocatable :: v
    integer, dimension(11) :: iParam, iPntr
    integer  :: ncv, lworkl
    integer  :: info, ido
    integer  :: i1, i2, j1, j2
    character(2) :: which
    character(1) :: bMat

    logical :: rVec, op_err
    logical, dimension(:), allocatable :: select


    !! which = part of the spectrum to be computed
    !!
    select case(eig%which)
    case(eig_which_sm)
       which = 'SM'
       if (eig%mode==eig_shift_invert) which = 'LM'

    case(eig_which_lm)
       which = 'LM'
       if (eig%mode==eig_shift_invert) which = 'SM'

    case default
       call quit("eigen_mod: arpack_ds: 'eig%which' = unknown")
    end select

    !! mode
    !!
    select case(eig%mode)
    case(eig_regular)
       iparam(7) = 1           

    case(eig_inverse)
       iparam(7) = 2           

    case(eig_shift_invert)
       iparam(7) = 3      

    case default
       call quit("eigen_mod: arpack_ds: 'eig%mode' = unknown")
    end select

    !! standard / generalised problem
    bmat = 'I'
    if (.NOT.eig%standard) then
       bmat = 'G'
    end if

    !! INITIALISATION
    !!
    ncv       = 2*eig%nev + 1     ! size for v
    if (ncv > eig%size) ncv = eig%size
    iparam(1) = 1                 ! ishifts = 1
    iparam(3) = eig%itMax         ! Max iter number
    iparam(4) = 1                 ! ???
    lworkl    = ncv*(ncv+8)
    call allocmem(workl, lworkl)
    call allocmem(resid, eig%size)
    call allocmem(workd, 3*eig%size)
    call allocmem(v, eig%size, ncv)
    !!
    resid = 0.0_rp
    workl = 0.0_rp
    workd = 0.0_rp
    v     = 0.0_rp

    !! ARPACK EIGEN-MODE COMPUTATION LOOP
    !!
    info      = 0
    resid     = 0.0_rp
    ido       = 0
    !!
    solve_loop: do while(.true.)

       call dsaupd( ido, bmat, eig%size, which, eig%nev, eig%tol, resid,&
            & ncv, v, eig%size, iparam, ipntr, workd, workl,&
            & lworkl, info )

       !! Checks the output parameter 'info'
       !! that must be equal to 0
       !!
       if (info /= 0) then
          if (info < 0 ) then
             write(*,*) "  info (see ARPACK 'dsaupd' doc) =", info
             call quit("eigen_mod: arpack_ds:&
                  & ARPACK 'dsaupd' error:,&
                  & info < 0, check parameters")
             
          else if (info == 1) then
             if (eig%verb>0) write(*,*) "  maximum number of iterations reached"
             call warning("eigen_mod: arpack_ds:&
                  & 'dsaupd' not converges", eig%verb)
             eig%ierr = 1
             exit solve_loop
             
          else if (info > 1 ) then
             write(*,*) "  info (see ARPACK 'dsaupd' doc) =", info
             call quit("eigen_mod: arpack_ds:&
                  & ARPACK 'dsaupd' error")
          end if
       end if

       !! Checks the output parameter 'ido'
       !! 
       if (ido .eq. -1 .or. ido .eq. 1) then
          
          if (eig%standard) then
             
             i1 = ipntr(1)
             i2 = ipntr(2)
             
             j1 = i1 + eig%size -1
             j2 = i2 + eig%size -1
             
             select case(eig%mode)
             case(eig_regular)
                
                call a(workd(i2:j2), workd(i1:j1))     
                eig%AEval = eig%AEval + 1
                
             case(eig_shift_invert)
                
                call op(workd(i2:j2), op_err, workd(i1:j1))     
                eig%OpEval = eig%OpEval + 1
                if (op_err) call warning(&
                     & "eigen_mod: arpack_ds:&
                     & check operation 'y = Op x' ", eig%verb-1)                

             end select
             
          else  !! GENERALISED EIGENPROBLEM 

             select case(eig%mode)
             case(eig_inverse)
                i1 = ipntr(1)
                i2 = ipntr(2)
             
                j1 = i1 + eig%size -1
                j2 = i2 + eig%size -1
                
                call a(workd(i2:j2), workd(i1:j1))     
                eig%AEval = eig%AEval + 1
                workd(i1:j1) = workd(i2:j2)
                
                call op(workd(i2:j2), op_err, workd(i1:j1))     
                eig%OpEval = eig%OpEval + 1
                if (op_err) call warning(&
                     & "eigen_mod: arpack_ds:&
                     & check operation 'y = Op x' ", eig%verb-1)

             case(eig_shift_invert)
                if (ido .eq. -1 ) then

                   i1 = ipntr(1)
                   i2 = ipntr(2)
             
                   j1 = i1 + eig%size -1
                   j2 = i2 + eig%size -1
                
                   call b(workd(i2:j2), workd(i1:j1))     
                   eig%BEval = eig%BEval + 1
                   workd(i1:j1) = workd(i2:j2)
                
                   call op(workd(i2:j2), op_err, workd(i1:j1))     
                   eig%OpEval = eig%OpEval + 1
                   if (op_err) call warning(&
                        & "eigen_mod: arpack_ds:&
                        & check operation 'y = Op x' ", eig%verb-1)

                else  !! ido = 1 here

                   i1 = ipntr(3)
                   i2 = ipntr(2)
             
                   j1 = i1 + eig%size -1
                   j2 = i2 + eig%size -1
                   
                   call op(workd(i2:j2), op_err, workd(i1:j1))     
                   eig%OpEval = eig%OpEval + 1
                   if (op_err) call warning(&
                        & "eigen_mod: arpack_ds:&
                        & check operation 'y = Op x' ", eig%verb-1)
                   
                end if

             end select
                
          end if

       else if (ido == 2) then

          i1 = ipntr(1)
          i2 = ipntr(2)
          
          j1 = i1 + eig%size -1
          j2 = i2 + eig%size -1
          
          call b(workd(i2:j2), workd(i1:j1))     
          eig%BEval = eig%BEval + 1

       else if (ido == 99) then
          exit solve_loop

       else
          write(*,*) "  ido (see ARPACK 'dsaupd' doc)  =", ido
          call quit("eigen_mod: arpack_ds:&
               & ARPACK 'dsaupd' error")
          
       end if
       
    end do solve_loop

    !! retrieve the number of performed Arnoldi's iterations 
    !!
    eig%iter = iparam(3)

    !! check convergence
    !!
    if (iparam(5) < eig%nev) then
          if (eig%verb>0) write(*,*) &
               &"  numer of converged eigen modes =", iparam(5)
          call warning("eigen_mod: arpack_ds:&
               & not converged", eig%verb)
          eig%ierr = 1
    end if

    !! POST-PROCESSING
    !!
    rvec = .true.
    call allocmem(select, ncv)
    call dseupd ( .true., 'A', select, eig%lambda, v, eig%size, &
         &        eig%shift, bmat, eig%size, which, eig%nev, &
         &        eig%tol, resid, ncv, v, eig%size, & 
         &        iparam, ipntr, workd, workl, lworkl, info )
    
    if ( info .ne. 0) then
          write(*,*) "  info (see dseupd arpack doc)   =", info
          call quit("eigen_mod: arpack_ds:&
               & ARPACK 'dseupd' error")

    end if

    !! RETRIEVE THE EIGEN-VECTORS
    !!
    eig%v(1:eig%size,1:eig%nev) = v(1:eig%size, 1:eig%nev)

  end subroutine arpack_ds


#endif

end module eigen_mod