BasicOps

Fortran module BasicOps: October 2016 (dj)

Setting kind of floating point numbers, complex unit and initialization of the seed for the random number generator. Contains as well sorting algorithms.

Authors

• 13. 12. Wall

• 4. Jaschke

Details

The following public subroutines / functions are defined for real / complex arrays if applicable.

procedure	include.f90	mpi.f90
ascending_hsort	X
descending_hsort	X
findtagindex	X
get_random_number	X
par_finalize		X
par_init		X
prop_error	X
raise_error	X
randomize	X
seed_init	X

BasicOps_f90.seed_init()

fortran-subroutine - March 2018 (dj) Set random seed using either intrinsic method or preparing for LAPACK’s DLARNV.

Arguments

rankINTEGER, OPTIONAL, in

Initialize seed with rank.

qtidINTEGER, OPTIONAL, in

Initialize seed with ID of quantum trajectory.

Source Code

show / hide f90 code

  subroutine seed_init(rank, qtid, distr)
      integer, intent(in), optional :: rank, qtid, distr

      if(use_dlarnv) then
          call seed_init_dlarnv(rank, qtid)
      else
          call seed_init_intrinsic(rank, qtid)
      end if
  end subroutine seed_init

BasicOps_f90.seed_init_dlarnv()

fortran-subroutine - March 2018 (dj) Set random seed for LAPACK’s DLARNV.

Arguments

rankINTEGER, OPTIONAL, in

Initialize seed with rank.

qtidINTEGER, OPTIONAL, in

Initialize seed with ID of quantum trajectory.

Source Code

show / hide f90 code

  subroutine seed_init_dlarnv(rank, qtid)
      integer, intent(in), optional :: rank, qtid

      ! Local variables
      ! ---------------

      ! Clock time
      integer :: clock

      ! Set default seed
      lapackseed = [11, 13, 17, 19]

      if(present(rank) .and. present(qtid)) then
          ! User specific seed - no clock time
          lapackseed(1) = rank
          lapackseed(2) = qtid
      elseif(present(qtid)) then
          ! Use qtid - no clock time
          lapackseed(2) = qtid
      elseif(present(rank)) then
          ! Use rank - no clock time
          lapackseed(1) = rank
      else
          ! Use clock time
          call system_clock(count=clock)
          lapackseed(3) = clock
      end if

  end subroutine seed_init_dlarnv

BasicOps_f90.seed_init_intrinsic()

fortran-subroutine - August 2016 (updated dj) Randomly seed the random number generator using the intrinsic system clock. Optional argument rank takes in the rank of a given processor to ensure no correlation of the random numbers on different processors.

Arguments

rankINTEGER, OPTIONAL, in

For parallel implementations when random number initialization should be different for each process in openmp or mpi or for quantum trajectories (QT). If rank is present, random numbers are reproducable and do not depend on system time!

qtidINTEGER, OPTIONAL, in

Specifies ID for quantum trajectory and allows to get a set of reproducable trajectories for the Lindblad master equation.

Source Code

show / hide f90 code

  subroutine seed_init_intrinsic(rank, qtid)
      integer, intent(in), optional :: rank, qtid

      ! Local variables
      ! ---------------

      ! Size of seed-returned by intrinsic procedure
      integer :: seedsize

      ! Clock time
      integer :: clock

      ! Seed passed to intrinsic procedure
      integer, dimension(:), allocatable :: seed

      ! dummy integer
      integer :: ii

      ! Ask for the size of the seed
      call random_seed(size=seedsize)

      ! Allocate the seed using the given seed size
      allocate(seed(seedsize))

      if(present(rank) .and. present(qtid)) then
          ! User specified seed - do not interfere with clock time
          do ii = 1, seedsize, 2
              seed(ii) = (rank + 1) * 17 * ii
          end do

          do ii = 2, seedsize, 2
              seed(ii) = (qtid + 1) * 19 * ii
          end do

      elseif(present(qtid)) then
          ! User specified seed - do not interfere with clock time
          do ii = 1, seedsize
              seed(ii) = (qtid + 1) * 19 * ii
          end do

      elseif(present(rank)) then
          ! User specified seed - do not interfere with clock time
          do ii = 1, seedsize
              seed(ii) = (rank + 1) * 17 * ii
          end do

      else
          ! Use clock time to initialize random seed
          call system_clock(count=clock)

          do ii = 1, seedsize
              seed(ii) = clock + 17 * (ii-1)
          end do

          ! 12 element seed for ubuntu system
          !seed = [1, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]
      end if

      ! Seed the random number generator using the given seed
      call random_seed(put=seed)

      deallocate(seed)

  end subroutine seed_init_intrinsic

