SLEPCDiagonalizer_module.F90

Go to the documentation of this file.

! Copyright 2019
!
! For a comprehensive list of the developers that contributed to these codes
! see the UK-AMOR website.
!
! This file is part of UKRmol-in (UKRmol+ suite).
!
!     UKRmol-in is free software: you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation, either version 3 of the License, or
!     (at your option) any later version.
!
!     UKRmol-in is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with  UKRmol-in (in source/COPYING). Alternatively, you can also visit
!     <https://www.gnu.org/licenses/>.
 
!> \brief   Diagonalizer type using SLEPc backend
!> \authors A Al-Refaie
!> \date    2017
!>
!> This module will be built only if SLEPc is available.
!>
!> \note 16/01/2019 - Jakub Benda: Unifom coding style and expanded documentation.
!>
module slepcdiagonalizer_module
 
    use blas_lapack_gbl,           only: blasint
    use const_gbl,                 only: stdout
    use consts_mpi_ci,             only: pass_to_cdenprop
    use cdenprop_defs,             only: civect
    use mpi_gbl,                   only: master, mpi_xermsg
    use precisn,                   only: longint, wp
 
    use basematrix_module,         only: basematrix
    use baseintegral_module,       only: baseintegral
    use diagonalizer_module,       only: basediagonalizer
    use diagonalizerresult_module, only: diagonalizerresult
    use options_module,            only: options
    use parallelization_module,    only: grid => process_grid
    use timing_module,             only: master_timer
    use slepcmatrix_module,        only: slepcmatrix
 
    use petsc
    use slepceps
 
#include <finclude/petsc.h>
#include <finclude/slepceps.h>
 
    implicit none
 
    public slepcdiagonalizer
 
    private
 
    type, extends(basediagonalizer) :: slepcdiagonalizer
    contains
        procedure, public :: diagonalize => diagonalize_slepc
       !procedure, public :: DoNonSLEPCMat
        procedure, public :: useslepcmat
        procedure, public :: selecteps
    end type
 
contains
 
    subroutine useslepcmat (this, matrix_elements, mat_out)
        class(slepcdiagonalizer)          :: this
        class(slepcmatrix), intent(in)    :: matrix_elements
        mat, pointer,       intent(inout) :: mat_out
 
        write(stdout,"('--------Optimized SLEPC Matrix Format chosen------')")
 
        mat_out => matrix_elements % get_PETSC_matrix()
 
        if (.not. associated(mat_out)) then
            stop "Matrix doesn;t exist for SLEPC, are you sure you created it?"
        end if
 
    end subroutine useslepcmat
 
 
    subroutine selecteps (this, eps)
        class(slepcdiagonalizer) :: this
        eps, intent(inout)       :: eps
        petscerrorcode           :: ierr
 
        write(stdout,"('--------KRYLOVSCHUR used as Diagonalizer------')")
 
        !     ** Create eigensolver context
        call epssettype(eps,epskrylovschur,ierr)
        !call EPSSetTrueResidual(eps,PETSC_TRUE,ierr);
        !call EPSLanczosSetReorthog(eps,EPS_LANCZOS_REORTHOG_SELECTIVE,ierr)
        if (ierr /= 0) then
            call mpi_xermsg('SLEPCDiagonalizer_module', 'SelectEPS', 'Failed to select eigensolver', int(ierr), 1)
        end if
 
    end subroutine selecteps
 
 
    subroutine diagonalize_slepc (this, matrix_elements, num_eigenpair, dresult, all_procs, option, integrals)
        class(slepcdiagonalizer)        :: this
        class(diagonalizerresult)       :: dresult
        class(basematrix),   intent(in) :: matrix_elements
        class(baseintegral), intent(in) :: integrals
        type(options),       intent(in) :: option
        integer,             intent(in) :: num_eigenpair
        logical,             intent(in) :: all_procs
 
        integer               :: mat_dimen, number_solved
        integer(blasint)      :: grid_master_context
        real(wp), allocatable :: diagonal_temp(:), eigen(:)
 
        type(civect) :: eigenvector
 
        !-----SLEPC values
        mat, pointer         :: a
        eps                  :: eps
        epstype              :: tname
        st                   :: st_eps
        epsconvergedreason   :: conv_reason
        petscreal            :: maxtol,tol, error
        petscscalar          :: kr, ki, current_shift
        petscscalar, pointer :: xx_v(:)
        vec                  :: xr, xi, f_xr
        petscint             :: nev, maxit, its, nconv, matrix_size, ido, n
        petscerrorcode       :: ierr
        vecscatter           :: vsctx
        ksp                  :: ksp
        pc                   :: pc
 
        maxit = option % max_iterations
        maxtol = option % max_tolerance
 
        if (maxit < 0) maxit = petsc_default_integer
 
        matrix_size = matrix_elements % get_matrix_size()
        mat_dimen = matrix_size
 
        write (stdout, "()")
        write (stdout, "('-----------WARNING----------------------- ')")
        write (stdout, "('When dealing with symmetries that are degenerate SLEPC may give incorrect values')")
        write (stdout, "('YOU HAVE BEEN WARNED')")
        write (stdout, "('---------------------------------------- ')")
        write (stdout, "('N: ',i8)") matrix_size
        write (stdout, "('Requested # of eigenpairs',i8)") num_eigenpair
        write (stdout, "('Maximum iterations',i8)") maxit
        write (stdout, "('Maximum tolerance',es9.2)") maxtol
        write (stdout, "()")
 
        select type(matrix_elements)
            type is (slepcmatrix)
                call this % UseSLEPCMat(matrix_elements, a)
            class is (basematrix)
                stop "Matrix format not yet implemented for SLEPC"
            !class is (BaseMatrix)
                !call this%DoNonSLEPCMat(matrix_elements,num_eigenpair,eigen,vecs,maxit,maxtol)
        end select
 
        call matcreatevecs(a, xr, xi, ierr)
 
        if (all_procs) then
            call vecscattercreatetoall(xr, vsctx, f_xr, ierr)
        else
            call vecscattercreatetozero(xr, vsctx, f_xr, ierr)
        end if
 
        call epscreate(int(grid % gcomm, kind(petsc_comm_world)), eps, ierr)
        call this % SelectEPS(eps)
        call epssetoperators(eps, a, petsc_null_mat, ierr)
        call epssetproblemtype(eps, eps_hep, ierr)
 
