A routine to invert the field matrix locally
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
class(supercell_type), | intent(inout) | :: | self |
subroutine sc_invert_local(self)
use mat_inv, only: invert_serial
use mp, only: broadcast_sub
use file_utils, only: error_unit
use gs2_time, only: code_dt
implicit none
class(supercell_type), intent(inout) :: self
complex, dimension(:,:), allocatable :: tmp_arr
real :: condition_number
real, parameter :: condition_number_warning_limit = sqrt(1.0/epsilon(condition_number))
character(len = 100) :: condition_number_message
!Exit early if possible
if(self%is_empty) return
if(.not.self%is_local) return
!Allocate temporary array
allocate(tmp_arr(self%ncol,self%nrow))
!First have to pull all row level data
call self%pull_rows_to_arr(tmp_arr)
!It may be better to hook into the GK_VERBOSITY level system
!rather than relying on the debug flag.
if (self%is_head) then
condition_number = l2norm(tmp_arr)
end if
!Now invert in place
if (self%is_head) then
call invert_serial(tmp_arr, self%nrow)
tmp_arr=transpose(tmp_arr)
end if
!It would be quite nice if we could do the following in a non-blocking manner
!but this would require us to either keep the potentially large tmp_arr allocated
!for all supercells or for us to be able to do this "in-place" (i.e. we use something
!like non-blocking scatterv rather than broadcast).
if (.not. self%is_all_local) then
if (self%sc_sub_pd%nproc > 0) call broadcast_sub(tmp_arr, &
self%head_iproc_pd, self%sc_sub_pd%id)
end if
if (self%is_head) then
condition_number = condition_number * l2norm(tmp_arr)
write(condition_number_message,'("Condition number for response with ik = ",I0,", is = ",I0," and dt = ",E14.6E3," is ",E14.6E3)') self%ik_ind, self%is_ind, code_dt, condition_number
if (condition_number > condition_number_warning_limit) then
write(error_unit(), '("Warning: ",A)') trim(condition_number_message)
end if
!In debug mode always display the condition number message to screen
if (debug) then
write(*,'(A)') trim(condition_number_message)
end if
end if
!Now push back to row level
call self%push_arr_to_rows(tmp_arr)
!Free memory
deallocate(tmp_arr)
end subroutine sc_invert_local