BasicOps_f90.get_random_number_real()

fortran-subroutine - March 2018 (dj) Generate a single real-valued random number.

Arguments

scREAL, out

On exit, random number generated.

distrINTEGER, OPTIONAL, in

Choose distribution, only used for DLARNV: 1: U(0, 1), 2: U(-1, 1), 3: N(0, 1). Default to U(0, 1).

Source Code

show / hide f90 code

  subroutine get_random_number_real(sc, distr, errst)
      real(KIND=rKind), intent(out) :: sc
      integer, intent(in), optional :: distr
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! Duplet for distribution (optional argument)
      integer :: dis

      !if(present(errst)) errst = 0

      if(use_dlarnv) then
          dis = 1
          if(present(distr)) dis = distr

          call dlarnv(dis, lapackseed, 1, sc)
      else
          call random_number(sc)

          if(present(distr)) then
              if(distr == 1) then
                  ! Pass - default is U(0, 1)
              elseif(distr == 2) then
                  ! Convert to U(-1, 1)
                  sc = 2.0_rKind * sc - 1.0_rKind
              else
                  errst = raise_error('get_random_number_real: '//&
                                      ' distr argument not valid.', 99, &
                                      'BasicOps_include.f90:232', errst=errst)
                  return
              end if
          end if
      end if

  end subroutine get_random_number_real

BasicOps_f90.get_random_number_array()

fortran-subroutine - March 2018 (dj) Generate a vector of real-valued random numbers.

Arguments

vecREAL, out

On exit, vector with random number generated.

distrINTEGER, OPTIONAL, in

Choose distribution, only used for DLARNV: 1: U(0, 1), 2: U(-1, 1), 3: N(0, 1). Default to U(0, 1).

Source Code

show / hide f90 code

  subroutine get_random_number_array(vec, distr, errst)
      real(KIND=rKind), dimension(:), intent(out) :: vec
      integer, intent(in), optional :: distr
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! Duplet for distribution (optional argument)
      integer :: dis

      !if(present(errst)) errst = 0

      if(use_dlarnv) then
          dis = 1
          if(present(distr)) dis = distr

          call dlarnv(dis, lapackseed, size(vec, 1), vec)
      else
          call random_number(vec)

          if(present(distr)) then
              if(distr == 1) then
                  ! Pass - default is U(0, 1)
              elseif(distr == 2) then
                  ! Convert to U(-1, 1)
                  vec = 2.0_rKind * vec - 1.0_rKind
              else
                  errst = raise_error('get_random_number_array: '//&
                                      ' distr argument not valid.', 99, &
                                      'BasicOps_include.f90:294', errst=errst)
                  return
              end if
          end if
      end if

  end subroutine get_random_number_array

BasicOps_f90.randomize_real()

fortran-subroutine - Randomize the elements of an array in the range [-1, 1]. March 2016 (update dj)

Arguments

arrREAL(*), out

Fill with random numbers from -1 to 1.

dimINTEGER, in

Dimension of the array arr.

Source Code

show / hide f90 code

  subroutine randomize_real(arr, dim)
      integer, intent(in) :: dim
      real(KIND=rKind), dimension(dim), intent(inout) :: arr

      call get_random_number(arr)
      arr = 2.0_rKind * arr - 1.0_rKind

  end subroutine randomize_real

BasicOps_f90.randomize_complex()

fortran-subroutine - Randomize the elements of an array in the range [-1-i, 1+i]. March 2016 (update dj)

Arguments

