"""
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**
* M. L. Wall
* D. 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 | |
+-------------------+-------------+---------+
"""
[docs]def seed_init():
"""
fortran-subroutine - March 2018 (dj)
Set random seed using either intrinsic method or preparing for
LAPACK's DLARNV.
**Arguments**
rank : INTEGER, OPTIONAL, in
Initialize seed with rank.
qtid : INTEGER, OPTIONAL, in
Initialize seed with ID of quantum trajectory.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def seed_init_dlarnv():
"""
fortran-subroutine - March 2018 (dj)
Set random seed for LAPACK's DLARNV.
**Arguments**
rank : INTEGER, OPTIONAL, in
Initialize seed with rank.
qtid : INTEGER, OPTIONAL, in
Initialize seed with ID of quantum trajectory.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def 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**
rank : INTEGER, 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!
qtid : INTEGER, OPTIONAL, in
Specifies ID for quantum trajectory and allows to get a set of
reproducable trajectories for the Lindblad master equation.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def get_random_number_real():
"""
fortran-subroutine - March 2018 (dj)
Generate a single real-valued random number.
**Arguments**
sc : REAL, out
On exit, random number generated.
distr : INTEGER, 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**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def get_random_number_array():
"""
fortran-subroutine - March 2018 (dj)
Generate a vector of real-valued random numbers.
**Arguments**
vec : REAL, out
On exit, vector with random number generated.
distr : INTEGER, 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**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def randomize_real():
"""
fortran-subroutine - Randomize the elements of an array in the range [-1, 1].
March 2016 (update dj)
**Arguments**
arr : REAL(\*), out
Fill with random numbers from -1 to 1.
dim : INTEGER, in
Dimension of the array `arr`.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def randomize_complex():
"""
fortran-subroutine - Randomize the elements of an array in the range [-1-i, 1+i].
March 2016 (update dj)
**Arguments**
arr : COMPLEX(\*), out
Fill with random numbers from -1-i to 1+i.
dim : INTEGER, in
Dimension of the array `arr`.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def descending_hsort():
"""
fortran-subroutine - ?? (mlw)
Sort array in descending order using the heapsort algorithm.
**Arguments**
arr : REAL(\*), in
Array to be sorted. Unchanged on exit.
indexout : INTEGER(\*), out
Contains the sorted indexing. e.g the biggest element is accessed
as arr(indexout(1)).
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def ascending_hsort_without_degeneracy():
"""
fortran-subroutine - ?? (mlw)
Sort A in ascending order using the heapsort algorithm
**Arguments**
arr : REAL(\*), in
Array to be sorted. Unchanged on exit.
indexout : INTEGER(\*), out
Contains the sorted indexing. e.g the smalest element is accessed
as arr(indexout(1)).
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def 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**
hash : INTEGER(*), ???
List of degenerated hashes to be sorted.
littlehash : INTEGER(*), out
length is SIZE(hash). Containing the list of hashes without
degeneracy. Not completly set (see `nunique`)
ind : INTEGER(*), out
length is SIZE(hash). i-th entry points to (one of) the i-th
biggest value of hash.
nunique : INTEGER, out
Number of different hashes, therefore `littlehash` is set from
elements 1 to `nunique`.
deg : INTEGER(*), 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**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def 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**
tag : REAL, in
Element to be located in list.
hash : REAL(\*), in
Sorted list searching for an entry matching tag. Must be sorted
ascending.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def 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**
tag : INTEGER, in
Element to be located in list.
hash : INTEGER(\*), in
Sorted list searching for an entry matching tag. Must be sorted
ascending.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def 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**
msg : CHARACTER(\*), in
Print this message in case of an error in order facilitate tracking
the error.
status : INTEGER, in
error status to be set.
file_line : CHARACTER(\*), in
String containing the file name and the line number where the error
is raised.
errst : INTEGER, OPTIONAL, in
Status of the error distinguishes between propagating error
when errst is present. If not present, stop is raised directly.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return
[docs]def prop_error():
"""
fortran-function - October 2016 (dj)
Handle posting of error messages in regard to the error status
from a previous call.
**Arguments**
msg : CHARACTER(\*), in
Print this message in case of an error in order facilitate tracking
the error.
file_line : CHARACTER(\*), in
String containing the file name and the line number where the error
is propagated.
errst : INTEGER, OPTIONAL, in
Status of the error as set from the last call.
If not present, it defaults to the no-error case.
**Source Code**
.. hidden-code-block:: fortran
:label: 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
"""
return