Contracted_Symbolic_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/>.
 
module contracted_symbolic_module
 
    use const_gbl,              only: stdout
    use consts_mpi_ci,          only: default_integral_threshold
    use containers,             only: bstree
    use integer_packing,        only: pack8ints, unpack8ints
    use iso_c_binding,          only: c_loc, c_ptr, c_f_pointer
    use mpi_gbl,                only: myrank, mpiint, mpi_mod_allgather, mpi_mod_rotate_cfp_arrays_around_ring, &
                                      mpi_mod_rotate_int_arrays_around_ring, mpi_xermsg, mpi_reduceall_max
    use precisn,                only: longint, wp
    use memorymanager_module,   only: master_memory
    use parallelization_module, only: grid => process_grid
    use symbolic_module,        only: symbolicelementvector
    use timing_module,          only: master_timer
 
    implicit none
 
    public contractedsymbolicelementvector
 
    private
 
    type, extends(bstree) :: contractedsymbolicelementvector
        private
 
        integer(longint), allocatable :: electron_integral(:,:)
        real(wp),         allocatable :: coeffs(:,:,:)
 
        integer  :: n = 0
        logical  :: constructed = .false.   
        real(wp) :: threshold  = 0.0
        integer  :: num_states_1 = 0
        integer  :: num_states_2 = 0
        integer  :: max_capacity = 0
        integer  :: expand_size = 100
    contains
        ! type-bound procedures used by the bstree interface
        procedure, public  :: compare => bstree_compare
 
        ! own type-bound procedures
        procedure, public  :: construct
        procedure, public  :: insert_ijklm_symbol
        procedure, public  :: insert_symbol
        procedure, public  :: remove_symbol_at
        procedure, public  :: is_empty
        procedure, public  :: clear
        procedure, public  :: check_same_integral
        procedure, public  :: get_integral_label
        procedure, public  :: synchronize_symbols
        procedure, public  :: estimate_synchronize_cost
        procedure, public  :: modify_coeff
        procedure, public  :: get_coefficient
        procedure, public  :: get_coeff_and_integral
        procedure, public  :: get_size
        procedure, public  :: get_num_targets_sym1
        procedure, public  :: get_num_targets_sym2
        procedure, public  :: add_symbols
        procedure, public  :: reduce_symbols
        procedure, public  :: print => print_symbols
       !procedure, public  :: prune_threshold
       !procedure, public  :: sort_symbols
        procedure, private :: expand_array
        procedure, private :: check_bounds
        procedure, public  :: destroy
        procedure, private :: check_constructed
        procedure, private :: synchronize_symbols_ii
       !procedure, public  :: count_occurance
       !procedure, public  :: construct
    end type  contractedsymbolicelementvector
 
    class(contractedsymbolicelementvector), pointer :: to_be_synched
 