arrCOMPLEX(*), out

Fill with random numbers from -1-i to 1+i.

dimINTEGER, in

Dimension of the array arr.

Source Code

show / hide f90 code

  subroutine randomize_complex(arr, dim)
      integer, intent(in) :: dim
      complex(KIND=rKind), dimension(dim), intent(inout) :: arr

      ! Local variables
      ! ---------------

      ! array for storing real random numbers
      real(KIND=rKind), dimension(:), allocatable :: rtmp

      allocate(rtmp(dim))

      call get_random_number(rtmp)
      arr = 2.0_rKind * rtmp - 1.0_rKind

      call get_random_number(rtmp)
      arr = arr + eye * (2.0_rKind * rtmp - 1.0_rKind)

      deallocate(rtmp)

  end subroutine randomize_complex

BasicOps_f90.descending_hsort()

fortran-subroutine - ?? (mlw) Sort array in descending order using the heapsort algorithm.

Arguments

arrREAL(*), in

Array to be sorted. Unchanged on exit.

indexoutINTEGER(*), out

Contains the sorted indexing. e.g the biggest element is accessed as arr(indexout(1)).

Source Code

show / hide f90 code

  subroutine descending_hsort(arr, indexout)
      real(KIND=rKind), intent(in) :: arr(:)
      integer, intent(out) :: indexout(:)

      ! Local variables
      ! ---------------

      ! for looping / indexing
      integer :: ii, jj, mm, ir, ix

      ! dimension of array to be sorted
      integer :: nn

      ! array to store [1, 2, ..., size(arr)], then updated when sorting
      integer :: indx(size(arr))

      ! temporary real
      real(KIND=rKind) :: temp

      ii = 1
      nn = size(arr)
      indx = (/ (ii, ii = 1, nn) /)
      mm = nn / 2 + 1
      ir = nn

      ! Fast return
      if(nn < 2) then
         indexout(:nn) = indx
         return
      end if

      do while(ir /= 1)
          if (mm > 1) then
              mm = mm - 1
              ix = indx(mm)
              temp = arr(ix)
          else
              ix = indx(ir)
              temp = arr(ix)
              indx(ir) = indx(1)
              ir = ir - 1
              if(ir == 1) then
                  indx(1) = ix
                  exit
              end if
          end if
          ii = mm
          jj = mm + mm

          do while(jj .le. ir)
              if(jj .le. ir) then
                  if(jj < ir) then
                      if(arr(indx(jj)) .lt. arr(indx(jj+1))) then
                          jj = jj + 1
                      end if
                  end if
                  if(temp < arr(indx(jj))) then
                      indx(ii) = indx(jj)
                      ii = jj
                      jj = jj + jj
                  else
                      jj = ir + 1
                  end if
                  cycle
              end if
              exit
          end do
          indx(ii) = ix
      end do

      forall(ii=1:size(indx))
          indexout(ii) = indx(size(indx)-ii+1)
      end forall

  end subroutine descending_hsort

BasicOps_f90.ascending_hsort_without_degeneracy()

fortran-subroutine - ?? (mlw) Sort A in ascending order using the heapsort algorithm

Arguments

arrREAL(*), in

Array to be sorted. Unchanged on exit.

indexoutINTEGER(*), out

Contains the sorted indexing. e.g the smalest element is accessed as arr(indexout(1)).

Source Code

