sc_invert_local Subroutine

private subroutine sc_invert_local(self)

A routine to invert the field matrix locally

Type Bound

supercell_type

Arguments

Type IntentOptional Attributes Name
class(supercell_type), intent(inout) :: self

Contents

Source Code


Source Code

  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