From d12c644a4cfa916c22ea64e58102e04c954e36ea Mon Sep 17 00:00:00 2001 From: vehre Date: Sat, 20 Jun 2015 15:05:23 +0200 Subject: [PATCH] Added files and readme. --- ArrayUtils.f90 | 41 + FileUtils.f90 | 1507 ++++++++++++++++++++++++++++++ IniObjects.f90 | 927 +++++++++++++++++++ Makefile | 74 ++ MatrixUtils.f90 | 2317 +++++++++++++++++++++++++++++++++++++++++++++++ MiscUtils.f90 | 205 +++++ MpiUtils.f90 | 172 ++++ ObjectLists.f90 | 960 ++++++++++++++++++++ RandUtils.f90 | 377 ++++++++ RangeUtils.f90 | 459 ++++++++++ Readme | 47 + StringUtils.f90 | 358 ++++++++ 12 files changed, 7444 insertions(+) create mode 100644 ArrayUtils.f90 create mode 100644 FileUtils.f90 create mode 100644 IniObjects.f90 create mode 100644 Makefile create mode 100644 MatrixUtils.f90 create mode 100644 MiscUtils.f90 create mode 100644 MpiUtils.f90 create mode 100644 ObjectLists.f90 create mode 100644 RandUtils.f90 create mode 100644 RangeUtils.f90 create mode 100644 Readme create mode 100644 StringUtils.f90 diff --git a/ArrayUtils.f90 b/ArrayUtils.f90 new file mode 100644 index 0000000..79c5720 --- /dev/null +++ b/ArrayUtils.f90 @@ -0,0 +1,41 @@ + module ArrayUtils + implicit none + + contains + + function IndexOf(aval,arr, n) + integer, intent(in) :: n, arr(n), aval + integer IndexOf, i + + do i=1,n + if (arr(i)==aval) then + IndexOf= i + return + end if + end do + IndexOf = 0 + + end function IndexOf + + function MaxIndex(arr, n) + integer, intent(in) :: n + real, intent(in) :: arr(n) + integer locs(1:1), MaxIndex + + locs = maxloc(arr(1:n)) + MaxIndex = locs(1) + + end function MaxIndex + + + function MinIndex(arr, n) + integer, intent(in) :: n + real, intent(in) :: arr(n) + integer locs(1:1), MinIndex + + locs = minloc(arr(1:n)) + MinIndex = locs(1) + + end function MinIndex + + end module ArrayUtils diff --git a/FileUtils.f90 b/FileUtils.f90 new file mode 100644 index 0000000..f9a0bc6 --- /dev/null +++ b/FileUtils.f90 @@ -0,0 +1,1507 @@ + module FileUtils + use MpiUtils + use MiscUtils + use StringUtils + use, intrinsic :: ISO_FORTRAN_ENV, only : INT64 + implicit none + !Utils using F2008 features + private + + integer, parameter :: File_size_int = KIND(INT64) + + Type :: TFileStream + integer :: unit = 0 + logical :: ReadOnly = .true. + character(LEN=:), allocatable :: mode + character(LEN=:), allocatable :: access !sequential (fortran or text) or stream (binary, with file position) + character(LEN=:), allocatable :: FileName + contains + procedure :: SetDefaultModes + procedure :: Open => TFileStream_Open + procedure :: OpenFile => TFileStream_OpenFile + procedure :: Flush => TFileStream_Flush + procedure :: CreateFile + procedure :: CreateOpenFile + procedure :: Position => TFileStream_Position + procedure :: Size => TFileStream_Size + procedure :: Close => TFileStream_Close + procedure :: Rewind => TFileStream_Rewind + procedure :: Error + procedure :: Opened + procedure :: WriteTrim + procedure :: CheckOpen + procedure, private :: ReadStringFunc + procedure, private :: ReadStringSub + procedure, private :: ReadItemSub + procedure, private :: ReadItemFunc + procedure, private :: ReadItems + procedure, private :: ReadArray + procedure, private :: ReadArray2 + procedure, private :: ReadArray2Func + procedure, private :: ReadArrayFunc + procedure, private :: WriteItemSub + procedure, private :: WriteArray + procedure, private :: WriteOneAndArray + procedure, private :: WriteArray2 + procedure, private :: WriteItems + procedure, private :: WriteSizedArray1 + procedure, private :: WriteSizedArray2 + procedure, private :: ReadSizedArray_R + procedure, private :: ReadSizedArray_D + procedure, private :: ReadSizedArray_I + procedure, private :: ReadSizedArray2_I + procedure, private :: ReadSizedArray2_D + procedure, private :: ReadSizedArray2_R + generic :: Write => WriteItems, WriteArray, WriteArray2, WriteOneAndArray + generic :: Read => ReadItems, ReadArray, ReadArray2 + generic :: ReadItem => ReadItemFunc, ReadArrayFunc, ReadArray2Func + generic :: ReadString => ReadStringSub + generic :: ReadStringItem => ReadStringFunc + generic :: WriteSizedArray => WriteSizedArray1,WriteSizedArray2 + generic :: ReadSizedArray => ReadSizedArray_R,ReadSizedArray_D,ReadSizedArray_I, & + & ReadSizedArray2_I, ReadSizedArray2_R, ReadSizedArray2_D + final :: TFileStream_Free + end type + + Type, extends(TFileStream) :: TBinaryFile + end type + + Type, extends(TFileStream) :: TSequentialFile + contains + procedure :: SetDefaultModes => TSequentialFile_SetDefaultModes + end type TSequentialFile + + Type, extends(TSequentialFile) :: TTextFile + character(LEN=20) :: RealFormat = '(*(E17.7))' + character(LEN=20) :: IntegerFormat = '(*(I10))' + logical :: AdvanceDefault = .true. + contains + procedure :: SetDefaultModes => TTextFile_SetDefaultModes + procedure :: ReadLine + procedure :: ReadLineSkipEmptyAndComments + procedure :: ReadNextContentLine !as above, but opens and closes file + procedure :: NewLine + procedure :: SkipLines + procedure :: Lines + procedure :: Columns + procedure :: WriteLeftAligned + procedure :: WriteItemTxt + procedure :: WriteArrayTxt + procedure :: WriteInLineItems + procedure :: WriteInLineTrim + procedure :: WriteFormat + procedure :: WriteTrim => WriteTrimTxt + procedure, private :: ReadStringSub => ReadStringTxt + procedure, private :: ReadItemTxt + procedure, private :: ReadArrayTxt + procedure, private :: WriteInLineItem + procedure, private :: WriteInLineArray + procedure, private :: DefaultAdvance + procedure, private :: WriteItem => TTextFile_WriteItem + procedure, private :: WriteArray => TTextFile_WriteArray + procedure, private :: WriteOneAndArray => TTextFile_WriteOneAndArray + procedure, private :: WriteArray2 => WriteArray2Txt + procedure, private :: ReadItemSub => ReadItemTxt + procedure, private :: ReadArray => ReadArrayTxt + procedure, private :: ReadArray2 => ReadArray2Txt + procedure, private :: WriteItems => WriteItemsTxt + generic :: WriteInLine => WriteInLineItem, WriteInLineArray + end type + + !Functions on filenames and text + !File instance below acts as a namespace + Type TFile + contains + procedure, nopass :: TxtNumberColumns + procedure, nopass :: TxtColumns + procedure, nopass :: TxtFileColumns + procedure, nopass :: TxtFileLines + procedure, nopass :: LastTopComment + procedure, nopass :: TopCommentLine + procedure, nopass :: LastLine => LastFileLine + procedure, nopass :: Size => FileSize + procedure, nopass :: Exists => FileExists + procedure, nopass :: ExtractName => ExtractFileName + procedure, nopass :: ExtractPath => ExtractFilePath + procedure, nopass :: Join => File_Join + procedure, nopass :: ChangeExt => ChangeFileExt + procedure, nopass :: CheckTrailingSlash + procedure, nopass :: IsFullPath + procedure, nopass :: ExtractExt => ExtractFileExt + procedure, nopass :: Delete => DeleteFile + procedure, nopass :: ReadTextMatrix + procedure, nopass :: ReadTextVector + procedure, nopass :: WriteTextMatrix + procedure, nopass :: WriteTextVector + procedure, nopass :: LoadTxt_1D + procedure, nopass :: LoadTxt_2D + procedure, nopass :: OpenTextFile + procedure, nopass :: CreateTextFile + procedure, nopass :: CharIsSlash + generic :: LoadTxt => LoadTxt_2D, LoadTxt_1D + generic :: SaveTxt => WriteTextMatrix, WriteTextVector + end type + + type(TFile), public, save :: File + + public TFileStream, TBinaryFile, TTextFile, TFile, File_size_int + contains + + function FileExists(aname) + character(LEN=*), intent(IN) :: aname + logical FileExists + + inquire(file=aname, exist = FileExists) + + end function FileExists + + function FileSize(name) result(fsize) + integer(file_size_int) fsize + character(LEN=*), intent(in)::name + + inquire(file=name, size=fsize) + + end function FileSize + + function TxtNumberColumns(InLine) result(n) + character(LEN=*) :: InLine + integer n,i + logical isNum + + n=0 + isNum=.false. + do i=1, len_trim(InLIne) + if (verify(InLine(i:i),'-+eE.0123456789') == 0) then + if (.not. IsNum) n=n+1 + IsNum=.true. + else + IsNum=.false. + end if + end do + + end function TxtNumberColumns + + function TxtColumns(InLine) result(n) + character(LEN=*) :: InLine + integer n,i + logical isNum + + n=0 + isNum=.false. + do i=1, len_trim(InLine) + if (InLine(i:i) > char(32)) then + if (.not. IsNum) n=n+1 + IsNum=.true. + else + IsNum=.false. + end if + end do + + end function TxtColumns + + + subroutine SetDefaultModes(this) + class(TFileStream) :: this + + if (.not. allocated(this%mode)) this%mode = 'unformatted' + if (.not. allocated(this%access)) this%access = 'stream' + + end subroutine SetDefaultModes + + subroutine TSequentialFile_SetDefaultModes(this) + class(TSequentialFile) :: this + + if (.not. allocated(this%access)) this%access= 'sequential' + call this%TFileStream%SetDefaultModes + + end subroutine TSequentialFile_SetDefaultModes + + subroutine TTextFile_SetDefaultModes(this) + class(TTextFile) :: this + + if (.not. allocated(this%mode)) this%mode = 'formatted' + call this%TSequentialFile%SetDefaultModes + + end subroutine TTextFile_SetDefaultModes + + + subroutine CheckOpen(this, forWrite) + class(TFileStream) :: this + logical, intent(in), optional :: forWrite + + if (this%unit/=0) return + if (DefaultFalse(forWrite) .and. this%ReadOnly) then + call this%Error('File not open for write') + else + call this%Error('File not opened') + end if + + end subroutine + + function Opened(this) + class(TFileStream) :: this + logical Opened + Opened = this%unit /=0 + end function Opened + + subroutine Error(this, msg, errormsg) + class(TFileStream) :: this + character(LEN=*), intent(IN) :: msg + character(LEN=*), intent(IN), optional :: errormsg + character(LEN=:), allocatable :: Filename + + if (.not. allocated(this%FileName)) then + FileName = '(no filename set)' + else + FileName = this%FileName + end if + if (present(errormsg)) then + call MpiStop(trim(errormsg)//' : '//FileName ) + else + call MpiStop(trim(msg)//' : '// FileName ) + end if + + end subroutine + + subroutine TFileStream_Close(this) + class(TFileStream), intent(inout) :: this + if (this%unit/=0) close(this%unit) + this%unit=0 + end subroutine TFileStream_Close + +#ifdef __GFORTRAN__ + impure elemental subroutine TFileStream_Free(this) +#else + subroutine TFileStream_Free(this) +#endif + Type(TFileStream), intent(inout) :: this + + call this%Close() + end subroutine TFileStream_Free + + subroutine TFileStream_Flush(this) + class(TFileStream) :: this + + call this%CheckOpen + flush(this%unit) + + end subroutine TFileStream_Flush + + subroutine TFileStream_Rewind(this) + class(TFileStream) :: this + + if (this%Opened()) rewind(this%unit) + + end subroutine TFileStream_Rewind + + function TFileStream_Position(this) + class(TFileStream) :: this + integer(file_size_int) TFileStream_Position + + call this%CheckOpen + if (this%access /= 'stream') call this%Error('Position requires access=stream') + inquire(this%unit, pos=TFileStream_Position) + + end function TFileStream_Position + + + function TFileStream_Size(this) + class(TFileStream) :: this + integer(file_size_int) TFileStream_Size + + if (this%Opened()) then + inquire(this%unit, size=TFileStream_Size) + else if (allocated(this%FileName)) then + TFileStream_Size = File%Size(this%FileName) + else + call this%Error('File not defined for size') + end if + + end function TFileStream_Size + + subroutine TFileStream_Open(this, aname, errormsg, status) + class(TFileStream) :: this + character(LEN=*), intent(IN) :: aname + character(LEN=*), intent(IN), optional :: errormsg + integer, intent(out), optional :: status + + call this%OpenFile(aname,errormsg=errormsg, status=status) + + end subroutine TFileStream_Open + + subroutine TFileStream_OpenFile(this, aname, mode, errormsg, forwrite, append, status) + class(TFileStream) :: this + character(LEN=*), intent(IN) :: aname + character(LEN=*), intent(IN), optional :: mode + character(LEN=*), intent(IN), optional :: errormsg + integer, intent(out), optional :: status + logical, intent(in), optional :: forwrite, append + character(LEN=:), allocatable :: amode, state, action, pos + integer out_status + + call this%Close() + call this%SetDefaultModes() + + this%FileName = trim(aname) + + amode = PresentDefault(this%mode, mode) + if (DefaultFalse(forwrite)) then + state = 'replace' + action = 'readwrite' + this%ReadOnly = .false. + else + state = 'old' + action = 'read' + this%ReadOnly = .true. + end if + if (DefaultFalse(append) .and. FileExists(aname)) then + pos = 'append' + state = 'old' + action = 'readwrite' + this%ReadOnly = .false. + else + pos='asis' + end if + + open(file=aname,form=amode,status=state, action=action, newunit=this%unit, & + & iostat=out_status, position =pos, access=this%access) + if (present(status)) then + status=out_status + if (out_status/=0) this%unit = 0 + else + if (out_status/=0) then + if (state == 'replace') then + call this%Error('Error creating file', errormsg) + else + call this%Error('File not found', errormsg) + end if + this%unit = 0 + end if + end if + end subroutine TFileStream_OpenFile + + subroutine CreateFile(this,aname, errormsg) + class(TFileStream) :: this + character(LEN=*), intent(IN) :: aname + character(LEN=*), intent(IN), optional :: errormsg + + call this%OpenFile(aname, errormsg=errormsg, forwrite =.true.) + + end subroutine CreateFile + + subroutine CreateOpenFile(this, aname, append, errormsg) + class(TFileStream) :: this + character(LEN=*), intent(IN) :: aname + logical, optional, intent(in) :: append + character(LEN=*), intent(IN), optional :: errormsg + + call this%OpenFile(aname, forwrite =.true., append=append, errormsg=errormsg) + + end subroutine CreateOpenFile + + + subroutine ReadStringSub(this, S, OK) + class(TFileStream) :: this + character(LEN=:), allocatable :: S + logical, optional :: OK + logical isOK + integer i, status + + call this%CheckOpen + read(this%unit, iostat=status) i + isOK = status==0 + if (isOK) then + if (allocated(S)) deallocate(S) + allocate(character(LEN=i)::S) + read(this%unit, iostat=status) S + isOK = status==0 + end if + if (present(OK)) OK = isOK + end subroutine ReadStringSub + + function ReadStringFunc(this) result(S) + class(TFileStream) :: this + character(LEN=:), allocatable :: S + call this%ReadString(S) + end function ReadStringFunc + + subroutine ReadItemSub(this, R, OK) + class(TFileStream) :: this + class(*), intent(out) :: R + logical, optional :: OK + logical res + integer status + + call this%CheckOpen + select type(R) + type is (real) + read(this%unit, iostat=status) R + type is (double precision) + read(this%unit, iostat=status) R + type is (integer) + read(this%unit, iostat=status) R + type is (logical) + read(this%unit, iostat=status) R + type is (character(LEN=*)) + read(this%unit, iostat=status) R + class default + call this%Error('Unknown type to read') + end select + + res = status==0 + if (status/=0 .and. (.not. IS_IOSTAT_END(status) .or. .not. present(OK))) & + & call this%Error('Error reading item') + if (present(OK)) OK = res + end subroutine ReadItemSub + + + function ReadItemFunc(this, R) result(res) + class(TFileStream) :: this + class(*), intent(out) :: R + logical :: res + + call this%ReadItemSub(R, res) + end function ReadItemFunc + + subroutine ReadArray(this, R, n, OK) + class(TFileStream) :: this + class(*) :: R(1:) + integer, intent(in), optional :: n + logical, optional :: OK + integer status + + call this%CheckOpen + select type(R) + type is (real) + read(this%unit, iostat=status) R(1:PresentDefault(size(R),n)) + type is (double precision) + read(this%unit, iostat=status) R(1:PresentDefault(size(R),n)) + type is (integer) + read(this%unit, iostat=status) R(1:PresentDefault(size(R),n)) + type is (logical) + read(this%unit, iostat=status) R(1:PresentDefault(size(R),n)) + class default + call this%Error('Unknown type to read') + end select + if (status/=0 .and. (.not. IS_IOSTAT_END(status) .or. .not. present(OK))) & + & call this%Error('Error reading item') + if (present(OK)) OK = status==0 + end subroutine ReadArray + + + function ReadArrayFunc(this, R, n) result(res) + class(TFileStream) :: this + class(*) :: R(1:) + integer, intent(in), optional :: n + logical :: res + + call this%ReadArray(R,n,res) + + end function ReadArrayFunc + + + function ReadArray2Func(this, R) result(res) + class(TFileStream) :: this + class(*) :: R(:,:) + logical :: res + call this%ReadArray2(R,res) + end function ReadArray2Func + + + subroutine ReadArray2(this, R, OK) + class(TFileStream) :: this + class(*) :: R(:,:) + logical, optional :: OK + integer status + + call this%CheckOpen + select type(R) + type is (real) + read(this%unit, iostat=status) R + type is (double precision) + read(this%unit, iostat=status) R + type is (integer) + read(this%unit, iostat=status) R + type is (logical) + read(this%unit, iostat=status) R + class default + call this%Error('Unknown type to read') + end select + if (status/=0 .and. (.not. IS_IOSTAT_END(status) .or. .not. present(OK))) & + & call this%Error('Error reading item') + if (present(OK)) OK = status==0 + end subroutine ReadArray2 + + subroutine ReadSizedArray_R(this, R) + class(TFileStream) :: this + Real, allocatable :: R(:) + integer sz + + call this%Read(sz) + if (allocated(R)) deallocate(R) + allocate(R(sz)) + call this%ReadArray(R) + + end subroutine ReadSizedArray_R + + + subroutine ReadSizedArray_D(this, R) + class(TFileStream) :: this + double precision, allocatable :: R(:) + integer sz + + call this%Read(sz) + if (allocated(R)) deallocate(R) + allocate(R(sz)) + call this%ReadArray(R) + + end subroutine ReadSizedArray_D + + subroutine ReadSizedArray_I(this, R) + class(TFileStream) :: this + integer, allocatable :: R(:) + integer sz + + call this%Read(sz) + if (allocated(R)) deallocate(R) + allocate(R(sz)) + call this%ReadArray(R) + + end subroutine ReadSizedArray_I + + subroutine ReadSizedArray2_R(this, R) + class(TFileStream) :: this + real, allocatable :: R(:,:) + integer sz1, sz2 + + call this%Read(sz1) + call this%Read(sz2) + if (allocated(R)) deallocate(R) + allocate(R(sz1,sz2)) + call this%ReadArray2(R) + + end subroutine ReadSizedArray2_R + + subroutine ReadSizedArray2_D(this, R) + class(TFileStream) :: this + double precision, allocatable :: R(:,:) + integer sz1, sz2 + + call this%Read(sz1) + call this%Read(sz2) + if (allocated(R)) deallocate(R) + allocate(R(sz1,sz2)) + call this%ReadArray2(R) + + end subroutine ReadSizedArray2_D + + subroutine ReadSizedArray2_I(this, R) + class(TFileStream) :: this + integer, allocatable :: R(:,:) + integer sz1, sz2 + + call this%Read(sz1) + call this%Read(sz2) + if (allocated(R)) deallocate(R) + allocate(R(sz1,sz2)) + call this%ReadArray2(R) + + end subroutine ReadSizedArray2_I + + subroutine WriteItemSub(this, R) + class(TFileStream) :: this + class(*), intent(in) :: R + + call this%CheckOpen(.true.) + select type(R) + type is (real) + Write(this%unit) R + type is (double precision) + Write(this%unit) R + type is (integer) + Write(this%unit) R + type is (logical) + Write(this%unit) R + type is (character(LEN=*)) + Write(this%unit) len(R) + Write(this%unit) R + class default + call this%Error('Unknown type to Write') + end select + + end subroutine WriteItemSub + + subroutine WriteTrim(this, S) + class(TFileStream) :: this + character(LEN=*), intent(in) :: S + call this%WriteItemSub(trim(S)) + end subroutine WriteTrim + + subroutine WriteSizedArray1(this, R, n) + class(TFileStream) :: this + class(*), intent(in) :: R(1:) + integer, intent(in), optional :: n + integer sz + + sz=PresentDefault(size(R),n) + call this%Write(sz) + call this%WriteArray(R,n) + + end subroutine WriteSizedArray1 + + subroutine WriteSizedArray2(this, R) + class(TFileStream) :: this + class(*), intent(in) :: R(:,:) + + call this%Write(size(R,dim=1)) + call this%Write(size(R,dim=2)) + call this%WriteArray2(R) + + end subroutine WriteSizedArray2 + + subroutine WriteOneAndArray(this, I, R) + class(TFileStream) :: this + class(*), intent(in) :: I + class(*), intent(in) :: R(:) + + call this%WriteItemSub(I) + call this%WriteArray(R) + + end subroutine WriteOneAndArray + + subroutine WriteArray(this, R, n) + class(TFileStream) :: this + class(*), intent(in) :: R(1:) + integer, intent(in), optional :: n + integer sz + + sz=PresentDefault(size(R),n) + call this%CheckOpen(.true.) + select type(R) + type is (real) + Write(this%unit) R(1:sz) + type is (double precision) + Write(this%unit) R(1:sz) + type is (integer) + Write(this%unit) R(1:sz) + type is (logical) + Write(this%unit) R(1:sz) + class default + call this%Error('Unknown type to Write') + end select + + end subroutine WriteArray + + subroutine WriteArray2(this, R) + class(TFileStream) :: this + class(*), intent(in) :: R(:,:) + + call this%CheckOpen(.true.) + select type(R) + type is (real) + Write(this%unit) R + type is (double precision) + Write(this%unit) R + type is (integer) + Write(this%unit) R + type is (logical) + Write(this%unit) R + class default + call this%Error('Unknown type to Write') + end select + + end subroutine WriteArray2 + + subroutine WriteItems(this, S1, S2,S3,S4,S5,S6) + class(TFileStream) :: this + class(*), intent(in) :: S1 + class(*), intent(in), optional :: S2, S3,S4,S5,S6 + + call this%WriteItemSub(S1) + if (present(S2)) call this%WriteItemSub(S2) + if (present(S3)) call this%WriteItemSub(S3) + if (present(S4)) call this%WriteItemSub(S4) + if (present(S5)) call this%WriteItemSub(S5) + if (present(S6)) call this%WriteItemSub(S6) + + end subroutine WriteItems + + subroutine ReadItems(this, S1, S2,S3,S4,S5,S6,OK) + class(TFileStream) :: this + class(*) S1 + class(*), optional :: S2,S3,S4,S5,S6 + logical, optional :: OK + + call this%ReadItemSub(S1,OK) + if (present(S2) .and. DefaultTrue(OK)) call this%ReadItemSub(S2,OK) + if (present(S3) .and. DefaultTrue(OK)) call this%ReadItemSub(S3,OK) + if (present(S4) .and. DefaultTrue(OK)) call this%ReadItemSub(S4,OK) + if (present(S5) .and. DefaultTrue(OK)) call this%ReadItemSub(S5,OK) + if (present(S6) .and. DefaultTrue(OK)) call this%ReadItemSub(S6,OK) + + end subroutine ReadItems + + !Text unformatted files + + function ReadLine(this, InLine, trimmed) result(OK) + class(TTextFile) :: this + character(LEN=:), allocatable, optional :: InLine + logical, intent(in), optional :: trimmed + integer, parameter :: line_buf_len= 1024*4 + character(LEN=line_buf_len) :: InS + logical :: OK, set + integer status, size + + call this%CheckOpen + OK = .false. + set = .true. + do + read (this%unit,'(a)',advance='NO',iostat=status, size=size) InS + OK = .not. IS_IOSTAT_END(status) + if (.not. OK) return + if (present(InLine)) then + if (set) then + InLine = InS(1:size) + set=.false. + else + InLine = InLine // InS(1:size) + end if + end if + if (IS_IOSTAT_EOR(status)) exit + end do + if (present(trimmed) .and. present(InLine)) then + if (trimmed) InLine = trim(adjustl(InLine)) + end if + + end function ReadLine + + function ReadLineSkipEmptyAndComments(this, InLine, comment) result(OK) + class(TTextFile) :: this + logical :: OK + character(LEN=:), allocatable :: InLine + character(LEN=:), allocatable, optional, intent(inout) :: comment + + if (present(comment)) then + if (.not. allocated(comment)) comment='' + end if + do + OK = this%ReadLine(InLine, trimmed=.true.) + if (.not. OK) return + if (InLine=='') cycle + if (InLine(1:1)/='#') then + return + else + if (present(comment)) comment = trim(InLine(2:)) + end if + end do + + end function ReadLineSkipEmptyAndComments + + function ReadNextContentLine(this, filename, InLine) result(OK) + class(TTextFile) :: this + character(LEN=*), intent(in) :: filename + character(LEN=:), intent(out), allocatable :: InLine + logical :: OK + + if (.not. this%Opened()) call this%Open(filename) + OK = this%ReadLineSkipEmptyAndComments(InLine) + if (.not. OK) call this%Close() + + end function ReadNextContentLine + + function SkipLines(this, n) result(OK) + class(TTextFile) :: this + integer, intent(in) :: n + logical OK + integer ix + + do ix = 1, n + if (.not. this%ReadLine()) then + OK = .false. + return + end if + end do + OK = .true. + + end function SkipLines + + function Columns(this) result(n) + class(TTextFile) :: this + integer n + character(LEN=:), allocatable :: InLine + + if (this%ReadLineSkipEmptyAndComments(InLine)) then + n = File%TxtNumberColumns(InLine) + else + n=0 + end if + call this%Rewind() + + end function Columns + + function Lines(this, nocomments) result(n) + class(TTextFile) :: this + logical, intent(in), optional :: nocomments + integer n + character(LEN=:), allocatable :: InLine + + n=0 + if (DefaultTrue(nocomments)) then + do while (this%ReadLineSkipEmptyAndComments(InLine)) + n = n+1 + end do + else + do while (this%ReadLine()) + n = n+1 + end do + end if + call this%Rewind() + + end function Lines + + function DefaultAdvance(this,advance) + class(TTextFile) :: this + logical, intent(in), optional :: advance + character(3) :: DefaultAdvance + + if (PresentDefault(this%AdvanceDefault,advance)) then + DefaultAdvance='YES' + else + DefaultAdvance='NO' + end if + + end function + + subroutine NewLine(this) + class(TTextFile) :: this + call this%WriteItemTxt('') + end subroutine + + subroutine WriteLeftAligned(this, Form, str) + class(TTextFile) :: this + character(LEN=*) str, Form + Character(LEN=max(len(str),128)) tmp + + call this%CheckOpen(.true.) + tmp = str + write(this%unit,form, advance='NO') tmp + + end subroutine WriteLeftAligned + + + subroutine WriteInLineItems(this, S1, S2,S3,S4,S5,S6) + class(TTextFile) :: this + class(*), intent(in) :: S1 + class(*), intent(in), optional :: S2, S3,S4,S5,S6 + + call this%WriteInLine(S1) + if (present(S2)) call this%WriteInLine(S2) + if (present(S3)) call this%WriteInLine(S3) + if (present(S4)) call this%WriteInLine(S4) + if (present(S5)) call this%WriteInLine(S5) + if (present(S6)) call this%WriteInLine(S6) + + end subroutine WriteInLineItems + + subroutine WriteItemsTxt(this, S1, S2,S3,S4,S5,S6) + class(TTextFile) :: this + class(*), intent(in) :: S1 + class(*), intent(in), optional :: S2, S3,S4,S5,S6 + + call this%WriteInLineItems(S1, S2,S3,S4,S5,S6) + call this%NewLine() + + end subroutine WriteItemsTxt + + subroutine WriteInLineItem(this,str,form) + class(TTextFile) :: this + class(*), intent(in) :: str + character(LEN=*), intent(in), optional :: form + call this%WriteItemTxt(str,form,.false.) + end subroutine WriteInLineItem + + subroutine WriteInLineArray(this,str,form,n) + class(TTextFile) :: this + class(*), intent(in) :: str(:) + character(LEN=*), intent(in), optional :: form + integer, intent(in), optional :: n + + call this%WriteArrayTxt(str,form,.false.,number=n) + end subroutine WriteInLineArray + + subroutine TTextFile_WriteItem(this,R) + class(TTextFile) :: this + class(*), intent(in) :: R + call this%WriteItemTxt(R,advance=.true.) + end subroutine TTextFile_WriteItem + + subroutine TTextFile_WriteArray(this,R,n) + class(TTextFile) :: this + class(*), intent(in) :: R(:) + integer, intent(in), optional :: n + call this%WriteArrayTxt(R,number=n,advance=.true.) + end subroutine TTextFile_WriteArray + + + subroutine TTextFile_WriteOneAndArray(this, I, R) + class(TTextFile) :: this + class(*), intent(in) :: I + class(*), intent(in) :: R(:) + + call this%WriteInlineItem(I) + call this%WriteArrayTxt(R) + + end subroutine TTextFile_WriteOneAndArray + + subroutine WriteItemTxt(this, str, form, advance) + class(TTextFile) :: this + class(*), intent(in) :: str + character(LEN=*), intent(in), optional :: form + logical, intent(in), optional :: advance + character(LEN=3) :: Ad + + call this%CheckOpen(.true.) + Ad = this%DefaultAdvance(advance) + select type(str) + type is (character(LEN=*)) + if (Ad=='YES') then + write(this%unit,PresentDefault('(a)',form)) trim(str) + else + write(this%unit,PresentDefault('(a)',form), advance=Ad) str + end if + type is (real) + write(this%unit,PresentDefault(this%RealFormat,form), advance=Ad) str + type is (double precision) + write(this%unit,PresentDefault(this%RealFormat,form), advance=Ad) str + type is (integer) + write(this%unit,PresentDefault(this%IntegerFormat,form), advance=Ad) str + type is (logical) + write(this%unit,PresentDefault('(L2)',form), advance=Ad) str + class default + call this%Error('unknown type to write') + end select + + end subroutine WriteItemTxt + + subroutine WriteArrayTxt(this, str, form, advance, number) + class(TTextFile) :: this + class(*), intent(in) :: str(:) + character(LEN=*), intent(in), optional :: form + logical, intent(in), optional :: advance + integer, intent(in), optional :: number + integer n + character(LEN=3) :: Ad + + Ad = this%DefaultAdvance(advance) + n = PresentDefault(size(str),number) + call this%CheckOpen(.true.) + select type(str) + type is (character(LEN=*)) + write(this%unit,PresentDefault('(a)',form), advance=Ad) str(1:n) + type is (real) + write(this%unit,PresentDefault(this%RealFormat,form), advance=Ad) str(1:n) + type is (double precision) + write(this%unit,PresentDefault(this%RealFormat,form), advance=Ad) str(1:n) + type is (integer) + write(this%unit,PresentDefault(this%IntegerFormat,form), advance=Ad) str(1:n) + type is (logical) + write(this%unit,PresentDefault('(*(L2))',form), advance=Ad) str(1:n) + class default + call this%Error('unknown type to write') + end select + + end subroutine WriteArrayTxt + + subroutine WriteArray2Txt(this, R) + class(TTextFile) :: this + class(*), intent(in) :: R(:,:) + integer i + + do i=1, size(R,1) + call this%WriteArrayTxt(R(i,:)) + end do + + end subroutine WriteArray2Txt + + + subroutine WriteTrimTxt(this, S) + class(TTextFile) :: this + character(LEN=*), intent(in) :: S + + call this%WriteItemTxt(S, advance=.true.) + + end subroutine WriteTrimTxt + + subroutine WriteInLineTrim(this, string) + class(TTextFile) :: this + character(LEN=*), intent(in) :: string + + call this%WriteItemTxt(trim(string), advance=.false.) + end subroutine WriteInLineTrim + + + subroutine ReadArray2Txt(this, R, OK) + class(TTextFile) :: this + class(*) :: R(:,:) + logical, optional :: OK + integer i + + do i=1, size(R,1) + call this%ReadArrayTxt(R(i,:), OK=OK) + if (present(OK)) then + if (.not. OK) return + end if + end do + + end subroutine ReadArray2Txt + + + subroutine ReadItemTxt(this, R, OK) + class(TTextFile) :: this + class(*), intent(out) :: R + logical, optional :: OK + integer status + + call this%CheckOpen + select type(R) + type is (character(LEN=*)) + Read(this%unit,'(a)', iostat=status) R + type is (real) + Read(this%unit,*, iostat=status) R + type is (double precision) + Read(this%unit,*, iostat=status) R + type is (integer) + Read(this%unit,*, iostat=status) R + class default + call this%Error('unknown type to Read') + end select + if (status/=0 .and. (.not. IS_IOSTAT_END(status) .or. .not. present(OK))) & + & call this%Error('Error reading item') + if (present(OK)) OK = status==0 + + end subroutine ReadItemTxt + + subroutine ReadArrayTxt(this, R, n, OK) + class(TTextFile) :: this + class(*) :: R(1:) + integer, intent(in), optional :: n + logical, optional :: OK + integer status + + call this%CheckOpen + select type(R) + type is (real) + Read(this%unit,*, iostat=status) R(1:PresentDefault(size(R),n)) + type is (double precision) + Read(this%unit,*, iostat=status) R(1:PresentDefault(size(R),n)) + type is (integer) + Read(this%unit,*, iostat=status) R(1:PresentDefault(size(R),n)) + class default + call this%Error('unknown type to Read') + end select + + if (status/=0 .and. (.not. IS_IOSTAT_END(status) .or. .not. present(OK))) & + & call this%Error('Error reading item') + if (present(OK)) OK = status==0 + + end subroutine ReadArrayTxt + + subroutine ReadStringTxt(this, S, OK) + class(TTextFile) :: this + character(LEN=:), allocatable :: S + logical, optional :: OK + logical isOK + + isOK = this%ReadLine(S) + if (present(OK)) OK = isOK + + end subroutine ReadStringTxt + + subroutine WriteFormat(this, formatst, i1,i2,i3,i4,i5,i6,i7,i8) + class(TTextFile) :: this + character(LEN=*), intent(in) :: formatst + class(*), intent(in) :: i1 + class(*), intent(in),optional :: i2,i3,i4,i5,i6,i7,i8 + + call this%CheckOpen(.true.) + write(this%unit,'(a)') FormatString(formatst,i1,i2,i3,i4,i5,i6,i7,i8) + + end subroutine WriteFormat + + !Misc functions + + function TopCommentLine(aname) result(res) + !Get top comment line in file, including # + character(LEN=*), intent(IN) :: aname + character(LEN=:), allocatable :: res + Type(TTextFile) :: F + + call F%Open(aname) + res='' + do while (res == '') + if (.not. F%ReadLine(res)) exit + end do + If (res(1:1)/='#') then + res = '' + end if + call F%Close() + + end function TopCommentLine + + + function LastTopComment(aname) result(res) + !Get content of last commented line at the top of file (e.g. column header), without # + character(LEN=*), intent(IN) :: aname + character(LEN=:), allocatable :: res, InLine + Type(TTextFile) :: F + + call F%Open(aname) + res='' + do while (F%ReadLine(InLine)) + if (trim(InLine)=='') cycle + if (InLine(1:1)=='#') then + res = trim(adjustl(InLine(2:))) + else + exit + end if + end do + call F%Close() + + end function LastTopComment + + + function TxtFileColumns(aname) result(n) + character(LEN=*), intent(IN) :: aname + integer n + Type(TTextFile) :: F + + call F%Open(aname) + n = F%Columns() + call F%Close() + + end function TxtFileColumns + + function TxtFileLines(aname) result(n) + character(LEN=*), intent(IN) :: aname + integer n + Type(TTextFile) :: F + + call F%Open(aname) + n = F%Lines() + call F%Close() + + end function TxtFileLines + + + function LastFileLine(aname) + character(LEN=*), intent(IN) :: aname + character(LEN=:), allocatable :: LastFileLine + Type(TTextFile) :: F + + LastFileLine = '' + call F%Open(aname) + do while (F%ReadLine(LastFileLine)) + end do + call F%Close() + end function LastFileLine + + + function CharIsSlash(C) + character, intent(in) :: C + logical CharIsSlash + character, parameter :: win_slash = char(92) + + CharIsSlash = C == win_slash .or. C == '/' + + end function CharIsSlash + + function IsFullPath(aname) + character(LEN=*), intent(IN) :: aname + logical IsFullPath + + IsFullPath = aname/='' + if (.not. IsFullPath) return + IsFullPath = CharIsSlash(aname(1:1)) + + end function IsFullpath + + function ExtractFilePath(aname) + character(LEN=*), intent(IN) :: aname + character(LEN=:), allocatable :: ExtractFilePath + integer i + + do i = len_trim(aname), 1, -1 + if (CharIsSlash(aname(i:i))) then + ExtractFilePath = aname(1:i) + return + end if + end do + ExtractFilePath = '' + + end function ExtractFilePath + + function ExtractFileExt(aname) + character(LEN=*), intent(IN) :: aname + character(LEN=:), allocatable :: ExtractFileExt + integer len, i + + len = len_trim(aname) + do i = len, 1, -1 + if (CharIsSlash(aname(i:i))) then + ExtractFileExt = '' + return + else if (aname(i:i)=='.') then + ExtractFileExt= aname(i:len) + return + end if + end do + ExtractFileExt = '' + + end function ExtractFileExt + + + function ExtractFileName(aname, no_ext, all_ext) + character(LEN=*), intent(IN) :: aname + character(LEN=:), allocatable :: ExtractFileName + logical, intent(in), optional :: no_ext, all_ext + integer alen, i + + alen = len_trim(aname) + do i = alen, 1, -1 + if (CharIsSlash(aname(i:i))) then + ExtractFileName = aname(i+1:alen) + exit + end if + end do + if (.not. allocated(ExtractFileName)) ExtractFileName = trim(aname) + if (DefaultFalse(no_ext)) then + do i = len(ExtractFileName), 1, -1 + if (ExtractFileName(i:i)=='.') then + ExtractFileName = ExtractFileName(1:i-1) + if (.not. DefaultFalse(all_ext)) exit + end if + end do + end if + + end function ExtractFileName + + function ChangeFileExt(aname,ext) + character(LEN=*), intent(IN) :: aname,ext + character(LEN=:), allocatable :: ChangeFileExt + integer len, i + + len = len_trim(aname) + do i = len, 1, -1 + if (aname(i:i)=='.') then + ChangeFileExt = aname(1:i) // trim(ext) + return + end if + end do + ChangeFileExt = trim(aname) // '.' // trim(ext) + + end function ChangeFileExt + + function CheckTrailingSlash(aname) + character(LEN=*), intent(in) :: aname + character(LEN=:), allocatable :: CheckTrailingSlash + integer alen + + alen = len_trim(aname) + if (CharIsSlash(aname(alen:alen))) then + CheckTrailingSlash = aname(1:alen) + else + CheckTrailingSlash = trim(aname)//'/' + end if + + end function CheckTrailingSlash + + + function File_Join(path, aname) + character(LEN=*), intent(in) :: path, aname + character(LEN=:), allocatable :: File_Join + + File_Join = CheckTrailingSlash(path)//trim(aname) + + end function File_Join + + + subroutine DeleteFile(aname) + character(LEN=*), intent(IN) :: aname + integer file_id, status + + if (FileExists(aname)) then + open(newunit = file_id, file = aname, iostat=status) + if (status/=0) return + close(unit = file_id, status = 'DELETE') + end if + + end subroutine DeleteFile + + + subroutine ReadTextVector(aname, vec, n) + character(LEN=*), intent(IN) :: aname + integer, intent(in) :: n + class(*), intent(out) :: vec(n) + integer j + Type(TTextFile) :: F + + call F%Open(aname) + do j=1,n + call F%Read(vec(j)) + end do + call F%Close() + + end subroutine ReadTextVector + + subroutine WriteTextVector(aname, vec, n, fmt) + character(LEN=*), intent(IN) :: aname + integer, intent(in), optional :: n + class(*), intent(in) :: vec(:) + character(LEN=*), intent(in), optional :: fmt + integer j + Type(TTextFile) :: F + + call F%CreateFile(aname) + if (present(fmt)) then + if (isFLoat(vec)) then + F%RealFormat = fmt + else + F%IntegerFormat = fmt + end if + end if + do j=1, PresentDefault(size(vec),n) + call F%Write(vec(j)) + end do + call F%Close() + + end subroutine WriteTextVector + + subroutine WriteTextMatrix(aname, mat, m, n, fmt) + character(LEN=*), intent(IN) :: aname + integer, intent(in), optional :: m, n + character(LEN=*), intent(in), optional :: fmt + class(*), intent(in) :: mat(:,:) + Type(TTextFile) :: F + + call F%CreateFile(aname) + if (present(fmt)) then + if (isFLoat(mat)) then + F%RealFormat = fmt + else + F%IntegerFormat = fmt + end if + end if + if (present(m) .or. present(n)) then + call F%Write(mat(1:PresentDefault(size(mat,1),m),1:PresentDefault(size(mat,2),n))) + else + call F%Write(mat) + end if + call F%Close() + + end subroutine WriteTextMatrix + + + subroutine ReadTextMatrix(aname, mat, inm,inn) + character(LEN=*), intent(IN) :: aname + integer, intent(in), optional :: inm,inn + real(kind(1.d0)), intent(out) :: mat(:,:) + character(LEN=:), allocatable :: InLine + integer j,k, status, n,m + Type(TTextFile) :: F + + m = PresentDefault(size(mat,dim=1),inm) + n = PresentDefault(size(mat,dim=2),inn) + call F%Open(aname) + + do j=1,m + status = 1 + if (.not. F%ReadLineSkipEmptyAndComments(InLine)) exit + read (InLine,*, iostat=status) mat(j,1:n) + if (status/=0) exit + end do + if (status/=0) then + call F%Rewind() !Try other possible format + do j=1,m + do k=1,n + status = 1 + if (.not. F%ReadLineSkipEmptyAndComments(InLine)) exit + read (InLine,*, iostat=status) mat(j,k) + if (status/=0) call F%Error( 'matrix file is the wrong size') + end do + end do + end if + + if (F%ReadLineSkipEmptyAndComments(InLine)) call F%Error( 'matrix file is the wrong size (too big)') + + call F%Close() + + end subroutine ReadTextMatrix + + subroutine LoadTxt_2D(aname, mat, m, n, comment) + character(LEN=*), intent(IN) :: aname + real(kind(1.d0)), allocatable :: mat(:,:) + integer mm, nn, j + Type(TTextFile) :: F + character(LEN=:), allocatable :: InLine + integer, optional, intent(out) :: m, n + character(LEN=:), allocatable, optional, intent(out) :: comment + integer status + + call F%Open(aname) + nn = F%Columns() + mm = F%Lines() + allocate(mat(mm,nn)) + j=1 + do while (F%ReadLineSkipEmptyAndComments(InLine, comment = comment)) + read (InLine,*, iostat=status) mat(j,1:nn) + if (status/=0) call F%Error( 'LoadTxt: error reading line:' //trim(InLine)) + j = j+1 + end do + call F%Close() + if (present(m)) m = mm + if (present(n)) n = nn + + end subroutine LoadTxt_2D + + subroutine LoadTxt_1D(aname, vec, n, comment) + character(LEN=*), intent(IN) :: aname + real(kind(1.d0)), allocatable :: vec(:) + integer nn, j + Type(TTextFile) :: F + character(LEN=:), allocatable :: InLine + integer, optional, intent(out) :: n + character(LEN=:), allocatable, optional, intent(out) :: comment + integer status + + call F%Open(aname) + nn = F%Lines() + allocate(vec(nn)) + j=1 + do while (F%ReadLineSkipEmptyAndComments(InLine, comment = comment)) + read (InLine,*, iostat=status) vec(j) + if (status/=0) call F%Error( 'LoadTxt: error reading line:' //trim(InLine)) + j = j+1 + end do + call F%Close() + if (present(n)) n = nn + + end subroutine LoadTxt_1D + + function CreateTextFile(fname) + character(LEN=*), intent(in) :: fname + Type(TTextFile) :: CreateTextFile + + call CreateTextFile%CreateFile(fname) + + end function CreateTextFile + + function OpenTextFile(fname) + character(LEN=*), intent(in) :: fname + Type(TTextFile) :: OpenTextFile + + call OpenTextFile%OpenFile(fname) + + end function OpenTextFile + + + end module FileUtils diff --git a/IniObjects.f90 b/IniObjects.f90 new file mode 100644 index 0000000..2b174bc --- /dev/null +++ b/IniObjects.f90 @@ -0,0 +1,927 @@ +!Module to read in name/value pairs from a file, with each line of the form line 'name = value' +!Should correctly interpret FITS headers +!Antony Lewis (http://cosmologist.info/). Released to the public domain. +!2014 using Fortran 2003/2008 + +module IniObjects + use FileUtils + use StringUtils + use MiscUtils + implicit none + private + + character(LEN=0), target :: Empty_String = '' + integer, parameter :: Ini_Enumeration_Len = 64 + + type TNameValue + character(LEN=:), allocatable :: Name + character(LEN=:), allocatable :: Value + end type TNameValue + + type TNameValue_pointer + Type(TNameValue), pointer :: P + end type TNameValue_pointer + + Type TNameValueList + integer :: Count = 0 + integer :: Delta = 128 + integer :: Capacity = 0 + logical :: ignoreDuplicates = .false. + logical :: AllowDuplicateKeys = .false. + !Use pointer so arrays can be moved around quickly without deep copies + type(TNameValue_pointer), dimension(:), allocatable :: Items + contains + procedure :: Init => TNameValueList_Init + procedure :: Clear => TNameValueList_Clear + procedure :: ValueOf => TNameValueList_ValueOf + procedure :: IndexOf => TNameValueList_IndexOf + procedure :: HasKey => TNameValueList_HasKey + procedure :: SetCapacity => TNameValueList_SetCapacity + procedure :: Delete => TNameValueList_Delete + procedure :: Error => TNameValueList_Error + procedure :: FailStop => TNameValueList_FailStop + procedure :: AddString => TNameValueList_Add + procedure :: Name => TNameValueList_Name + procedure :: Value => TNameValueList_Value + procedure :: Override => TNameValueList_Override + procedure :: TNameValueList_AddDouble + procedure :: TNameValueList_AddReal + procedure :: TNameValueList_AddInt + procedure :: TNameValueList_AddLogical + generic :: Add => AddString, TNameValueList_AddDouble, & + & TNameValueList_AddReal, TNameValueList_AddInt,& + & TNameValueList_AddLogical + final :: TNameValueList_Free + end Type TNameValueList + + Type, extends(TNameValueList) :: TIniFile + logical :: SlashComments = .false. + logical :: Echo_Read = .false. + logical :: Fail_on_not_found = .false. + logical :: ExpandEnvironmentVariables = .true. + character(LEN=:), allocatable :: Original_filename + Type(TNameValueList) :: ReadValues + contains + procedure :: Open => Ini_Open + procedure :: Open_Fromlines => Ini_Open_Fromlines + procedure :: Close => Ini_Close + procedure :: Read_String => Ini_Read_String + procedure :: Read_String_Default => Ini_Read_String_Default + procedure :: Read_String_Array => Ini_Read_String_Array + procedure :: Read_Int_Array => Ini_Read_Int_Array + procedure :: Read_Int => Ini_Read_Int + procedure :: Read_Double => Ini_Read_Double + procedure :: Read_Double_Array => Ini_Read_Double_Array + procedure :: Read_Real => Ini_Read_Real + procedure :: Read_Real_Array => Ini_Read_Real_Array + procedure :: Read_Logical => Ini_Read_Logical + procedure :: Read_Enumeration => Ini_Read_Enumeration + procedure :: Read_Enumeration_List => Ini_Read_Enumeration_List + procedure :: TestEqual => Ini_TestEqual + procedure :: SaveReadValues => Ini_SaveReadValues + procedure :: Key_To_Arraykey => Ini_Key_To_Arraykey + procedure :: EnumerationValue => Ini_EnumerationValue + procedure :: ResolveLinkedFile => Ini_ResolveLinkedFile + procedure :: ExpandEnvironment => Ini_ExpandEnvironment + procedure, private :: NameValue_AddLine => Ini_NameValue_AddLine + procedure, private :: EmptyCheckDefault => Ini_EmptyCheckDefault + procedure, private :: Ini_Read_Real_Change + procedure, private :: Ini_Read_Double_Change + procedure, private :: Ini_Read_Int_Change + procedure, private :: Ini_Read_Logical_Change + procedure, private :: Ini_Read_String_Change + generic :: Read => Ini_Read_Real_Change, Ini_Read_Double_Change, Ini_Read_Int_Change, & + & Ini_Read_Logical_Change, Ini_Read_String_Change + end Type TIniFile + + + public TNameValueList, TIniFile, Ini_Enumeration_Len +contains + + subroutine TNameValueList_Init(this, ignoreDuplicates) + class(TNameValueList) :: this + logical, intent(in), optional :: ignoreDuplicates + + this%Count = 0 + this%Capacity = 0 + this%Delta = 128 + this%ignoreDuplicates = DefaultFalse(ignoreDuplicates) + + end subroutine TNameValueList_Init + + subroutine TNameValueList_Clear(this) + class(TNameValueList) :: this + integer i, status + + do i=this%count,1,-1 + deallocate (this%Items(i)%P, stat = status) + end do + deallocate (this%Items, stat = status) + call this%Init() + + end subroutine TNameValueList_Clear + + subroutine TNameValueList_Free(this) + Type(TNameValueList) :: this + call this%Clear() + end subroutine TNameValueList_Free + + function TNameValueList_ValueOf(this, AName) result(AValue) + class(TNameValueList), intent(in) :: this + character(LEN=*), intent(in) :: AName + character(LEN=:), pointer :: AValue + integer i + + i = this%IndexOf(AName) + if (i/=-1) then + AValue => this%Items(i)%P%Value + else + AValue => Empty_String + end if + + end function TNameValueList_ValueOf + + function TNameValueList_Name(this, i) result(Name) + class(TNameValueList), intent(in) :: this + integer i + character(LEN=:), pointer :: Name + + Name => this%Items(i)%P%Name + + end function TNameValueList_Name + + function TNameValueList_Value(this, i) result(Value) + class(TNameValueList), intent(in) :: this + integer i + character(LEN=:), pointer :: Value + + Value => this%Items(i)%P%Value + + end function TNameValueList_Value + + + function TNameValueList_IndexOf(this, AName) result (AValue) + class(TNameValueList), intent(in) :: this + character(LEN=*), intent(in) :: AName + integer :: AValue + integer i + + do i=1, this%Count + if (this%Items(i)%P%Name == AName) then + AValue = i + return + end if + end do + AValue = -1 + + end function TNameValueList_IndexOf + + function TNameValueList_HasKey(this, AName) result (AValue) + class(TNameValueList), intent(in) :: this + character(LEN=*), intent(in) :: AName + logical :: AValue + + AValue = this%IndexOf(AName) /= -1 + + end function TNameValueList_HasKey + + subroutine TNameValueList_Override(this, Settings, only_if_exists) + class(TNameValueList) :: this + class(TNameValueList), intent(in) :: Settings + logical, intent(in), optional :: only_if_exists + integer i, ix + character(LEN=:), pointer :: name + + do i=1, Settings%Count + name => Settings%Name(i) + ix = this%IndexOf(name) + if (ix/=-1) then + this%Items(ix)%P%Value = Settings%Value(i) + elseif (.not. DefaultFalse(only_if_exists)) then + call this%Add(name, Settings%Value(i)) + end if + end do + + end subroutine TNameValueList_Override + + subroutine TNameValueList_Add(this, AName, AValue, only_if_undefined) + class(TNameValueList) :: this + character(LEN=*), intent(in) :: AName, AValue + logical, optional, intent(in) :: only_if_undefined + logical isDefault + + isDefault = DefaultTrue(only_if_undefined) + + if ((.not. this%AllowDuplicateKeys .or. isDefault) .and. this%HasKey(AName)) then + if (this%ignoreDuplicates .or. isDefault) return + call this%Error('duplicate key name',AName) + end if + if (this%Count == this%Capacity) call this%SetCapacity(this%Capacity + this%Delta) + this%Count = this%Count + 1 + allocate(this%Items(this%Count)%P) + this%Items(this%Count)%P%Name = trim(adjustl(AName)) + this%Items(this%Count)%P%Value = trim(adjustl(AValue)) + + end subroutine TNameValueList_Add + + subroutine TNameValueList_AddReal(this, AName, AValue) + class(TNameValueList) :: this + character(LEN=*), intent(in) :: AName + real, intent(in) :: AValue + character(LEN=32) tmp + + write(tmp,*) AValue + call this%AddString(AName, Tmp) + + end subroutine TNameValueList_AddReal + + subroutine TNameValueList_AddDouble(this, AName, AValue) + class(TNameValueList) :: this + character(LEN=*), intent(in) :: AName + double precision, intent(in) :: AValue + character(LEN=32) tmp + + write(tmp,*) AValue + call this%AddString(AName, Tmp) + + end subroutine TNameValueList_AddDouble + + + subroutine TNameValueList_AddInt(this, AName, AValue) + class(TNameValueList) :: this + character(LEN=*), intent(in) :: AName + integer, intent(in) :: AValue + character(LEN=32) tmp + + write(tmp,*) AValue + call this%AddString(AName, Tmp) + + end subroutine TNameValueList_AddInt + + + subroutine TNameValueList_AddLogical(this, AName, AValue) + class(TNameValueList) :: this + character(LEN=*), intent(in) :: AName + logical, intent(in) :: AValue + character(LEN=32) tmp + + write(tmp,*) AValue + call this%AddString(AName, Tmp) + + end subroutine TNameValueList_AddLogical + + + subroutine TNameValueList_SetCapacity(this, C) + class(TNameValueList) :: this + integer C + type(TNameValue_pointer), dimension(:), allocatable :: TmpItems + + if (this%Count > 0) then + if (C < this%Count) call this%Error('TNameValueList_SetCapacity, smaller than Count') + allocate(TmpItems(C)) + TmpItems(:this%Count) = this%Items(:this%Count) + call move_alloc(TmpItems, this%Items) + else + allocate(this%Items(C)) + end if + this%Capacity = C + + end subroutine TNameValueList_SetCapacity + + subroutine TNameValueList_Delete(this, i) + class(TNameValueList) :: this + integer, intent(in) :: i + + deallocate(this%Items(i)%P) + if (this%Count > 1) this%Items(i:this%Count-1) = this%Items(i+1:this%Count) + this%Count = this%Count -1 + + end subroutine TNameValueList_Delete + + subroutine TNameValueList_FailStop(this) + class(TNameValueList) :: this + + stop + end subroutine TNameValueList_FailStop + + subroutine TNameValueList_Error(this, Msg, Key) + class(TNameValueList) :: this + character(LEN=*) :: Msg + character(LEN=*), optional :: Key + + if (present(Key)) then + write(*,*) 'Error for key "'//trim(Key)//'" : '//Msg + else + write(*,*) 'Error :'//Msg + end if + call this%FailStop() + + end subroutine TNameValueList_Error + + subroutine Ini_EmptyCheckDefault(this, Key, Default) + class(TIniFile) :: this + character(LEN=*), intent(in) :: Key + class(*), optional, intent(in) :: Default + + if (.not. present(Default)) call this%Error('missing key',Key) + + end subroutine Ini_EmptyCheckDefault + + function Ini_ExpandEnvironment(this, InValue) result(res) + !expand $(PLACEHOLDER) with environment variables, as Makefile + class(TIniFile) :: this + character(LEN=:), allocatable, intent(in) :: InValue + character(LEN=:), allocatable :: res, S + integer i + + i = index(InValue,'$(') + if (i > 0) then + res = InValue(:i-1) + S = InValue(i:) + i=1 + do while (i <= len(S)) + if (S(i:i)=='$') then + if (S(i+1:i+1)=='$') then + res = res // '$' + i = i + 1 + else if (S(i+1:i+1)=='(') then + S = S(i+2:) + i = index(S,')') + if (i==0) then + call this%Error('bad environment placeholder: '//InValue) + end if + res = res //GetEnvironmentVariable(S(:i-1)) + end if + else + res = res // S(i:i) + end if + i = i + 1 + end do + else + res = InValue + end if + + end function Ini_ExpandEnvironment + + subroutine Ini_NameValue_AddLine(this,AInLine,only_if_undefined) + class(TIniFile) :: this + character (LEN=*), intent(IN) :: AInLine + integer EqPos, slashpos, lastpos + logical, optional, intent(in) :: only_if_undefined + logical isDefault + character (LEN=len(AInLine)) :: AName, InLine + character(LEN=:), allocatable :: val + + isDefault = DefaultFalse(only_if_undefined) + + InLine=trim(adjustl(AInLine)) + EqPos = scan(InLine,'=') + if (EqPos/=0 .and. InLine(1:1)/='#' .and. .not. StringStarts(InLine,'COMMENT')) then + AName = trim(InLine(1:EqPos-1)) + + val = adjustl(InLine(EqPos+1:)) + if (this%SlashComments) then + slashpos=scan(val,'/') + if (slashpos /= 0) then + val = val(1:slashpos-1) + end if + end if + if (this%ExpandEnvironmentVariables) then + val = this%ExpandEnvironment(val) + end if + lastpos=len_trim(val) + if (lastpos>1) then + if (val(1:1)=='''' .and. val(lastpos:lastpos)=='''') then + val = val(2:lastpos-1) + end if + end if + call this%Add(AName, val, only_if_undefined = isDefault) + end if + + end subroutine Ini_NameValue_AddLine + + function Ini_ResolveLinkedFile(this, name, thisfilename) result(IncludeFile) + class(TIniFile) :: this + character(LEN=*), intent(in) :: name, thisfilename + character(LEN=:), allocatable :: IncludeFile + + if (.not. File%IsFullPath(name)) then + IncludeFile= File%ExtractPath(thisfilename)//trim(name) + if (File%Exists(IncludeFile)) then + if (File%Exists(name) .and. name/=IncludeFile) & + call this%Error(trim(thisfilename)// & + ' , ambiguous multiple matches to include file: '//trim(name)) + else + IncludeFile= name + end if + else + IncludeFile= name + end if + if (.not. File%Exists(IncludeFile)) then + call this%Error(trim(thisfilename)//' , include file not found: '//trim(name)) + end if + + end function Ini_ResolveLinkedFile + + recursive subroutine Ini_Open(this, filename, error, slash_comments, append,only_if_undefined) + class(TIniFile) :: this + character (LEN=*), intent(IN) :: filename + logical, intent(OUT), optional :: error + logical, optional, intent(IN) :: slash_comments + logical, optional, intent(in) :: append, only_if_undefined + character (LEN=:), allocatable :: IncludeFile + character(LEN=:), allocatable :: InLine + integer lastpos, i, status + Type(TNameValueList) IncudeFiles, DefaultValueFiles + logical isDefault + Type(TTextFile) :: F + + + isDefault = DefaultFalse(only_if_undefined) + + if (.not. DefaultFalse(append)) then + call this%TNameValueList%Init() + call this%ReadValues%Init(.true.) + this%Original_filename = filename + end if + + this%SlashComments = DefaultFalse(slash_comments) + + call F%Open(filename, status=status) + if (status/=0) then + if (present(error)) then + error=.true. + return + else + call this%Error('Ini_Open, file not found: '//trim(filename)) + end if + end if + + call IncudeFiles%Init() + call DefaultValueFiles%Init() + + do + if (.not. F%ReadLineSkipEmptyAndComments(InLine)) exit + if (InLine == 'END') exit + if (StringStarts(InLine,'INCLUDE(')) then + lastpos = scan(InLine,')') + if (lastpos/=0) then + call IncudeFiles%Add(InLine(9:lastpos-1),'') + else + call this%Error('Ini_Open, error in INCLUDE line: '//trim(filename)) + end if + elseif (StringStarts(InLine,'DEFAULT(')) then + !Settings to read in as defaults, overridden by subsequent re-definitions + lastpos = scan(InLine,')') + if (lastpos/=0) then + call DefaultValueFiles%Add(InLine(9:lastpos-1),'') + else + call this%Error('Ini_Open, error in DEFAULT line: '//trim(filename)) + end if + elseif (InLine /= '') then + call this%NameValue_AddLine(InLine, only_if_undefined=isDefault) + end if + end do + + call F%Close() + if (present(error)) error=.false. + + do i=1, IncudeFiles%Count + if (DefaultFalse(error)) exit + IncludeFile = this%ResolveLinkedFile(IncudeFiles%Items(i)%P%Name, filename) + call this%Open(IncludeFile, error, slash_comments, append=.true.,only_if_undefined=isDefault) + end do + do i=1, DefaultValueFiles%Count + if (DefaultFalse(error)) exit + IncludeFile = this%ResolveLinkedFile(DefaultValueFiles%Items(i)%P%Name, filename) + call this%Open(IncludeFile, error, slash_comments, append=.true., only_if_undefined=.true.) + end do + call IncudeFiles%Clear() + call DefaultValueFiles%Clear() + + end subroutine Ini_Open + + subroutine Ini_Open_Fromlines(this, Lines, NumLines, slash_comments) + class(TIniFile) :: this + integer, intent(IN) :: NumLines + character (LEN=*), dimension(NumLines), intent(IN) :: Lines + logical, intent(IN), optional :: slash_comments + integer i + + call this%TNameValueList%Init() + call this%ReadValues%Init(.true.) + + if (present(slash_comments)) then + this%SlashComments = slash_comments + end if + + do i=1,NumLines + call this%NameValue_AddLine(Lines(i)) + end do + + end subroutine Ini_Open_Fromlines + + + subroutine Ini_Close(this) + class(TIniFile) :: this + + call this%TNameValueList%Clear() + call this%ReadValues%Clear() + + end subroutine Ini_Close + + function Ini_Read_String(this, Key, NotFoundFail) result(AValue) + class(TIniFile) :: this + character (LEN=*), intent(IN) :: Key + logical, optional, intent(IN) :: NotFoundFail + character(LEN=:), pointer :: AValue + + AValue => this%ValueOf(Key) + + if (AValue/='') then + call this%ReadValues%Add(Key, AValue) + if (this%Echo_Read) write (*,*) trim(Key)//' = ',trim(AValue) + return + end if + if (PresentDefault(this%fail_on_not_found, NotFoundFail)) then + call this%Error('key not found',Key) + end if + + end function Ini_Read_String + + function Ini_Read_String_Default(this, Key, Default, AllowBlank, EnvDefault) result(AValue) + !Returns from this first; if not found tries environment variable if EnvDefault is true, + !If not found returns Default (unless not present, in which case gives error) + !AllBlank specifies whether an empty string counts as a valid result to give instead of Default + class(TIniFile) :: this + character (LEN=*), intent(IN) :: Key + character (LEN=*), intent(IN), optional ::Default + logical, intent(in), optional :: AllowBlank + logical, optional, intent(IN) :: EnvDefault + character(LEN=:), allocatable :: AValue + logical is_present + + if (this%HasKey(Key)) then + AValue = this%Read_String(Key, .false.) + if (AValue/='' .or. DefaultFalse(AllowBlank)) return + else + AValue='' + end if + if (DefaultFalse(EnvDefault)) then + AValue = GetEnvironmentVariable(Key,is_present) + if (DefaultFalse(AllowBlank) .and. is_present) return + end if + if (AValue=='') then + if (.not. present(Default)) call this%Error('key not found',Key) + AValue = Default + end if + call this%ReadValues%Add(Key, AValue) + if (this%Echo_Read) write (*,*) trim(Key)//' = ',trim(AValue) + + end function Ini_Read_String_Default + + function Ini_TestEqual(this, Key, value, EmptyOK) result(OK) + class(TIniFile) :: this + character (LEN=*), intent(IN) :: Key, value + logical, intent(in), optional :: EmptyOK + character(LEN=:), pointer :: val + logical :: OK + + val => this%ValueOf(Key) + OK = val == value + if (.not. OK .and. present(EmptyOK)) then + OK = EmptyOK .and. val=='' + end if + + end function Ini_TestEqual + + function Ini_Key_To_Arraykey(this,Key, index) result(AValue) + class(TIniFile), intent(in) :: this + character (LEN=*), intent(IN) :: Key + integer, intent(in) :: index + character(LEN=:), allocatable :: AValue + character(LEN=32) :: numstr + + write (numstr,*) index + numstr=adjustl(numstr) + AValue = trim(Key) // '(' // trim(numStr) // ')' + + end function Ini_Key_To_Arraykey + + function Ini_Read_String_Array(this, Key, index, NotFoundFail) result(AValue) + class(TIniFile) :: this + integer, intent(in) :: index + character (LEN=*), intent(IN) :: Key + logical, optional, intent(IN) :: NotFoundFail + character(LEN=:), pointer :: AValue + character(LEN=:), allocatable :: ArrayKey + + ArrayKey = this%Key_To_Arraykey(Key,index) + AValue => this%Read_String(ArrayKey, NotFoundFail) + + end function Ini_Read_String_Array + + function Ini_Read_Int_Array(this,Key, index, Default) + !Reads Key(1), Key(2), etc. + class(TIniFile) :: this + integer Ini_Read_Int_Array + integer, optional, intent(IN) :: Default + integer, intent(in) :: index + character(LEN=*), intent(IN) :: Key + character(LEN=:), allocatable :: ArrayKey + + ArrayKey = this%Key_To_Arraykey(Key,index) + Ini_Read_Int_Array = this%Read_Int(ArrayKey, Default) + + end function Ini_Read_Int_Array + + function Ini_Read_Int(this, Key, Default, min, max, OK) + class(TIniFile) :: this + integer Ini_Read_Int + integer, optional, intent(IN) :: Default, min, max + logical, intent(out), optional :: OK + character(LEN=*), intent(IN) :: Key + character(LEN=:), pointer :: S + integer status + + S => this%Read_String(Key,.not. present(Default)) + if (S == '') then + call this%EmptyCheckDefault(Key,Default) + Ini_Read_Int = Default + call this%ReadValues%Add(Key, Default) + else + if (verify(trim(S),'-+0123456789') /= 0) then + status=1 + if (present(OK)) then + Ini_Read_Int=-1 + OK = .false. + return + end if + else + read (S,*, iostat=status) Ini_Read_Int + end if + if (status/=0) call this%Error('error reading integer',Key) + if (present(max)) then + if (Ini_Read_Int > max) call this%Error('value > max',Key) + end if + if (present(min)) then + if (Ini_Read_Int < min) call this%Error('value < min',Key) + end if + end if + if (present(OK)) OK = .true. + + end function Ini_Read_Int + + integer function Ini_EnumerationValue(this, S, Names) + class(TIniFile) :: this + character(LEN=*) :: S + character(LEN=Ini_Enumeration_Len), intent(in) :: Names(:) + integer i + + do i=1, size(Names) + if (S==Names(i)) then + Ini_EnumerationValue = i + return + end if + end do + Ini_EnumerationValue = -1 + + end function Ini_EnumerationValue + + function Ini_Read_Enumeration(this, Key, Names, Default) + class(TIniFile) :: this + character(LEN=*), intent(in) :: Key + character(LEN=Ini_Enumeration_Len), intent(in) :: Names(:) + integer, optional, intent(in) :: Default + integer Ini_Read_Enumeration + character(LEN=:), pointer :: S + logical OK + + Ini_Read_Enumeration = this%Read_Int(Key, Default, OK = OK) + if (OK) then + if (Ini_Read_Enumeration<1 .or. Ini_Read_Enumeration> size(Names)) & + & call this%Error('enumeration value not valid',Key) + else + S => this%ValueOf(Key) + Ini_Read_Enumeration = this%EnumerationValue(S, Names) + if (Ini_Read_Enumeration<0) call this%Error('"'//S//'" enumeration name not recognised',Key) + end if + + end function Ini_Read_Enumeration + + subroutine Ini_Read_Enumeration_List(this, Key, Names, Enums, nvalues, max_enums, Default) + class(TIniFile) :: this + character(LEN=*), intent(in) :: Key + character(LEN=Ini_Enumeration_Len), intent(in) :: Names(:) + integer, allocatable, intent(out) :: Enums(:) + character(LEN=*), optional, intent(in) :: Default + integer, intent(in), optional :: nvalues, max_enums + integer, parameter :: defmax = 128 + integer :: maxvalues + integer, allocatable :: Values(:) + character(LEN=Ini_Enumeration_Len+8) part + character(LEN=:), allocatable :: InLine + integer i, slen, pos, n, status + + InLine = this%Read_String_Default(Key, Default) + pos = 1 + slen = len_trim(InLine) + n=0 + maxvalues = PresentDefault(PresentDefault(defmax, max_enums), nvalues) + allocate(Values(maxvalues)) + do + do while (pos <= slen) + if (IsWhiteSpace(InLine(pos:pos))) then + pos = pos+1 + else + exit + endif + end do + read(InLine(pos:), *, iostat=status) part + if (status/=0) exit + pos = pos + len_trim(part) + i= this%EnumerationValue(part, Names) + if (i<0) call this%Error('"'//part//'" enumeration name not recognised',Key) + n=n+1 + if (n > defmax) call this%Error('More than maximum enumeration values', Key) + values(n) = i + if (n == maxvalues) exit + end do + if (present(nvalues)) then + if (n==1) then + !Fill whole array with same value + allocate(Enums(nvalues), source= values(1)) + return + elseif (n/=nvalues) then + call this%Error('Wrong number of enumeration values', Key) + end if + end if + allocate(Enums, source= values(:n)) + + end subroutine Ini_Read_Enumeration_List + + function Ini_Read_Double(this,Key, Default, min, max) + class(TIniFile) :: this + double precision Ini_Read_Double + double precision, optional, intent(IN) :: Default, min,max + character (LEN=*), intent(IN) :: Key + character(LEN=:), pointer :: S + integer status + + S => this%Read_String(Key,.not. present(Default)) + if (S == '') then + call this%EmptyCheckDefault(Key,Default) + Ini_Read_Double = Default + call this%ReadValues%Add(Key, Default) + else + read (S,*, iostat=status) Ini_Read_Double + if (status/=0) call this%Error('error reading double',Key) + end if + if (present(max)) then + if (Ini_Read_Double > max) call this%Error('value > max',Key) + end if + if (present(min)) then + if (Ini_Read_Double < min) call this%Error('value < min',Key) + end if + + end function Ini_Read_Double + + function Ini_Read_Double_Array(this,Key, index, Default, min, max) + !Reads Key(1), Key(2), etc. + class(TIniFile) :: this + double precision Ini_Read_Double_Array + double precision, optional, intent(IN) :: Default, min, max + integer, intent(in) :: index + character(LEN=*), intent(IN) :: Key + character(LEN=:), allocatable :: ArrayKey + + ArrayKey = this%Key_To_Arraykey(Key,index) + Ini_Read_Double_Array = this%Read_Double(ArrayKey, Default, min, max) + + end function Ini_Read_Double_Array + + + function Ini_Read_Real(this,Key, Default, min, max) + class(TIniFile) :: this + real Ini_Read_Real + real, optional, intent(IN) :: Default, min, max + character(LEN=*), intent(IN) :: Key + character(LEN=:), pointer :: S + integer status + + S => this%Read_String(Key,.not. present(Default)) + if (S == '') then + call this%EmptyCheckDefault(Key,Default) + Ini_Read_Real = Default + call this%ReadValues%Add(Key, Default) + else + read (S,*, iostat=status) Ini_Read_Real + if (status/=0) call this%Error('error reading real',Key) + end if + if (present(max)) then + if (Ini_Read_Real > max) call this%Error('value > max',Key) + end if + if (present(min)) then + if (Ini_Read_Real < min) call this%Error('value < min',Key) + end if + + end function Ini_Read_Real + + + function Ini_Read_Real_Array(this,Key, index, Default, min, max) + !Reads Key(1), Key(2), etc. + class(TIniFile) :: this + real Ini_Read_Real_Array + real, optional, intent(IN) :: Default, min, max + integer, intent(in) :: index + character(LEN=*), intent(IN) :: Key + character(LEN=:), allocatable :: ArrayKey + + ArrayKey = this%Key_To_Arraykey(Key,index) + Ini_Read_Real_Array = this%Read_Real(ArrayKey, Default,min,max) + + end function Ini_Read_Real_Array + + function Ini_Read_Logical(this, Key, Default) + class(TIniFile) :: this + logical Ini_Read_Logical + logical, optional, intent(IN) :: Default + character(LEN=*), intent(IN) :: Key + character(LEN=:), pointer :: S + integer status + + S => this%Read_String(Key,.not. present(Default)) + if (S == '') then + call this%EmptyCheckDefault(Key,Default) + Ini_Read_Logical = Default + call this%ReadValues%Add(Key, Default) + else + if (verify(trim(S),'10TF') /= 0) then + status=1 + else + read (S,*, iostat=status) Ini_Read_Logical + end if + if (status/=0) call this%Error('error reading logical',Key) + end if + + end function Ini_Read_Logical + + + subroutine Ini_SaveReadValues(this, afile) + class(TIniFile) :: this + character(LEN=*), intent(in) :: afile + integer i + Type(TTextFile) :: F + + call F%CreateFile(afile) + + do i=1, this%ReadValues%Count + call F%Write(this%ReadValues%Items(i)%P%Name, ' = ', this%ReadValues%Items(i)%P%Value) + end do + + call F%Close() + + end subroutine Ini_SaveReadValues + + + subroutine Ini_Read_Int_Change(this, Key, Current, min, max) + class(TIniFile) :: this + character(LEN=*), intent(IN) :: Key + integer, intent(inout) :: Current + integer, optional, intent(IN) :: min, max + Current = this%Read_Int(Key,Current,min,max) + end subroutine Ini_Read_Int_Change + + subroutine Ini_Read_Double_Change(this,Key, Current, min, max) + class(TIniFile) :: this + character(LEN=*), intent(IN) :: Key + double precision, intent(inout) :: Current + double precision, optional, intent(IN) :: min,max + Current = this%Read_Double(Key, Current, min, max) + end subroutine Ini_Read_Double_Change + + subroutine Ini_Read_Real_Change(this,Key, Current, min, max) + class(TIniFile) :: this + character(LEN=*), intent(IN) :: Key + real, intent(inout) :: Current + real, optional, intent(IN) :: min,max + Current = this%Read_Real(Key, Current, min, max) + end subroutine Ini_Read_Real_Change + + subroutine Ini_Read_Logical_Change(this,Key, Current) + class(TIniFile) :: this + character(LEN=*), intent(IN) :: Key + logical, intent(inout) :: Current + Current = this%Read_Logical(Key, Current) + end subroutine Ini_Read_Logical_Change + + subroutine Ini_Read_String_Change(this,Key, Current) + class(TIniFile) :: this + character(LEN=*), intent(IN) :: Key + character(LEN=:), allocatable, intent(inout) :: Current + + Current = this%Read_String_Default(Key, Current) + + end subroutine Ini_Read_String_Change + +end module IniObjects diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..d0e96e6 --- /dev/null +++ b/Makefile @@ -0,0 +1,74 @@ +# + +all: Release Debug ReleaseMPI DebugMPI + +MPIF90C ?= mpif90 + +# There are ___ kinds of F90 flags defined: +# F90COMMONFLAGS: These flags are specified for each kind (Debug, Release) of compilation +# F90DEBUGFLAGS: These flags are only used for creating Debug artefacts +# F90RELEASEFLAGS: These flags are only used for creating Release artefacts + +# For standalone compiling set the compiler +ifortErr = $(shell which ifort >/dev/null 2>&1; echo $$?) +ifeq "$(ifortErr)" "0" + +#Intel compiler +F90C ?= ifort +F90COMMONFLAGS ?= -fpp -W0 -WB -gen_dep=$$*.d +F90DEBUGFLAGS ?= -g +F90RELEASEFLAGS ?= -fast +# Intel has a special archiver for libraries. +AR ?= xiar + +else + +F90C ?= gfortran +F90COMMONFLAGS ?= -cpp -ffree-line-length-none -fmax-errors=4 -MMD +F90DEBUGFLAGS ?= -g -O0 +F90RELEASEFLAGS ?= -O3 -ffast-math + +endif + +# When no library archiver is set yet, use ar. +AR ?= ar + +SRCS = MiscUtils.f90 StringUtils.f90 ArrayUtils.f90 MpiUtils.f90 FileUtils.f90 \ + IniObjects.f90 RandUtils.f90 ObjectLists.f90 MatrixUtils.f90 RangeUtils.f90 +OBJS = $(patsubst %.f90,%.o,$(SRCS)) + +Release: + $(MAKE) OUTPUT_DIR=Release F90FLAGS="$(F90RELEASEFLAGS)" directories + +ReleaseMPI: + $(MAKE) OUTPUT_DIR=ReleaseMPI F90C=$(MPIF90C) F90FLAGS="$(F90RELEASEFLAGS) -DMPI" directories + +Debug: + $(MAKE) OUTPUT_DIR=Debug F90FLAGS="$(F90DEBUGFLAGS)" directories + +DebugMPI: + $(MAKE) OUTPUT_DIR=DebugMPI F90C=$(MPIF90C) F90FLAGS="$(F90DEBUGFLAGS) -DMPI" directories + +%.o: $(SRC_DIR)/%.f90 + $(F90C) $(F90COMMONFLAGS) $(F90FLAGS) -o $*.o -c $< + +libutils.a: $(OBJS) + $(AR) -r libutils.a $(OBJS) + +clean: + -rm -fr Debug* Release* + +# export all variables to the sub-make below +export + +# Build in output dir to get file dependencies correctly and make use +# of delta builds, where only changed files and their dependents are +# rebuild. +directories: + mkdir -p $(OUTPUT_DIR) + $(MAKE) -C $(OUTPUT_DIR) -f../Makefile SRC_DIR=.. libutils.a + +# Include dependency files generated by a previous run of make +-include $(SRCS:.f90=.d) + +.PHONY: directories clean Release Debug libutils.a diff --git a/MatrixUtils.f90 b/MatrixUtils.f90 new file mode 100644 index 0000000..8ca3f5e --- /dev/null +++ b/MatrixUtils.f90 @@ -0,0 +1,2317 @@ + !Matrix utility routines. Uses BLAS/LAPACK. Mostly wrapper routines. + !Generally (but not always) assumes that all matrix arrays are defined at exactly correct size + !Not complete + !Antony Lewis May 2003-2014 + !http://cosmologist.info/utils/ + + + module MatrixUtils + use MpiUtils + use FileUtils + use StringUtils + implicit none + + logical, parameter :: Matrix_runmsgs = .false. +#ifdef MATRIX_SINGLE + integer, parameter :: dm = KIND(1.0) +#else + integer, parameter :: dm = KIND(1.d0) +#endif + !Precision of matrix operators + !If changing also need to change prefix on LAPACK routine names + integer, parameter :: Mat_F90=1, Mat_Norm=2, Mat_DC =3 !Normal, basic BLAS/LAPACK or divide and conquer + integer, parameter :: matrix_method = mat_DC + + real Matrix_StartTime + + Type TMatrixType + real(dm), dimension(:,:), pointer :: M + end Type TMatrixType + + complex(dm), parameter :: COne = (1._dm,0._dm), CZero = (0._dm,0._dm) + real(dm), parameter :: ROne = 1._dm, RZero = 0._dm + real, parameter :: SOne = 1., SZero = 0. + + character(LEN=*), parameter :: Matrix_IO_fmt= '(*(1E17.7))' + + contains + + + function GetMatrixTime() + real GetMatrixTime + real atime + + call cpu_time(atime) + + GetMatrixTime = atime + + + end function GetMatrixTime + + subroutine Matrix_start(Name) + character(LEN=*), intent(in) :: Name + + if (Matrix_runmsgs) then + Matrix_StartTime = GetMatrixTime() + Write(*,*) 'Matrix_'//trim(Name) //' start' + end if + end subroutine Matrix_start + + subroutine Matrix_end(Name) + character(LEN=*), intent(in) :: Name + + if (Matrix_runmsgs) then + Write(*,*) 'Matrix_'//trim(Name) //' end: ', GetMatrixTime() - Matrix_StartTime + end if + end subroutine Matrix_end + + subroutine Matrix_WriteFileRow(aunit, vec,n) + integer, intent(in) :: aunit + integer, intent(in) :: n + real(dm) :: vec(n) + + write (aunit, Matrix_IO_fmt) vec(1:n) + + end subroutine Matrix_WriteFileRow + + subroutine Matrix_Write(aname, mat, forcetable, commentline) + character(LEN=*), intent(in) :: aname + character(LEN=*), intent(in), optional :: commentline + real(dm), intent(in) :: mat(:,:) + logical, intent(in), optional :: forcetable + integer i,k + integer shp(2) + logical WriteTab + Type(TTextFile) :: F + + shp = shape(mat) + WriteTab = shp(2)<=50 + if (present(forcetable)) then + if (forcetable) WriteTab = .true. + end if + call F%CreateFile(aname) + if (present(commentline)) then + call F%Write('#'//commentline) + end if + do i=1, shp(1) + if (.not. WriteTab) then + do k=1, shp(2) + call F%Write(mat(i,k)) + end do + else + write(F%unit,F%RealFormat) mat(i,1:shp(2)) + !workaround for ifort 14 bug + !call F%Write(mat(i,1:shp(2))) + end if + end do + call F%Close() + + end subroutine Matrix_Write + + subroutine Matrix_Write_double(aname, mat, forcetable) + character(LEN=*), intent(in) :: aname + double precision, intent(in) :: mat(:,:) + logical, intent(in), optional :: forcetable + integer i,k + integer shp(2) + logical WriteTab + Type(TTextFile) F + + shp = shape(mat) + WriteTab = shp(2)<=50 + if (present(forcetable)) then + if (forcetable) WriteTab = .true. + end if + call F%CreateFile(aname) + do i=1, shp(1) + if (.not. WriteTab) then + do k=1, shp(2) + call F%Write(mat(i,k)) + end do + else + call F%Write(mat(i,1:shp(2))) + end if + end do + call F%Close() + + end subroutine Matrix_Write_double + + + subroutine Matrix_Write_Binary(aname, mat) + character(LEN=*), intent(in) :: aname + real(dm), intent(in) :: mat(:,:) + Type(TBinaryFile) F + + call F%CreateFile(aname) + Write(F%unit) mat + call F%Close() + + end subroutine Matrix_Write_Binary + + + subroutine MatrixSym_Write_Binary(aname, mat) + character(LEN=*), intent(in) :: aname + real(dm), intent(in) :: mat(:,:) + integer i + integer shp(2) + Type(TBinaryFile) F + + shp = shape(mat) + if (shp(1) /= shp(2)) call MpiStop('MatrixSym_Write_Binary: Not square matrix') + if (shp(1) == 0) return + + call F%CreateFile(aname) + do i=1,shp(1) + write(F%unit) mat(i:shp(2),i) + end do + call F%Close() + + end subroutine MatrixSym_Write_Binary + + subroutine MatrixSym_Write_Binary_Single(aname, mat) + character(LEN=*), intent(in) :: aname + real(dm), intent(in) :: mat(:,:) + integer i + integer shp(2) + Type(TBinaryFile) F + + shp = shape(mat) + if (shp(1) /= shp(2)) call MpiStop('MatrixSym_Write_Binary_Single: Not square matrix') + if (shp(1) == 0) return + + call F%CreateFile(aname) + do i=1,shp(1) + write(F%unit) real(mat(i:shp(2),i), kind(1.0)) + end do + call F%Close() + + end subroutine MatrixSym_Write_Binary_Single + + + + subroutine Matrix_WriteVec(aname, vec) + character(LEN=*), intent(in) :: aname + real(dm), intent(in) :: vec(:) + integer i + Type(TTextFile) F + + call F%CreateFile(aname) + do i=1, size(vec) + call F%Write(vec(i)) + end do + call F%Close() + + end subroutine Matrix_WriteVec + + + subroutine Matrix_Read_Binary(aname, mat) + character(LEN=*), intent(in) :: aname + real(dm), intent(out) :: mat(:,:) + Type(TBinaryFile) F + + call F%Open(aname) + read (F%unit) mat + call F%Close() + + end subroutine Matrix_Read_Binary + + + subroutine MatrixSym_Read_Binary(aname, mat) + character(LEN=*), intent(in) :: aname + real(dm), intent(out) :: mat(:,:) + integer i + integer shp(2) + Type(TBinaryFile) F + + shp = shape(mat) + if (shp(1) /= shp(2)) call MpiStop( 'MatrixSym_Read_Binary: Not square matrix') + if (shp(1) == 0) return + + call F%Open(aname) + do i=1,shp(1) + read (F%unit) mat(i:shp(1),i) + mat(i,i:shp(1)) = mat(i:shp(1),i) + end do + call F%Close() + + end subroutine MatrixSym_Read_Binary + + subroutine MatrixSym_Read_Binary_Single(aname, mat) + character(LEN=*), intent(in) :: aname + real, intent(out) :: mat(:,:) + integer i + integer shp(2) + Type(TBinaryFile) F + + shp = shape(mat) + if (shp(1) /= shp(2)) call MpiStop( 'MatrixSym_Read_Binary: Not square matrix') + if (shp(1) == 0) return + + call F%Open(aname) + do i=1,shp(1) + read (F%unit) mat(i:shp(1),i) + mat(i,i:shp(1)) = mat(i:shp(1),i) + end do + call F%Close() + + end subroutine MatrixSym_Read_Binary_Single + + + + subroutine Matrix_Read(aname, mat) + character(LEN=*), intent(IN) :: aname + real(dm), intent(out) :: mat(:,:) + integer j,k + integer shp(2) + real(dm) tmp + Type(TTextFile) :: F + + shp = shape(mat) + + call F%Open(aname) + + do j=1,shp(1) + read (F%unit,*, end = 200, err=100) mat(j,1:shp(2)) + end do + goto 120 + +100 rewind(F%unit) !Try other possible format + do j=1,shp(1) + do k=1,shp(2) + read (F%unit,*, end = 200) mat(j,k) + end do + end do + +120 read (F%unit,*, err = 150, end =150) tmp + goto 200 + +150 call F%Close() + return + +200 call MpiStop('Matrix_Read: file '//trim(aname)//' is the wrong size') + + + end subroutine Matrix_Read + + subroutine Matrix_ReadSingle(aname, mat) + character(LEN=*), intent(IN) :: aname + real, intent(out) :: mat(:,:) + integer j,k + integer shp(2) + real tmp + Type(TTextFile) :: F + + shp = shape(mat) + + call F%Open(aname) + + do j=1,shp(1) + read (F%unit,*, end = 200, err=100) mat(j,1:shp(2)) + end do + goto 120 + +100 rewind(F%unit) !Try other possible format + do j=1,shp(1) + do k=1,shp(2) + read (F%unit,*, end = 200) mat(j,k) + end do + end do + +120 read (F%unit,*, err = 150, end =150) tmp + goto 200 + +150 call F%Close() + return + +200 call MpiStop('Matrix_Read:Single file '//trim(aname)//' is the wrong size') + + + end subroutine Matrix_ReadSingle + + + function Matrix_Diag(M, n) + integer, intent(in) :: n + real(dm), intent(in) :: M(:,:) + real(dm) Matrix_Diag(n) + integer i + + do i=1,n + Matrix_Diag(i) = M(i,i) + end do + + end function Matrix_Diag + + function ILAENV_wrap(i,S1,S2,a,b,c,d) + integer ILAENV_wrap + integer, intent(in) :: i,a,b,c,d + character(LEN=*), intent(in) :: S1, S2 + integer, external :: ILAENV + + !If you don't have ILAENV in math library, change routine to return some positive integer + !that is a guess at the blocksize +#ifdef MATRIX_SINGLE + ILAENV_wrap = 16 +#else + ILAENV_wrap = ILAENV(i,S1,S2,a,b,c,d) +#endif + !!!IFC + end function ILAENV_wrap + + + subroutine Matrix_Diagonalize(M, diag, n) + !Does m = U diag U^T, returning U in M + integer, intent(in) :: n + real(dm), intent(inout):: m(n,n) + real(dm), intent(out) :: diag(n) + integer ierr, tmpsize + real(dm), allocatable, dimension(:) :: tmp + + call Matrix_Start('Diagonalize') +#ifdef MATRIX_SINGLE + tmpsize = max( (ILAENV_wrap(1,'SSYTRD','U',n,n,n,n)+2)*N,max(1,3*n-1)) !3*n**2 + allocate(tmp(tmpsize)); + call SSYEV('V','U',n,m,n,diag,tmp,tmpsize,ierr) !evalues and vectors of symmetric matrix +#else + tmpsize = max( (ILAENV_wrap(1,'DSYTRD','U',n,n,n,n)+2)*N,max(1,3*n-1)) !3*n**2 + allocate(tmp(tmpsize)); + call DSYEV('V','U',n,m,n,diag,tmp,tmpsize,ierr) !evalues and vectors of symmetric matrix +#endif + if (ierr /= 0) call MpiStop('Error in Matrix_Diagonalize') + deallocate(tmp) + call Matrix_End('Diagonalize') + + end subroutine Matrix_Diagonalize + + subroutine Matrix_Diagonalize_DC(M, diag, n) + !Complex version. Does m = U diag U^dag, returning U in M + integer, intent(in) :: n + real(dm), intent(inout):: m(n,n) + real(dm), intent(out) :: diag(n) + integer ierr, tmpsize ,isize + real(dm), allocatable, dimension(:) :: tmp + integer, allocatable,dimension(:):: iwork + + call Matrix_Start('Diagonalize') + + if (matrix_method == Mat_DC) then + !Divide and conquer + tmpsize = 1 + 6*N + 2*N**2 + isize = 3+5*N + allocate(tmp(tmpsize)) + allocate(iwork(isize)) +#ifdef MATRIX_SINGLE + call SSYEVD('V','U',n,M,n,diag,tmp,tmpsize,iwork,isize,ierr) !evalues and vectors of hermitian matrix +#else + call DSYEVD('V','U',n,M,n,diag,tmp,tmpsize,iwork,isize,ierr) !evalues and vectors of hermitian matrix +#endif + deallocate(iwork) + deallocate(tmp) + else + call Matrix_Diagonalize(M, diag, n) + end if + + if (ierr /= 0) call MpiStop('Error in Matrix_Diagonalize') + + call Matrix_End('Diagonalize') + + end subroutine Matrix_Diagonalize_DC + + + + subroutine Matrix_Root(M, n, pow) + !Does M**pow for symmetric M using U D**pow U^T + !Not optimized for large matrices + integer, intent(in) :: n + real(dm), intent(inout):: M(n,n) + real(dm) :: Tmp(n,n) + real(dm), intent(in) :: pow + + real(dm) :: diag(n) + integer i + + call Matrix_Diagonalize(M, diag, n) + Tmp = M + diag = diag**pow + do i = 1, n + M(:,i) = M(:,i)*diag(i) + end do + M = matmul(M,transpose(Tmp)) + + end subroutine Matrix_Root + + + subroutine Matrix_Diagonalize_Partial(M, diag, n, emin,emax, nfound) + !Real version. Does m = U diag U^dag, returning U in M + !Assumes up to nfound values will be found. nfound set to true value on exit + integer, intent(in) :: n + real(dm), intent(inout):: m(:,:) + real(dm), intent(out) :: diag(:) + real(dm), intent(in) :: emin,emax + integer, intent(inout) :: nfound + integer ierr, worksize, LIWork + real(dm), allocatable, dimension(:) :: work + real(dm), allocatable, dimension(:,:) :: tmp + integer, allocatable,dimension(:):: supp,iwork + real(dm) wsize(1) + real(dm) atol + integer ISize(1) + + atol = 1d-9 + call Matrix_Start('Matrix_Diagonalize_Partial') + allocate(tmp(n,nfound)) + allocate(Supp(n)) + !Query + WorkSize = -1 + LIWork = -1 +#ifdef MATRIX_SINGLE + call SSYEVR('V','V','U',n,M,Size(M,DIM=1),emin,emax,0,0,atol,nfound,diag,tmp,Size(TMP,DIM=1),& + Supp,WSize,WorkSize,ISize,LIWork,ierr ) +#else + call DSYEVR('V','V','U',n,M,Size(M,DIM=1),emin,emax,0,0,atol,nfound,diag,tmp,Size(TMP,DIM=1),& + Supp,WSize,WorkSize,ISize,LIWork,ierr ) +#endif + WorkSize = Real(WSize(1)) + LIWork = ISize(1) + allocate(Work(WorkSize),IWork(LIWork)) +#ifdef MATRIX_SINGLE + call SSYEVR('V','V','U',n,M,Size(M,DIM=1),emin,emax,0,0,atol,nfound,diag,tmp,Size(TMP,DIM=1),& + Supp,Work,WorkSize,IWork,LIWork,ierr ) +#else + call DSYEVR('V','V','U',n,M,Size(M,DIM=1),emin,emax,0,0,atol,nfound,diag,tmp,Size(TMP,DIM=1),& + Supp,Work,WorkSize,IWork,LIWork,ierr ) +#endif + deallocate(Supp,Work,IWork) + if (ierr /= 0) call MpiStop('Matrix_Diagonalize_Partial: Error') + M(1:n,1:nfound) = tmp(1:n,1:nfound) !nfound now different + deallocate(tmp) + call Matrix_End('Matrix_Diagonalize_Partial') + + end subroutine Matrix_Diagonalize_Partial + + + subroutine Matrix_CDiagonalize_Partial(M, diag, n, emin,emax, nfound) + !Complex version. Does m = U diag U^dag, returning U in M + !Assumes up to nfound values will be found. nfound set to true value on exit + integer, intent(in) :: n + complex(dm), intent(inout):: m(:,:) + real(dm), intent(out) :: diag(:) + real(dm), intent(in) :: emin,emax + integer, intent(inout) :: nfound + integer ierr, worksize, LRWork, LIWork + real(dm), allocatable, dimension(:) :: Rwork + complex(dm), allocatable, dimension(:) :: work + complex(dm), allocatable, dimension(:,:) :: tmp + integer, allocatable,dimension(:):: supp,iwork + complex(dm) wsize(1) + real(dm) Rsize(1), atol + integer ISize(1) + + atol = 1d-9 + call Matrix_Start('Matrix_CDiagonalize_Partial') + allocate(tmp(n,nfound)) + allocate(Supp(n)) + !Query + WorkSize = -1 + LRWork = -1 + LIWork = -1 +#ifdef MATRIX_SINGLE + call CHEEVR('V','V','U',n,M,Size(M,DIM=1),emin,emax,0,0,atol,nfound,diag,tmp,Size(TMP,DIM=1),& + Supp,WSize,WorkSize,RSize,LRWork,ISize,LIWork,ierr ) +#else + call ZHEEVR('V','V','U',n,M,Size(M,DIM=1),emin,emax,0,0,atol,nfound,diag,tmp,Size(TMP,DIM=1),& + Supp,WSize,WorkSize,RSize,LRWork,ISize,LIWork,ierr ) +#endif + WorkSize = Real(WSize(1)) + LRWork = RSize(1) + LIWork = ISize(1) + allocate(Work(WorkSize),RWork(LRWork),IWork(LIWork)) +#ifdef MATRIX_SINGLE + call CHEEVR('V','V','U',n,M,Size(M,DIM=1),emin,emax,0,0,atol,nfound,diag,tmp,Size(TMP,DIM=1),& + Supp,Work,WorkSize,RWork,LRWork,IWork,LIWork,ierr ) +#else + call ZHEEVR('V','V','U',n,M,Size(M,DIM=1),emin,emax,0,0,atol,nfound,diag,tmp,Size(TMP,DIM=1),& + Supp,Work,WorkSize,RWork,LRWork,IWork,LIWork,ierr ) +#endif + deallocate(Supp,Work,RWork,IWork) + if (ierr /= 0) call MpiStop('Matrix_CDiagonalize_Partial: Error') + M(1:n,1:nfound) = tmp(1:n,1:nfound) !nfound now different + deallocate(tmp) + call Matrix_End('Matrix_CDiagonalize_Partial') + + + end subroutine + + + subroutine Matrix_CDiagonalize(M, diag, n) + !Complex version. Does m = U diag U^dag, returning U in M + integer, intent(in) :: n + complex(dm), intent(inout):: m(n,n) + real(dm), intent(out) :: diag(n) + integer ierr, tmpsize ,isize, rworksize + real(dm), allocatable, dimension(:) :: Rwork + complex(dm), allocatable, dimension(:) :: tmp + integer, allocatable,dimension(:):: iwork + + call Matrix_Start('CDiagonalize') + + if (matrix_method == Mat_DC) then + !Divide and conquer + tmpsize = 2*N + N**2 + rworksize = 1 + 4*N + 2*N*int(log(real(N))/log(2.)+1) + 3*N**2 + isize = (2 + 5*N)*4 + allocate(tmp(tmpsize),rwork(rworksize)) + allocate(iwork(isize)) +#ifdef MATRIX_SINGLE + call CHEEVD('V','U',n,M,n,diag,tmp,tmpsize,Rwork,rworksize,iwork,isize,ierr) !evalues and vectors of hermitian matrix +#else + call ZHEEVD('V','U',n,M,n,diag,tmp,tmpsize,Rwork,rworksize,iwork,isize,ierr) !evalues and vectors of hermitian matrix +#endif + deallocate(iwork) + else + rworksize = max(1, 3*n-2) +#ifdef MATRIX_SINGLE + tmpsize = max( (ILAENV_wrap(1,'CHETRD','U',n,n,n,n)+1)*N,max(1,2*n-1)) ! 3*n**2 + allocate(tmp(tmpsize),rwork(rworksize)); + call CHEEV('V','U',n,m,n,diag,tmp,tmpsize,Rwork,ierr) !evalues and vectors of hermitian matrix +#else + tmpsize = max( (ILAENV_wrap(1,'ZHETRD','U',n,n,n,n)+1)*N,max(1,2*n-1)) ! 3*n**2 + allocate(tmp(tmpsize),rwork(rworksize)); + call ZHEEV('V','U',n,m,n,diag,tmp,tmpsize,Rwork,ierr) !evalues and vectors of hermitian matrix +#endif + end if + + if (ierr /= 0) call MpiStop('Error in Matrix_CDiagonalize') + deallocate(tmp,rwork) + + call Matrix_End('CDiagonalize') + + end subroutine Matrix_CDiagonalize + + function Matrix_CTrace(M) + complex(dm), intent(in) :: M(:,:) + complex(dm) tmp,Matrix_CTrace + integer i + + if (size(M,dim=1) /= size(M,dim=2)) call MpiStop('Matrix_CTrace: non-square matrix') + tmp =0 + do i=1,size(M,dim=1) + tmp = tmp + M(i,i) + end do + Matrix_CTrace = tmp + + end function Matrix_CTrace + + function Matrix_Trace(M) + real(dm), intent(in) :: M(:,:) + real(dm) tmp,Matrix_Trace + integer i + + if (size(M,dim=1) /= size(M,dim=2)) call mpiStop('Matrix_Trace: non-square matrix') + tmp =0 + do i=1,size(M,dim=1) + tmp = tmp + M(i,i) + end do + Matrix_Trace = tmp + + end function Matrix_Trace + + + function MatrixSym_LogDet(mat) result (logDet) + real(dm), intent(in) :: mat(:,:) + real(dm) logDet + real(dm) Tmp(size(mat,dim=1),size(mat,dim=1)) + integer i + + if (size(mat,dim=1) /= size(mat,dim=2)) call mpiStop('MatrixSym_LogDet: non-square matrix') + Tmp = mat + call Matrix_Cholesky(tmp) + logDet =0 + do i=1, size(mat,dim=1) + logDet = logDet + log(tmp(i,i)) + end do + logDet = 2._dm*logDet + + end function MatrixSym_LogDet + + + subroutine Matrix_CRotateSymm(Mat,U,m,Out,triangular) + !Gets U^dag Mat U + integer, intent(in) ::m + complex(dm), intent(in) :: Mat(:,:),U(:,:) + complex(dm) Out(:,:) + complex(dm), dimension(:,:), allocatable :: C + integer n + logical, intent(in), optional :: triangular + logical :: triang + + call Matrix_Start('CRotateSymm') + + if (present(triangular)) then + triang=triangular + else + triang=.false. + end if + + n = Size(Mat,DIM=1) + if (n /= Size(Mat,DIM=2)) call mpiStop('Matrix_CRotateSymm: Need square matrix') + if (n /= Size(U,DIM=1)) call MpiStop('Matrix_CRotateSymm: Matrix size mismatch') + if (Size(Out,DIM=1) < m .or. Size(Out,DIM=2) < m) & + call MpiStop('Matrix_CRotateSymm: Wrong output size') + + if (matrix_method == Mat_F90) then + Out = matmul(matmul(transpose(conjg(U(1:n,1:m))),Mat),U(1:n,1:m)) + else +#ifdef MATRIX_SINGLE + if (triang) then + if (m/=n) call MpiStop('Matrix_CRotateSymm: Matrices must be same size') + call CHEMM('L','U',n,n,COne,Mat,Size(Mat,DIM=1),U,Size(U,DIM=1),CZero,Out,Size(Out,Dim=1)) + call CTRMM('Left','Upper','Complex-Transpose','Not-unit',n,n,COne,U,Size(U,DIM=1),Out,Size(Out,Dim=1)) + else + allocate(C(n,m)) + call CHEMM('L','U',n,m,COne,Mat,Size(Mat,DIM=1),U,Size(U,DIM=1),CZero,C,n) + call CGEMM('C','N',m,m,n,COne,U,Size(U,DIM=1),C,n,CZero,Out,Size(Out,Dim=1)) + deallocate(C) + end if +#else + if (triang) then + if (m/=n) call MpiStop('Matrix_CRotateSymm: Matrices must be same size') + call ZHEMM('L','U',n,n,COne,Mat,Size(Mat,DIM=1),U,Size(U,DIM=1),CZero,Out,Size(Out,Dim=1)) + call ZTRMM('Left','Upper','Complex-Transpose','Not-unit',n,n,COne,U,Size(U,DIM=1),Out,Size(Out,Dim=1)) + else + allocate(C(n,m)) + call ZHEMM('L','U',n,m,COne,Mat,Size(Mat,DIM=1),U,Size(U,DIM=1),CZero,C,n) + call ZGEMM('C','N',m,m,n,COne,U,Size(U,DIM=1),C,n,CZero,Out,Size(Out,Dim=1)) + deallocate(C) + end if +#endif + end if + call Matrix_End('CRotateSymm') + + + end subroutine Matrix_CRotateSymm + + subroutine Matrix_RotateSymm(Mat,U,m,Out, triangular) + !Gets U^T Mat U + !If triangular U = Upper triangular (U^T lower triangular) + integer, intent(in) ::m + real(dm), intent(in) :: Mat(:,:),U(:,:) + real(dm) Out(:,:) + real(dm), dimension(:,:), allocatable :: C + logical, intent(in), optional :: triangular + logical triang + integer n + + call Matrix_Start('RotateSymm') + + if (present(triangular)) then + triang=triangular + else + triang=.false. + end if + + n = Size(Mat,DIM=1) + if (n /= Size(Mat,DIM=2)) call MpiStop('Matrix_RotateSymm: Need square matrix') + if (n /= Size(U,DIM=1)) call MpiStop('Matrix_RotateSymm: Matrix size mismatch') + if (Size(Out,DIM=1) < m .or. Size(Out,DIM=2) < m) & + call MpiStop('Matrix_RotateSymm: Wrong output size') + + if (matrix_method == Mat_F90) then + Out = matmul(matmul(transpose(U(1:n,1:m)),Mat),U(1:n,1:m)) + else +#ifdef MATRIX_SINGLE + if (triang) then + if (m/=n) call MpiStop('Matrix_RotateSymm: Matrices must be same size') + call SSYMM('L','U',n,n,ROne,Mat,Size(Mat,DIM=1),U,Size(U,DIM=1),RZero,Out,Size(Out,Dim=1)) + call STRMM('Left','Upper','Transpose','Not-unit',n,n,ROne,U,Size(U,DIM=1),Out,Size(Out,Dim=1)) + else + allocate(C(n,m)) + call SSYMM('L','U',n,m,ROne,Mat,Size(Mat,DIM=1),U,Size(U,DIM=1),RZero,C,n) + call SGEMM('T','N',m,m,n,ROne,U,Size(U,DIM=1),C,n,RZero,Out,Size(Out,Dim=1)) + deallocate(C) + end if +#else + if (triang) then + if (m/=n) call MpiStop('Matrix_RotateSymm: Matrices must be same size') + call DSYMM('L','U',n,n,ROne,Mat,Size(Mat,DIM=1),U,Size(U,DIM=1),RZero,Out,Size(Out,Dim=1)) + call DTRMM('Left','Upper','Transpose','Not-unit',n,n,ROne,U,Size(U,DIM=1),Out,Size(Out,Dim=1)) + else + allocate(C(n,m)) + call DSYMM('L','U',n,m,ROne,Mat,Size(Mat,DIM=1),U,Size(U,DIM=1),RZero,C,n) + call DGEMM('T','N',m,m,n,ROne,U,Size(U,DIM=1),C,n,RZero,Out,Size(Out,Dim=1)) + deallocate(C) + end if +#endif + end if + call Matrix_End('RotateSymm') + + + end subroutine Matrix_RotateSymm + + + subroutine Matrix_RotateAntiSymm(Mat,U,m,Out) + !Gets U^T Mat U + !Where Mat = -Mat^T + integer, intent(in) ::m + real(dm), intent(in) :: Mat(:,:),U(:,:) + real(dm) Out(:,:) + real(dm), dimension(:,:), allocatable :: C + integer i,j,n + + call Matrix_Start('RotateAntiSymm') + + n = Size(Mat,DIM=1) + if (n /= Size(Mat,DIM=2)) call MpiStop('Matrix_RotateAntiSymm: Need square matrix') + if (n /= Size(U,DIM=1)) call MpiStop('Matrix_RotateAntiSymm: Matrix size mismatch') + if (Size(Out,DIM=1) < m .or. Size(Out,DIM=2) < m) & + call MpiStop('Matrix_RotateAntiSymm: Wrong output size') + + if (matrix_method == Mat_F90) then + Out = matmul(matmul(transpose(U(1:n,1:m)),Mat),U(1:n,1:m)) + else + allocate(C(n,m)) + C = U(1:n,1:m) +#ifdef MATRIX_SINGLE + call STRMM('Left','Lower','Not-Transpose','Not-unit',n,m,ROne,Mat,Size(Mat,DIM=1),C,Size(C,Dim=1)) + call SGEMM('T','N',m,m,n,ROne,U,Size(U,DIM=1),C,n,RZero,Out,Size(Out,Dim=1)) +#else + call DTRMM('Left','Lower','Not-Transpose','Not-unit',n,m,ROne,Mat,Size(Mat,DIM=1),C,Size(C,Dim=1)) + call DGEMM('T','N',m,m,n,ROne,U,Size(U,DIM=1),C,n,RZero,Out,Size(Out,Dim=1)) +#endif + deallocate(C) + end if + + do i=1, m + do j=1,i + Out(j,i) = Out(j,i) - Out(i,j) + out(i,j) = -Out(j,i) + end do + end do + + call Matrix_End('RotateAntiSymm') + + end subroutine Matrix_RotateAntiSymm + + subroutine Matrix_CMult_SymmRight(Mat,U,Out,a,b) + complex(dm), intent(in) :: Mat(:,:),U(:,:) + complex(dm) Out(:,:) + complex(dm), intent(in), optional :: a,b + complex(dm) mult, beta + integer n,m + + call Matrix_Start('CMult_SymmRight') + + m = Size(Mat,DIM=1) + n = Size(U,DIM=2) + if (n /= Size(Mat,DIM=2) .or. n/=Size(U,DIM=1)) & + call MpiStop('Matrix_CMult_SymmRight: Size mismatch') + if (present(a)) then + mult = a + else + mult = COne + end if + if (present(b)) then + beta = b + else + beta = CZero + end if + if (matrix_method == Mat_F90) then + if (beta /= CZero) then + out = a*MatMul(Mat,U) + beta*Out + else + out = MatMul(Mat,U) + if (mult /= COne) Out = Out*mult + end if + else +#ifdef MATRIX_SINGLE + call CHEMM('R','U',m,n,mult,U,Size(U,DIM=1),Mat,Size(Mat,DIM=1),beta,Out,Size(Out,DIM=1)) +#else + call ZHEMM('R','U',m,n,mult,U,Size(U,DIM=1),Mat,Size(Mat,DIM=1),beta,Out,Size(Out,DIM=1)) +#endif + end if + + call Matrix_End('CMult_SymmRight') + + end subroutine Matrix_CMult_SymmRight + + + subroutine Matrix_CMult_SymmLeft(Mat,U,Out,a,b) + complex(dm), intent(in) :: Mat(:,:),U(:,:) + complex(dm) Out(:,:) + complex(dm), intent(in), optional :: a,b + complex(dm) mult, beta + integer n,m + + call Matrix_Start('CMult_SymmLeft') + + m = Size(Mat,DIM=1) + n = Size(U,DIM=2) + if (m /= Size(U,DIM=1) .or. m/=Size(Mat,DIM=2)) & + call MpiStop('Matrix_CMult_SymmLeft: Size mismatch') + if (present(a)) then + mult = a + else + mult = COne + end if + if (present(b)) then + beta = b + else + beta = CZero + end if + if (matrix_method == Mat_F90) then + if (beta /= CZero) then + out = a*MatMul(Mat,U) + beta*Out + else + out = MatMul(Mat,U) + if (mult /= COne) Out = Out*mult + end if + else +#ifdef MATRIX_SINGLE + call CHEMM('L','U',m,n,mult,Mat,Size(Mat,DIM=1),U,Size(U,DIM=1),beta,Out,Size(Out,DIM=1)) +#else + call ZHEMM('L','U',m,n,mult,Mat,Size(Mat,DIM=1),U,Size(U,DIM=1),beta,Out,Size(Out,DIM=1)) +#endif + end if + + call Matrix_End('CMult_SymmLeft') + + end subroutine Matrix_CMult_SymmLeft + + + subroutine Matrix_CMult(Mat,U,Out,a,b) + ! Out = a*Mat U + b*out + complex(dm), intent(in) :: Mat(:,:),U(:,:) + complex(dm) Out(:,:) + complex(dm), intent(in), optional :: a,b + complex(dm) mult, beta + integer m,n,k + + call Matrix_Start('CMult') + + m = Size(Mat,DIM=1) + n = Size(U,DIM=2) + k = Size(Mat,DIM=2) + if (k /= Size(U,DIM=1)) call MpiStop('Matrix_Mult: Matrix size mismatch') + if (present(a)) then + mult = a + else + mult = COne + end if + if (present(b)) then + beta = b + else + beta = CZero + end if + + if (matrix_method == Mat_F90) then + if (beta /= CZero) then + out = a*MatMul(Mat,U) + beta*Out + else + out = MatMul(Mat,U) + if (mult /= COne) Out = Out*mult + end if + else +#ifdef MATRIX_SINGLE + call CGEMM('N','N',m,n,k,mult,Mat,m,U,k,beta,Out,Size(Out,Dim=1)) +#else + call ZGEMM('N','N',m,n,k,mult,Mat,m,U,k,beta,Out,Size(Out,Dim=1)) +#endif + end if + call Matrix_End('CMult') + + + end subroutine Matrix_CMult + + + subroutine Matrix_MultSq_RepRight(Mat,U,a) + !U = a*Mat*U + real(dm), intent(in) :: Mat(:,:) + real(dm), intent(inout) ::U(:,:) + real(dm), intent(in), optional :: a + real(dm) aa + integer m,n + real(dm), dimension(:,:), allocatable :: tmp + + + m = Size(Mat,DIM=1) + n = Size(Mat,DIM=2) + if (m /= n) call MpiStop('Matrix_MultSq: Matrix size mismatch') + m = Size(U,DIM=1) + n = Size(U,DIM=2) + if (m /= n) call MpiStop('Matrix_MultSq: Matrix size mismatch') + + allocate(tmp(n,n)) + if (present(a)) then + aa=a + else + aa=ROne + end if + + call Matrix_Mult(Mat,U,tmp,aa) + U = tmp + deallocate(tmp) + + end subroutine Matrix_MultSq_RepRight + + subroutine Matrix_MultTri(Mat,L, side) + ! Mat -> L Mat or Mat L where L is lower triangular + real(dm), intent(inout) :: Mat(:,:) + real(dm), intent(in) :: L(:,:) + character(LEN=*), intent(in) :: side + integer m,n + + call Matrix_Start('Matrix_MultTri') + + m = Size(Mat,DIM=1) + n = Size(Mat,DIM=2) + + if (side(1:1)=='L') then + if (Size(L,DIM=2) /= m) call MpiStop('Matrix_MultTri: Matrix size mismatch') + else + if (Size(L,DIM=1) /= n) call MpiStop('Matrix_MultTri: Matrix size mismatch') + end if +#ifdef MATRIX_SINGLE + call STRMM(side,'Lower','Not-Transpose','Not-unit',m,n,ROne,L,Size(L,DIM=1),Mat,Size(Mat,Dim=1)) +#else + call DTRMM(side,'Lower','Not-Transpose','Not-unit',m,n,ROne,L,Size(L,DIM=1),Mat,Size(Mat,Dim=1)) +#endif + call Matrix_End('Matrix_MultTri') + + end subroutine Matrix_MultTri + + + + subroutine Matrix_Mult(Mat,U,Out,a,b) + ! Out = a*Mat U + b*out + real(dm), intent(in) :: Mat(:,:),U(:,:) + real(dm) :: Out(:,:) + real(dm), intent(in), optional :: a,b + real(dm) mult, beta + integer m,n,k + + call Matrix_Start('Mult') + + + m = Size(Mat,DIM=1) + n = Size(U,DIM=2) + k = Size(Mat,DIM=2) + if (k /= Size(U,DIM=1)) call MpiStop('Matrix_Mult: Matrix size mismatch') + + + if (present(a)) then + mult = a + else + mult = ROne + end if + if (present(b)) then + beta = b + else + beta = RZero + end if + + if (matrix_method == Mat_F90) then + if (beta /= RZero) then + out = a*MatMul(Mat,U) + beta*Out + else + out = MatMul(Mat,U) + if (mult /= ROne) Out = Out*mult + end if + else +#ifdef MATRIX_SINGLE + call SGEMM('N','N',m,n,k,mult,Mat,m,U,k,beta,Out,Size(Out,Dim=1)) +#else + call DGEMM('N','N',m,n,k,mult,Mat,m,U,k,beta,Out,Size(Out,Dim=1)) +#endif + end if + call Matrix_End('Mult') + + + end subroutine Matrix_Mult + + subroutine Matrix_Mult_SymmLeft(Mat,U,Out,a,b) + real(dm), intent(in) :: Mat(:,:),U(:,:) + real(dm) Out(:,:) + real(dm), intent(in), optional :: a,b + real(dm) mult, beta + integer n,m + + call Matrix_Start('Mult_SymmLeft') + + m = Size(Mat,DIM=1) + n = Size(U,DIM=2) + if (m /= Size(U,DIM=1) .or. m/=Size(Mat,DIM=2)) & + call MpiStop('Matrix_Mult_SymmLeft: Size mismatch') + if (present(a)) then + mult = a + else + mult = ROne + end if + if (present(b)) then + beta = b + else + beta = RZero + end if + if (matrix_method == Mat_F90) then + if (beta /= RZero) then + out = a*MatMul(Mat,U) + beta*Out + else + out = MatMul(Mat,U) + if (mult /= COne) Out = Out*mult + end if + else +#ifdef MATRIX_SINGLE + call SSYMM('L','U',m,n,mult,Mat,Size(Mat,DIM=1),U,Size(U,DIM=1),beta,Out,Size(Out,DIM=1)) +#else + call DSYMM('L','U',m,n,mult,Mat,Size(Mat,DIM=1),U,Size(U,DIM=1),beta,Out,Size(Out,DIM=1)) +#endif + end if + + call Matrix_End('Mult_SymmLeft') + + end subroutine Matrix_Mult_SymmLeft + + + subroutine Matrix_Mult_SymmRight(Mat,U,Out,a,b) + ! Out = a*Mat U + b*out + real(dm), intent(in) :: Mat(:,:),U(:,:) + real(dm) Out(:,:) + real(dm), intent(in), optional :: a,b + real(dm) mult, beta + integer n,m + + call Matrix_Start('Mult_SymmRight') + + m = Size(Mat,DIM=1) + n = Size(U,DIM=2) + if (n /= Size(Mat,DIM=2) .or. n/=Size(U,DIM=1)) & + call MpiStop('Matrix_Mult_SymmRight: Size mismatch') + if (present(a)) then + mult = a + else + mult = ROne + end if + if (present(b)) then + beta = b + else + beta = RZero + end if + if (matrix_method == Mat_F90) then + if (beta /= RZero) then + out = a*MatMul(Mat,U) + beta*Out + else + out = MatMul(Mat,U) + if (mult /= ROne) Out = Out*mult + end if + else +#ifdef MATRIX_SINGLE + call SSYMM('R','U',m,n,mult,U,Size(U,DIM=1),Mat,Size(Mat,DIM=1),beta,Out,Size(Out,DIM=1)) +#else + call DSYMM('R','U',m,n,mult,U,Size(U,DIM=1),Mat,Size(Mat,DIM=1),beta,Out,Size(Out,DIM=1)) +#endif + end if + + call Matrix_End('Mult_SymmRight') + + end subroutine Matrix_Mult_SymmRight + + + subroutine Matrix_CMultGen(Mat,m,k,U,n,Out) + ! out(1:m,1:n) = MatMul(Mat(1:m,1:k),U(1:k,1:n)) + integer, intent(in) :: m,k,n + complex(dm), intent(in) :: Mat(:,:),U(:,:) + complex(dm) Out(:,:) + + call Matrix_Start('CMultGen') + + if (SIZE(Out,DIM=1) m. ') + + call Matrix_Start('SVD_VT') + + if (present(U) .and. Matrix_method == Mat_DC) then + !Use divide and conquer + allocate(IWork(8*MIN(M,N))) + WorkSize= -1 !3*min(M,N)*min(M,N) +max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)) +#ifdef MATRIX_SINGLE + call SGESDD('O',m,n, Mat, m ,D,U,m,Mat,n,OptWk,WorkSize,IWork,ierr) +#else + call DGESDD('O',m,n, Mat, m ,D,U,m,Mat,n,OptWk,WorkSize,IWork,ierr) +#endif + WorkSize = nint(OptWk) + allocate(rv1(WorkSize)) +#ifdef MATRIX_SINGLE + call SGESDD('O',m,n, Mat, m ,D,U,m,Mat,n,rv1,WorkSize,IWork,ierr) +#else + call DGESDD('O',m,n, Mat, m ,D,U,m,Mat,n,rv1,WorkSize,IWork,ierr) +#endif + deallocate(IWOrk) + else + call MpiStop('Matrix_SVD_VT Not no-U non-DC case') + end if + + if (ierr/=0) call MpiStop('error in Matrix_SVD_VT') + deallocate(rv1) + + call Matrix_End('SVD_VT') + + end subroutine Matrix_SVD_VT + + + subroutine Matrix_CSVD_VT(Mat,m, n, D, U) + !Do singular value decomposition of m x n matrix Mat + !Mat = U D V^dag + !returns V^dag in Mat, vector D of diagonal elements of, unitary matrix U + integer, intent(in) :: m,n + complex(dm), intent(inout) :: Mat(m,n) + complex(dm), intent(out),optional :: U(m,m) + real(dm), intent(out) :: D(*) + + integer WorkSize, ierr + integer,allocatable, dimension (:) :: IWork + complex(dm), allocatable, dimension (:) :: rv1 + real(dm), allocatable, dimension (:) :: rwork + + if (n<=m) call MpiStop('Matrix_CSVD_VT assumed n>m. If equal use SVD_U.') + + call Matrix_Start('CSVD_VT') + + + if (present(U) .and. Matrix_method == Mat_DC) then + !Use divide and conquer + WorkSize= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N) + 5*N !Add on 5N.. + allocate(rv1(WorkSize)) + allocate(rwork(5*min(M,N)*min(M,N) + 5*min(M,N) )) + allocate(IWork(8*MIN(M,N))) +#ifdef MATRIX_SINGLE + call CGESDD('O',m,n, Mat, m ,D,U,m,Mat,n,rv1,WorkSize,rwork,IWork,ierr) +#else + call ZGESDD('O',m,n, Mat, m ,D,U,m,Mat,n,rv1,WorkSize,rwork,IWork,ierr) +#endif + deallocate(IWOrk) + else + allocate(rwork((max(3*min(m,n),5*min(m,n)-4)))) + WorkSize= 3*max(m,n)**2 + allocate(rv1(WorkSize), STAT = ierr) + if (ierr /=0) then + WorkSize= MAX(3*MIN(M,N)+MAX(M,N),5*MIN(M,N)) + allocate(rv1(WorkSize)) + end if +#ifdef MATRIX_SINGLE + if (present(U)) then + call CGESVD('S','O',m,n, Mat, m , D,U,m,Mat,n,rv1,WorkSize,rwork,ierr) + else + call CGESVD('N','O',m,n, Mat, m , D,Mat,m,Mat,n,rv1,WorkSize,rwork,ierr) + end if +#else + if (present(U)) then + call ZGESVD('S','O',m,n, Mat, m , D,U,m,Mat,n,rv1,WorkSize,rwork,ierr) + else + call ZGESVD('N','O',m,n, Mat, m , D,Mat,m,Mat,n,rv1,WorkSize,rwork,ierr) + end if +#endif + end if + + if (ierr/=0) call MpiStop('error in Matrix_SVD_VT') + deallocate(rv1) + + call Matrix_End('SVD_VT') + + end subroutine Matrix_CSVD_VT + + subroutine Matrix_CSVD_U(Mat,m, n, D, VT) + !Do singular value decomposition of m x n matrix Mat + !Mat = U D VT + !returns U in Mat, vector D of diagonal elements of, unitary matrix V + integer, intent(in) :: m,n + complex(dm), intent(inout) :: Mat(m,n) + complex(dm), intent(out),optional :: VT(n,n) + real(dm), intent(out) :: D(*) + integer WorkSize, ierr + integer,allocatable, dimension (:) :: IWork + complex(dm), allocatable, dimension (:) :: rv1 + real(dm), allocatable, dimension (:) :: rwork + + call Matrix_Start('CSVD_U') + + if (m=n') + + if (present(VT) .and. Matrix_method == Mat_DC) then + WorkSize= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N) + 5*N !Add on 5N.. + allocate(rv1(WorkSize)) + allocate(rwork(5*min(M,N)*min(M,N) + 5*min(M,N) )) + allocate(IWork(8*MIN(M,N))) +#ifdef MATRIX_SINGLE + call CGESDD('O',m,n, Mat, m ,D,Mat,m,VT,n,rv1,WorkSize,rwork,IWork,ierr) +#else + call ZGESDD('O',m,n, Mat, m ,D,Mat,m,VT,n,rv1,WorkSize,rwork,IWork,ierr) +#endif + deallocate(IWOrk) + else + allocate(rwork((max(3*min(m,n),5*min(m,n)-4)))) + + WorkSize= 3*max(m,n)**2 + allocate(rv1(WorkSize), STAT = ierr) + if (ierr /=0) then + WorkSize= MAX(3*MIN(M,N)+MAX(M,N),5*MIN(M,N)) + allocate(rv1(WorkSize)) + end if +#ifdef MATRIX_SINGLE + if (present(VT)) then + call CGESVD('O','S',m,n, Mat, m , D,Mat,m,VT,n,rv1,WorkSize,rwork,ierr) + else + call CGESVD('O','N',m,n, Mat, m , D,Mat,m,Mat,n,rv1,WorkSize,rwork,ierr) + end if +#else + if (present(VT)) then + call ZGESVD('O','S',m,n, Mat, m , D,Mat,m,VT,n,rv1,WorkSize,rwork,ierr) + else + call ZGESVD('O','N',m,n, Mat, m , D,Mat,m,Mat,n,rv1,WorkSize,rwork,ierr) + end if +#endif + end if + if (ierr/=0) call MpiStop('error in Matrix_SVD_U') + call Matrix_End('CSVD_U') + + deallocate(rv1,rwork) + + end subroutine Matrix_CSVD_U + + + subroutine Matrix_CSVD_AllVT(Mat,m, n, D, VT) + !n>m + !Do singular value decomposition of m x n matrix Mat + !Mat = U D V^dag + !returns all nxn V^dag in VT, vector D of diagonal elements of + integer, intent(in) :: m,n + complex(dm), intent(inout) :: Mat(m,n) + complex(dm), intent(out):: VT(n,n) + real(dm), intent(out) :: D(m) + + integer WorkSize, ierr + complex(dm), allocatable, dimension (:) :: rv1 + complex(dm), allocatable, dimension (:,:) :: U + integer, allocatable, dimension(:) :: IWork + real(dm), allocatable, dimension (:) :: rwork + + + call Matrix_Start('CSVD_AllVT') + + if (Matrix_method == Mat_DC) then + !Divide and conquer doesn't seem to provide outputs we want here + WorkSize= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N) + 5*N !Add on 5N.. + allocate(rv1(WorkSize)) + allocate(rwork(5*min(M,N)*min(M,N) + 5*min(M,N) )) + allocate(IWork(8*MIN(M,N))) + allocate(U(m,m)) +#ifdef MATRIX_SINGLE + call CGESDD('A',m,n, Mat, m ,D,U,m,VT,n,rv1,WorkSize,rwork,IWork,ierr) +#else + call ZGESDD('A',m,n, Mat, m ,D,U,m,VT,n,rv1,WorkSize,rwork,IWork,ierr) +#endif + deallocate(U) + deallocate(IWork) + else + WorkSize= 2*m*n + 2*max(n,m) + allocate(rwork(5*max(m,n))) + allocate(rv1(WorkSize), STAT = ierr) + if (ierr /=0) then + WorkSize= MAX(3*MIN(M,N)+MAX(M,N),5*MIN(M,N)) + allocate(rv1(WorkSize)) + end if +#ifdef MATRIX_SINGLE + call CGESVD('N','A',m,n, Mat, m , D,Mat,m,VT,n,rv1,WorkSize,rwork,ierr) +#else + call ZGESVD('N','A',m,n, Mat, m , D,Mat,m,VT,n,rv1,WorkSize,rwork,ierr) +#endif + end if + + if (ierr/=0) call MpiStop('error in Matrix_SVD_AllVT') + deallocate(rv1,rwork) + + call Matrix_End('CSVD_AllVT') + + + end subroutine Matrix_CSVD_allVT + + + subroutine Matrix_DiagPreMul(D,M) + ! M -> matmul(diag(D),M) + real(dm), intent(inout) :: M(:,:) + real(dm), intent(in) :: D(:) + integer i + + if (Size(D) /= SiZE(M,DIM=1)) call MpiStop('Matrix_DiagPreMul: Wrong size') + do i = 1, size(D) + M(i,:) = M(i,:)*D(i) + end do + + end subroutine Matrix_DiagPreMul + + + subroutine Matrix_SolveSymm(M,a,soln) + real(dm), intent(out) :: soln(:) + real(dm), intent(in):: M(:,:),a(:) + integer IPIV(size(a)),info + real(dm), dimension(:,:), allocatable :: tmp + real(dm), dimension(:), allocatable :: work + integer n, WorkSize + + n=Size(M,DIM=1) + if (n<=1) return + if (Size(M,DIM=2)/=n) call MpiStop('Matrix_SolveSq: non-square matrix') + call Matrix_Start('SolveSymm') + + + WorkSize = n**2 + allocate(work(WorkSize)) + allocate(tmp(n,n)) + tmp = M +#ifdef MATRIX_SINGLE + call SSYTRF('U',n,tmp,n,IPIV, work,WorkSize,info) +#else + call DSYTRF('U',n,tmp,n,IPIV, work,WorkSize,info) +#endif + deallocate(work) + if (info/=0) call MpiStop('error in SolveSymm') + soln(1:n) = a(1:n) +#ifdef MATRIX_SINGLE + call SSYTRS('U',n,1,tmp,n,IPIV,soln,n,info) +#else + call DSYTRS('U',n,1,tmp,n,IPIV,soln,n,info) +#endif + if (info/=0) call MpiStop('error (2) in SolveSymm') + deallocate(tmp) + + call Matrix_End('SolveSymm') + + + end subroutine Matrix_SolveSymm + + + subroutine Matrix_SolveASymm(M,a,soln) + real(dm), intent(out) :: soln(:) + real(dm), intent(in):: M(:,:),a(:) + integer IPIV(size(a)),info + real(dm), dimension(:,:), allocatable :: tmp + integer n + + n=Size(M,DIM=1) + if (n<=1) return + if (Size(M,DIM=2)/=n) call MpiStop('Matrix_SolveSq: non-square matrix') + + call Matrix_Start('SolveASymm') + + allocate(tmp(n,n)) + tmp = M +#ifdef MATRIX_SINGLE + call SGETRF(n,n,tmp,n,IPIV, info) +#else + call DGETRF(n,n,tmp,n,IPIV, info) +#endif + if (info/=0) call MpiStop('error in SolveASymm') + soln(1:n) = a(1:n) +#ifdef MATRIX_SINGLE + call SGETRS('N',n,1,tmp,n,IPIV,Soln,n,info) +#else + call DGETRS('N',n,1,tmp,n,IPIV,Soln,n,info) +#endif + if (info/=0) call MpiStop('error (2) in SolveASymm') + deallocate(tmp) + + call Matrix_End('SolveASymm') + + end subroutine Matrix_SolveASymm + + function Matrix_vecdot(vec1,vec2) + real(dm) vec1(:),vec2(:) + real(dm) Matrix_vecdot + integer n +#ifdef MATRIX_SINGLE + real(dm) sdot + external sdot +#else + real(dm) ddot + external ddot +#endif + n=size(vec1) + if (n/=size(vec2)) call MpiStop('Matrix_vecdot: size mismatch') +#ifdef MATRIX_SINGLE + Matrix_vecdot = sdot(n, vec1, 1, vec2, 1) +#else + Matrix_vecdot = ddot(n, vec1, 1, vec2, 1) +#endif + end function Matrix_vecdot + + function Matrix_QuadForm(Mat,vec) + !Get vec^T*Mat*vec where Mat is symmetric + real(dm) Matrix_QuadForm + real(dm) vec(:) + real(dm) Mat(:,:) + real(dm), dimension(:), allocatable :: out + integer n + + n=size(vec) + allocate(out(n)) + call Matrix_MulVecSymm(Mat,vec,out) + Matrix_QuadForm = Matrix_vecdot(vec, out) + deallocate(out) + + end function Matrix_QuadForm + + subroutine Matrix_MulVec(Mat,vec,Out,a,b) + ! Out = a*Mat*vec + b*out + real(dm), intent(in) :: Mat(:,:) + real(dm) vec(:) + real(dm) Out(:) + real(dm), intent(in), optional :: a,b + real(dm) mult, beta + integer m,n + + call Matrix_Start('MulVec') + + m = Size(Mat,DIM=1) + n = Size(Vec) + if (Size(Mat,DIM=2) /= n) call MpiStop('Matrix_MulVec: size mismatch') + if (present(a)) then + mult = a + else + mult = ROne + end if + if (present(b)) then + beta = b + else + beta = RZero + end if + + if (matrix_method == Mat_F90) then + if (beta /= RZero) then + out = a*MatMul(Mat,Vec) + beta*Out + else + out = MatMul(Mat,Vec) + if (mult /= ROne) Out = Out*mult + end if + else +#ifdef MATRIX_SINGLE + call SGEMV('N',m,n,mult,Mat,m,vec, 1,beta, Out,1) +#else + call DGEMV('N',m,n,mult,Mat,m,vec, 1,beta, Out,1) +#endif + end if + call Matrix_End('MulVec') + + end subroutine Matrix_MulVec + + subroutine Matrix_MulVecSingle(Mat,vec,Out,a,b) + ! Out = a*Mat*vec + b*out + real, intent(in) :: Mat(:,:) + real vec(:) + real Out(:) + real, intent(in), optional :: a,b + real mult, beta + integer m,n + + call Matrix_Start('MulVecSingle') + + m = Size(Mat,DIM=1) + n = Size(Vec) + if (Size(Mat,DIM=2) /= n) call MpiStop('Matrix_MulVecSingle: size mismatch') + if (present(a)) then + mult = a + else + mult = SOne + end if + if (present(b)) then + beta = b + else + beta = SZero + end if + + if (matrix_method == Mat_F90) then + if (beta /= SZero) then + out = a*MatMul(Mat,Vec) + beta*Out + else + out = MatMul(Mat,Vec) + if (mult /= SOne) Out = Out*mult + end if + else + call SGEMV('N',m,n,mult,Mat,m,vec, 1,beta, Out,1) + end if + call Matrix_End('MulVecSingle') + + end subroutine Matrix_MulVecSingle + + + + + subroutine Matrix_MulVecSymm(Mat,vec,Out,a,b) + ! Out = a*Mat*vec + b*out + real(dm), intent(in) :: Mat(:,:) + real(dm) vec(:) + real(dm) Out(:) + real(dm), intent(in), optional :: a,b + real(dm) mult, beta + integer m,n + + call Matrix_Start('MulVecSymm') + + m = Size(Mat,DIM=1) + n = Size(Vec) + if (m /= n) call MpiStop('Matrix_MulVecSymm: size mismatch') + if (present(a)) then + mult = a + else + mult = ROne + end if + if (present(b)) then + beta = b + else + beta = RZero + end if + + if (matrix_method == Mat_F90) then + if (beta /= RZero) then + out = a*MatMul(Mat,Vec) + beta*Out + else + out = MatMul(Mat,Vec) + if (mult /= ROne) Out = Out*mult + end if + else +#ifdef MATRIX_SINGLE + call SSYMV('U',m,mult,Mat,m,vec, 1,beta, Out,1) +#else + call DSYMV('U',m,mult,Mat,m,vec, 1,beta, Out,1) +#endif + end if + call Matrix_End('MulVecSymm') + + end subroutine Matrix_MulVecSymm + + subroutine Matrix_MulVecSymmSingle(Mat,vec,Out,a,b) + ! Out = a*Mat*vec + b*out + real, intent(in) :: Mat(:,:) + real vec(:) + real Out(:) + real, intent(in), optional :: a,b + real mult, beta + integer m,n + + call Matrix_Start('MulVecSymm') + + m = Size(Mat,DIM=1) + n = Size(Vec) + if (m /= n) call MpiStop('Matrix_MulVecSymm: size mismatch') + if (present(a)) then + mult = a + else + mult = SOne + end if + if (present(b)) then + beta = b + else + beta = SZero + end if + + if (matrix_method == Mat_F90) then + if (beta /= RZero) then + out = a*MatMul(Mat,Vec) + beta*Out + else + out = MatMul(Mat,Vec) + if (mult /= ROne) Out = Out*mult + end if + else + call SSYMV('U',m,mult,Mat,m,vec, 1,beta, Out,1) + end if + call Matrix_End('MulVecSymmSingle') + + end subroutine Matrix_MulVecSymmSingle + + function Matrix_vecdotSingle(vec1,vec2) + real vec1(:),vec2(:) + real Matrix_vecdotSingle + integer n + real sdot + external sdot + + n=size(vec1) + if (n/=size(vec2)) call MpiStop('Matrix_vecdotSingle: size mismatch') + Matrix_vecdotSingle = sdot(n, vec1, 1, vec2, 1) + + end function Matrix_vecdotSingle + + + subroutine Matrix_InverseArrayMPI(Arr,nmat) + !Invert array of matrices by sending each to separate CPU + integer, intent(in) :: nmat +#ifdef __GFORTRAN__ + Type(TMatrixType), target :: Arr(:) +#else + Type(TMatrixType), target :: Arr(*) +#endif + Type(TMatrixType), pointer :: AM + integer n + integer i,MpiID, MpiSize + integer sz +#ifdef MPI + integer j, ierr, sid + Type(TMatrixType), target :: tmp +#endif + + call MpiStat(MpiID, MpiSize) + if (MpiId==0) then + n=nmat + sz = Size(Arr(1)%M,DIM=1) + end if + ! if (MpiID==0) then + ! do i=1,nmat + ! print *,'inverting',i + ! call Matrix_inverse(Arr(i)%M) + ! end do + ! end if + ! return +#ifdef MPI + if (MpiID==0) print *, 'MatrixInverseArray: starting' + call MPI_BCAST(n,1,MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(sz,1,MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if (MpiID/=0) then + allocate(tmp%M(sz,sz)) + AM => tmp + end if +#endif + + do i= 1,n + if (MpiID==0) AM => Arr(i) +#ifdef MPI + if (mod(i,MpiSize)/=MpiID) then + !Do nothing + if (MpiId==0) then + j=mod(i,MpiSize) + call MPI_SEND(AM%M,size(AM%M),MPI_DOUBLE_PRECISION, j, 1, MPI_COMM_WORLD, ierr) + end if + else + if (MpiId/=0) then + !Get from main thread + call MPI_RECV(AM%M,size(AM%M),MPI_DOUBLE_PRECISION, 0, 1, MPI_COMM_WORLD,MPI_STATUS_IGNORE, ierr) + end if +#endif + call Matrix_Inverse(AM%M) + +#ifdef MPI + if (MpiID==0) then + do j = max(1,i-MpiSize+1),i-1 + sid = mod(j,MpiSize) + call MPI_RECV(Arr(j)%M,size(Arr(j)%M),MPI_DOUBLE_PRECISION, sid, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + end do + else + call MPI_SEND(AM%M,size(AM%M),MPI_DOUBLE_PRECISION, 0, 1, MPI_COMM_WORLD, ierr) + end if + end if +#endif + end do + + +#ifdef MPI + if (MpiID==0) then + do j=n - mod(n,MpiSize) +1 ,n + sid= mod(j,MpiSize) + call MPI_RECV(ARr(j)%M,Size(ARr(j)%M),MPI_DOUBLE_PRECISION, sid, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + end do + else + deallocate(tmp%M) + end if +#endif + if (MpiID==0) print *, 'MatrixInverseArray: Done' + + + end subroutine Matrix_InverseArrayMPI + + + + end module MatrixUtils diff --git a/MiscUtils.f90 b/MiscUtils.f90 new file mode 100644 index 0000000..af276ae --- /dev/null +++ b/MiscUtils.f90 @@ -0,0 +1,205 @@ + module MiscUtils + implicit none + + INTERFACE PresentDefault + module procedure PresentDefault_S, PresentDefault_L, PresentDefault_I, PresentDefault_R, PresentDefault_D + END INTERFACE PresentDefault + + INTERFACE IfThenElse + module procedure IfThenElse_S, IfThenElse_L, IfThenElse_I, IfThenElse_R, IfThenElse_D + END INTERFACE IfThenElse + + INTERFACE IsFloat + module procedure IsFloat0, IsFloat1, IsFloat2 + END INTERFACE IsFloat + + contains + + function DefaultFalse(S) result(Sout) + logical, intent(in), optional :: S + logical :: Sout + + if (present(S)) then + SOut = S + else + SOut = .false. + end if + + end function DefaultFalse + + function DefaultTrue(S) result(Sout) + logical, intent(in), optional :: S + logical :: Sout + + if (present(S)) then + SOut = S + else + SOut = .true. + end if + + end function DefaultTrue + + + function PresentDefault_S(default, S) result(Sout) + character(LEN=*), intent(in), target :: default + character(LEN=*), intent(in), target, optional :: S + character(LEN=:), pointer :: Sout + + if (present(S)) then + SOut => S + else + SOut => default + end if + end function PresentDefault_S + + function PresentDefault_L(default, S) result(Sout) + logical, intent(in) :: default + logical, intent(in), optional :: S + logical :: Sout + + if (present(S)) then + SOut = S + else + SOut = default + end if + end function PresentDefault_L + + function PresentDefault_I(default, S) result(Sout) + integer, intent(in) :: default + integer, intent(in), optional :: S + integer :: Sout + + if (present(S)) then + SOut = S + else + SOut = default + end if + end function PresentDefault_I + + function PresentDefault_R(default, S) result(Sout) + real, intent(in) :: default + real, intent(in), optional :: S + real:: Sout + + if (present(S)) then + SOut = S + else + SOut = default + end if + end function PresentDefault_R + + function PresentDefault_D(default, S) result(Sout) + double precision, intent(in) :: default + double precision, intent(in), optional :: S + double precision:: Sout + + if (present(S)) then + SOut = S + else + SOut = default + end if + end function PresentDefault_D + + + function IfThenElse_S(flag, either, or) result(IfThenElse) + logical, intent(in) :: flag + character(LEN=:), pointer :: IfThenElse + character(LEN=*), target :: either, or + + if (flag) then + IfThenElse => either + else + IfThenElse => or + end if + + end function + + function IfThenElse_L(flag, either, or) result(IfThenElse) + logical, intent(in) :: flag + logical :: IfThenElse, either, or + + if (flag) then + IfThenElse = either + else + IfThenElse = or + end if + + end function + + function IfThenElse_I(flag, either, or) result(IfThenElse) + logical, intent(in) :: flag + integer :: IfThenElse, either, or + + if (flag) then + IfThenElse = either + else + IfThenElse = or + end if + + end function + + function IfThenElse_R(flag, either, or) result(IfThenElse) + logical, intent(in) :: flag + real :: IfThenElse, either, or + + if (flag) then + IfThenElse = either + else + IfThenElse = or + end if + + end function + + function IfThenElse_D(flag, either, or) result(IfThenElse) + logical, intent(in) :: flag + double precision :: IfThenElse, either, or + + if (flag) then + IfThenElse = either + else + IfThenElse = or + end if + + end function + + logical function isFloat0(R) + class(*), intent(in) :: R + + select type(R) + type is (real) + isFloat0 = .true. + type is (double precision) + isFloat0 = .true. + class default + isFloat0 = .false. + end select + end function isFloat0 + + logical function isFloat1(R) + class(*), intent(in) :: R(:) + + isFloat1 = isFloat0(R(LBOUND(R,1))) + + end function isFloat1 + + logical function isFloat2(R) + class(*), intent(in) :: R(:,:) + + isFloat2 = isFloat0(R(LBOUND(R,1),LBOUND(R,2))) + + end function isFloat2 + + function IndexOf(aval,arr, n) + integer, intent(in) :: n, arr(n), aval + integer :: IndexOf, i + + do i=1,n + if (arr(i)==aval) then + IndexOf= i + return + end if + end do + IndexOf = 0 + + end function IndexOf +end module MiscUtils diff --git a/MpiUtils.f90 b/MpiUtils.f90 new file mode 100644 index 0000000..a0f5d10 --- /dev/null +++ b/MpiUtils.f90 @@ -0,0 +1,172 @@ + module MpiUtils + implicit none + +#ifdef MPI + include "mpif.h" +#endif + + integer, parameter :: TTimer_dp = Kind(1.d0) + + Type TTimer + real(TTimer_dp) start_time + contains + procedure :: Start => TTimer_Start + procedure :: Time => TTimer_Time + procedure :: WriteTime => TTimer_WriteTime + end type TTimer + + contains + + function GetMpiRank() + integer GetMpiRank +#ifdef MPI + integer ierror + call mpi_comm_rank(mpi_comm_world,GetMPIrank,ierror) + if (ierror/=MPI_SUCCESS) call MpiStop('MPI fail') +#else + GetMpiRank=0 +#endif + + end function GetMpiRank + + function IsMainMPI() + logical IsMainMPI + + IsMainMPI = GetMpiRank() == 0 + + end function IsMainMPI + + subroutine MpiStop(Msg) + character(LEN=*), intent(in), optional :: Msg + integer i +#ifdef MPI + integer ierror, MpiRank +#endif + + if (present(Msg)) write(*,*) trim(Msg) + +#ifdef MPI + call mpi_comm_rank(mpi_comm_world,MPIrank,ierror) + write (*,*) 'MpiStop: ', MpiRank + call MPI_ABORT(MPI_COMM_WORLD,i) +#endif + i=1 !put breakpoint on this line to debug + stop + + end subroutine MpiStop + + subroutine MpiStat(MpiID, MpiSize) + implicit none + integer MpiID,MpiSize +#ifdef MPI + integer ierror + call mpi_comm_rank(mpi_comm_world,MpiID,ierror) + if (ierror/=MPI_SUCCESS) call MpiStop('MpiStat: MPI rank') + call mpi_comm_size(mpi_comm_world,MpiSize,ierror) +#else + MpiID=0 + MpiSize=1 +#endif + end subroutine MpiStat + + subroutine MpiQuietWait + !Set MPI thread to sleep, e.g. so can run openmp on cpu instead +#ifdef MPI + integer flag, ierr, STATUS(MPI_STATUS_SIZE) + integer i, MpiId, MpiSize + + call MpiStat(MpiID, MpiSize) + if (MpiID/=0) then + do + call MPI_IPROBE(0,0,MPI_COMM_WORLD,flag, MPI_STATUS_IGNORE,ierr) + if (flag/=0) then + call MPI_RECV(i,1,MPI_INTEGER, 0,0,MPI_COMM_WORLD,status,ierr) + exit + end if + call sleep(1) + end do + end if +#endif + end subroutine + + subroutine MpiWakeQuietWait +#ifdef MPI + integer j, MpiId, MpiSize, ierr,r + + call MpiStat(MpiID, MpiSize) + if (MpiID==0) then + do j=1, MpiSize-1 + call MPI_ISSEND(MpiId,1,MPI_INTEGER, j,0,MPI_COMM_WORLD,r,ierr) + end do + end if +#endif + end subroutine MpiWakeQuietWait + + subroutine MpiShareString(S, from) + character(LEN=:), allocatable :: S + integer from +#ifdef MPI + integer inlen, rank, ierror + + rank = GetMpiRank() + + if (rank==from) inlen=len(S) + + CALL MPI_Bcast(inlen, 1, MPI_INTEGER, from, MPI_COMM_WORLD, ierror) + if (ierror/=MPI_SUCCESS) call MpiStop('MpiShareString: fail') + + if (rank /= from ) allocate(character(inlen)::S) + CALL MPI_Bcast(S, LEN(S), MPI_CHARACTER, from, MPI_COMM_WORLD, ierror) +#endif + end subroutine MpiShareString + + + function TimerTime() + real(TTimer_dp) time + real(TTimer_dp) :: TimerTime +#ifdef MPI + TimerTime = MPI_WTime() +#else + call cpu_time(time) + TimerTime= time +#endif + end function TimerTime + + subroutine TTimer_Start(this) + class(TTimer) :: this + this%start_time = TimerTime() + end subroutine TTimer_Start + + real(TTimer_dp) function TTimer_Time(this) + class(TTimer) :: this + TTimer_Time = TimerTime() - this%start_time + end function TTimer_Time + + subroutine TTimer_WriteTime(this,Msg, start) + class(TTimer) :: this + character(LEN=*), intent(in), optional :: Msg + real(TTimer_dp), optional :: start + real(TTimer_dp) T, DeltaT + character(LEN=:), allocatable :: tmp + + if (present(start)) then + T=start + else + T=this%start_time + end if + + DeltaT = TimerTime() - T + if (present(Msg)) then + tmp = trim(Msg)//': ' + if (DeltaT > 0.00002 .and. DeltaT < 1000 .and. len_trim(tmp)<24) then + write (*,'(a25,f10.5)') tmp, DeltaT + else + write (*,*) trim(Msg)//': ', DeltaT + end if + end if + if (.not. present(start)) this%start_time = TimerTime() + + end subroutine TTimer_WriteTime + + + end module MpiUtils diff --git a/ObjectLists.f90 b/ObjectLists.f90 new file mode 100644 index 0000000..0b42ec8 --- /dev/null +++ b/ObjectLists.f90 @@ -0,0 +1,960 @@ + + module ObjectLists + !Implement lists of arbitrary objects + !AL Feb 2014 + use FileUtils + implicit none + + private + + type Object_pointer + class(*), pointer :: p => null() + class(*), pointer :: Object => null() + end type Object_pointer + + type Object_array_pointer + class(*), pointer :: p(:) => null() + end type Object_array_pointer + +#ifdef SINGLE + integer, parameter :: list_prec = Kind(1.0) +#else + integer, parameter :: list_prec = Kind(1.d0) +#endif + + Type, abstract :: TSaveLoadStateObject + contains + procedure :: SaveState + procedure :: LoadState + end Type TSaveLoadStateObject + + Type, extends(TSaveLoadStateObject) :: TObjectList + integer :: Count =0 + integer :: Delta = 32 + integer :: DeltaScale = 10 + !expanding expanding Items, expand array by Delta + Count/DeltaScale + integer :: Capacity = 0 + logical :: OwnsObjects = .true. + Type(Object_pointer), allocatable :: Items(:) + contains + procedure :: AddArray + procedure :: AddItem + procedure :: AddItemPointer + procedure :: AddCopy + procedure :: AssignPointers + procedure :: DeleteItem + procedure :: DeleteRange + procedure :: FreeItem + procedure :: SetCapacity + procedure :: ArrayItem + procedure :: ArrayItemIndex + procedure :: SaveBinary + procedure :: ReadBinary + procedure :: Thin + procedure :: Sort + procedure :: SortArr + procedure :: Swap + procedure :: Compare + procedure :: Clear + procedure :: DeltaSize + procedure :: QuickSort + procedure :: QuickSortArr + procedure :: RemoveDuplicates + procedure :: CheckIndex + procedure :: Error + procedure :: SaveState => TObjectList_SaveState + procedure :: LoadState => TObjectList_LoadState + procedure :: Object => TObjectList_Object !Secondary object associated with index + FINAL :: finalize + generic :: Add => AddItem, AddArray + end Type TObjectList + + Type, extends(TObjectList) :: TOwnedIntrinsicList + contains + procedure :: AddItem => TOwnedIntrinsicList_AddItem + procedure :: LoadState => TOwnedIntrinsicList_LoadState + procedure :: SaveState => TOwnedIntrinsicList_SaveState + end type + + Type, extends(TOwnedIntrinsicList):: TRealCompareList + contains + procedure :: Compare => CompareReal + end Type TRealCompareList + + Type, extends(TRealCompareList):: TRealList + contains + procedure :: TRealList_Item + procedure :: AddArrayItems => TRealList_AddArrayItems + procedure :: AsArray =>TRealList_AsArray + generic :: Item => TRealList_Item + !could add read and save state here + end Type TRealList + + Type, extends(TRealCompareList):: TRealArrayList + contains + procedure :: Value => TRealArrayList_Value + procedure :: RealArrItem => TRealArrayList_Item + generic :: Item => Value, RealArrItem + end Type TRealArrayList + + Type, extends(TObjectList):: TIntegerArrayList + contains + procedure :: Value => TIntegerArrayList_Value + procedure :: IntegerArrItem => TIntegerArrayList_Item + generic :: Item => Value, IntegerArrItem + end Type TIntegerArrayList + + + Type, extends(TOwnedIntrinsicList) :: TStringList + contains + procedure :: CharAt => TStringList_CharAt + procedure :: Compare => TStringList_Compare + procedure :: StringItem => TStringList_Item + procedure :: SetFromString => TStringList_SetFromString + procedure :: ReadColumnsGetArray => TStringList_ReadColumnsGetArray + procedure :: IndexOf => TStringList_IndexOf + procedure :: ValueOf => TStringList_ValueOf + procedure :: WriteItems + generic :: Item => StringItem + end Type TStringList + + + public list_prec, TSaveLoadStateObject, TObjectList, TRealArrayList, TRealList, TIntegerArrayList, TStringList + contains + + subroutine LoadState(this,F) + class(TSaveLoadStateObject) :: this + class(TFileStream) :: F + end subroutine LoadState + + subroutine SaveState(this,F) + class(TSaveLoadStateObject) :: this + class(TFileStream) :: F + end subroutine SaveState + + + subroutine Clear(this, itemsOnly) + Class(TObjectList) :: this + integer i + logical, intent(in), optional :: itemsOnly + logical eachItem + + if (allocated(this%Items)) then + eachItem = .true. + if (present(itemsOnly)) eachItem=.not. itemsOnly + if (eachItem) then + do i=1,this%count + call this%FreeItem(i) + end do + end if + deallocate (this%Items) + end if + this%Count = 0 + this%Capacity = 0 + + end subroutine Clear + + subroutine finalize(this) + Type(TObjectList) :: this + call this%Clear() + end subroutine finalize + + subroutine Error(this, msg) + Class(TObjectList) :: this + character(LEN=*), intent(in) :: msg + + write(*,*) msg + stop + + end subroutine Error + + subroutine AddItem(this, C, Object) + Class(TObjectList) :: this + class(*), intent(in), target :: C + class(*), intent(in), target, optional :: Object + class(*), pointer :: CP + class(*), pointer :: ObjectP + + !This all looks a bit unneccessary, just trying to avoid ifort bugs, c.f. + !http://software.intel.com/en-us/forums/topic/390944 + !This subroutine does *not* help directly, but derived types can use AddItemPointer which is OK. + + CP=> C + if (present(Object)) then + ObjectP => Object + else + nullify(ObjectP) + end if + call this%AddItemPointer(CP, ObjectP) + + end subroutine AddItem + + + subroutine AddItemPointer(this, C, Object) + Class(TObjectList) :: this + class(*), intent(in), pointer :: C + class(*), intent(in), pointer, optional :: Object + + if (this%Count == this%Capacity) call this%SetCapacity(this%Capacity + this%DeltaSize()) + this%Count = this%Count + 1 + this%Items(this%Count)%P=>C + if (present(Object)) this%Items(this%Count)%Object=> Object + + end subroutine AddItemPointer + + subroutine AddCopy(this, C, Object) + Class(TObjectList) :: this + class(*), intent(in) :: C + class(*), intent(in), optional :: Object + class(*), pointer :: P + class(*), pointer :: PO + + nullify(PO) + if (this%OwnsObjects) then + allocate(P, source=C) + if (present(Object)) allocate(PO, source=Object) + else + call this%Error('ObjectLists: Cannot add copy to un-owned list') + end if + call this%AddItemPointer(P, PO) + + end subroutine AddCopy + + subroutine AssignPointers(this, L2, ixmin, ixmax) + Class(TObjectList) :: this, L2 + integer, intent(in), optional :: ixmin, ixmax + integer i1,i2 + + call this%Clear() + i1=1 + i2=L2%Count + if (present(ixmin)) i1=ixmin + if (present(ixmax)) i2=ixmax + call this%SetCapacity(i2-i1+1) + this%Items = L2%Items(i1:i2) + this%Count = i2-i1+1 + this%OwnsObjects = .false. + + end subroutine AssignPointers + + integer function DeltaSize(this) + Class(TObjectList) :: this + + DeltaSize= this%Delta + this%Count/this%DeltaScale + + end function + + subroutine SetCapacity(this, C) + Class(TObjectList) :: this + integer C + Type(Object_pointer), dimension(:), allocatable :: TmpItems + + if (this%Count > 0) then + if (C < this%Count) call this%Error('ObjectLists: SetCapacity: smaller than Count') + allocate(TmpItems(C)) + TmpItems(:this%Count) = this%Items(:this%Count) + call move_alloc(TmpItems, this%Items) + else + allocate(this%Items(C)) + end if + this%Capacity = C + end subroutine SetCapacity + + subroutine FreeItem(this, i) + Class(TObjectList) :: this + integer, intent(in) :: i + logical want_Dealloc + + if (associated(this%Items(i)%P)) then + want_Dealloc = this%OwnsObjects + select type (point => this%Items(i)%P) + class is (object_array_pointer) + if (this%OwnsObjects .and. associated(Point%P)) deallocate(Point%P) + want_Dealloc = .true. + !type is (real(kind=list_prec)) + ! want_Dealloc = .false. + end select + if (want_Dealloc) deallocate(this%Items(i)%P) + this%Items(i)%P=> null() + end if + + if (associated(this%Items(i)%Object) .and. this%OwnsObjects) deallocate(this%Items(i)%Object) + this%Items(i)%Object=> null() + + end subroutine FreeItem + + subroutine DeleteItem(this, i) + Class(TObjectList) :: this + integer, intent(in) :: i + + call this%FreeItem(i) + if (this%Count > 1) this%Items(i:this%Count-1) = this%Items(i+1:this%Count) + this%Items(this%Count)%P => null() + this%Items(this%Count)%Object => null() + this%Count = this%Count -1 + + end subroutine DeleteItem + + subroutine DeleteRange(this, i1,i2) + Class(TObjectList) :: this + integer, intent(in) :: i1,i2 + integer i, dN + + do i=i1,i2 + call this%FreeItem(i) + end do + dN= i2-i1+1 + if (i2 null() + this%Items(i)%Object => null() + end do + this%Count = this%Count - dN + + end subroutine DeleteRange + + subroutine AddArray(this, P) + Class (TObjectList) :: this + class(*), target, intent(in) :: P(:) + class(*), pointer :: Pt + + allocate(object_array_pointer::Pt) + call this%AddItemPointer(Pt) + select type (Pt) + class is (object_array_pointer) + if (this%ownsObjects) then + allocate(Pt%P, source= P) + else + Pt%P => P + end if + end select + end subroutine AddArray + + subroutine CheckIndex(this,i) + Class(TObjectList) :: this + integer, intent(in) :: i + + if (i>this%Count .or. i<1) call this%Error('Item out of range') + + end subroutine CheckIndex + + !why this crashes in ifort 13 I do not know.. + !subroutine AddArray(this, P) + !Class (TObjectList) :: this + !class(*), target, intent(in) :: P(:) + !Type(object_array_pointer), pointer :: Pt + ! + !allocate(Pt) + !call this%AddItem(Pt) + !if (this%ownsObjects) then + ! allocate(Pt%P(1:SIZE(P)), source= P) + !else + ! Pt%P => P + !end if + !end subroutine AddArray + + function ArrayItem(this, i) result(P) + Class(TObjectList) :: this + integer, intent(in) :: i + Class(*), pointer :: P(:) + + call this%CheckIndex(i) + select type (Point=> this%Items(i)%P) + class is (object_array_pointer) + P => Point%P + class default + call this%Error('ObjectLists: item is not array item') + end select + + end function ArrayItem + + function ArrayItemIndex(this, i, j) result(P) + Class(TObjectList) :: this + integer, intent(in) :: i, j + Class(*), pointer :: P + Class(*), pointer :: Arr(:) + + Arr => this%ArrayItem(i) + P => Arr(j) + + end function ArrayItemIndex + + function TObjectList_Object(this, i) + class(TObjectList) :: this + integer, intent(in) :: i + class(*), pointer :: TObjectList_Object + + call this%CheckIndex(i) + TObjectList_Object => this%Items(i)%Object + + end function TObjectList_Object + + subroutine SaveBinary(this,fid) + Class(TObjectList) :: this + integer, intent(in) :: fid + integer i,k + class(*), pointer :: P(:) + + write (fid) this%Count + do i=1,this%Count + select type (Item=> this%Items(i)%P) + class is (object_array_pointer) + P => this%ArrayItem(i) + select type (Point=> P) + Type is (real) + k=1 + write(fid) size(P),k + write(fid) Point + Type is (double precision) + k=2 + write(fid) size(P),k + write(fid) Point + Type is (integer) + k=3 + write(fid) size(P),k + write(fid) Point + Type is (logical) + k=4 + write(fid) size(P),k + write(fid) Point + class default + call this%Error('TObjectList: Unknown type to save') + end select + Type is (character(LEN=*)) + k=5 + write(fid) len(Item), k + write(fid) Item + class default + call this%Error('TObjectList: not implemented non-array save') + end select + end do + + end subroutine SaveBinary + + subroutine ReadBinary(this,fid) + Class(TObjectList) :: this + integer, intent(in) :: fid + integer num,i,sz, k + real, pointer :: ArrR(:) + double precision, pointer :: ArrD(:) + integer, pointer :: ArrI(:) + logical, pointer :: ArrL(:) + class(*), pointer :: St + + call this%Clear() + this%OwnsObjects = .false. + read (fid) num + call this%SetCapacity(num) + do i=1,num + read(fid) sz, k + if (k==1) then + allocate(ArrR(sz)) + read(fid) ArrR + call this%AddArray(ArrR) + else if (k==2) then + allocate(ArrD(sz)) + read(fid) ArrD + call this%AddArray(ArrD) + else if (k==3) then + allocate(ArrI(sz)) + read(fid) ArrI + call this%AddArray(ArrI) + else if (k==4) then + allocate(ArrL(sz)) + read(fid) ArrL + call this%AddArray(ArrL) + else if (k==5) then + allocate(character(sz)::St) !Ifort required class(*) pointer + select type (St) + type is (character(LEN=*)) + read(fid) St + end select + call this%AddItemPointer(St) + else + call this%Error('TObjectList ReadBinary - unknown object type') + end if + end do + this%Count = num + this%OwnsObjects = .true. + + end subroutine ReadBinary + + subroutine Thin(this, i) + Class(TObjectList):: this + integer, intent(in) :: i + integer newCount, j + Type(Object_pointer), dimension(:), pointer :: TmpItems + + if (this%Count > 1) then + newCount = (this%Count-1)/i+1 + allocate(TmpItems(newCount)) + TmpItems= this%Items(1:this%Count:i) + if (this%OwnsObjects) then + do j=1,this%count + if (mod(j-1,i)/=0) call this%FreeItem(j) + end do + end if + deallocate(this%Items) + this%Capacity = newCount + allocate(this%Items, source = TmpItems) + this%Count = newCount + deallocate(TmpItems) + end if + + end subroutine Thin + + subroutine Swap(this, i, j) + Class(TObjectList) :: this + integer, intent(in) :: i, j + type(Object_pointer) :: temp + + temp = this%Items(i) + this%Items(i) = this%Items(j) + this%Items(j) = temp + + end subroutine Swap + + recursive subroutine QuickSortArr(this, Lin, R, index) + !Sorts an array of pointers by the value of the index'th entry + Class(TObjectList) :: this + integer, intent(in) :: Lin, R, index + integer I, J, L + class(*), pointer :: P + + L = Lin + do + I = L + J = R + P => this%ArrayItemIndex((L + R)/2, index) + do + do while (this%Compare(this%ArrayItemIndex(I, Index),P) < 0) + I = I + 1 + end do + + do while (this%Compare(this%ArrayItemIndex(J,Index), P) > 0) + J = J - 1 + end do + + if (I <= J) then + call this%Swap(I,J) + I = I + 1 + J = J - 1 + end if + if (I > J) exit + + end do + if (L < J) call this%QuickSortArr(L, J, index) + L = I + if (I >= R) exit + end do + + end subroutine QuickSortArr + + subroutine SortArr(this, index) + Class(TObjectList) :: this + integer, intent(in) :: index + + if (this%Count>1) call this%QuickSortArr(1, this%Count, index) + + end subroutine SortArr + + + recursive subroutine QuickSort(this, Lin, R) + Class(TObjectList) :: this + integer, intent(in) :: Lin, R + integer I, J, L + class(*), pointer :: P + + L= Lin + do + I = L + J = R + P => this%Items((L + R)/2)%P + do + do while (this%Compare(this%Items(I)%P,P) < 0) + I = I + 1 + end do + + do while (this%Compare(this%Items(J)%P, P) > 0) + J = J - 1 + end do + + if (I <= J) then + call this%Swap(I,J) + I = I + 1 + J = J - 1 + end if + if (I > J) exit + end do + if (L < J) call this%QuickSort(L, J) + L = I + if (I >= R) exit + end do + + end subroutine QuickSort + + subroutine Sort(this) + Class(TObjectList) :: this + + if (this%Count>1) call this%QuickSort(1, this%Count) + + end subroutine Sort + + integer function Compare(this, R1, R2) result(comp) + Class(TObjectList) :: this + class(*) R1,R2 + + comp=0 !equality + call this%Error('TObjectList: Compare must be defined for derived type') + + end function Compare + + subroutine RemoveDuplicates(this) + Class(TObjectList) :: this + integer i + + do i=this%Count-1, 1, -1 + if (this%Compare(this%Items(i+1)%P, this%Items(i)%P)==0) call this%DeleteItem(i+1) + end do + + end subroutine RemoveDuplicates + + subroutine TObjectList_SaveState(this,F) + class(TObjectList) :: this + class(TFileStream) :: F + integer i + + call F%Write(this%Count) + do i=1,this%Count + select type (item => this%Items(i)%P) + class is (TSaveLoadStateObject) + call item%SaveState(F) + class default + call this%Error('TObjectList_SaveState: List contains non-TSaveLoadStateObject item') + end select + end do + + end subroutine TObjectList_SaveState + + subroutine TObjectList_LoadState(this,F) + class(TObjectList) :: this + class(TFileStream) :: F + integer i, count + + if (.not. F%ReadItem(count) .or. count/=this%Count) & + & call this%Error('TObjectList_LoadState count mismatch (objects must exist before load)') + do i=1,this%Count + select type (item => this%Items(i)%P) + class is (TSaveLoadStateObject) + call item%LoadState(F) + class default + call this%Error('List contains non-TSaveLoadStateObject item') + end select + end do + + end subroutine TObjectList_LoadState + + + !TOwnedIntrinsicList + + subroutine TOwnedIntrinsicList_LoadState(this,F) + class(TOwnedIntrinsicList) :: this + class(TFileStream) :: F + call this%ReadBinary(F%unit) + end subroutine TOwnedIntrinsicList_LoadState + + subroutine TOwnedIntrinsicList_SaveState(this,F) + class(TOwnedIntrinsicList) :: this + class(TFileStream) :: F + call this%SaveBinary(F%unit) + end subroutine TOwnedIntrinsicList_SaveState + + + subroutine TOwnedIntrinsicList_AddItem(this, C, Object) + Class(TOwnedIntrinsicList) :: this + class(*), intent(in), target :: C + class(*), intent(in), target, optional :: Object + + call this%TObjectList%AddCopy(C, Object) + + end subroutine TOwnedIntrinsicList_AddItem + + !TRealCompareList + integer function CompareReal(this, R1, R2) result(comp) + Class(TRealCompareList) :: this + class(*) R1,R2 + real(list_prec) R + + select type (RR1 => R1) + type is (real(list_prec)) + select type (RR2 => R2) + type is (real(list_prec)) + R = RR1-RR2 + if (R< 0) then + comp =-1 + elseif (R>0) then + comp = 1 + else + comp = 0 + end if + return + end select + class default + call this%Error('TRealList: Compare not defined for this type') + end select + + end function CompareReal + + + !TRealList: List of reals + function TRealList_Item(this,i) result(R) + Class(TRealList) :: this + integer, intent(in) :: i + real(list_prec) R + + call this%CheckIndex(i) + select type (pt=>this%Items(i)%P) + type is (real(kind=list_prec)) + R = pt + class default + call this%Error('TRealList: object of wrong type') + end select + + end function TRealList_Item + + subroutine TRealList_AddArrayItems(this, A) + Class(TRealList) :: this + real(kind=list_prec), intent(in) :: A(:) + integer i + + do i=1, size(A) + call this%AddItem(A(i)) + end do + + end subroutine TRealList_AddArrayItems + + function TRealList_AsArray(this) result(A) + Class(TRealList) :: this + real(kind=list_prec):: A(this%Count) + integer i + + do i=1, size(A) + A(i) = this%Item(i) + end do + + end function TRealList_AsArray + + !TRealArrayList: List of arrays of reals + + function TRealArrayList_Item(this, i) result(P) + Class(TRealArrayList) :: this + integer, intent(in) :: i + real(list_prec), pointer :: P(:) + class(*), pointer :: Item(:) + + Item => this%ArrayItem(i) + select type (pt=>Item) + type is (real(kind=list_prec)) + P=> pt + class default + call this%Error('TRealArrayList: object of wrong type') + end select + + end function TRealArrayList_Item + + function TRealArrayList_Value(this, i, j) result(P) + Class(TRealArrayList) :: this + integer, intent(in) :: i, j + real(list_prec) :: P + class(*), pointer :: C + + C => this%ArrayItemIndex(i,j) + select type (Arr=> C) + Type is (real(list_prec)) + P = Arr + end select + + end function TRealArrayList_Value + + !TIntegerArrayList: List of arrays of reals + + function TIntegerArrayList_Item(this, i) result(P) + Class(TIntegerArrayList) :: this + integer, intent(in) :: i + integer, pointer :: P(:) + class(*), pointer :: Item(:) + + Item => this%ArrayItem(i) + select type (pt=>Item) + type is (integer) + P=> pt + class default + call this%Error('TIntegerArrayList: object of wrong type') + end select + + end function TIntegerArrayList_Item + + function TIntegerArrayList_Value(this, i, j) result(P) + Class(TIntegerArrayList) :: this + integer, intent(in) :: i, j + Integer :: P + class(*), pointer :: C + + C => this%ArrayItemIndex(i,j) + select type (Arr=> C) + Type is (integer) + P = Arr + end select + + end function TIntegerArrayList_Value + + !!! TStringList + + function TStringList_Item(this,i) result(S) + Class(TStringList) :: this + integer, intent(in) :: i + character(LEN=:), pointer :: S + + call this%CheckIndex(i) + select type (pt=>this%Items(i)%P) + type is (character(LEN=*)) + S => pt + class default + call this%Error('TStringList: object of wrong type') + end select + + end function TStringList_Item + + function TStringList_ValueOf(this, key) + class(TStringList) :: this + character(LEN=*), intent(in) :: key + integer :: i + character(LEN=:), pointer :: TStringList_ValueOf + + i = this%IndexOf(key) + if (i==-1) call this%Error('TStringList_ValueOf key not found:'//key) + select type (value=>this%Items(i)%Object) + type is (character(LEN=*)) + TStringList_ValueOf=> value + class default + call this%Error('TStringList_ValueOf Object is not a string') + end select + + end function TStringList_ValueOf + + + subroutine TStringList_SetFromString(this, S, valid_chars_in) + class(TStringList) :: this + character(Len=*), intent(in) :: S + character(Len=*), intent(in), optional :: valid_chars_in + character(LEN=:), allocatable :: item + integer i,j + character(LEN=:), allocatable :: valid_chars + + if (present(valid_chars_in)) then + valid_chars = trim(valid_chars_in) + else + valid_chars='abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_-.' + endif + + call this%Clear() + allocate(item, source=S) + j=0 + do i=1, len_trim(S) + if (verify(S(i:i),valid_chars) == 0) then + j=j+1 + item(j:j) = S(i:i) + else + if (trim(S(i:i))/='') then + write (*,*) 'Invalid character in: '//trim(S) + end if + if (j>0) call this%Add(item(1:j)) + j=0 + end if + end do + if (j>0) call this%Add(item(1:j)) + + end subroutine TStringList_SetFromString + + + subroutine TStringList_ReadColumnsGetArray(this, filename, array) + class(TStringList) :: this + character(LEN=*), intent(in) :: filename + real(list_prec), intent(out), allocatable :: array(:,:) + character(LEN=:), allocatable :: comment + + call File%LoadTxt(filename, array,comment = comment) + call this%SetFromString(comment) + if (this%Count /= size(array,1)) & + call this%Error('Column header does not match number of columns: ' //trim(filename)) + + end subroutine TStringList_ReadColumnsGetArray + + function TStringList_IndexOf(this, S) result(index) + class(TStringList) :: this + character(LEN=*), intent(in) :: S + integer index, i + + do i=1,this%Count + if (this%Item(i)==S) then + index = i + return + end if + end do + index=-1 + + end function TStringList_IndexOf + + function TStringList_CharAt(this, i, j) result(C) + Class(TStringList) :: this + integer, intent(in) :: i, j + character :: C + character(LEN=:), pointer :: P + + call this%CheckIndex(i) + P => this%Item(i) + C = P(j:j) + + end function TStringList_CharAt + + integer function TStringList_Compare(this, R1, R2) result(comp) + Class(TStringList) :: this + class(*) R1,R2 + + select type (RR1 => R1) + type is (character(LEN=*)) + select type (RR2 => R2) + type is (character(LEN=*)) + if (RR1 < RR2) then + comp =-1 + elseif (RR1>RR2) then + comp = 1 + else + comp = 0 + end if + return + end select + class default + call this%Error('TStringList_Compare: not defined for this type') + end select + + end function TStringList_Compare + + subroutine WriteItems(this, unit) + use, intrinsic :: iso_fortran_env, only : output_unit + Class(TStringList) :: this + integer, optional, intent(in) :: unit + integer i, aunit + + if (present(unit)) then + aunit = unit + else + aunit = output_unit + end if + do i=1, this%Count + write(aunit,*) this%Item(i) + end do + + end subroutine WriteItems + + end module ObjectLists diff --git a/RandUtils.f90 b/RandUtils.f90 new file mode 100644 index 0000000..f10544d --- /dev/null +++ b/RandUtils.f90 @@ -0,0 +1,377 @@ + + module RandUtils + use MpiUtils + implicit none + + integer :: rand_inst = 0 + integer, parameter :: krand = KIND(1.d0) + integer :: Rand_Feedback = 1 + + interface RandRotation + module procedure RandRotationS, RandRotationD + end interface + + contains + + + !subroutine init_random_seed() + !integer, allocatable :: seed(:) + !integer :: i, n, un, istat, dt(8), pid, t(2), s + !integer(8) :: count, tms + ! + !call random_seed(size = n) + !allocate(seed(n)) + !! First try if the OS provides a random number generator + !! Fallback to XOR:ing the current time and pid. The PID is + !! useful in case one launches multiple instances of the same + !! program in parallel. + !call system_clock(count) + !if (count /= 0) then + ! t = transfer(count, t) + !else + ! call date_and_time(values=dt) + ! tms = (dt(1) - 1970) * 365_8 * 24 * 60 * 60 * 1000 & + ! + dt(2) * 31_8 * 24 * 60 * 60 * 1000 & + ! + dt(3) * 24 * 60 * 60 * 60 * 1000 & + ! + dt(5) * 60 * 60 * 1000 & + ! + dt(6) * 60 * 1000 + dt(7) * 1000 & + ! + dt(8) + ! t = transfer(tms, t) + !end if + !s = ieor(t(1), t(2)) + !pid = getpid() + 1099279 ! Add a prime + !s = ieor(s, pid) + !if (n >= 3) then + ! seed(1) = t(1) + 36269 + ! seed(2) = t(2) + 72551 + ! seed(3) = pid + ! if (n > 3) then + ! seed(4:) = s + 37 * (/ (i, i = 0, n - 4) /) + ! end if + !else + ! seed = s + 37 * (/ (i, i = 0, n - 1 ) /) + !end if + !call random_seed(put=seed) + !end subroutine init_random_seed + + + + subroutine initRandom(i, i2) + implicit none + integer, optional, intent(in) :: i + integer, optional, intent(in) :: i2 + integer seed_in,kl,ij + character(len=10) :: fred + real(krand) :: klr + + if (present(i)) then + seed_in = i + else + seed_in = -1 + end if + if (seed_in /=-1) then + if (present(i2)) then + kl=i2 + if (i2 > 30081) call MpiStop('initRandom:second seed too large') + else + kl = 9373 + end if + ij = i + else + call system_clock(count=ij) + ij = mod(ij + rand_inst*100, 31328) + call date_and_time(time=fred) + read (fred,'(e10.3)') klr + kl = mod(int(klr*1000), 30081) + end if + + if (Rand_Feedback > 0 ) write(*,'(" Random seeds:",1I6,",",1I6," rand_inst:",1I4)') ij,kl,rand_inst + call rmarin(ij,kl) + + end subroutine initRandom + + subroutine RandIndices(indices, nmax, n) + integer, intent(in) :: nmax, n + integer indices(n),i, ix + integer tmp(nmax) + + if (n> nmax) call MpiStop('Error in RandIndices, n > nmax') + do i=1, nmax + tmp(i)=i + end do + do i=1, n + ix = int(ranmar()*(nmax +1 -i)) + 1 + indices(i) = tmp(ix) + tmp(ix) = tmp(nmax+1-i) + end do + + end subroutine RandIndices + + subroutine RandRotationS(R, N) + !this is most certainly not the world's most efficient or robust random rotation generator + integer, intent(in) :: N + real R(N,N), vec(N), norm + integer i,j + + do j = 1, N + do + do i = 1, N + vec(i) = real(Gaussian1()) + end do + do i = 1, j-1 + vec = vec - sum(vec*R(i,:))*R(i,:) + end do + norm = sum(vec**2) + if (norm > 1e-3) exit + end do + R(j,:) = vec / sqrt(norm) + end do + + end subroutine RandRotationS + + + subroutine RandRotationD(R, N) + !this is most certainly not the world's most efficient or robust random rotation generator + integer, intent(in) :: N + double precision R(N,N), vec(N), norm + integer i,j + + do j = 1, N + do + do i = 1, N + vec(i) = Gaussian1() + end do + do i = 1, j-1 + vec = vec - sum(vec*R(i,:))*R(i,:) + end do + norm = sum(vec**2) + if (norm > 1e-3) exit + end do + R(j,:) = vec / sqrt(norm) + end do + + end subroutine RandRotationD + + + double precision function GAUSSIAN1() + implicit none + double precision R, V1, V2, FAC + integer, save :: iset = 0 + double precision, save :: gset + + !Box muller + if (ISET==0) then + R=2 + do while (R >= 1.d0) + V1=2.d0*ranmar()-1.d0 + V2=2.d0*ranmar()-1.d0 + R=V1**2+V2**2 + end do + FAC=sqrt(-2.d0*log(R)/R) + GSET=V1*FAC + GAUSSIAN1=V2*FAC + ISET=1 + else + GAUSSIAN1=GSET + ISET=0 + endif + end function GAUSSIAN1 + + + double precision function CAUCHY1() + implicit none + + Cauchy1 = Gaussian1()/max(1d-15,abs(Gaussian1())) + + end function CAUCHY1 + + + real function RANDEXP1() + ! + ! Random-number generator for the exponential distribution + ! Algorithm EA from J. H. Ahrens and U. Dieter, + ! Communications of the ACM, 31 (1988) 1330--1337. + ! Coded by K. G. Hamilton, December 1996, with corrections. + ! + real u, up, g, y + + real, parameter :: alog2= 0.6931471805599453 + real, parameter :: a = 5.7133631526454228 + real, parameter :: b = 3.4142135623730950 + real, parameter :: c = -1.6734053240284925 + real, parameter :: p = 0.9802581434685472 + real, parameter :: aa = 5.6005707569738080 + real, parameter :: bb = 3.3468106480569850 + real, parameter :: hh = 0.0026106723602095 + real, parameter :: dd = 0.0857864376269050 + + u = ranmar() + do while (u<=0) ! Comment out this block + u = ranmar() ! if your RNG can never + enddo ! return exact zero + g = c + u = u+u + do while (u<1.0) + g = g + alog2 + u = u+u + enddo + u = u-1.0 + if (u<=p) then + randexp1 = g + aa/(bb-u) + return + endif + do + u = ranmar() + y = a/(b-u) + up = ranmar() + if ((up*hh+dd)*(b-u)**2 <= exp(-(y+c))) then + randexp1 = g+y + return + endif + enddo + + end function RANDEXP1 + + + ! This random number generator originally appeared in ''Toward a Universal + ! Random Number Generator'' by George Marsaglia and Arif Zaman. + ! Florida State University Report: FSU-SCRI-87-50 (1987) + ! + ! It was later modified by F. James and published in ''A Review of Pseudo- + ! random Number Generators'' + ! + ! THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE. + ! (However, a newly discovered technique can yield + ! a period of 10^600. But that is still in the development stage.) + ! + ! It passes ALL of the tests for random number generators and has a period + ! of 2^144, is completely portable (gives bit identical results on all + ! machines with at least 24-bit mantissas in the floating point + ! representation). + ! + ! The algorithm is a combination of a Fibonacci sequence (with lags of 97 + ! and 33, and operation "subtraction plus one, modulo one") and an + ! "arithmetic sequence" (using subtraction). + ! + ! On a Vax 11/780, this random number generator can produce a number in + ! 13 microseconds. + !======================================================================== + ! + ! PROGRAM TstRAN + ! INTEGER IJ, KL, I + ! Thee are the seeds needed to produce the test case results + ! IJ = 1802 + ! KL = 9373 + ! + ! + ! Do the initialization + ! call rmarin(ij,kl) + ! + ! Generate 20000 random numbers + ! do 10 I = 1, 20000 + ! x = RANMAR() + !10 continue + ! + ! If the random number generator is working properly, the next six random + ! numbers should be: + ! 6533892.0 14220222.0 7275067.0 + ! 6172232.0 8354498.0 10633180.0 + ! + ! + ! + ! write(6,20) (4096.0*4096.0*RANMAR(), I=1,6) + !20 format (3f12.1) + ! end + ! + subroutine RMARIN(IJ,KL) + ! This is the initialization routine for the random number generator RANMAR() + ! NOTE: The seed variables can have values between: 0 <= IJ <= 31328 + ! 0 <= KL <= 30081 + !The random number sequences created by these two seeds are of sufficient + ! length to complete an entire calculation with. For example, if sveral + ! different groups are working on different parts of the same calculation, + ! each group could be assigned its own IJ seed. This would leave each group + ! with 30000 choices for the second seed. That is to say, this random + ! number generator can create 900 million different subsequences -- with + ! each subsequence having a length of approximately 10^30. + ! + ! Use IJ = 1802 & KL = 9373 to test the random number generator. The + ! subroutine RANMAR should be used to generate 20000 random numbers. + ! Then display the next six random numbers generated multiplied by 4096*4096 + ! If the random number generator is working properly, the random numbers + ! should be: + ! 6533892.0 14220222.0 7275067.0 + ! 6172232.0 8354498.0 10633180.0 + double precision U(97), C, CD, CM, S, T + integer I97, J97,i,j,k,l,m + integer ij,kl + integer ii,jj + + + ! INTEGER IRM(103) + + common /RASET1/ U, C, CD, CM, I97, J97 + if( IJ < 0 .or. IJ > 31328 .or. & + KL < 0 .or. KL > 30081 ) then + print '(A)', ' The first random number seed must have a value between 0 and 31328' + print '(A)',' The second seed must have a value between 0 and 30081' + stop + endif + I = mod(IJ/177, 177) + 2 + J = mod(IJ , 177) + 2 + K = mod(KL/169, 178) + 1 + L = mod(KL, 169) + do II = 1, 97 + S = 0.0 + T = 0.5 + do JJ = 1, 24 + M = mod(mod(I*J, 179)*K, 179) + I = J + J = K + K = M + L = mod(53*L+1, 169) + if (mod(L*M, 64) >= 32) then + S = S + T + endif + T = 0.5 * T +3 continue + end do + U(II) = S +2 continue + end do + C = 362436.0 / 16777216.0 + CD = 7654321.0 / 16777216.0 + CM = 16777213.0 /16777216.0 + I97 = 97 + J97 = 33 + + end subroutine RMARIN + + double precision function RANMAR() + ! This is the random number generator proposed by George Marsaglia in + ! Florida State University Report: FSU-SCRI-87-50 + ! It was slightly modified by F. James to produce an array of pseudorandom + ! numbers. + double precision U(97), C, CD, CM + integer I97, J97 + double precision uni + + common /RASET1/ U, C, CD, CM, I97, J97 + ! INTEGER IVEC + UNI = U(I97) - U(J97) + if( UNI < 0.0 ) UNI = UNI + 1.0 + U(I97) = UNI + I97 = I97 - 1 + if(I97 == 0) I97 = 97 + J97 = J97 - 1 + if(J97 == 0) J97 = 97 + C = C - CD + if( C < 0.d0 ) C = C + CM + UNI = UNI - C + if( UNI < 0.d0 ) UNI = UNI + 1.0 ! bug? + RANMAR = UNI + + end function RANMAR + + + end module RandUtils diff --git a/RangeUtils.f90 b/RangeUtils.f90 new file mode 100644 index 0000000..adac6b8 --- /dev/null +++ b/RangeUtils.f90 @@ -0,0 +1,459 @@ + !A collection of ranges, consisting of sections of minimum step size + !Useful for getting piecewise uniform point samplings for numerical integration + !(with possible mix of log and linear spacing) + !Antony Lewis, http://cosmologist.info/ + + module RangeUtils + use MiscUtils + implicit none + private + + type, public :: TRange + integer :: start_index + integer :: steps + logical :: IsLog + double precision :: Low, High + double precision :: delta + double precision :: delta_max, delta_min !for log spacing, the non-log max and min step size + contains + + end type TRange + + type, public :: TRanges + integer :: count = 0 + integer :: npoints = 0 + double precision :: Lowest, Highest + type(TRange), allocatable :: R(:) + logical :: has_dpoints = .false. + double precision, dimension(:), allocatable :: points, dpoints + !dpoints is (points(i+1)-points(i-1))/2 + double precision :: RangeTol = 0.1d0 + !fraction of bin width we are prepared for merged bin widths to increase by + contains + procedure :: Init => TRanges_Free + procedure :: Free => TRanges_Free + procedure :: IndexOf => TRanges_IndexOf + procedure :: GetArray => TRanges_GetArray + procedure :: Getdpoints => TRanges_Getdpoints + procedure :: Add_delta => TRanges_Add_delta + procedure :: Add => TRanges_Add + procedure :: Write => TRanges_Write + end type TRanges + + contains + +#ifdef __GFORTRAN__ + impure elemental subroutine TRanges_Free(this) +#else + subroutine TRanges_Free(this) +#endif + class(TRanges), intent(inout) :: this + + if (allocated(this%R)) deallocate(this%R) + if (allocated(this%points)) deallocate(this%points) + if (allocated(this%dpoints)) deallocate(this%dpoints) + this%count = 0 + this%npoints = 0 + this%has_dpoints = .false. + end subroutine TRanges_Free + + + function TRanges_IndexOf(this, tau) result(pointstep) + class(TRanges), intent(in) :: this + double precision, intent(in) :: tau + integer :: pointstep, i + + pointstep=0 + do i=1, this%count + associate(AReg => this%R(i)) + if (tau < AReg%High .and. tau >= AReg%Low) then + if (AReg%IsLog) then + pointstep = AReg%start_index + int(log(tau / AReg%Low) / AReg%delta) + else + pointstep = AReg%start_index + int((tau - AReg%Low) / AReg%delta) + end if + return + end if + end associate + end do + + if (tau >= this%Highest) then + pointstep = this%npoints + else + stop 'TRanges_IndexOf: value out of range' + end if + end function TRanges_IndexOf + + + subroutine TRanges_GetArray(this, want_dpoints) + class(TRanges), intent(inout) :: this + logical, intent(in), optional :: want_dpoints + integer :: i,j,ix + + this%has_dpoints = DefaultTrue(want_dpoints) + + ! Dealloc/Realloc only, when not enough space present. + if (allocated(this%points) .and. size(this%points) < this%npoints) & + deallocate(this%points) + if (.not. allocated(this%points)) allocate(this%points(this%npoints)) + + ix=0 + do i=1, this%count + associate (AReg => this%R(i)) + do j = 0, AReg%steps-1 + ix=ix+1 + if (AReg%IsLog) then + this%points(ix) = AReg%Low*exp(j*AReg%delta) + else + this%points(ix) = AReg%Low + AReg%delta*j + end if + end do + end associate + end do + ix =ix+1 + this%points(ix) = this%Highest + if (ix /= this%npoints) stop 'TRanges_GetArray: ERROR' + + if (this%has_dpoints) call this%Getdpoints() + end subroutine TRanges_GetArray + + + subroutine TRanges_Getdpoints(this, half_ends) + class(TRanges), intent(inout) :: this + logical, intent(in), optional :: half_ends + integer :: i + logical :: halfs + + halfs = DefaultTrue(half_ends) + + ! Dealloc/Realloc only, when not enough space present. + if (allocated(this%dpoints) .and. size(this%dpoints) < this%npoints) & + deallocate(this%dpoints) + if (.not. allocated(this%dpoints)) allocate(this%dpoints(this%npoints)) + + do i=2, this%npoints-1 + this%dpoints(i) = (this%points(i+1) - this%points(i-1))/2 + end do + if (halfs) then + this%dpoints(1) = (this%points(2) - this%points(1))/2 + this%dpoints(this%npoints) = (this%points(this%npoints) - this%points(this%npoints-1))/2 + else + this%dpoints(1) = (this%points(2) - this%points(1)) + this%dpoints(this%npoints) = (this%points(this%npoints) - this%points(this%npoints-1)) + end if + end subroutine TRanges_Getdpoints + + + subroutine TRanges_Add_delta(this, t_start, t_end, t_approx_delta, IsLog) + class(TRanges), intent(inout) :: this + logical, intent(in), optional :: IsLog + double precision, intent(in) :: t_start, t_end, t_approx_delta + integer :: n + logical :: WantLog + + WantLog = DefaultFalse(IsLog) + + if (t_end <= t_start) & + stop 'TRanges_Add_delta: end must be larger than start' + if (t_approx_delta <=0) stop 'TRanges_Add_delta: delta must be > 0' + + if (WantLog) then + n = max(1,int(log(t_end/t_start)/t_approx_delta + 1.d0 - this%RangeTol)) + else + n = max(1,int((t_end-t_start)/t_approx_delta + 1.d0 - this%RangeTol)) + end if + call this%Add(t_start, t_end, n, WantLog) + end subroutine TRanges_Add_delta + + + subroutine TRanges_Add(this, t_start, t_end, nstep, IsLog) + class(TRanges), intent(inout) :: this + logical, intent(in), optional :: IsLog + double precision, intent(in) :: t_start, t_end + integer, intent(in) :: nstep + type(TRange), allocatable, target :: NewRanges(:) + double precision, allocatable :: EndPoints(:), RequestDelta(:) + integer :: ixin, nreg, ix, i,j, nsteps + double precision :: min_request, max_request, min_log_step, max_log_step + double precision :: delta, diff, max_delta + logical :: WantLog + + WantLog = DefaultFalse(IsLog) + + if (WantLog) then + delta = log(t_end / t_start) / nstep + else + delta = (t_end - t_start) / nstep + end if + + if (t_end <= t_start) stop 'TRanges_Add: end must be larger than start' + if (nstep <= 0) stop 'TRanges_Add: nstep must be > 0' + + nreg = this%count + 1 + allocate(NewRanges(nreg)) + if (allocated(this%R) .and. this%count > 0) then +#if defined(__IBMCPP__) || defined(__xlC__) || defined(__xlc__) + ! avoid IBM compiler bug, from Angel de Vicente + ! detection of IBM compiler by preprocessors symbols taken from boost + do i= 1, this%count + NewRanges(i) = this%R(i) + end do +#else + NewRanges(1:this%count) = this%R(1:this%count) +#endif + end if + + allocate(EndPoints(0:nreg * 2)) + associate (AReg => NewRanges(nreg)) + AReg%Low = t_start + AReg%High = t_end + AReg%delta = delta + AReg%steps = nstep + AReg%IsLog = WantLog + end associate + + !Get end point in order + ix = 0 + do i=1, nreg + associate (AReg => NewRanges(i)) + if (ix==0) then + EndPoints(1) = AReg%Low + EndPoints(2) = AReg%High + ix = 2 + else + ixin = ix + do j=1,ixin + if (AReg%Low < EndPoints(j)) then + EndPoints(j+1:ix+1) = EndPoints(j:ix) + EndPoints(j) = AReg%Low + ix=ix+1 + exit + end if + end do + if (ixin == ix) then + ix = ix+1 + EndPoints(ix) = AReg%Low + ix = ix+1 + EndPoints(ix) = AReg%High + else + ixin = ix + do j=1,ixin + if (AReg%High < EndPoints(j)) then + EndPoints(j+1:ix+1) = EndPoints(j:ix) + EndPoints(j) = AReg%High + ix=ix+1 + exit + end if + end do + if (ixin == ix) then + ix = ix+1 + EndPoints(ix) = AReg%High + end if + end if + end if + end associate + end do + + !remove duplicate points + ixin = ix + ix = 1 + do i=2, ixin + if (EndPoints(i) /= EndPoints(ix)) then + ix=ix+1 + EndPoints(ix) = EndPoints(i) + end if + end do + + + !ix is the number of end points + this%Lowest = EndPoints(1) + this%Highest = EndPoints(ix) + this%count = 0 + + max_delta = this%Highest - this%Lowest + + if (.not. allocated(this%R) .or. size(this%R) < ix + 1) then + if (allocated (this%R)) deallocate (this%R) + allocate (this%R(ix + 1)) + end if + + allocate(RequestDelta(ix)) + + do i=1, ix - 1 + associate (AReg => this%R(i)) + AReg%Low = EndPoints(i) + AReg%High = EndPoints(i+1) + + ! max_delta = EndPoints(i+1) - EndPoints(i) + delta = max_delta + AReg%IsLog = .false. + + do j=1, nreg + if (AReg%Low >= NewRanges(j)%Low .and. AReg%Low < NewRanges(j)%High) then + if (NewRanges(j)%IsLog) then + if (AReg%IsLog) then + delta = min(delta,NewRanges(j)%delta) + else + min_log_step = AReg%Low*(exp(NewRanges(j)%delta)-1) + if (min_log_step < delta) then + max_log_step = AReg%High*(1-exp(-NewRanges(j)%delta)) + if (delta < max_log_step) then + delta = min_log_step + else + AReg%IsLog = .true. + delta = NewRanges(j)%delta + end if + end if + end if + else !New Range is not log + if (AReg%IsLog) then + max_log_step = AReg%High*(1-exp(-delta)) + if (NewRanges(j)%delta < max_log_step) then + min_log_step = AReg%Low*(exp(delta)-1) + if (min_log_step < NewRanges(j)%delta) then + AReg%IsLog = .false. + delta = min_log_step + else + delta = - log(1- NewRanges(j)%delta/AReg%High) + end if + end if + else + delta = min(delta, NewRanges(j)%delta) + end if + end if + end if + end do + + if (AReg%IsLog) then + Diff = log(AReg%High/AReg%Low) + else + Diff = AReg%High - AReg%Low + endif + if (delta >= Diff) then + AReg%delta = Diff + AReg%steps = 1 + else + AReg%steps = max(1,int(Diff/delta + 1.d0 - this%RangeTol)) + AReg%delta = Diff / AReg%steps + end if + + this%count = this%count + 1 + RequestDelta(this%count) = delta + + if (AReg%IsLog) then + if (AReg%steps ==1) then + AReg%Delta_min = AReg%High - AReg%Low + AReg%Delta_max = AReg%Delta_min + else + AReg%Delta_min = AReg%Low*(exp(AReg%delta)-1) + AReg%Delta_max = AReg%High*(1-exp(-AReg%delta)) + end if + else + AReg%Delta_max = AReg%delta + AReg%Delta_min = AReg%delta + end if + end associate + end do + + + !Get rid of tiny TRanges + ix = this%count + do i=ix, 1, -1 + associate (AReg => this%R(i)) + if (AReg%steps ==1) then + Diff = AReg%High - AReg%Low + if (AReg%IsLog) then + min_request = AReg%Low*(exp(RequestDelta(i))-1) + max_request = AReg%High*(1-exp(-RequestDelta(i))) + else + min_request = RequestDelta(i) + max_request = min_request + end if + if (i/= this%count) then !from i/= ix Mar08 + associate (LastReg => this%R(i+1)) + if (RequestDelta(i) >= AReg%delta .and. Diff <= LastReg%Delta_min & + .and. LastReg%Delta_min <= max_request) then + + LastReg%Low = AReg%Low + if (Diff > LastReg%Delta_min*this%RangeTol) then + LastReg%steps = LastReg%steps + 1 + end if + if (LastReg%IsLog) then + LastReg%delta = log(LastReg%High/LastReg%Low) / LastReg%steps + else + LastReg%delta = (LastReg%High -LastReg%Low) / LastReg%steps + end if + this%R(i:this%Count-1) = this%R(i+1:this%Count) + this%Count = this%Count -1 + cycle + end if + end associate + end if + if (i/=1) then + associate (LastReg => this%R(i-1)) + if (RequestDelta(i) >= AReg%delta .and. Diff <= LastReg%Delta_max & + .and. LastReg%Delta_max <= min_request) then + LastReg%High = AReg%High + !AlMat08 LastReg%Low = AReg%Low + if (Diff > LastReg%Delta_max*this%RangeTol) then + LastReg%steps = LastReg%steps + 1 + end if + if (LastReg%IsLog) then + LastReg%delta = log(LastReg%High/LastReg%Low) / LastReg%steps + else + LastReg%delta = (LastReg%High -LastReg%Low) / LastReg%steps + end if + this%R(i:this%Count-1) = this%R(i+1:this%Count) + this%Count = this%Count -1 + end if + end associate + end if + end if + end associate + end do + + + !Set up start indices and get total number of steps + nsteps = 1 + do i = 1, this%Count + associate (AReg => this%R(i)) + AReg%Start_index = nsteps + nsteps = nsteps + AReg%steps + if (AReg%IsLog) then + if (AReg%steps ==1) then + AReg%Delta_min = AReg%High - AReg%Low + AReg%Delta_max = AReg%Delta_min + else + AReg%Delta_min = AReg%Low*(exp(AReg%delta)-1) + AReg%Delta_max = AReg%High*(1-exp(-AReg%delta)) + end if + else + AReg%Delta_max = AReg%delta + AReg%Delta_min = AReg%delta + end if + end associate + end do + + this%npoints = nsteps + + deallocate(NewRanges, EndPoints, RequestDelta) + end subroutine TRanges_Add + + + subroutine TRanges_Write(this) + class(TRanges), intent(in), target :: this + integer :: i + + do i=1,this%count + associate (AReg => this%R(i)) + if (AReg%IsLog) then + Write (*,'("Range ",I3,":", 3E14.4," log")') i, AReg%Low, AReg%High, AReg%delta + else + Write (*,'("Range ",I3,":", 3E14.4," linear")') i, AReg%Low, AReg%High, AReg%delta + end if + end associate + end do + end subroutine TRanges_Write + + + end module RangeUtils diff --git a/Readme b/Readme new file mode 100644 index 0000000..3e5ed6f --- /dev/null +++ b/Readme @@ -0,0 +1,47 @@ +=================== +ForUtils +=================== +:ForUtils: Various fortran 2003/2008 utility classes. +:Version: 0.1 +:Author: Antony Lewis +:Homepage: https://github.com/cmbant/forutils + + +Description +============ + +ForUtils is a package of Fortran classes and convenience functions for +Fortran 2003/2008 programs. The utilities comprise these areas: + +* ArrayUtils - Find (minimal/maximal) index of an element in an array. +* FileUtils - Class for handling file access. +* IniObjects - Read/Write key/value style configuration files with array and +default value support. +* MatrixUtils - Read/Write matrices and interface to some BLAS/LAPACK routines. +* MiscUtils - Utility functions for optional arguments. +* MpiUtils - Wrappers for mpi-routines to compile without mpi library. +* ObjectLists - Lists of arbitrary objects including specializations for vectors. +* RandUtils - Some functions to generate random numbers. +* RangeUtils - Maintain sets of intervals. +* StringUtils - Utilities for strings, like concat of distinct types a.s.o. + + +Getting Started +================ + +Clone this git repository: + + $ git clone https://github.com/cmbant/forutils + +Compile: + + $ make all + + +Dependencies +============= +* Fortran 2008 compatible compiler - E.g., ifort 14+, gfortran 5.2+ (current +trunk will do). +* MPI library - Only when you want the MpiUtils fully functional. Without an +MPI library MpiUtils compile, but the functions are merely no-ops and the +makefile target DebugMPI and ReleaseMPI can not be build. diff --git a/StringUtils.f90 b/StringUtils.f90 new file mode 100644 index 0000000..61d5e4d --- /dev/null +++ b/StringUtils.f90 @@ -0,0 +1,358 @@ + module StringUtils + use MiscUtils + implicit none + + + INTERFACE RealToStr + module procedure SingleToStr, DoubleToStr + END INTERFACE RealToStr + + contains + + function IsWhiteSpace(C) + character, intent(in) :: C + logical IsWhiteSpace + + IsWhiteSpace = (C==' ') .or. (C==char(9)) + + end function IsWhiteSpace + + function GetParamCount() + integer GetParamCount + + GetParamCount = command_argument_count() + + end function GetParamCount + + function GetParam(i) + character(LEN=:), allocatable :: GetParam + integer, intent(in) :: i + character(LEN=:), allocatable :: tmp + integer l + + if (GetParamCount() < i) then + GetParam = '' + else + call get_command_argument(i,length=l) + allocate(character(l)::tmp) + call get_command_argument(i,value=tmp) + GetParam = trim(tmp) + end if + + end function GetParam + + function GetEnvironmentVariable(name, is_present) result(value) + character(LEN=*), intent(in) :: name + character(LEN=:), allocatable :: value + logical, intent(out), optional :: is_present + integer L, status + + call get_environment_variable(name, length=L, status=status) + if (present(is_present)) is_present = status==0 + if (status==0) then + allocate(character(L)::value) + call get_environment_variable(name, value, status=status) + end if + if (status/=0) value='' + + end function GetEnvironmentVariable + + + function StringStarts(S, substring, index) result(OK) + character(LEN=*), intent(in) :: S, substring + integer, intent(in), optional :: index + logical OK + integer start + + if (present(index)) then + start = index + else + start =1 + end if + + OK = S(start:min(len(S),start+len_trim(substring)-1))==substring + + end function StringStarts + + function StringTrimmed(S, trimmed) result(newS) + character(LEN=*), intent(in) :: S + logical, intent(in), optional :: trimmed + character(LEN=:), allocatable :: newS + + if (DefaultFalse(trimmed)) then + newS = trim(S) + else + newS = S + end if + + end function StringTrimmed + + subroutine StringReplace(FindS, RepS, S) + character(LEN=*), intent(in) :: FindS, RepS + character(LEN=:), allocatable, intent(inout) :: S + integer i + + i = index(S,FindS) + if (i>0) then + S = S(1:i-1)//trim(RepS)//S(i+len_trim(FindS):len_trim(S)) + end if + + end subroutine StringReplace + + function StringEscape(S, C, escape) result(newS) + character(LEN=*), intent(in) :: S + character, intent(in) :: C + character, intent(in), optional :: escape + character(LEN=:), allocatable :: newS, esc + integer i + character, parameter :: backslash = char(92) + + if (present(escape)) then + esc= escape + else + esc = backslash + end if + newS = '' + do i=1, len_trim(S) + if (S(i:i)==C) then + newS = newS //esc// C + else + newS = newS //S(i:i) + end if + end do + + end function StringEscape + + function Join(separator, S, S1,S2,S3,S4,S5,S6, trimmed) result(newS) + character(LEN=*), intent(in) :: Separator, S, S1 + character(LEN=*), optional :: S2,S3,S4,S5,S6 + character(LEN=:), allocatable :: newS + logical, intent(in), optional :: trimmed + + newS = StringTrimmed(S,trimmed)//Separator//StringTrimmed(S1,trimmed) + if (present(S2)) newS = newS //Separator //StringTrimmed(S2,trimmed) + if (present(S3)) newS = newS //Separator //StringTrimmed(S3,trimmed) + if (present(S4)) newS = newS //Separator //StringTrimmed(S4,trimmed) + if (present(S5)) newS = newS //Separator //StringTrimmed(S5,trimmed) + if (present(S6)) newS = newS //Separator //StringTrimmed(S6,trimmed) + + end function Join + + function numcat(S, num) + character(LEN=*) S + character(LEN=:), allocatable :: numcat + integer num + + numcat = concat(S,num) + end function numcat + + function IntToStr(I, minlen) + integer , intent(in) :: I + character(LEN=:), allocatable :: IntToStr + integer, intent(in), optional :: minlen + integer n + character (LEN=128) :: form, tmp + + if (present(minlen)) then + n = minlen + if (I<0) n=n+1 + form = concat('(I',n,'.',minlen,')') + write (tmp,form) i + IntToStr = trim(tmp) + else + write (tmp,*) i + IntToStr = trim(adjustl(tmp)) + end if + + end function IntToStr + + function StrToInt(S) + integer :: StrToInt + character(LEN=*), intent(in) :: S + + read(S,*) StrToInt + + end function StrToInt + + subroutine StringAppend(S,X) + character(LEN=:), allocatable, intent(inout) :: S + class(*), intent(in) :: X + + if (.not. allocated(S)) S='' + select type (X) + type is (character(LEN=*)) + S = S // trim(X) + type is (integer) + S = S // IntToStr(X) + type is (real) + S = S // RealToStr(X) + type is (double precision) + S=S //RealToStr(X) + class default + stop 'StringAppend: Unknown type' + end select + end subroutine + + function concat(S1,S2,S3,S4,S5,S6,S7,S8) result(outstr) + character(LEN=*), intent(in) :: S1 + class(*), intent(in) :: S2 + class(*), intent(in), optional :: S3, S4, S5, S6,S7,S8 + character(LEN=:), allocatable :: outstr + + outstr=S1 + call StringAppend(outstr,S2) + if (present(S3)) then + call StringAppend(outstr,S3) + if (present(S4)) then + call StringAppend(outstr,S4) + if (present(S5)) then + call StringAppend(outstr,S5) + if (present(S6)) then + call StringAppend(outstr,S6) + if (present(S7)) then + call StringAppend(outstr,S7) + if (present(S8)) then + call StringAppend(outstr,S8) + end if + end if + end if + end if + end if + end if + + end function concat + + + function DoubleToStr(R, figs) + double precision, intent(in) :: R + integer, intent(in), optional :: figs + character(LEN=:), allocatable :: DoubleToStr + + DoubleToStr = SingleToStr(real(R),figs) + + end function DoubleToStr + + function SingleToStr(R, figs) + real, intent(in) :: R + integer, intent(in), optional :: figs + character(LEN=:), allocatable :: SingleToStr + character(LEN=30) tmp + + if (abs(R)>=0.001 .or. R==0.) then + write (tmp,'(f12.6)') R + + tmp = adjustl(tmp) + if (present(figs)) then + SingleToStr = tmp(1:figs) + else + if (abs(R)>10000) then + write(tmp,*) R + SingleToStr = trim(adjustl(tmp)) + else + SingleToStr = tmp(1:6) + end if + end if + else + if (present(figs)) then + write (tmp,trim(numcat('(E',figs))//'.2)') R + else + write (tmp,'(G9.2)') R + end if + SingleToStr = trim(adjustl(tmp)) + end if + + end function SingleToStr + + function SubNextFormat(S, X) result(OK) + character(LEN=:), allocatable :: S + class(*) X + logical OK + integer ix, P, n + character c + character(LEN=:), allocatable :: form, fform, rep + + P=1 + do + ix=scan(S(P:),'%') + OK = ix/=0 .and. ix < len(S) + if (.not. OK) return + c = S(ix+P:ix+P) + if (c=='%') then + P=P+Ix+1 + else + exit + end if + end do + form = '' + do while( verify(c,'0123456789') == 0) + form = form // c + P=P+1 + c= S(ix+P:ix+P) + end do + select type (X) + type is (integer) + if (len(form)>0) then + n= StrToInt(form) + fform = 'I'//IntToStr(n) + if (form(1:1)=='0') fform=fform//'.'//IntToStr(n) + allocate(character(n)::rep) + write(rep,'('//fform//')') X + else + rep = IntToStr(X) + end if + if (c=='d' .or. c=='u') then + call StringReplace('%'//form//c, rep, S) + else + write(*,*) 'Wrong format for type: '//trim(S) + stop + end if + type is (Character(LEN=*)) + if (c/='s') then + write(*,*) 'Wrong format for type: '//trim(S) + stop + end if + call StringReplace('%s', X, S) + type is (double precision) + if (c/='f') then + write(*,*) 'Wrong format for type: '//trim(S) + stop + end if + call StringReplace('%f', RealToStr(X),S) + type is (real) + if (c/='f') then + write(*,*) 'Wrong format for type: '//trim(S) + stop + end if + call StringReplace('%f', RealToStr(X),S) + class default + stop 'Unsupported format type' + end select + + end function SubNextFormat + + function FormatString(formatst, i1,i2,i3,i4,i5,i6,i7,i8, allow_unused) result(S) + character(LEN=*), intent(in) :: formatst + class(*), intent(in),optional :: i1,i2,i3,i4,i5,i6,i7,i8 + logical, optional, intent(in) :: allow_unused + character(LEN=:), allocatable :: S + logical OK + !Note that this routine is incomplete and very simple (so buggy in complex cases) + !(should not substitute on the previously substituted string, etc, etc..) + !Can do things like FormatString('case %d, ans = %03d%%',i,percent) + S = formatst + OK = .true. + if (present(i1)) OK = SubNextFormat(S, i1) + if (OK .and. present(i2)) OK = SubNextFormat(S, i2) + if (OK .and. present(i3)) OK = SubNextFormat(S, i3) + if (OK .and. present(i4)) OK = SubNextFormat(S, i4) + if (OK .and. present(i5)) OK = SubNextFormat(S, i5) + if (OK .and. present(i6)) OK = SubNextFormat(S, i6) + if (OK .and. present(i7)) OK = SubNextFormat(S, i7) + if (OK .and. present(i8)) OK = SubNextFormat(S, i8) + if (.not. OK .and. .not. DefaultFalse(allow_unused)) stop 'FormatString: Wrong number or kind of formats in string' + call StringReplace('%%', '%', S) + + end function FormatString + + + end module StringUtils