SLEPC_utilities.f90

Go to the documentation of this file.

!------------------------------------------------------------------------------!
!------------------------------------------------------------------------------!
!------------------------------------------------------------------------------!
module slepc_utilities
#include <PB3D_macros.h>
#include <wrappers.h>
#include <slepc/finclude/slepceps.h>
! for slepc 3.6.0:
!#include <slepc/finclude/slepcepsdef.h>
! for slepc 3.5.3:
!#include <finclude/slepcepsdef.h>
    use str_utilities
    use output_ops
    use messages
    use slepceps
    use num_vars, only: iu, dp, max_str_ln
 
    implicit none
    private
    public insert_block_mat
#if ldebug
    public debug_insert_block_mat
#endif
    
    ! global variables
#if ldebug
    logical :: debug_insert_block_mat = .false.                                 
#endif
    
contains
    integer function insert_block_mat(mds,block,mat,r_id,ind,n_r,transp,&
        &overwrite,ind_insert) result(ierr)
        use x_vars, only: modes_type
        use num_vars, only: use_pol_flux_f
#if ldebug
        use num_vars, only: rank
        use input_utilities, only: pause_prog
#endif
        
        character(*), parameter :: rout_name = 'insert_block_mat'
        
        ! input / output
        type(modes_type), intent(in), target :: mds
        petscscalar, intent(in) :: block(:,:)                                   
        mat, intent(inout) :: mat                                               
        petscint, intent(in) :: r_id                                            
        petscint, intent(in) :: ind(2)                                          
        petscint, intent(in) :: n_r                                             
        petscbool, intent(in), optional :: transp                               
        petscbool, intent(in), optional :: overwrite                            
        petscbool, intent(in), optional :: ind_insert                           
        
        ! local variables
        petscint, pointer :: nm_x(:,:)                                          ! m (pol. flux) or n (tor. flux)
        petscbool :: transp_loc                                                 ! local copy of transp
        petscbool :: overwrite_loc                                              ! local copy of overwrite
        character(len=max_str_ln) :: err_msg                                    ! error message
        petscint :: operation                                                   ! either ADD_VALUES or INSERT_VALUES
        petscscalar, allocatable :: block_loc(:,:)                              ! local block, possibly shifted from block
        petscbool :: ind_insert_loc                                             ! local ind_insert
        
        ! initialize ierr
        ierr = 0
        
        ! set up local transp, overwrite and ind_insert
        transp_loc = .false.
        if (present(transp)) transp_loc = transp
        overwrite_loc = .false.
        if (present(overwrite)) overwrite_loc = overwrite
        ind_insert_loc = .false.
        if (present(ind_insert)) ind_insert_loc = ind_insert
        
        ! set operation
        if (overwrite_loc) then
            operation = insert_values
        else
            operation = add_values
        end if
        
        ! set nm_X
        if (use_pol_flux_f) then
            nm_x => mds%m
        else
            nm_x => mds%n
        end if
        
#if ldebug
        ! user output
        if (debug_insert_block_mat) then
            if (.not.ind_insert_loc) call sleep(rank)
            call writo('>>> rank '//trim(i2str(rank))//&
                &': at (k,m) = '//trim(i2str(r_id))//' + ('//&
                &trim(i2str(ind(1)))//','//trim(i2str(ind(2)))//'):',&
                &persistent=.true.)
            call lvl_ud(1)
        end if
#endif
        
        ! only set values if within matrix range
        if (minval(r_id+ind).ge.0 .and. maxval(r_id+ind).lt.n_r) then
            ! set error message
            err_msg = 'Couldn''t add values to matrix'
            
            ! initialize
            allocate(block_loc(size(block,1),size(block,2)))
            
            ! untransposed block
            call setup_local_block(nm_x,r_id,ind(1:2),block,block_loc)
            
            ! set values
            call matsetvaluesblocked(mat,1,r_id+ind(1),1,r_id+ind(2),&
                &transpose(block_loc),operation,ierr)
            chckerr(err_msg)
            
            ! untransposed block
            if (transp_loc) then
#if ldebug
                if (debug_insert_block_mat) then
                    call writo('Also at the transposed place:')
                end if
#endif
                
                ! transposed block
                call setup_local_block(nm_x,r_id,ind(2:1:-1),&
                    &transpose(conjg(block)),block_loc)
                
                call matsetvaluesblocked(mat,1,r_id+ind(2),1,r_id+ind(1),&
                    &transpose(block_loc),operation,ierr)
                chckerr(err_msg)
            end if
            
            ! clean up
            deallocate(block_loc)
        else
#if ldebug
            if (debug_insert_block_mat) &
                &call writo('Due to out of range no local block is added',&
                &persistent=.true.)
#endif
        end if
        
        ! clean up
        nullify(nm_x)
        
#if ldebug
        if (debug_insert_block_mat) then
            call pause_prog(ind_insert_loc)
            call lvl_ud(-1)
        endif
#endif
    contains
        subroutine setup_local_block(nm_X,r_id,ind,block,block_loc)
            use x_vars, only: n_mod_x
            
            ! input / output
            petscint, intent(in) :: nm_x(:,:)                                   
            petscint, intent(in) :: r_id                                        
            petscint, intent(in) :: ind(2)                                      
            petscscalar, intent(in) :: block(:,:)                               
            petscscalar, intent(inout) :: block_loc(:,:)                        
            
            ! local variables
            petscint :: k, m                                                   ! counters
            petscint :: k_loc, m_loc                                            ! local k and m
            
            block_loc = 0._dp
            do m = 1,n_mod_x
                do k = 1,n_mod_x
                    k_loc = k + nm_x(r_id+1,k) - &
                        &nm_x(min(r_id+1+ind(1),size(nm_x,1)),k)
                    m_loc = m + nm_x(r_id+1,m) - &
                        &nm_x(min(r_id+1+ind(2),size(nm_x,1)),m)
                    if (k_loc.ge.1 .and. m_loc.ge.1 .and. &
                        &k_loc.le.n_mod_x .and. m_loc.le.n_mod_x) &
                        block_loc(k_loc,m_loc) = block(k,m)
                end do
            end do
            
#if ldebug
            if (debug_insert_block_mat) then
                call writo('following local block is going to be added by &
                    &rank '//trim(i2str(rank))//':',persistent=.true.)
                call writo('Re =',persistent=.true.)
                call print_ar_2(rp(block_loc))
                call writo('Im =',persistent=.true.)
                call print_ar_2(ip(block_loc))
                if (overwrite_loc) then
                    call writo('(with operation INSERT_VALUES)',&
                        &persistent=.true.)
                else
                    call writo('(with operation ADD_VALUES)',&
                        &persistent=.true.)
                end if
            endif
#endif
        end subroutine setup_local_block
    end function insert_block_mat
end module slepc_utilities