!#if defined(PETSC_HAVE_MUMPS)
!            if(num_eigenpair > 20) then
!                call EPSGetST(eps,st_eps,ierr);
!                call STSetType(st_eps,STSINVERT,ierr)
!                call EPSSetST(eps,st_eps,ierr)
!                call EPSSetWhichEigenpairs(eps,EPS_ALL,ierr);
!                   call EPSSetInterval(eps,-90,-80,ierr);
!                call STGetKSP(st_eps,ksp,ierr);
!                call KSPSetType(ksp,KSPPREONLY,ierr);
!                call KSPGetPC(ksp,pc,ierr);
!
!                call PCSetType(pc,PCCHOLESKY,ierr);
!
!                call EPSKrylovSchurSetDetectZeros(eps,PETSC_TRUE,ierr)
!                call PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS,ierr)
!
!                 call PetscOptionsInsertString("-mat_mumps_icntl_13 1 -mat_mumps_icntl_24 1 -mat_mumps_cntl_3 1e-12",ierr);
!
!            else
!                call EPSSetWhichEigenpairs(eps,    EPS_SMALLEST_REAL ,ierr);
!            endif
!#else
!                call EPSSetWhichEigenpairs(eps,    EPS_SMALLEST_REAL ,ierr);
!#endif
 
        !call MatView(A,    PETSC_VIEWER_DRAW_WORLD,ierr)
 
        !** Set operators. In this case, it is a standard eigenvalue problem
 
        !Create the ST
        n = num_eigenpair
        call epssetdimensions(eps, n, petsc_decide, petsc_decide, ierr)
        !  call EPSSetType(eps,EPSJD,ierr)
        call epssettolerances(eps, maxtol, maxit, ierr)
        call master_timer % start_timer("SLEPC solve")
        !call EPSView(eps,    PETSC_VIEWER_STDOUT_WORLD,ierr)
        call epssolve(eps, ierr)
        if (ierr /= 0) then
            call mpi_xermsg('SLEPCDiagonalizer_module', 'SelectEPS', 'Failed to select eigensolver', int(ierr), 1)
        end if
        call master_timer % stop_timer("SLEPC solve")
        !call EPSGetConvergedReason(eps,conv_reason,ierr)
 
        ! write(stdout,"('Converged reason ',i8)") conv_reason
        call epsgetiterationnumber(eps, its, ierr)
        call epsgettolerances(eps, tol, maxit, ierr)
 
        if (grid % grank == master) then
            write (stdout, *)
            write (stdout, "(('/ Stopping condition: tol=',1P,E11.4,', maxit=',I0,', stopped at it=',I0))") tol, maxit, its
        end if
 
        call epsgetconverged(eps, nconv, ierr)
 
        if (nconv < num_eigenpair) then
            write (stdout, *) 'Not all requested eigenpairs have converged!!!!'
            write (stdout, *) 'Only ', nconv, ' have converged against ', num_eigenpair, ' requested'
            stop "Diagoanlization not completed"
        end if
 
        allocate(eigen(num_eigenpair))
 
        do ido = 0, min(nconv, n) - 1
            ! ** Get converged eigenpairs: i-th eigenvalue is stored in kr
            ! ** (real part) and ki (imaginary part)
            call epsgeteigenpair(eps, ido, kr, ki, xr, xi, ierr)
            eigen(ido + 1) = petscrealpart(kr)
        end do
 
        ! store Hamiltonian information into the CDENPROP vector
        if (iand(option % vector_storage_method, pass_to_cdenprop) /= 0) then
            call dresult % export_header(option, integrals)
        end if
 
        ! write Hamiltonian information into the output disk file (the current context master will do this)
        if (grid % grank == master) then
            call dresult % write_header(option, integrals)
        end if
 
        ! save eigenvalues by master of the context (to disk), or by everyone (to internal arrays)
        if (all_procs .or. grid % grank == master) then
            call dresult % handle_eigenvalues(eigen, matrix_elements % diagonal, num_eigenpair, mat_dimen)
        end if
 
        ! when the post-processing stage is required, set up the needed ScaLAPACK matrices
        if (iand(dresult % vector_storage, pass_to_cdenprop) /= 0) then
 
            ! split off a 1-by-1 (i.e. master only) sub-grid
            grid_master_context = grid % gcntxt