show / hide f90 code

  subroutine ascending_hsort_without_degeneracy(arr, indexout)
      real(KIND=rKind), intent(in) :: arr(:)
      integer, intent(out) :: indexout(:)

      ! Local variables
      ! ---------------

      ! for looping / indexing
      integer :: ii, jj, mm, ir, ix

      ! dimension of array to be sorted
      integer :: nn

      ! array to store [1, 2, ..., size(arr)], then updated when sorting
      integer :: indx(size(indexout))

      ! temporary real
      real(KIND=rKind) :: temp

      ii = 1
      nn = size(indexout)
      indx = (/ (ii, ii = 1, nn) /)
      mm = nn / 2 + 1
      ir = nn

      ! Fast return
      if(nn < 2) then
          indexout(:nn) = indx
          return
      end if

      do while(ir /= 1)
          if(mm > 1) then
              mm = mm - 1
              ix = indx(mm)
              temp = arr(ix)
          else
              ix = indx(ir)
              temp = arr(ix)
              indx(ir) = indx(1)
              ir = ir - 1
              if(ir == 1) then
                  indx(1) = ix
                  exit
              end if
          end if
          ii = mm
          jj = mm + mm

          do while(jj .le. ir)
              if(jj .le. ir) then
                  if(jj < ir) then
                      if(arr(indx(jj)) < arr(indx(jj+1))) then
                          jj = jj + 1
                      end if
                  end if
                  if(temp < arr(indx(jj))) then
                      indx(ii) = indx(jj)
                      ii = jj
                      jj = jj + jj
                  else
                      jj = ir + 1
                  end if
                  cycle
              end if
              exit
          end do
          indx(ii) = ix
      end do

      forall(ii = 1:nn)
          indexout(ii) = indx(ii)
      end forall

  end subroutine ascending_hsort_without_degeneracy

BasicOps_f90.ascending_hsort_with_degeneracy()

fortran-subroutine - September 2017 (dj, updated) Given a hash which is possibly degnerate, construct a hash table “littlehash” of the unique values, and return the number of unique values nunique along with their degeneracies and positions in the full hash table. i.e. the value littlehash(j) occurs for hash(ind(deg(j):deg(j+1))). Assumes littlehash and ind are SIZE(hash) in length and deg is SIZE(hash)+1 in length.

Arguments

hashINTEGER(*), ???

List of degenerated hashes to be sorted.

littlehashINTEGER(*), out

length is SIZE(hash). Containing the list of hashes without degeneracy. Not completly set (see nunique)

indINTEGER(*), out

length is SIZE(hash). i-th entry points to (one of) the i-th biggest value of hash.

nuniqueINTEGER, out

Number of different hashes, therefore littlehash is set from elements 1 to nunique.

degINTEGER(*), out

deg has length SIZE(hash) + 1 The first index of a degenerate i-th hash can be found with ind(deg(i) + 1) up to the last with ind(deg(i+1)) (Other description is most probably wrong, since deg(1) = 0)

Source Code

show / hide f90 code

  subroutine ascending_hsort_with_degeneracy(hash, littlehash, ind, &
                                             nunique, deg, errst)

      REAL(KIND=rKind) :: hash(:),littlehash(:)
      INTEGER :: ind(:), deg(:)
      INTEGER :: nunique
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! size of array to be sorted
      integer :: nn

      !if(present(errst)) errst = 0

      nn = size(hash)
      call ascending_hsort(hash, ind)

      ! Fast return for empty array
      if(nn == 0) then
          nunique = 0
          return
      end if

      nunique = 1
      littlehash(nunique) = hash(ind(1))
      deg(nunique) = 0

      do ii = 2, nn
          if(hash(ind(ii)) /= hash(ind(ii-1))) then
              nunique = nunique + 1
              littlehash(nunique) = hash(ind(ii))
              deg(nunique) = ii - 1
          end if
      end do

      !if(nunique + 1 > size(deg, 1)) then
      !    errst = raise_error('ascending_hsort_with_degeneracy'//&
      !                        ' : deg out of bounds.', 99, &
      !                        errst=errst)
      !    return
      !end if

      deg(nunique + 1) = nn

  end subroutine ascending_hsort_with_degeneracy

BasicOps_f90.FindTagIndex_real()

fortran-function - ?? (mlw) Find the index of element tag in the sorted list (or hash table) hash using binary search. If the element is not found then return -1.

Arguments

tagREAL, in

Element to be located in list.

hashREAL(*), in

Sorted list searching for an entry matching tag. Must be sorted ascending.

Source Code