contains
 
    integer function bstree_compare (this, i, j, data) result (verdict)
 
        class(contractedsymbolicelementvector), intent(in) :: this
        integer,                      intent(in) :: i, j
        type(c_ptr),        optional, intent(in) :: data
 
        integer(longint)          :: ii(2), jj(2)
        integer(longint), pointer :: kk(:)
 
        ! get the default index set to compare with (if provided)
        if (present(data)) call c_f_pointer(data, kk, (/ 2 /))
 
        ! obtain pointers to the index sets
        if (i <= 0) then; ii = kk; else; ii = this % electron_integral(:,i); end if
        if (j <= 0) then; jj = kk; else; jj = this % electron_integral(:,j); end if
 
        ! compare the index sets using lexicographical comparison
        if (ii(1) < jj(1)) then
            verdict = -1
        else if (ii(1) > jj(1)) then
            verdict = +1
        else if (ii(2) < jj(2)) then
            verdict = -1
        else if (ii(2) > jj(2)) then
            verdict = +1
        else
            verdict = 0
        end if
 
    end function bstree_compare
 
 
    subroutine check_constructed (this)
        class(contractedsymbolicelementvector) :: this
        if (.not. this % constructed) then
            write (stdout, "('Vector::constructed - Vector is not constructed')")
            stop "Vector::constructed - Vector is not constructed"
        end if
 
    end subroutine check_constructed
 
 
    subroutine construct (this, n1, n2, threshold, initial_size)
        class(contractedsymbolicelementvector) :: this
        real(wp), optional, intent(in)         :: threshold
        integer,  optional, intent(in)         :: initial_size
 
        integer :: err, n1, n2
 
        if (present(threshold)) then
            this % threshold = threshold
        else
            this % threshold = default_integral_threshold
        end if
 
        this % expand_size = 100
        this % max_capacity = this % expand_size
        this % num_states_1 = n1
        this % num_states_2 = n2
 
        !Allocate the vectors
        allocate(this % electron_integral(2, this % max_capacity), this % coeffs(n1, n2, this % max_capacity), stat = err)
        call master_memory % track_memory (kind(this % electron_integral), &
                                           size(this % electron_integral), err, 'CONSYMBOLIC::ELECTRONINT')
        call master_memory % track_memory(kind(this % coeffs), &
                                          size(this % coeffs), err, 'CONSYMBOLIC::ELECTRONCOEFF')
        if (err /= 0) then
            call mpi_xermsg('Contracted_Symbolic_module', 'construct', 'SymbolicVector arrays not allocated', 1, 1)
        end if
 
        !Do an intial clear
        call this % clear
 
        !This has been succesfully constructed
        this % constructed = .true.
 
    end subroutine construct
 
 
    logical function check_same_integral (this, integral, idx) result (found)
 
        class(contractedsymbolicelementvector), intent(in)  :: this
        integer(longint),  target,              intent(in)  :: integral(2)
        integer,                                intent(out) :: idx
 
        call this % check_constructed
 
        idx   = this % locate(-1, c_loc(integral))
        found = idx > 0
 
    end function check_same_integral
 
 
    subroutine insert_ijklm_symbol (this, i, j, k, l, m, coeffs, check_same_)
        class(contractedsymbolicelementvector) :: this
        integer,  intent(in)            :: i, j, k, l, m
        real(wp), intent(in)            :: coeffs(this % num_states_1, this % num_states_2)
        logical,  intent(in), optional  :: check_same_
 
        integer(longint) :: integral_label(2)
        logical          :: check_same
 
        if (present(check_same_)) then
            check_same = check_same_
        else
            check_same = .true.
        end if
 
        call pack8ints(i, j, k, l, m, 0, 0, 0, integral_label)
        call this%insert_symbol(integral_label, coeffs, check_same)
 
    end subroutine insert_ijklm_symbol
 
 
    subroutine insert_symbol (this, integral_label, coeffs, check_same_)
        class(contractedsymbolicelementvector) :: this
        integer(longint), intent(in)           :: integral_label(2)
        real(wp),         intent(in)           :: coeffs(this % num_states_1, this % num_states_2)
        logical,          intent(in), optional :: check_same_
 
        integer :: idx = 0, check_idx
        logical :: check_same
 
        if (present(check_same_)) then
            check_same = check_same_
        else
            check_same = .true.
        end if
 
        ! Filter small coefficients
        !if(abs(coeff) < this%threshold) return
 
        if (.not. this % check_same_integral(integral_label, idx)) then
                !If not then update number of elements
                this % n = this % n + 1
                !The idx is the last element
                idx = this % n
                !If we've exceeded capacity then expand
                if (this % n > this % max_capacity) then
                    call this % expand_array()
                end if
                !Place into the array
                this % electron_integral(:,idx) = integral_label(:)
                !Zero the coefficient since we can guareentee the compiler will do so
                this % coeffs(:,:,idx) = 0.0_wp
                !Insert it into the binary tree
                call this % insert(idx)
        end if
 
        !Now we have an index, lets add in the coeffcient
        this % coeffs(:,:,idx) = this % coeffs(:,:,idx) + coeffs(:,:)
 
    end subroutine insert_symbol
 
 
    subroutine expand_array (this)
        class(contractedsymbolicelementvector) :: this
        integer(longint), allocatable :: temp_integral(:,:)
        real(wp),         allocatable :: temp_coeffs(:,:,:)
        integer :: temp_capacity, ido, err
 
        !write(stdout,"('Size of arrays = ',5i12)") this%max_capacity,this%num_states_1,this%num_states_2,size(this%electron_integral),size(this%coeffs)
 
        !Of course lets check if we are constructed
        call this % check_constructed
 
        !Now the tempcapacity is the current max capacity
        temp_capacity = this % max_capacity
 
        allocate(temp_integral(2, this % max_capacity), &
                 temp_coeffs(this % num_states_1, this % num_states_2, this % max_capacity), stat = err)
        call master_memory % track_memory(kind(temp_integral), size(temp_integral), err, 'CONSYMBOLIC::EXP::TEMPINT')
        call master_memory % track_memory(kind(temp_coeffs), size(temp_coeffs), err, 'CONSYMBOLIC::EXP::TEMPCOEFF')
 
        !Copy over the elements into the temporary array
        do ido = 1, temp_capacity
            temp_integral(1,ido) = this % electron_integral(1,ido)
            temp_integral(2,ido) = this % electron_integral(2,ido)
            temp_coeffs(:,:,ido) = this % coeffs(:,:,ido)
        end do
 
        call master_memory % free_memory(kind(this % electron_integral), size(this % electron_integral))
        call master_memory % free_memory(kind(this % coeffs), size(this % coeffs))
 
        !deallocate our old values
        deallocate(this % electron_integral)
        deallocate(this % coeffs)
 
        !Increase the max capacity
        this % max_capacity = this % max_capacity + this % expand_size
 
        !Allocate our new array with the new max capacity
        allocate(this % electron_integral(2, this % max_capacity), &
                 this % coeffs(this % num_states_1, this % num_states_2, this % max_capacity), stat = err)
        call master_memory % track_memory (kind(this % electron_integral), &
                                           size(this % electron_integral), err,'CONSYMBOLIC::EXP::ELECTRONINT_EXP')
        call master_memory % track_memory(kind(this % coeffs), &
                                          size(this % coeffs), err, 'CONSYMBOLIC::EXP::ELECTRONCOEFF_EXP')
        !Zero the elements
        this % electron_integral = -1
        this % coeffs = 0.0_wp
 
        do ido = 1, temp_capacity
            !Copy our old values back
            this % electron_integral(1,ido) = temp_integral(1,ido)
            this % electron_integral(2,ido) = temp_integral(2,ido)
            this % coeffs(:,:,ido)          = temp_coeffs(:,:,ido)
        end do
 
        call master_memory % free_memory(kind(temp_integral), size(temp_integral))
        call master_memory % free_memory(kind(temp_coeffs), size(temp_coeffs))
 
        deallocate(temp_integral, temp_coeffs)
 
    end subroutine expand_array
 
 
    subroutine    add_symbols (this, rhs, alpha)
        class(contractedsymbolicelementvector)   :: this
        class(symbolicelementvector), intent(in) :: rhs
        real(wp)         :: alpha(this % num_states_1, this % num_states_2), int_coeff
        integer          :: ido
        integer(longint) :: label(2)
 
        !Loop through each rhs symbolic and add it to me scaled by alpha
        do ido = 1, rhs % get_size()
            call rhs % get_coeff_and_integral(ido, int_coeff, label)
            call this % insert_symbol(label, int_coeff * alpha)
        end do
 
    end subroutine add_symbols
 
 
    subroutine reduce_symbols (this, rhs)
        class(contractedsymbolicelementvector)             :: this
        class(contractedsymbolicelementvector), intent(in) :: rhs
        integer          :: ido
        integer(longint) :: label(2)
 
        !Loop through each rhs symbolic and add it to me scaled by alpha
        do ido = 1, rhs % get_size()
            call this % insert_symbol(rhs % electron_integral(:,ido), rhs % coeffs(:,:,ido))
        end do
 
    end subroutine reduce_symbols
 
 
    subroutine synchronize_symbols (this)
        class(contractedsymbolicelementvector) :: this
 
        integer :: my_num_symbols, largest_num_symbols, ierr
        integer :: ido, proc_id, jdo
 
        integer(longint), allocatable, target :: labels(:,:)
        integer(longint), pointer             :: label_ptr(:)
        integer(longint)                      :: procs_num_of_symbols_int, procs_num_of_symbols_coeff, size_labels, size_coeffs
 
        real(wp), allocatable, target :: coeffs(:,:,:)
        real(wp), pointer             :: coeffs_ptr(:)
 
        if (grid % gprocs <= 1) then
            return
        end if
 
        call master_timer % start_timer("Symbol Synchronize")
 
        my_num_symbols = this % get_size()
        call mpi_reduceall_max(my_num_symbols, largest_num_symbols, grid % lcomm)
 
        if (largest_num_symbols == 0) then
            call master_timer % stop_timer("Symbol Synchronize")
            return
        end if
 
        call master_memory % track_memory (kind(labels), largest_num_symbols * 2, 0, 'CONSYMBOLIC::MPISYNCH::LABELS')
        call master_memory % track_memory (kind(coeffs), largest_num_symbols * this % num_states_1 * this % num_states_2, &
                                            0, 'CONSYMBOLIC::MPISYNCH::COEFFS')
 
        allocate(labels(2, largest_num_symbols), coeffs(this % num_states_1, this % num_states_2, largest_num_symbols), stat = ierr)
 
        size_labels = size(labels, kind = longint)
        size_coeffs = size(coeffs, kind = longint)
 
        label_ptr(1:size_labels) => labels(:,:)
        coeffs_ptr(1:size_coeffs) => coeffs(:,:,:)
 
        labels = 0
        coeffs = 0
 
        if (my_num_symbols > 0) then
            labels(1,1:my_num_symbols) = this % electron_integral(1,1:my_num_symbols)
            labels(2,1:my_num_symbols) = this % electron_integral(2,1:my_num_symbols)
            do ido = 1, my_num_symbols
                coeffs(1:this%num_states_1,1:this%num_states_2,ido) = this % coeffs(1:this%num_states_1,1:this%num_states_2,ido)
            end do
        end if
 
        procs_num_of_symbols_int = my_num_symbols * 2
        procs_num_of_symbols_coeff = my_num_symbols * this % num_states_1 * this % num_states_2
        do proc_id = 1, grid % lprocs - 1
            call mpi_mod_rotate_int_arrays_around_ring(procs_num_of_symbols_int, label_ptr, size_labels, grid % lcomm)
            call mpi_mod_rotate_cfp_arrays_around_ring(procs_num_of_symbols_coeff, coeffs_ptr, size_coeffs, grid % lcomm)
 
            do ido = 1, procs_num_of_symbols_int / 2
                call this % insert_symbol(labels(1:2,ido), coeffs(1:this%num_states_1,1:this%num_states_2,ido))
            end do
        end do
 
        call master_timer % stop_timer("Symbol Synchronize")
        call master_memory % free_memory(kind(labels), size(labels))
        call master_memory % free_memory(kind(coeffs), size(coeffs))
 
        deallocate(labels, coeffs)
 
    end subroutine synchronize_symbols
 
 
    subroutine synchronize_symbols_ii (this)
        class(contractedsymbolicelementvector), target :: this
 
        integer(longint), allocatable, target :: labels(:,:)
        integer(longint), pointer             :: label_ptr(:)
 
        real(wp), allocatable, target :: coeffs(:,:,:)
        real(wp), pointer             :: coeffs_ptr(:)
 
        integer :: my_num_symbols, largest_num_symbols, procs_num_of_symbols_int, procs_num_of_symbols_coeff, ido, proc_id, jdo
        integer :: ierr, local_communicator, global_communicator, odd_even
 
        integer(kind=mpiint) :: master_comm, global_rank, global_nprocs, error, proc, tag = 1
 
        if (grid % gprocs <= 1) then
            return
        end if
 
    end subroutine synchronize_symbols_ii
 
 
    logical function check_bounds (this, i)
        class(contractedsymbolicelementvector) :: this
        integer, intent(in)                    :: i
 
        if (i <= 0 .or. i > this % n) then
            write (stdout, "('SymbolicVector::check_bounds - Out of Bounds access')")
            stop "SymbolicVector::check_bounds - Out of Bounds access"
            check_bounds = .false.
        else
            check_bounds = .true.
        end if
 
    end function check_bounds
 
 
    subroutine remove_symbol_at (this, idx)
        class(contractedsymbolicelementvector) :: this
        integer, intent(in) :: idx
        integer             :: ido
 
        if (this % check_bounds(idx)) then
            do ido = idx + 1, this % n
                this % electron_integral(1, ido - 1) = this % electron_integral(1, ido)
                this % electron_integral(2, ido - 1) = this % electron_integral(2, ido)
                !this%coeffs(ido-1)        = this%coeffs(ido)
            end do
            this % electron_integral(:, this % n) = 0
           !this % coeffs(this % n) = 0
            this % n = this % n - 1
        end if
 
    end subroutine remove_symbol_at
 
 
    logical function is_empty (this)
 
        class(contractedsymbolicelementvector) :: this
 
        is_empty = (this % n == 0)
 
    end function is_empty
 
 
    subroutine clear (this)
 
        class(contractedsymbolicelementvector) :: this
 
        this % n = 0
        call this % bstree % destroy
 
    end subroutine clear
 
 
    subroutine modify_coeff (this, idx, coeff)
        class(contractedsymbolicelementvector) :: this
        integer, intent(in) :: idx
        real(wp) :: coeff
        !this%coeffs(idx) = this%coeffs(idx)+ coeff
    end subroutine modify_coeff
 
 
    function get_integral_label (this, idx)
        class(contractedsymbolicelementvector) :: this
        integer, intent(in)                    :: idx
        integer(longint), dimension(2)         :: get_integral_label
 
        if (this % check_bounds(idx)) get_integral_label(:) = this % electron_integral(:,idx)
 
    end function get_integral_label
 
 
    real(wp) function get_coefficient(this, idx, n1, n2)
        class(contractedsymbolicelementvector) :: this
        integer, intent(in) :: idx, n1, n2
 
        !if(this%check_bounds(idx)==.true.) get_coefficient=this%coeffs(idx)
        get_coefficient = this % coeffs(n1, n2, idx)
 
    end function get_coefficient
 
 
    subroutine get_coeff_and_integral (this, idx, coeff, label)
        class(contractedsymbolicelementvector), intent(in) :: this
        integer,          intent(in)  :: idx
        real(wp),         intent(out) :: coeff(this % num_states_1, this % num_states_2)
        integer(longint), intent(out) :: label(2)
 
        if (this % check_bounds(idx)) then
            label(1:2) = this % electron_integral(1:2,idx)
            coeff(:,:) = this % coeffs(:,:,idx)
        end if
 
    end subroutine get_coeff_and_integral
 
 
    integer function get_size (this)
        class(contractedsymbolicelementvector) :: this
        get_size = this % n
    end function get_size
 
 
    integer function get_num_targets_sym1 (this)
        class(contractedsymbolicelementvector) :: this
        get_num_targets_sym1 = this % num_states_1
    end function get_num_targets_sym1
 
 
    integer function get_num_targets_sym2 (this)
        class(contractedsymbolicelementvector) :: this
        get_num_targets_sym2 = this % num_states_2
    end function get_num_targets_sym2
 
 
    subroutine destroy (this)
        class(contractedsymbolicelementvector) :: this
 
        if (allocated(this % electron_integral)) then
            call master_memory % free_memory(kind(this % electron_integral), size(this % electron_integral))
            deallocate(this % electron_integral)
        end if
        if(allocated(this % coeffs)) then
            call master_memory % free_memory(kind(this % coeffs), size(this % coeffs))
            deallocate(this % coeffs)
        end if
 
        call this % bstree % destroy
 
        this % constructed = .false.
 
    end subroutine destroy
 
 
    subroutine print_symbols (this)
        class(contractedsymbolicelementvector) :: this
        integer  :: labels(8), ido
        real(wp) :: coeff
 
        if (.not. this % is_empty()) write (stdout, "('Outputting symbolic elements....')")
 
        do ido = 1, this % n
            call unpack8ints(this % electron_integral(1,ido), labels)
            write (1923 + myrank, "(5i4,' -- ',5es14.3)") labels(1:5), this % coeffs(:,:,ido)
        end do
 
    end subroutine print_symbols
 
 
    logical function estimate_synchronize_cost (this)
        class(contractedsymbolicelementvector) :: this
        integer(longint) :: memory_cost
        integer          :: my_num_symbols, largest_num_symbols
        integer          :: gathered_bool, global_bool(grid % gprocs)
 
        if (grid % gprocs <= 1) then
            estimate_synchronize_cost = .true.
            return
        end if
 
        my_num_symbols = this % get_size()
        gathered_bool = 0
        estimate_synchronize_cost = .true.
 
        call mpi_reduceall_max(my_num_symbols, largest_num_symbols, grid % gcomm)
 
        memory_cost = 2 * largest_num_symbols * (2 + this % num_states_1 * this % num_states_2) * 8
 
        if (memory_cost >= master_memory % get_scaled_available_memory(0.75_wp)) gathered_bool = 1
 
        global_bool = 0
 
        call mpi_mod_allgather(gathered_bool, global_bool, grid % gcomm)
 
        if (sum(global_bool) /= 0) then
            estimate_synchronize_cost = .false.
        end if
 
    end function estimate_synchronize_cost
 
end module contracted_symbolic_module