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