show / hide f90 code

  function FindTagIndex_real(tag, hash) result(tagidx)
      real(KIND=rKind), intent(in) :: tag
      real(KIND=rKind), intent(in) :: hash(:)
      integer :: tagidx

      ! Local variables
      ! ---------------

      ! start, mid, and end point of the search interval
      integer :: min, max, mid

      max = size(hash)

      ! Fast return (empty array)
      if(max == 0) then
          tagidx = -1
          return
      end if

      min = 1
      mid = (min + max) / 2

      ! Fast return
      if((tag < hash(1)) .or. (tag > hash(max))) then
          tagidx = -1
          return
      end if

      do while(min <= max)
          mid = (min + max) / 2
          if(tag > hash(mid)) then
              min = mid + 1
          else if(tag < hash(mid)) then
              max = mid - 1
          else
              tagidx = mid
              return
              exit
          end if
      end do
      tagidx = -1

  end function FindTagIndex_real

BasicOps_f90.FindTagIndex_int()

fortran-function - ?? (mlw) Find the index of element tag in the sorted list (or hash table) hash using binary search. If the element is not found then return -1.

Arguments

tagINTEGER, in

Element to be located in list.

hashINTEGER(*), in

Sorted list searching for an entry matching tag. Must be sorted ascending.

Source Code

show / hide f90 code

  function FindTagIndex_int(tag, hash) result(tagidx)
      integer, intent(in) :: tag
      integer, intent(in) :: hash(:)
      integer :: tagidx

      ! Local variables
      ! ---------------

      ! start, mid, and end point of the search interval
      integer :: min, max, mid

      max = size(hash)

      ! Fast return (empty array)
      if(max == 0) then
          tagidx = -1
          return
      end if

      min = 1
      mid = (min + max) / 2

      ! Fast return
      if((tag < hash(1)) .or. (tag > hash(max))) then
          tagidx = -1
          return
      end if

      do while(min <= max)
          mid = (min + max) / 2
          if(tag > hash(mid)) then
              min = mid + 1
          else if(tag < hash(mid)) then
              max = mid - 1
          else
              tagidx = mid
              return
              exit
          end if
      end do
      tagidx = -1

  end function FindTagIndex_int

BasicOps_f90.raise_error()

fortran-function - October 2016 (dj) Call this function in case an error occurs. It will either stop the program or propagate the error message through the code.

Arguments

msgCHARACTER(*), in

Print this message in case of an error in order facilitate tracking the error.

statusINTEGER, in

error status to be set.

file_lineCHARACTER(*), in

String containing the file name and the line number where the error is raised.

errstINTEGER, OPTIONAL, in

Status of the error distinguishes between propagating error when errst is present. If not present, stop is raised directly.

Source Code

show / hide f90 code

  function raise_error(msg, status, file_line, errst) result(res)
      character(len=*), intent(in) :: msg
      integer, intent(in) :: status
      character(len=*), intent(in), optional :: file_line
      integer, intent(in), optional :: errst
      integer :: res

      if(.not. present(errst)) then
          if(present(file_line)) then
              write(elog, *) file_line//': '//msg
          else
              write(elog, *) msg
          end if
          stop 'raise_error'
      end if

      if(present(file_line)) then
          write(elog, *) file_line//': '//msg
      else
          write(elog, *) msg
      end if
      res = status

  end function raise_error

BasicOps_f90.prop_error()

fortran-function - October 2016 (dj) Handle posting of error messages in regard to the error status from a previous call.

Arguments

msgCHARACTER(*), in

Print this message in case of an error in order facilitate tracking the error.

file_lineCHARACTER(*), in

String containing the file name and the line number where the error is propagated.

errstINTEGER, OPTIONAL, in

Status of the error as set from the last call. If not present, it defaults to the no-error case.

Source Code

show / hide f90 code

  function prop_error(msg, file_line, errst) result(has_error)
      character(len=*), intent(in) :: msg
      character(len=*), intent(in), optional :: file_line
      integer, intent(in), optional :: errst
      logical :: has_error

      if(present(errst)) then
          if(errst == 0) then
              has_error = .false.
          else
              if(present(file_line)) then
                  write(elog, *) file_line//': '//msg
              else
                  write(elog, *) msg
              end if
              has_error = .true.
          end if
      else
          has_error = .false.
      end if

  end function prop_error