#if defined(usempi) && defined(scalapack)
            call blacs_gridinit(grid_master_context, 'R', 1_blasint, 1_blasint)
#endif
 
            ! set up ScaLAPACK storage for a single eigenvector (master-only ScaLAPACK matrix)
            eigenvector % blacs_context = grid_master_context
            eigenvector % CV_is_scalapack = .true.
            call eigenvector % init_CV(mat_dimen, 1)
 
            ! set up ScaLAPACK storage for all eigenvectors, distributed among processes in the current BLACS grid
            dresult % ci_vec % blacs_context   = grid % gcntxt
            dresult % ci_vec % CV_is_scalapack = .true.
            call dresult % ci_vec % init_CV(mat_dimen, num_eigenpair)
 
        end if
 
        do ido = 0, min(nconv, n) - 1
            ! ** Get converged eigenpairs: i-th eigenvalue is stored in kr
            ! ** (real part) and ki (imaginary part)
            call epsgeteigenpair(eps, ido, kr, ki, xr, xi, ierr)
 
            ! ** Compute the relative error associated to each eigenpair
            !call EPSComputeError(eps,ido,EPS_ERROR_RELATIVE,error,ierr)
            !write(stdout,*) PetscRealPart(kr), error
            !160      format (1P,'   ',E12.4,'       ',E12.4)
            !eigen(ido+1) = PetscRealPart(kr)
            ! write(stdout,*)eigen(ido+1)
 
            !Grab the vectors
            call vecscatterbegin(vsctx, xr, f_xr, insert_values, scatter_forward, ierr)
            call vecscatterend(vsctx, xr, f_xr, insert_values, scatter_forward, ierr)
            if (all_procs .or. grid % grank == master) then
                call vecgetarrayread(f_xr, xx_v, ierr)
                call dresult % handle_eigenvector(xx_v(1:mat_dimen), mat_dimen)
            end if
            if (iand(dresult % vector_storage, pass_to_cdenprop) /= 0) then
                if (grid % grank == master) then
                    eigenvector % CV(1:mat_dimen, 1) = xx_v(1:mat_dimen)
                end if
                call dresult % ci_vec % redistribute(eigenvector, grid % gcntxt, mat_dimen, 1, 1, 1, 1, ido + 1)
            end if
            if (all_procs .or. grid % grank == master) then
                call vecrestorearrayread(f_xr, xx_v, ierr)
            end if
        end do
 
        if (iand(dresult % vector_storage, pass_to_cdenprop) /= 0) then
            call dresult % export_eigenvalues(eigen(1:num_eigenpair), matrix_elements % diagonal, num_eigenpair, &
                                              mat_dimen, dresult % ci_vec % ei, dresult % ci_vec % iphz)
            write (stdout, '(/," Eigensolutions exported into the CIVect structure.",/)')
        end if
 
        ! finish writing the solution file and let other processes access it
        if (grid % grank == master) then
            call dresult % finalize_solutions(option)
        end if
 
        call epsdestroy(eps, ierr)
        call vecdestroy(xr, ierr)
        call vecdestroy(xi, ierr)
 
        deallocate(eigen)
 
    end subroutine diagonalize_slepc
 
end module slepcdiagonalizer_module