diff --git a/io_factory.f90 b/io_factory.f90 new file mode 100644 index 0000000..36b37bc --- /dev/null +++ b/io_factory.f90 @@ -0,0 +1,37 @@ +module io_factory + use mpi_aux + use io_handler_base + use io_handler_ftn +#ifdef TROVE_USE_MPI_ + use io_handler_mpi +#endif + use errors + implicit none + + contains + subroutine allocateHandler(ioHandler) + class(ioHandlerBase), allocatable, intent(out) :: ioHandler +#ifdef TROVE_USE_MPI_ + if (blacs_size > 1) then + ! Only allocate MPI IO when compiled with MPI *and* + ! running with more than one MPI process + allocate(ioHandlerMPI::ioHandler) + else + allocate(ioHandlerFTN::ioHandler) + endif +#else + allocate(ioHandlerFTN::ioHandler) +#endif + end subroutine + + subroutine openFile(ioHandler, filename, err, action, position, status, form, access) + class(ioHandlerBase), allocatable, intent(out) :: ioHandler + character (len = *), intent(in) :: filename + type(ErrorType), intent(inout) :: err + character (len = *), intent(in) :: action + character (len = *), intent(in), optional :: position, status, form, access + + call allocateHandler(ioHandler) + call ioHandler%open(filename, err, action, position, status, form, access) + end subroutine +end module diff --git a/io_handler_base.f90 b/io_handler_base.f90 index eb4e9ff..ef78fca 100644 --- a/io_handler_base.f90 +++ b/io_handler_base.f90 @@ -7,16 +7,19 @@ module io_handler_base type, abstract :: ioHandlerBase contains procedure(open), deferred :: open + procedure(seek), deferred :: seek generic :: write => writeScalar, write1DArray, write2DArray, write2DArrayDistBlacs, write2DArrayDistColumn procedure(writeScalar), deferred :: writeScalar procedure(write1DArray), deferred :: write1DArray procedure(write2DArray), deferred :: write2DArray procedure(write2DArrayDistBlacs), deferred :: write2DArrayDistBlacs procedure(write2DArrayDistColumn), deferred :: write2DArrayDistColumn - generic :: read => readScalar, read1DArray, read2DArray + generic :: read => readScalar, read1DArray, read2DArray, read2DArrayDistBlacs, read2DArrayDistColumn procedure(readScalar), deferred :: readScalar procedure(read1DArray), deferred :: read1DArray procedure(read2DArray), deferred :: read2DArray + procedure(read2DArrayDistBlacs), deferred :: read2DArrayDistBlacs + procedure(read2DArrayDistColumn), deferred :: read2DArrayDistColumn end type ioHandlerBase abstract interface @@ -74,6 +77,25 @@ subroutine read2DArray(this, object) class(ioHandlerBase) :: this class(*), dimension(:,:), intent(out) :: object end subroutine + subroutine read2DArrayDistBlacs(this, object, descr, block_type) + import ioHandlerBase + import MPI_Datatype + class(ioHandlerBase) :: this + class(*), dimension(:,:), intent(out) :: object + integer, intent(in) :: descr(9) ! Description array outputted from co_block_type_init + type(MPI_Datatype), intent(in) :: block_type + end subroutine + subroutine read2DArrayDistColumn(this, object, dimen) + import ioHandlerBase + class(ioHandlerBase) :: this + class(*), intent(out) :: object(:,:) + integer, intent(in) :: dimen ! Dimension of entire distributed array + end subroutine + subroutine seek(this, offset) + import ioHandlerBase + class(ioHandlerBase) :: this + integer, intent(in) :: offset + end subroutine end interface private diff --git a/io_handler_ftn.f90 b/io_handler_ftn.f90 index fb48e26..ee7bc06 100644 --- a/io_handler_ftn.f90 +++ b/io_handler_ftn.f90 @@ -12,6 +12,7 @@ module io_handler_ftn integer :: iounit = 0 integer :: stat = 0 logical :: isOpen = .false. + character (len=20) :: accessVal contains procedure :: writeScalar => writeScalarFTN procedure :: write1DArray => write1DArrayFTN @@ -21,32 +22,20 @@ module io_handler_ftn procedure :: readScalar => readScalarFTN procedure :: read1DArray => read1DArrayFTN procedure :: read2DArray => read2DArrayFTN + procedure :: read2DArrayDistBlacs => read2DArrayDistBlacsFTN + procedure :: read2DArrayDistColumn => read2DArrayDistColumnFTN procedure :: open procedure :: close + procedure :: seek final :: destroyIoHandlerFTN end type ioHandlerFTN - ! Constructor - interface ioHandlerFTN - procedure :: newIoHandlerFTN - end interface ioHandlerFTN - private public :: ioHandlerFTN contains - type(ioHandlerFTN) function newIoHandlerFTN(fname, err, action, position, status, form, access) result(this) - ! writer FTN constructor - type(ErrorType), intent(inout) :: err - character (len = *), intent(in) :: fname - character (len = *), intent(in) :: action - character (len = *), intent(in), optional :: position, status, form, access - - call this%open(fname, err, action, position, status, form, access) - end function - subroutine destroyIoHandlerFTN(this) type(ioHandlerFTN) :: this call this%close() @@ -83,17 +72,21 @@ subroutine open(this, fname, err, action, position, status, form, access) end if if (present(access)) then - accessVal = access + this%accessVal = access else - accessVal = 'sequential' - end if + this%accessVal = 'sequential' + endif - print *, "FTN: Opening ", trim(fname), " with ", \ - trim(positionVal), " ", trim(statusVal), " ", trim(formVal), " ", trim(accessVal) + print *, "FTN: Opening ", trim(fname), " with ", & + trim(action), " ", & + trim(positionVal), " ", & + trim(statusVal), " ", & + trim(formVal), " ", & + trim(this%accessVal) - open(newunit=this%iounit, action=action,\ - form=formVal, position=positionVal, status=statusVal, file=fname,\ - iostat=this%stat) + open(newunit=this%iounit, action=action, & + form=formVal, position=positionVal, status=statusVal, & + access=this%accessVal, file=fname, iostat=this%stat) if (this%stat == 0) then this%isOpen = .true. @@ -110,10 +103,23 @@ subroutine close(this) endif end subroutine + subroutine seek(this, offset) + class(ioHandlerFTN) :: this + INTEGER, PARAMETER :: SEEK_SET = 0, SEEK_CUR = 1, SEEK_END = 2 + integer, intent(in) :: offset + integer :: total_offset + + if (trim(this%accessVal) == "sequential" .and. offset .ne. 0) then + ! Add two bookend offsets + total_offset = offset + 2*4 + endif + call fseek(this%iounit, total_offset, SEEK_CUR) + end subroutine + subroutine writeScalarFTN(this, object) class(ioHandlerFTN) :: this class(*), intent(in) :: object - print *, "writing scalar object with FTN IO" + select type(object) type is (integer) write(this%iounit) object @@ -122,16 +128,16 @@ subroutine writeScalarFTN(this, object) type is (complex) write(this%iounit) object type is (character(len=*)) - write(this%iounit) object + write(this%iounit) trim(object) class default - print *, "ERROR: Tried to write unsupported type" + stop "ioHandlerFTN%writeScalarFTN: Tried to write unsupported type" end select end subroutine subroutine write1DArrayFTN(this, object) class(ioHandlerFTN) :: this class(*), dimension(:), intent(in) :: object - print *, "writing 1D array with FTN IO" + select type(object) type is (integer) write(this%iounit) object @@ -140,14 +146,14 @@ subroutine write1DArrayFTN(this, object) type is (complex) write(this%iounit) object class default - print *, "ERROR: Tried to write unsupported type" + stop "ioHandlerFTN%write1DArrayFTN: Tried to write unsupported type" end select end subroutine subroutine write2DArrayFTN(this, object) class(ioHandlerFTN) :: this class(*), dimension(:,:), intent(in) :: object - print *, "writing 2D array with FTN IO" + select type(object) type is (integer(int32)) write(this%iounit) object @@ -162,7 +168,7 @@ subroutine write2DArrayFTN(this, object) type is (complex(kind=8)) write(this%iounit) object class default - print *, "ERROR: Tried to write unsupported type" + stop "ioHandlerFTN%write2DArrayFTN: Tried to write unsupported type" end select end subroutine @@ -192,48 +198,91 @@ subroutine write2DArrayDistColumnFTN(this, object, mdimen) subroutine readScalarFTN(this, object) class(ioHandlerFTN) :: this class(*), intent(out) :: object - print *, "reading object with FTN IO" + select type(object) - type is (integer) + type is (integer(int32)) read(this%iounit) object - type is (real) + type is (integer(int64)) read(this%iounit) object - type is (complex) + type is (real(real32)) + read(this%iounit) object + type is (real(real64)) + read(this%iounit) object + type is (complex(kind=4)) + read(this%iounit) object + type is (complex(kind=8)) + read(this%iounit) object + type is (character(len=*)) read(this%iounit) object class default - print *, "Unsupported type!" + stop "ioHandlerFTN%readScalarFTN: Tried to read unsupported type" end select end subroutine subroutine read1DArrayFTN(this, object) class(ioHandlerFTN) :: this class(*), dimension(:), intent(out) :: object - print *, "reading 1D array with FTN IO" + select type(object) - type is (integer) + type is (integer(int32)) read(this%iounit) object - type is (real) + type is (integer(int64)) read(this%iounit) object - type is (complex) + type is (real(real32)) + read(this%iounit) object + type is (real(real64)) + read(this%iounit) object + type is (complex(kind=4)) + read(this%iounit) object + type is (complex(kind=8)) read(this%iounit) object class default - print *, "Unsupported type!" + stop "ioHandlerFTN%read1DArrayFTN: Tried to read unsupported type" end select end subroutine subroutine read2DArrayFTN(this, object) class(ioHandlerFTN) :: this class(*), dimension(:,:), intent(out) :: object - print *, "reading 2D array with FTN IO" + select type(object) - type is (integer) + type is (integer(int32)) read(this%iounit) object - type is (real) + type is (integer(int64)) read(this%iounit) object - type is (complex) + type is (real(real32)) + read(this%iounit) object + type is (real(real64)) + read(this%iounit) object + type is (complex(kind=4)) + read(this%iounit) object + type is (complex(kind=8)) read(this%iounit) object class default - print *, "Unsupported type!" + stop "ioHandlerFTN%read2DArrayFTN: Tried to read unsupported type" end select end subroutine + + subroutine read2DArrayDistBlacsFTN(this, object, descr, block_type) + ! Read blacs-distributed arrays + + class(ioHandlerFTN) :: this + class(*), dimension(:,:), intent(out) :: object + integer, intent(in) :: descr(9) ! Description array outputted from co_block_type_init + type(MPI_Datatype), intent(in) :: block_type + + ! Using the fortran io_handler means array isn't distributed, just read normally + call this%read2DArray(object) + end subroutine + + subroutine read2DArrayDistColumnFTN(this, object, dimen) + ! read arrays distributed as columns using co_distr_data + + class(ioHandlerFTN) :: this + class(*), dimension(:,:), intent(out) :: object + integer, intent(in) :: dimen + + ! Using the fortran io_handler means array isn't distributed, just read normally + call this%read2DArray(object) + end subroutine end module diff --git a/io_handler_mpi.f90 b/io_handler_mpi.f90 index 4d12362..8b0cef3 100644 --- a/io_handler_mpi.f90 +++ b/io_handler_mpi.f90 @@ -17,7 +17,7 @@ type is (character(len=*)); \ call function(handle, obj, size, mytype, status, err); \ class default; \ - write(*,*) 'Not covered'; \ + stop 'MPI: tried to handle unsupported type'; \ end select @@ -31,10 +31,11 @@ module io_handler_mpi implicit none type, extends(ioHandlerBase) :: ioHandlerMPI - integer (kind=MPI_Offset_kind) :: offset = 0 - integer :: rank = 0 + integer (kind=MPI_Offset_kind) :: bookendBytes = 0 + integer :: rank = -1 type(MPI_File) :: fileh logical :: isOpen = .false. + character (len=20) :: accessVal contains procedure :: writeScalar => writeScalarMPI procedure :: write1DArray => write1DArrayMPI @@ -44,31 +45,20 @@ module io_handler_mpi procedure :: readScalar => readScalarMPI procedure :: read1DArray => read1DArrayMPI procedure :: read2DArray => read2DArrayMPI + procedure :: read2DArrayDistBlacs => read2DArrayDistBlacsMPI + procedure :: read2DArrayDistColumn => read2DArrayDistColumnMPI procedure :: open procedure :: close + procedure :: seek final :: destroyIoHandlerMPI end type ioHandlerMPI - interface ioHandlerMPI - procedure :: newIoHandlerMPI - end interface ioHandlerMPI - private public :: ioHandlerMPI contains - type(ioHandlerMPI) function newIoHandlerMPI(fname, err, action, position, status, form, access) result(this) - ! writer MPI constructor - character (len = *), intent(in) :: fname - type(ErrorType), intent(inout) :: err - character (len = *), intent(in) :: action - character (len = *), intent(in), optional :: position, status, form, access - - call this%open(fname, err, action, position, status, form, access) - end function - subroutine destroyIoHandlerMPI(this) type(ioHandlerMPI) :: this call this%close() @@ -81,7 +71,7 @@ subroutine open(this, fname, err, action, position, status, form, access) type(ErrorType), intent(inout) :: err character (len = *), intent(in) :: action character (len = *), intent(in), optional :: position, status, form, access - character (len = 20) :: positionVal, statusVal, formVal, accessVal + character (len = 20) :: positionVal, statusVal, formVal integer :: ierr if (this%isOpen) then @@ -107,22 +97,30 @@ subroutine open(this, fname, err, action, position, status, form, access) end if if (present(access)) then - accessVal = access + this%accessVal = access else - accessVal = 'sequential' - end if + this%accessVal = 'sequential' + endif - print *, "MPI: Opening ", trim(fname), " with ", \ - trim(positionVal), " ", trim(statusVal), " ", trim(formVal), " ", trim(accessVal) + if (trim(this%accessVal) == "sequential") then + this%bookendBytes = 4 + else + this%bookendBytes = 0 + endif ! FIXME use above flags to change open behaviour call MPI_Comm_rank(MPI_COMM_WORLD, this%rank, ierr) + if(this%rank == 0) then + print *, "MPI: Opening ", trim(fname), " with ", \ + trim(action), " ", trim(positionVal), " ", trim(statusVal), " ", trim(formVal), " ", trim(this%accessVal) + endif + ! FIXME is there a better way to set MPI_MODE_* flags? if(trim(action) == 'write') then call MPI_File_open(MPI_COMM_WORLD, fname, MPI_MODE_WRONLY + MPI_MODE_CREATE, MPI_INFO_NULL, this%fileh, ierr) - else + else if(trim(action) == 'read') then call MPI_File_open(MPI_COMM_WORLD, fname, MPI_MODE_RDONLY, MPI_INFO_NULL, this%fileh, ierr) endif @@ -140,6 +138,14 @@ subroutine close(this) ! FIXME handle error end subroutine close + subroutine seek(this, offset) + class(ioHandlerMPI) :: this + integer, intent(in) :: offset + integer(kind=MPI_OFFSET_KIND) :: total_offset + + call MPI_File_seek(this%fileh, offset + 2*this%bookendBytes, MPI_SEEK_CUR) + end subroutine + subroutine getMPIVarInfo(object, byteSize, mpiType) class(*), intent(in) :: object integer, intent(out) :: byteSize @@ -172,15 +178,38 @@ subroutine getMPIVarInfo(object, byteSize, mpiType) end select end subroutine + subroutine writeBookendBytes(this, bytes) + ! write first and last bookends containing object size in bytes + + class(ioHandlerMPI) :: this + integer :: bytes, ierr + integer(kind=MPI_OFFSET_KIND) :: offset + + if (this%rank == 0 .and. trim(this%accessVal)=="sequential") then + ! Get original position + call MPI_File_get_position(this%fileh, offset, ierr) + + ! Write bookends + call MPI_File_write(this%fileh, bytes, 1, MPI_INTEGER, & + MPI_STATUS_IGNORE, ierr) + call MPI_File_seek(this%fileh, int(bytes,MPI_OFFSET_KIND), MPI_SEEK_CUR, ierr) + call MPI_File_write(this%fileh, bytes, 1, MPI_INTEGER, & + MPI_STATUS_IGNORE, ierr) + + ! Reset position + call MPI_File_seek(this%fileh, offset, MPI_SEEK_SET, ierr) + endif + end subroutine + subroutine writeScalarMPI(this, object) class(ioHandlerMPI) :: this class(*), intent(in) :: object integer :: byteSize, ierr, length type(MPI_Datatype) :: mpiType + integer(kind = MPI_OFFSET_KIND) :: offset call getMPIVarInfo(object, byteSize, mpiType) - this%offset = this%offset + 4+byteSize+4 select type(object) type is (character(len=*)) @@ -189,22 +218,21 @@ subroutine writeScalarMPI(this, object) length = 1 end select - if (this%rank /= 0) then - return - end if + call writeBookendBytes(this, byteSize) - call MPI_File_write(this%fileh, byteSize, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, ierr) - MPI_WRAPPER(MPI_File_write, this%fileh, object, length, mpiType, & - MPI_STATUS_IGNORE, ierr) - call MPI_File_write(this%fileh, byteSize, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, ierr) + call MPI_File_seek(this%fileh, this%bookendBytes, MPI_SEEK_CUR, ierr) + if (this%rank == 0) then + MPI_WRAPPER(MPI_File_write, this%fileh, object, length, mpiType, MPI_STATUS_IGNORE, ierr) + else + call MPI_File_seek(this%fileh, int(byteSize,kind=MPI_OFFSET_KIND), MPI_SEEK_CUR, ierr) + end if + call MPI_File_seek(this%fileh, this%bookendBytes, MPI_SEEK_CUR, ierr) end subroutine subroutine write1DArrayMPI(this, object) class(ioHandlerMPI) :: this class(*), intent(in) :: object(:) - print *, "ERROR: 1D array saving not currently supported" + stop "ERROR: 1D array saving not currently supported by MPI writer" end subroutine subroutine write2DArrayMPI(this, object) @@ -212,27 +240,24 @@ subroutine write2DArrayMPI(this, object) class(*), intent(in) :: object(:,:) type(MPI_Datatype) :: mpiType - integer :: byteSize, globalSize, ierr, length + integer :: byteSize, globalSize, ierr, length, arrSizeBytes - integer(kind = MPI_OFFSET_KIND) :: arrSizeBytes + integer(kind = MPI_OFFSET_KIND) :: offset globalSize = size(object) call getMPIVarInfo(object(1,1), byteSize, mpiType) arrSizeBytes = globalSize*byteSize - this%offset = this%offset + 4 + arrSizeBytes + 4 + call writeBookendBytes(this, arrSizeBytes) - if (this%rank /= 0) then - return + call MPI_File_seek(this%fileh, this%bookendBytes, MPI_SEEK_CUR, ierr) + if (this%rank == 0) then + MPI_WRAPPER(MPI_File_write, this%fileh, object, globalSize, mpiType, MPI_STATUS_IGNORE, ierr) + else + call MPI_File_seek(this%fileh, int(arrSizeBytes,kind=MPI_OFFSET_KIND), MPI_SEEK_CUR, ierr) end if - - call MPI_File_write(this%fileh, arrSizeBytes, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, ierr) - MPI_WRAPPER(MPI_File_write, this%fileh, object, globalSize, mpiType, & - MPI_STATUS_IGNORE, ierr) - call MPI_File_write(this%fileh, arrSizeBytes, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, ierr) + call MPI_File_seek(this%fileh, this%bookendBytes, MPI_SEEK_CUR, ierr) end subroutine subroutine write2DArrayDistBlacsMPI(this, object, descr, block_type) @@ -242,8 +267,8 @@ subroutine write2DArrayDistBlacsMPI(this, object, descr, block_type) type(MPI_Datatype), intent(in) :: block_type ! subarray type outputed from co_block_type_init type(MPI_Datatype) :: mpiType - integer :: byteSize, globalSize, ierr - integer(kind = MPI_OFFSET_KIND) :: arrSizeBytes + integer :: byteSize, arrSizeBytes, globalSize, ierr + integer(kind = MPI_OFFSET_KIND) :: offset, disp integer :: dims(2) @@ -253,27 +278,21 @@ subroutine write2DArrayDistBlacsMPI(this, object, descr, block_type) call getMPIVarInfo(object(1,1), byteSize, mpiType) arrSizeBytes = globalSize*byteSize - if (this%rank == 0) then - ! write first and last bookends containing array size in bytes - call MPI_File_write(this%fileh, arrSizeBytes, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, ierr) - call MPI_File_seek(this%fileh, arrSizeBytes, MPI_SEEK_CUR, ierr) - call MPI_File_write(this%fileh, arrSizeBytes, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, ierr) - endif - ! Offset first bookend - this%offset = this%offset + 4 - ! Set file view including offset - call MPI_File_set_view(this%fileh, this%offset, mpiType, block_type, & + call writeBookendBytes(this, arrSizeBytes) + + call MPI_File_get_position(this%fileh, offset, ierr) + ! Get initial displacement in file + call MPI_File_get_byte_offset(this%fileh, offset, disp, ierr) + + ! Set file view including offsetting bookend + call MPI_File_set_view(this%fileh, disp+this%bookendBytes, mpiType, block_type, & 'native', MPI_INFO_NULL, ierr) ! Write array in parallel - MPI_WRAPPER(MPI_File_write_all, this%fileh, object, size(object), mpiType, & - MPI_STATUS_IGNORE, ierr) - ! Offset by size of array and end bookend integer - this%offset = this%offset + arrSizeBytes + 4 - ! Reset file view back to regular ol bytes - call MPI_File_set_view(this%fileh, this%offset, MPI_BYTE, MPI_BYTE, & + MPI_WRAPPER(MPI_File_write_all, this%fileh, object, size(object), mpiType, MPI_STATUS_IGNORE, ierr) + ! Reset file view + call MPI_File_set_view(this%fileh, int(0,MPI_OFFSET_KIND), MPI_BYTE, MPI_BYTE, & 'native', MPI_INFO_NULL, ierr) + call MPI_File_seek(this%fileh, disp+this%bookendBytes+arrSizeBytes+this%bookendBytes, MPI_SEEK_SET) end subroutine subroutine write2DArrayDistColumnMPI(this, object, mdimen) @@ -282,62 +301,137 @@ subroutine write2DArrayDistColumnMPI(this, object, mdimen) integer, intent(in) :: mdimen ! Dimension of entire distributed array type(MPI_Datatype) :: mpiType - integer :: byteSize, globalSize, ierr, writestat - integer(kind = MPI_OFFSET_KIND) :: arrSizeBytes + integer :: byteSize, globalSize, ierr, writestat, arrSizeBytes + integer(kind = MPI_OFFSET_KIND) :: offset, disp globalSize = mdimen**2 call getMPIVarInfo(object(1,1), byteSize, mpiType) arrSizeBytes = globalSize*byteSize - ! TODO what if format isn't sequential?? - if (this%rank == 0) then - ! write first and last bookends containing array size in bytes - call MPI_File_write(this%fileh, arrSizeBytes, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, ierr) - call MPI_File_seek(this%fileh, arrSizeBytes, MPI_SEEK_CUR, ierr) - call MPI_File_write(this%fileh, arrSizeBytes, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, ierr) - endif - ! Offset first bookend - this%offset = this%offset + 4 - ! Seek to byte after bookend - call MPI_File_seek_shared(this%fileh, this%offset, MPI_SEEK_SET, ierr) + call writeBookendBytes(this, arrSizeBytes) + + ! Get individual pointer offset + call MPI_File_get_position(this%fileh, offset, ierr) + ! Set shared pointer to individual pointer + bookend + call MPI_File_seek_shared(this%fileh, offset+this%bookendBytes, MPI_SEEK_SET, ierr) ! Write array in parallel MPI_WRAPPER(MPI_File_write_ordered,this%fileh,object,1,mpitype_column,MPI_STATUS_IGNORE,ierr) - ! Offset by size of array and end bookend integer - this%offset = this%offset + arrSizeBytes + 4 - ! Ensure all file pointers point to end of array - call MPI_File_seek(this%fileh, this%offset, MPI_SEEK_SET, ierr) + ! Skip over last bookend + call MPI_File_seek_shared(this%fileh, int(this%bookendBytes,MPI_OFFSET_KIND), MPI_SEEK_CUR, ierr) + + ! Set individual pointer to match shared + call MPI_File_get_position_shared(this%fileh, offset, ierr) + call MPI_File_seek(this%fileh, offset, MPI_SEEK_SET, ierr) end subroutine subroutine readScalarMPI(this, object) class(ioHandlerMPI) :: this class(*), intent(out) :: object - print *, "reading object to MPI IO" - ! Example object handling + + integer :: byteSize, ierr, length + type(MPI_Datatype) :: mpiType + + call getMPIVarInfo(object, byteSize, mpiType) + select type(object) - type is (integer) - print *, object - type is (real) - print *, object - type is (complex) - print *, object + type is (character(len=*)) + length = len(object) class default - print *, "Unsupported type!" + length = 1 end select + + call MPI_File_seek(this%fileh, int(this%bookendBytes,MPI_OFFSET_KIND), MPI_SEEK_CUR, ierr) + MPI_WRAPPER(MPI_File_read_all, this%fileh, object, length, mpiType, MPI_STATUS_IGNORE, ierr) + call MPI_File_seek(this%fileh, int(this%bookendBytes,MPI_OFFSET_KIND), MPI_SEEK_CUR, ierr) end subroutine subroutine read1DArrayMPI(this, object) class(ioHandlerMPI) :: this class(*), dimension(:), intent(out) :: object - print *, "reading 1D array to MPI IO" + stop "reading 1D array with MPI IO not supported" end subroutine subroutine read2DArrayMPI(this, object) class(ioHandlerMPI) :: this - class(*), dimension(:,:), intent(out) :: object - print *, "reading 2D array to MPI IO" + class(*), intent(out) :: object(:,:) + + type(MPI_Datatype) :: mpiType + integer :: byteSize, globalSize, ierr, length, arrSizeBytes + + integer(kind = MPI_OFFSET_KIND) :: offset + + globalSize = size(object) + + call getMPIVarInfo(object(1,1), byteSize, mpiType) + arrSizeBytes = globalSize*byteSize + + call MPI_File_seek(this%fileh, this%bookendBytes, MPI_SEEK_CUR, ierr) + MPI_WRAPPER(MPI_File_read_all, this%fileh, object, globalSize, mpiType, MPI_STATUS_IGNORE, ierr) + call MPI_File_seek(this%fileh, this%bookendBytes, MPI_SEEK_CUR, ierr) + end subroutine + + subroutine read2DArrayDistBlacsMPI(this, object, descr, block_type) + class(ioHandlerMPI) :: this + class(*), intent(out) :: object(:,:) + integer, intent(in) :: descr(9) ! Description array outputted from co_block_type_init + type(MPI_Datatype), intent(in) :: block_type ! subarray type outputed from co_block_type_init + type(MPI_Datatype) :: mpiType + integer :: byteSize, arrSizeBytes, globalSize, ierr + integer(kind = MPI_OFFSET_KIND) :: offset, disp + + integer :: dims(2) + + dims(:) = descr(3:4) + globalSize = dims(1)*dims(2) + + call getMPIVarInfo(object(1,1), byteSize, mpiType) + arrSizeBytes = globalSize*byteSize + + call MPI_File_get_position(this%fileh, offset, ierr) + ! Get initial displacement in file + call MPI_File_get_byte_offset(this%fileh, offset, disp, ierr) + ! Set file view including offsetting bookend + call MPI_File_set_view(this%fileh, disp+this%bookendBytes, mpiType, block_type, & + 'native', MPI_INFO_NULL, ierr) + ! Read array in parallel + MPI_WRAPPER(MPI_File_read_all, this%fileh, object, size(object), mpiType, MPI_STATUS_IGNORE, ierr) + ! Reset file view back to regular ol bytes, including bookends and array we've just written + call MPI_File_set_view(this%fileh, int(0,MPI_OFFSET_KIND), MPI_BYTE, MPI_BYTE, & + 'native', MPI_INFO_NULL, ierr) + call MPI_File_seek(this%fileh, disp+this%bookendBytes+arrSizeBytes+this%bookendBytes, MPI_SEEK_SET) + end subroutine + + subroutine read2DArrayDistColumnMPI(this, object, dimen) + class(ioHandlerMPI) :: this + class(*), intent(out) :: object(:,:) + integer, intent(in) :: dimen ! Dimension of entire distributed array + + type(MPI_Datatype) :: mpiType_ + integer :: byteSize, arrSizeBytes, globalSize, ierr + integer(kind = MPI_OFFSET_KIND) :: offset, disp + + call getMPIVarInfo(object(1,1), byteSize, mpiType_) + arrSizeBytes = dimen*dimen*byteSize + + ! Get individual pointer offset + call MPI_File_get_position(this%fileh, offset, ierr) + call MPI_File_get_byte_offset(this%fileh, offset, disp, ierr) + + ! Set shared pointer to individual pointer + bookend + call MPI_File_seek_shared(this%fileh, disp+this%bookendBytes, MPI_SEEK_SET, ierr) + + ! Read array in parallel + MPI_WRAPPER(MPI_File_read_ordered,this%fileh,object,1,mpitype_column,MPI_STATUS_IGNORE,ierr) + + ! Set shared pointer to point to byte after last bookend + call MPI_File_seek_shared(this%fileh, disp+this%bookendBytes+arrSizeBytes+this%bookendBytes, MPI_SEEK_SET) + + ! Set individual pointer to match shared + call MPI_File_get_position_shared(this%fileh, offset, ierr) + ! Ensure every process syncs the correct shared offset + call MPI_Barrier(MPI_COMM_WORLD) + call MPI_File_seek(this%fileh, offset, MPI_SEEK_SET, ierr) end subroutine end module diff --git a/makefile b/makefile index bff92c5..a44a80e 100644 --- a/makefile +++ b/makefile @@ -14,6 +14,12 @@ COMPILER ?= intel MODE ?= release USE_MPI ?= 0 +ifneq ($(USE_MPI),1) +ifneq ($(USE_MPI),0) +$(error USE_MPI "$(USE_MPI)" should be set to 1 to enable MPI or 0 (default) to disable MPI.) +endif +endif + # Intel ####### ifeq ($(strip $(COMPILER)),intel) @@ -29,7 +35,8 @@ ifeq ($(strip $(COMPILER)),intel) endif LAPACK = -mkl=parallel - ifneq ($(strip $(USE_MPI)),0) + ifeq ($(USE_MPI),1) + FC=mpiifort LAPACK += -lmkl_scalapack_lp64 -lmkl_blacs_intelmpi_lp64 endif @@ -37,7 +44,7 @@ ifeq ($(strip $(COMPILER)),intel) ########## else ifeq ($(strip $(COMPILER)),gfortran) FC = gfortran - FFLAGS = -cpp -std=gnu -fopenmp -march=native -ffree-line-length-512 -fcray-pointer -I$(OBJDIR) -J$(OBJDIR) + FFLAGS = -cpp -std=gnu -fopenmp -march=native -ffree-line-length-none -fcray-pointer -I$(OBJDIR) -J$(OBJDIR) GCC_VERSION_GT_10 := $(shell expr `gcc -dumpversion | cut -f1 -d.` \>= 10) ifeq "${GCC_VERSION_GT_10}" "1" @@ -54,8 +61,9 @@ else ifeq ($(strip $(COMPILER)),gfortran) endif LAPACK = -L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_gf_lp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm -ldl - ifneq ($(strip $(USE_MPI)),0) + ifneq ($(USE_MPI),0) # Assume we're using openmpi with gfortran + FC = mpif90 LAPACK += -lmkl_blacs_openmpi_lp64 -lmkl_scalapack_lp64 endif else @@ -64,8 +72,7 @@ endif CPPFLAGS = -D_EXTFIELD_DEBUG_ -ifneq ($(strip $(USE_MPI)),0) - FC = mpif90 +ifeq ($(USE_MPI),1) FFLAGS += -DTROVE_USE_MPI_ endif @@ -92,7 +99,7 @@ user_pot_dir=. TARGET=$(BINDIR)/$(EXE) MPI_SRCS = -ifneq ($(strip $(USE_MPI)),0) +ifeq ($(USE_MPI),1) MPI_SRCS += io_handler_mpi.f90 endif @@ -104,7 +111,7 @@ SRCS := timer.f90 accuracy.f90 diag.f90 dipole.f90 extfield.f90 fields.f90 fwigx pot_abcd.f90 pot_c2h4.f90 pot_c2h6.f90 pot_c3h6.f90 pot_ch3oh.f90 \ pot_xy2.f90 pot_xy3.f90 pot_xy4.f90 pot_zxy2.f90 pot_zxy3.f90 \ prop_xy2.f90 prop_xy2_quad.f90 prop_xy2_spinrot.f90 prop_xy2_spinspin.f90 \ - io_handler_base.f90 io_handler_ftn.f90 \ + io_handler_base.f90 io_handler_ftn.f90 io_factory.f90 \ refinement.f90 richmol_data.f90 rotme_cart_tens.f90 symmetry.f90 tran.f90 trove.f90 $(pot_user).f90 $(MPI_SRCS) OBJS := ${SRCS:.f90=.o} @@ -159,25 +166,29 @@ tarball: checkin: ci -l Makefile *.f90 +ifeq ($(USE_MPI),1) test: regression-tests unit-tests-nompi unit-tests-mpi +else +test: regression-tests unit-tests-nompi +endif regression-tests: $(TARGET) echo "Running regression tests" cd test/regression; ./run_regression_tests.sh -unit-tests-nompi: $(TARGET) +unit-tests-nompi: io_handler_ftn.o $(MAKE) -C test/unit LAPACK="$(LAPACK)" test_io echo "Running unit tests without MPI" test/unit/test_io -ifneq ($(strip $(USE_MPI)),0) -unit-tests-mpi: $(TARGET) +ifeq ($(USE_MPI),1) +unit-tests-mpi: io_handler_mpi.o $(MAKE) -C test/unit LAPACK="$(LAPACK)" test_mpi_io echo "Running unit tests with MPI" - mpirun -n 4 --mca opal_warn_on_missing_libcuda 0 test/unit/test_mpi_io + mpirun -n 4 test/unit/test_mpi_io else -unit-tests-mpi: $(TARGET) - echo "Skipping unit tests with MPI (USE_MPI not set)" +unit-tests-mpi: + $(error set USE_MPI=1 to compile & test with MPI) endif ################################################################################ @@ -218,7 +229,7 @@ mol_xy.o: mol_xy.f90 accuracy.o moltype.o mol_zxy2.o: mol_zxy2.f90 accuracy.o moltype.o mol_zxy3.o: mol_zxy3.f90 accuracy.o moltype.o lapack.o mpi_aux.o: mpi_aux.f90 accuracy.o timer.o -perturbation.o: perturbation.f90 accuracy.o molecules.o moltype.o lapack.o plasma.o fields.o timer.o symmetry.o me_numer.o diag.o mpi_aux.o io_handler_base.o io_handler_ftn.o $(MPI_OBJS) +perturbation.o: perturbation.f90 accuracy.o molecules.o moltype.o lapack.o plasma.o fields.o timer.o symmetry.o me_numer.o diag.o mpi_aux.o io_factory.o io_handler_base.o io_handler_ftn.o $(MPI_OBJS) plasma.o: plasma.f90 accuracy.o timer.o pot_abcd.o: pot_abcd.f90 accuracy.o moltype.o lapack.o pot_c2h4.o: pot_c2h4.f90 accuracy.o moltype.o @@ -239,8 +250,9 @@ richmol_data.o: richmol_data.f90 accuracy.o timer.o rotme_cart_tens.o: rotme_cart_tens.f90 accuracy.o timer.o fwigxjpf.o moltype.o accuracy.o symmetry.o: symmetry.f90 accuracy.o timer.o timer.o: timer.f90 accuracy.o -tran.o: tran.f90 accuracy.o timer.o me_numer.o molecules.o fields.o moltype.o symmetry.o perturbation.o mpi_aux.o +tran.o: tran.f90 accuracy.o timer.o me_numer.o molecules.o fields.o moltype.o symmetry.o perturbation.o mpi_aux.o io_factory.o io_handler_base.o io_handler_ftn.o trove.o: trove.f90 accuracy.o fields.o perturbation.o symmetry.o timer.o moltype.o dipole.o refinement.o tran.o extfield.o io_handler_base.o: io_handler_base.f90 errors.o mpi_aux.o io_handler_ftn.o: io_handler_ftn.f90 io_handler_base.o errors.o mpi_aux.o io_handler_mpi.o: io_handler_mpi.f90 io_handler_base.o errors.o mpi_aux.o +io_factory.o: io_factory.f90 io_handler_base.o io_handler_ftn.o mpi_aux.o $(MPI_OBJS) diff --git a/mpi_aux.f90 b/mpi_aux.f90 index a020bc1..33e2fb3 100644 --- a/mpi_aux.f90 +++ b/mpi_aux.f90 @@ -49,9 +49,10 @@ module mpi_aux integer, parameter :: MPI_OFFSET_KIND=8 #endif - integer,dimension(:),allocatable :: proc_sizes, proc_offsets, send_or_recv + integer,dimension(:),allocatable :: send_or_recv integer :: comm_size, mpi_rank - integer :: co_startdim, co_enddim + integer :: co_startdim, co_enddim, co_curr_dimen, & + co_blocksize, co_localsize logical :: comms_inited = .false., distr_inited=.false. type(MPI_Datatype) :: mpitype_column type(MPI_Datatype),dimension(:), allocatable :: mpi_blocktype @@ -211,69 +212,42 @@ subroutine co_init_distr(dimen, startdim, enddim, blocksize) integer,intent(out) :: startdim, enddim, blocksize integer :: ierr -#ifdef TROVE_USE_MPI_ - integer,dimension(:),allocatable :: starts, ends - integer :: localsize, proc_index, localsize_ + integer :: proc_index, localsize integer :: i, to_calc, ioslice_width, ioslice_maxwidth if (.not. comms_inited) stop "COMMS NOT INITIALISED" !if (distr_inited) stop "DISTRIBUTION ALREADY INITIALISED" - proc_index = mpi_rank+1 + if (dimen < comm_size) then + stop "co_init_distr: Cannot distribute matrix of dimension less than comm_size" + endif if (.not. distr_inited) then - allocate(proc_sizes(comm_size),proc_offsets(comm_size),send_or_recv(comm_size),starts(comm_size),ends(comm_size),stat=ierr) + allocate(send_or_recv(comm_size),stat=ierr) if (ierr .gt. 0) stop "CO_INIT_DISTR ALLOCATION FAILED" - else - allocate(starts(comm_size),ends(comm_size),stat=ierr) endif + co_curr_dimen = dimen ! set which dimension array we're currently distributing + if (comm_size .eq. 1) then - startdim = 1 - enddim = dimen co_startdim = 1 co_enddim = dimen - blocksize = dimen*dimen + co_blocksize = dimen*dimen + co_localsize = dimen send_or_recv(1) = 0 else - if (mpi_rank .eq. 0) then !root - - localsize = dimen/comm_size - localsize_ = int(1+real(dimen/comm_size)) - - starts(1) = 1 - ends(1) = localsize_ - proc_sizes(1) = localsize_*(comm_size*localsize_) - proc_offsets(1) = 0 - - do i=2,comm_size-1 - starts(i) = (i-1)*localsize_+1 - ends(i) = i*localsize_ - proc_sizes(i) = localsize_ * (comm_size*localsize_)!dimen - proc_offsets(i) = localsize_*(i-1)*(comm_size*localsize_)!dimen - end do - - starts(comm_size) = (i-1) * localsize_ + 1 - ends(comm_size) = dimen!comm_size*localsize_!dimen - proc_sizes(comm_size) = localsize_*comm_size*localsize_!dimen - - proc_offsets(comm_size) = (comm_size-1)*localsize_*(comm_size*localsize_)!dimen + localsize = 1+dimen/comm_size + co_startdim = mpi_rank*localsize + 1 + if (mpi_rank == comm_size-1) then + ! Last process gets full dimension + co_enddim = dimen + else + co_enddim = (mpi_rank+1)*localsize endif + co_blocksize = localsize*(comm_size*localsize) - call mpi_bcast(starts, comm_size, mpi_integer, 0, mpi_comm_world) - call mpi_bcast(ends, comm_size, mpi_integer, 0, mpi_comm_world) - call mpi_bcast(proc_sizes, comm_size, mpi_integer, 0, mpi_comm_world) - call mpi_bcast(proc_offsets, comm_size, mpi_integer, 0, mpi_comm_world) - - - - blocksize = proc_sizes(proc_index) - startdim = starts(proc_index) - enddim = ends(proc_index) - - co_startdim = startdim - co_enddim = enddim + co_localsize = localsize if(.not. distr_inited) then allocate(mpi_blocktype(comm_size)) @@ -292,47 +266,55 @@ subroutine co_init_distr(dimen, startdim, enddim, blocksize) endif endif + proc_index = mpi_rank+1 if (i.eq.proc_index) then send_or_recv(i) = 0 elseif ( ((i.gt.(proc_index - to_calc) .and. i.lt.proc_index)) .or. & ((proc_index-to_calc).lt.1 .and. (i-comm_size).gt.(proc_index-to_calc))) then send_or_recv(i) = 1 ! send - call co_create_type_subarray(int(1+real(dimen/comm_size)), blocksize, int(1+real(dimen/comm_size)), i, mpi_blocktype(i)) +#ifdef TROVE_USE_MPI_ + call co_create_type_subarray(co_localsize, co_blocksize, co_localsize, i, mpi_blocktype(i)) +#endif else send_or_recv(i) = -1 ! recv endif end do - endif - ioslice_width = enddim-startdim+1 - ioslice_maxwidth = (int(1+real(dimen/comm_size))) + startdim = co_startdim + enddim = co_enddim + blocksize = co_blocksize + + ioslice_width = co_enddim-co_startdim+1 +#ifdef TROVE_USE_MPI_ if (comm_size .eq. 1) then call co_create_type_column(dimen,dimen,dimen) else - call co_create_type_column(dimen, comm_size*ioslice_maxwidth, ioslice_width) - endif - - deallocate(starts,ends) - -#else - if (.not. comms_inited) stop "COMMS NOT INITIALISED" - if (.not. distr_inited) then - allocate(send_or_recv(1),stat=ierr) - if (ierr .gt. 0) stop "CO_INIT_DISTR ALLOCATION FAILED" + call co_create_type_column(ioslice_width, dimen, comm_size*co_localsize) endif - startdim = 1 - enddim = dimen - co_startdim = 1 - co_enddim = dimen - blocksize = dimen*dimen - send_or_recv(1) = 0 #endif distr_inited = .true. end subroutine co_init_distr + subroutine co_validate_dimensions(dimen) + integer, intent(in) :: dimen + + if (dimen .ne. co_curr_dimen) then + stop "Tried to use an array of different dimension to the current distributed setup" + endif + end subroutine + + subroutine co_create_distr_array(arr, dimen) + real(rk), pointer, intent(out) :: arr(:,:) + integer, intent(in) :: dimen + + call co_validate_dimensions(dimen) + + allocate(arr(comm_size*co_localsize, co_startdim:co_enddim)) + end subroutine + ! ! Distribute the contents of an array among processes. ! If only using one process or not using MPI, do nothing. @@ -464,7 +446,7 @@ subroutine co_write_matrix_distr(x, longdim, lb, ub, outfile) end subroutine co_write_matrix_distr #ifdef TROVE_USE_MPI_ - subroutine co_create_type_column(extent, blocksize, ncols) + subroutine co_create_type_column(ncols, extent, blocksize) integer, intent(in) :: extent, blocksize, ncols integer :: ierr,writecount diff --git a/perturbation.f90 b/perturbation.f90 index eb78c07..be212f8 100644 --- a/perturbation.f90 +++ b/perturbation.f90 @@ -15,6 +15,7 @@ module perturbation use symmetry , only : SymmetryInitialize,sym use me_numer use diag + use io_factory use io_handler_base use io_handler_ftn #ifdef TROVE_USE_MPI_ @@ -127,9 +128,7 @@ module perturbation ! type PTcontrME integer(ik) :: isize ! Number of expansion coeffs. - !real(rk),pointer :: icoeff(:,:) ! Expansion powers real(rk),pointer :: me(:,:) ! matrix elements - ! end type PTcontrME type PTcoeffT @@ -139,7 +138,6 @@ module perturbation real(rk),pointer :: coeff(:,:,:) integer(ik),pointer :: IndexQ(:,:) integer(ik),pointer :: ifromsparse(:) - !integer(ik),pointer :: itosparse(:) type(PTintcoeffs1dT),pointer :: icoeff(:) type(PTintcoeffs1dT),pointer :: iuniq(:) integer(ik),pointer :: ifield(:) @@ -155,7 +153,6 @@ module perturbation logical :: initialized = .false. type(PTcoeffs1dT),pointer :: abcissa(:) type(PTcoeffs1dT),pointer :: weight(:) - !type(PTcoeffs3dT),pointer :: deriv(:) real(rk),pointer :: poten(:) real(rk),pointer :: gvib(:,:,:) real(rk),pointer :: grot(:,:,:) @@ -170,7 +167,6 @@ module perturbation integer(ik),pointer :: nsize(:) ! number of points along each mode integer(hik) :: total_size ! total size of the dvr basis set real(ark),pointer :: drho(:) ! dvr integration step - end type PTdvrT @@ -7337,15 +7333,15 @@ subroutine PThamiltonian_contract(jrot) real(rk) :: zpe integer :: slevel,dimen_s,max_dim,iterm,jterm,total_roots,icontr,ierr ! - integer(ik) :: iunit,unitO,unitC,rec_len,irec_len,chkptIO + class(ioHandlerBase), allocatable :: kinetmatHandler + integer(ik) :: unitO,unitC,rec_len,irec_len,chkptIO integer(ik) :: ncontr,maxcontr,maxcontr0 character(len=cl) :: task character(len=4) :: jchar character(len=cl) :: unitfname,filename,statusf='old',symchar logical :: only_store = .false. logical :: no_diagonalization = .false. - !AT - type(MPI_File) :: mpiiofile + integer :: startdim,enddim,localrootsize ! ! A special case when the diagonlization is to be skipped @@ -7479,22 +7475,22 @@ subroutine PThamiltonian_contract(jrot) ! !----------Only allocate if we are putting vectors into memory---------------! if(trim(job%diagonalizer(1:13))/='READ-ENERGIES') then - matsize = int(dimen_s,hik)*int(job%nroots(isym),hik) - if (job%verbose>=4) write(out,"('Allocate array b',i7,'x',i7,' = ',i8)") dimen_s,job%nroots(isym),matsize - allocate (a(dimen_s,job%nroots(isym)),bterm(job%nroots(isym),2),stat=alloc) - ! - a = 0 - ! - call ArrayStart('PThamiltonian_contract:b',alloc,1,kind(a),matsize) - ! - bterm = 1 - ! - if (job%verbose>=4) call MemoryReport - ! - if (job%verbose>=1) then - write (out,"(//'Size of the symmetrized hamiltonian = ',i7,' Symmetry = ',a4)") dimen_s,sym%label(isym) - endif - endif + matsize = int(dimen_s,hik)*int(job%nroots(isym),hik) + if (job%verbose>=4) write(out,"('Allocate array b',i7,'x',i7,' = ',i8)") dimen_s,job%nroots(isym),matsize + allocate (a(dimen_s,job%nroots(isym)),bterm(job%nroots(isym),2),stat=alloc) + ! + a = 0 + ! + call ArrayStart('PThamiltonian_contract:b',alloc,1,kind(a),matsize) + ! + bterm = 1 + ! + if (job%verbose>=4) call MemoryReport + ! + if (job%verbose>=1) then + write (out,"(//'Size of the symmetrized hamiltonian = ',i7,' Symmetry = ',a4)") dimen_s,sym%label(isym) + endif + endif ! call diagonalization_contract(jrot,isym,dimen_s,a,zpe,rlevel,total_roots,bterm,k_row(isym,1:dimen_s)) ! @@ -7665,12 +7661,7 @@ subroutine PThamiltonian_contract(jrot) call TimerStart('Calculating the Hamiltonian matrix') ! task = 'top' - if (trim(job%kinetmat_format).eq.'MPIIO') then - call PTrestore_rot_kinetic_matrix_elements_mpi(jrot,task,mpiiofile,dimen,& - ncontr,maxcontr) - else - call PTrestore_rot_kinetic_matrix_elements(jrot,task,iunit,dimen,ncontr,maxcontr) - endif + call PTrestore_rot_kinetic_matrix_elements(jrot,task,kinetmatHandler,dimen,ncontr,maxcontr) ! ! We have two calculation options: fast and cheap and slow but expensive. ! @@ -7687,20 +7678,10 @@ subroutine PThamiltonian_contract(jrot) call TimerStart('Restoring KE matrix') ! task = 'rot' - if (trim(job%kinetmat_format).eq.'MPIIO') then - call PTrestore_rot_kinetic_matrix_elements_mpi(jrot,task,mpiiofile,dimen,& - ncontr,maxcontr) - else - call PTrestore_rot_kinetic_matrix_elements(jrot,task,iunit,dimen,ncontr,maxcontr) - endif + call PTrestore_rot_kinetic_matrix_elements(jrot,task,kinetmatHandler,dimen,ncontr,maxcontr) ! task = 'cor' - if (trim(job%kinetmat_format).eq.'MPIIO') then - call PTrestore_rot_kinetic_matrix_elements_mpi(jrot,task,mpiiofile,dimen,& - ncontr,maxcontr) - else - call PTrestore_rot_kinetic_matrix_elements(jrot,task,iunit,dimen,ncontr,maxcontr) - endif + call PTrestore_rot_kinetic_matrix_elements(jrot,task,kinetmatHandler,dimen,ncontr,maxcontr) ! call TimerStop('Restoring KE matrix') ! @@ -7715,13 +7696,8 @@ subroutine PThamiltonian_contract(jrot) call TimerStart('Restoring KE matrix') ! task = 'vib' - if (trim(job%kinetmat_format).eq.'MPIIO') then - call PTrestore_rot_kinetic_matrix_elements_mpi(jrot,task,mpiiofile,dimen,& - ncontr,maxcontr) - else - call PTrestore_rot_kinetic_matrix_elements(jrot,task,iunit,dimen,ncontr,maxcontr) - endif - + call PTrestore_rot_kinetic_matrix_elements(jrot,task,kinetmatHandler,dimen,ncontr,maxcontr) + call TimerStop('Restoring KE matrix') ! if (job%verbose>=5) write(out,"(/' ...done!')") @@ -7794,12 +7770,7 @@ subroutine PThamiltonian_contract(jrot) if (job%verbose>=4) write(out,"(/' Construct the Hamiltonian matrix...')") ! task = 'top-icontr' - if (trim(job%kinetmat_format).eq.'MPIIO') then - call PTrestore_rot_kinetic_matrix_elements_mpi(jrot,task,mpiiofile,dimen,& - ncontr,maxcontr) - else - call PTrestore_rot_kinetic_matrix_elements(jrot,task,iunit,dimen,ncontr,maxcontr) - endif + call PTrestore_rot_kinetic_matrix_elements(jrot,task,kinetmatHandler,dimen,ncontr,maxcontr) ! if (job%verbose>=5) write(out,"(' N Arrays of ',f12.5,'Gb each will be allocated (N is the number of processors)')") & real(sym%Nrepresen,rk)*real(PT%max_deg_size,rk)*real(max_dim,rk)*real(rk,rk)/1024.0_rk**3 @@ -7819,32 +7790,17 @@ subroutine PThamiltonian_contract(jrot) if (FLrotation.and.jrot/=0) then ! task = 'rot-icontr' - if (trim(job%kinetmat_format).eq.'MPIIO') then - call PTrestore_rot_kinetic_matrix_elements_mpi(jrot,task,mpiiofile,dimen,& - ncontr,maxcontr) - else - call PTrestore_rot_kinetic_matrix_elements(jrot,task,iunit,dimen,ncontr,maxcontr,icontr) - endif + call PTrestore_rot_kinetic_matrix_elements(jrot,task,kinetmatHandler,dimen,ncontr,maxcontr,icontr) ! task = 'cor-icontr' - if (trim(job%kinetmat_format).eq.'MPIIO') then - call PTrestore_rot_kinetic_matrix_elements_mpi(jrot,task,mpiiofile,dimen,& - ncontr,maxcontr) - else - call PTrestore_rot_kinetic_matrix_elements(jrot,task,iunit,dimen,ncontr,maxcontr,icontr) - endif + call PTrestore_rot_kinetic_matrix_elements(jrot,task,kinetmatHandler,dimen,ncontr,maxcontr,icontr) ! endif ! if ( PTvibrational_me_calc ) then ! task = 'vib-icontr' - if (trim(job%kinetmat_format).eq.'MPIIO') then - call PTrestore_rot_kinetic_matrix_elements_mpi(jrot,task,mpiiofile,dimen,& - ncontr,maxcontr) - else - call PTrestore_rot_kinetic_matrix_elements(jrot,task,iunit,dimen,ncontr,maxcontr,icontr) - endif + call PTrestore_rot_kinetic_matrix_elements(jrot,task,kinetmatHandler,dimen,ncontr,maxcontr,icontr) ! endif ! @@ -7915,7 +7871,6 @@ subroutine PThamiltonian_contract(jrot) if (FLrotation.and.jrot/=0.and..not.job%vib_rot_contr) then ! if (job%verbose>=4) write(out,"(' Rotational part...')") - ! ! AT - Initialise parallel distribution (if not done yet) ! TODO - find correct spot to place this w/ appropriate guard clauses ! NOTE don't know `ncontr` until after calling PTrestore.. with `task=='top'` @@ -7923,12 +7878,7 @@ subroutine PThamiltonian_contract(jrot) ! task = 'rot' ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call PTrestore_rot_kinetic_matrix_elements_mpi(jrot,task,mpiiofile,dimen,& - ncontr,maxcontr) - else - call PTrestore_rot_kinetic_matrix_elements(jrot,task,iunit,dimen,ncontr,maxcontr) - endif + call PTrestore_rot_kinetic_matrix_elements(jrot,task,kinetmatHandler,dimen,ncontr,maxcontr) ! ! $omp parallel private(mat_t,alloc_p) allocate (mat_t(sym%Nrepresen,PT%max_deg_size,max_dim),stat=alloc_p) @@ -7978,18 +7928,12 @@ subroutine PThamiltonian_contract(jrot) enddo ! deallocate(grot) - call Arraystop('grot-matrix') ! if (job%verbose>=4) write(out,"(' Coriolis part...')") ! task = 'cor' ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call PTrestore_rot_kinetic_matrix_elements_mpi(jrot,task,mpiiofile,dimen,& - ncontr,maxcontr) - else - call PTrestore_rot_kinetic_matrix_elements(jrot,task,iunit,dimen,ncontr,maxcontr) - endif + call PTrestore_rot_kinetic_matrix_elements(jrot,task,kinetmatHandler,dimen,ncontr,maxcontr) ! ! $omp parallel private(mat_t,alloc_p) allocate (mat_t(sym%Nrepresen,PT%max_deg_size,max_dim),stat=alloc_p) @@ -8039,7 +7983,6 @@ subroutine PThamiltonian_contract(jrot) enddo ! deallocate(gcor) - call Arraystop('gcor-matrix') ! endif ! @@ -8051,12 +7994,7 @@ subroutine PThamiltonian_contract(jrot) ! task = 'vib' ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call PTrestore_rot_kinetic_matrix_elements_mpi(jrot,task,mpiiofile,dimen,& - ncontr,maxcontr) - else - call PTrestore_rot_kinetic_matrix_elements(jrot,task,iunit,dimen,ncontr,maxcontr) - endif + call PTrestore_rot_kinetic_matrix_elements(jrot,task,kinetmatHandler,dimen,ncontr,maxcontr) ! ! $omp parallel private(mat_t,alloc_p) allocate (mat_t(sym%Nrepresen,PT%max_deg_size,max_dim),stat=alloc_p) @@ -8121,7 +8059,7 @@ subroutine PThamiltonian_contract(jrot) call TimerStop('Calculating the Hamiltonian matrix') ! if (job%verbose>=4) write(out,"('...done!')") - if (mpi_rank.eq.0) then!mpiio + if (mpi_rank.eq.0) then do isym = 1,sym%Nrepresen if (.not.job%select_gamma(isym)) cycle enddo @@ -8461,7 +8399,7 @@ subroutine PThamiltonian_contract(jrot) endif ! enddo - endif!mpiio + endif ! if ( trim(job%IOeigen_action)=='SAVE'.or.trim(job%IOeigen_action)=='APPEND' ) then ! @@ -8496,182 +8434,91 @@ subroutine PThamiltonian_contract(jrot) end subroutine PThamiltonian_contract ! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!! MPIIO !!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine open_chkptfile_mpi(fileh, filename, mode) -#ifdef TROVE_USE_MPI_ - use mpi_f08 -#endif - use mpi_aux - - type(MPI_File),intent(inout) :: fileh - character(len=*),intent(in) :: filename, mode - -#ifdef TROVE_USE_MPI_ - integer :: amode, ierr - - select case(mode) - case('read') - amode = mpi_mode_rdonly - case('write') - amode = mpi_mode_wronly+mpi_mode_create - end select - - call MPI_File_open(mpi_comm_world, filename, amode, mpi_info_null, fileh, ierr) - if (ierr.gt.0) then - if (mpi_rank .eq. 0) write(*,*) "Error opening MPI-IO-formatted Vib. kinetic checkpoint file. ", filename - stop "MPI_PTrestore_rot_kinetic_matrix_elements - Error opening MATELEM MPI-IO input file" - endif - - !File errors indicate big trouble, so we set errors to be fatal - not likely to be recoverable - call MPI_File_set_errhandler(fileh, MPI_ERRORS_ARE_FATAL) -#endif - end subroutine open_chkptfile_mpi - - subroutine close_chkptfile_mpi(fileh) -#ifdef TROVE_USE_MPI_ - use mpi_f08 -#endif - use mpi_aux - - type(MPI_File), intent(inout) :: fileh -#ifdef TROVE_USE_MPI_ - integer :: ierr - - call mpi_file_close(fileh, ierr) - if (ierr.gt.0) then - if (mpi_rank .eq. 0) write(*,*) "Error closing MPI-IO-formatted Vib. kinetic checkpoint file." - stop "MPI_PTrestore_rot_kinetic_matrix_elements - Error closing MATELEM MPI-IO input file" - endif -#endif - end subroutine close_chkptfile_mpi - - subroutine PTrestore_rot_kinetic_matrix_elements_mpi(jrot, task, fileh, dimen, & + subroutine PTrestore_rot_kinetic_matrix_elements(jrot, task, ioHandler, dimen, & ncontr, maxcontr, icontr) -#ifdef TROVE_USE_MPI_ - use mpi_f08 -#endif use mpi_aux integer(ik),intent(in) :: jrot character(len=cl),intent(in) :: task - type(MPI_File),intent(inout) :: fileh + class(ioHandlerBase), allocatable, intent(inout) :: ioHandler integer(ik),intent(in) :: dimen integer(ik),intent(inout),optional :: ncontr integer(ik),intent(inout),optional :: maxcontr integer(ik),intent(in),optional :: icontr -#ifdef TROVE_USE_MPI_ - type(MPI_File) :: fileh_slice + type(ErrorType) :: err + + class(ioHandlerBase), allocatable :: sliceHandler character(len=cl) :: job_id,filename,readbuf - integer(kind=MPI_Offset_kind) :: file_offset + integer :: file_offset integer :: ierr - integer(hik) :: rootsize,rootsize_,rootsize2,rootsize2_,nprocs,tid,icontr1,icontr2 + integer(hik) :: nprocs,tid,icontr1,icontr2 + integer :: OMP_GET_NUM_THREADS,omp_get_thread_num integer(ik),allocatable :: imat_t(:,:) real(rk),allocatable :: mat_t(:,:),mat_(:,:) real(rk) :: f_t + integer(ik) :: chkptIO, chkptIO_ + integer(ik) :: maxcontr0 + integer(hik) :: rootsize2_ + integer :: k1,k2,islice - !AT - integer :: localmatrix_x,localmatrix_y - integer(hik) :: localrootsize if ( trim(job%IOkinet_action)/='VIB_READ'.and.trim(job%IOkinet_action)/='READ'.and..not.job%contrci_me_fast) return if ( trim(job%IOswap_matelem)/='NONE') return - rootsize = int(ncontr*(ncontr+1)/2,hik) - rootsize_ = int(maxcontr*(maxcontr+1)/2,hik) - ! - rootsize2 = int(ncontr*ncontr,hik) - rootsize2_ = int(maxcontr*maxcontr,hik) - ! - !dimen = max(min(int(PT%Maxcontracts*job%compress),PT%Maxcontracts),1) - if (mpi_rank .eq. 0 .and. (job%verbose>=6.and.present(icontr)) ) write(out,"('icontr = ',i9)") icontr - !AT - determine task-local matrix dimensions - if (comm_size.eq.1) then - localmatrix_y = ncontr - else - localmatrix_y = int(1+real(ncontr/comm_size))*comm_size - endif - localmatrix_x = co_enddim-co_startdim+1 - localrootsize = int(localmatrix_x*localmatrix_y,hik) - select case (trim(task)) case('top') - job_id = '[MPI-IO] Vib. matrix elements of the rot. kinetic' + job_id = 'Vib. matrix elements of the rot. kinetic' !TODO - MPI-compatible IOStart !call IOStart(trim(job_id),fileh) - call open_chkptfile_mpi(fileh, job%kinetmat_file, 'read') + call openFile(ioHandler, job%kinetmat_file, err, action='read', & + position='rewind', status='old', form='unformatted') + HANDLE_ERROR(err) - !Collective read woo - call MPI_File_read_all(fileh, readbuf, 7, mpi_character, mpi_status_ignore, ierr) - - if (readbuf(1:7) /= '[MPIIO]') then - if (mpi_rank .eq. 0) write(*,*) "Invalid MPIIO identifier to MPI-IO-formatted Vib. kinetic checkpoint file." - call mpi_barrier(mpi_comm_world, ierr) - stop "MPI_PTrestore_rot_kinetic_matrix_elements - bogus file format header" - endif - - call MPI_File_read_all(fileh, readbuf, 18, mpi_character, mpi_status_ignore, ierr) + call ioHandler%read(readbuf(1:18)) if (readbuf(1:18) /= 'Start Kinetic part') then - if (mpi_rank .eq. 0) write(*,*) "Invalid header to MPI-IO-formatted Vib. kinetic checkpoint file." - call mpi_barrier(mpi_comm_world, ierr) - stop "MPI_PTrestore_rot_kinetic_matrix_elements - bogus file format" + if (mpi_rank .eq. 0) write (out,"(' Vib. kinetic checkpoint file ',a,' has bogus header: ',a)") & + job%kinetmat_file,readbuf(1:18) + stop 'PTrestore_rot_kinetic_matrix_elements - bogus file format' endif - call MPI_File_read_all(fileh, ncontr, 1, mpi_integer, mpi_status_ignore, ierr) + call ioHandler%read(ncontr) if (jrot==0.and.PT%Maxcontracts/=ncontr) then if (mpi_rank .eq. 0) write (out,"(' Vib. kinetic checkpoint file ',a)") job%kinetmat_file if (mpi_rank .eq. 0) write (out,"(' Actual and stored basis sizes at J=0 do not agree ',2i9)") PT%Maxcontracts,ncontr - call mpi_barrier(mpi_comm_world, ierr) - stop 'MPI_PTrestore_rot_kinetic_matrix_elements - in file - illegal nroots ' + stop 'PTrestore_rot_kinetic_matrix_elements - in file - illegal nroots ' end if - rootsize = int(ncontr*(ncontr+1)/2,hik) - rootsize2 = int(ncontr*ncontr,hik) - if (mpi_rank .eq. 0.and.job%verbose>=6) write(out,"(/'Restore_rot_kin...: Number of elements: ',i8)") PT%Maxcontracts ! Read the indexes of the J=0 contracted basis set. if (.not.FLrotation.or.jrot==0) then - - !allocate (imat_t(0:PT%Nclasses,ncontr),stat=ierr) - !call ArrayStart('mat_t',ierr,size(imat_t),kind(imat_t)) - - call MPI_File_read_all(fileh, readbuf, 10, mpi_character, mpi_status_ignore, ierr) + call ioHandler%read(readbuf(1:10)) if (readbuf(1:10)/='icontr_cnu') then if (mpi_rank .eq. 0) write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_cnu is missing ',a)") job%kinetmat_file,readbuf(1:10) - call mpi_barrier(mpi_comm_world, ierr) - stop 'MPI_PTrestore_rot_kinetic_matrix_elements - in file - icontr_cnu missing' + stop 'PTrestore_rot_kinetic_matrix_elements - in file - icontr_cnu missing' end if - file_offset = (PT%Nclasses+1)*int(ncontr,MPI_OFFSET_KIND)*mpi_int_size - call MPI_File_seek(fileh, file_offset, MPI_SEEK_CUR) + file_offset = (PT%Nclasses+1)*int(ncontr,MPI_OFFSET_KIND)*sizeof(0) + call ioHandler%seek(file_offset) - call MPI_File_read_all(fileh, readbuf, 11, mpi_character, mpi_status_ignore, ierr) + call ioHandler%read(readbuf(1:11)) if (readbuf(1:11)/='icontr_ideg') then if (mpi_rank .eq. 0) write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_ideg is missing ',a)") job%kinetmat_file,readbuf(1:11) - call mpi_barrier(mpi_comm_world, ierr) - stop 'MPI_PTrestore_rot_kinetic_matrix_elements - in file - icontr_ideg missing' + stop 'PTrestore_rot_kinetic_matrix_elements - in file - icontr_ideg missing' end if - file_offset = (PT%Nclasses+1)*int(ncontr,MPI_OFFSET_KIND)*mpi_int_size - call MPI_File_seek(fileh, file_offset, MPI_SEEK_CUR) - file_offset = file_offset+file_offset+mpi_int_size+46 - call MPI_File_seek_shared(fileh, file_offset, MPI_SEEK_SET) - - !deallocate(imat_t) - + file_offset = (PT%Nclasses+1)*int(ncontr,MPI_OFFSET_KIND)*sizeof(0) + call ioHandler%seek(file_offset) else allocate (PT%icontr_cnu(0:PT%Nclasses,ncontr),PT%icontr_ideg(0:PT%Nclasses,ncontr),stat=ierr) @@ -8684,36 +8531,33 @@ subroutine PTrestore_rot_kinetic_matrix_elements_mpi(jrot, task, fileh, dimen, & allocate (PT%icontr2icase(ncontr,2),stat=ierr) call ArrayStart('PT%contractive_space',ierr,size(PT%icontr2icase),kind(PT%icontr2icase)) - call MPI_File_read_all(fileh, readbuf, 10, mpi_character, mpi_status_ignore, ierr) + call ioHandler%read(readbuf(1:10)) if (readbuf(1:10)/='icontr_cnu') then if (mpi_rank .eq. 0) write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_cnu is missing ',a)") job%kinetmat_file,readbuf(1:10) - call mpi_barrier(mpi_comm_world, ierr) - stop 'MPI_PTrestore_rot_kinetic_matrix_elements - in file - icontr_cnu missing' + stop 'PTrestore_rot_kinetic_matrix_elements - in file - icontr_cnu missing' end if - call MPI_File_read_all(fileh, PT%icontr_cnu, (PT%Nclasses+1)*ncontr, mpi_integer, mpi_status_ignore, ierr) + call ioHandler%read(PT%icontr_cnu) - call MPI_File_read_all(fileh, readbuf, 11, mpi_character, mpi_status_ignore, ierr) + call ioHandler%read(readbuf(1:11)) if (readbuf(1:11)/='icontr_ideg') then if (mpi_rank .eq. 0) write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_ideg is missing ',a)") job%kinetmat_file,readbuf(1:11) - call mpi_barrier(mpi_comm_world, ierr) - stop 'MPI_PTrestore_rot_kinetic_matrix_elements - in file - icontr_ideg missing' + stop 'PTrestore_rot_kinetic_matrix_elements - in file - icontr_ideg missing' end if - call MPI_File_read_all(fileh, PT%icontr_ideg, (PT%Nclasses+1)*ncontr, mpi_integer, mpi_status_ignore, ierr) + call ioHandler%read(PT%icontr_ideg) endif if (job%vib_rot_contr.and.PTvibrational_me_calc) then - call MPI_File_read_all(fileh, readbuf, 7, mpi_character, mpi_status_ignore, ierr) + call ioHandler%read(readbuf(1:7)) if (readbuf(1:7)/='vib-rot') then if (mpi_rank .eq. 0) write (out,"(' Vib. kinetic checkpoint file ',a,': label vib-rot is missing ',a)") job%kinetmat_file,readbuf(1:7) - call mpi_barrier(mpi_comm_world, ierr) - stop 'MPI_PTrestore_rot_kinetic_matrix_elements - in file - vib-rot missing' + stop 'PTrestore_rot_kinetic_matrix_elements - in file - vib-rot missing' end if - call close_chkptfile_mpi(fileh) + deallocate(ioHandler) write(job_id,"('single swap_matrix')") @@ -8722,7 +8566,9 @@ subroutine PTrestore_rot_kinetic_matrix_elements_mpi(jrot, task, fileh, dimen, & filename = trim(job%matelem_suffix)//"0"//'.chk' - call open_chkptfile_mpi(fileh, filename, 'read') + call openFile(ioHandler, filename, err, action='read', & + status='old', form='unformatted') + HANDLE_ERROR(err) endif @@ -8738,8 +8584,7 @@ subroutine PTrestore_rot_kinetic_matrix_elements_mpi(jrot, task, fileh, dimen, & if (maxcontr>ncontr) then if (mpi_rank .eq. 0) write (out,"(' Actual and stored basis sizes at J=0 do not agree (maxcontr,ncontr) ',2i8)") maxcontr,ncontr - call mpi_barrier(mpi_comm_world, ierr) - stop 'MPI_PTrestore_rot_kinetic_matrix_elements - in file - illegal ncontr ' + stop 'PTrestore_rot_kinetic_matrix_elements - in file - illegal ncontr ' end if case('rot') ! @@ -8747,33 +8592,20 @@ subroutine PTrestore_rot_kinetic_matrix_elements_mpi(jrot, task, fileh, dimen, & ! ! Read the rotational part only ! - !allocate(mat_(maxcontr,maxcontr),stat=ierr) - !call ArrayStart('PThamiltonian_contract: mat_',ierr,1,kind(mat_),rootsize2_) - ! if (.not.job%IOmatelem_split) then ! - call MPI_File_read_all(fileh, readbuf, 5, mpi_character, mpi_status_ignore, ierr) + call ioHandler%read(readbuf(1:5)) if (readbuf(1:5)/='g_rot') then if(mpi_rank .eq. 0) write (out,"(' Vib. kinetic checkpoint file ',a,': g_rot is missing ',a)") job%kinetmat_file,readbuf(1:5) - call mpi_barrier(mpi_comm_world, ierr) - stop 'MPI_PTrestore_rot_kinetic_matrix_elements - in file - g_rot missing' - end if + stop 'PTrestore_rot_kinetic_matrix_elements - in file - g_rot missing' + endif ! allocate(grot(3,3),stat=ierr) ! - islice = 0 - ! do k1 = 1,3 - ! do k2 = 1,3 - ! - islice = islice + 1 - ! - allocate(grot(k1,k2)%me(localmatrix_y,co_startdim:co_enddim),stat=ierr) - call ArrayStart('grot-matrix',ierr,1,kind(f_t),localrootsize) - ! - call co_read_matrix_distr(grot(k1,k2)%me, ncontr, co_startdim, co_enddim, fileh) - ! + call co_create_distr_array(grot(k1,k2)%me, ncontr) + call ioHandler%read(grot(k1,k2)%me, ncontr) enddo enddo ! @@ -8784,25 +8616,19 @@ subroutine PTrestore_rot_kinetic_matrix_elements_mpi(jrot, task, fileh, dimen, & islice = 0 ! do k1 = 1,3 - ! do k2 = 1,3 ! islice = islice + 1 ! - call divided_slice_open_mpi(islice,fileh_slice,'g_rot',job%matelem_suffix) - ! - allocate(grot(k1,k2)%me(localmatrix_y,co_startdim:co_enddim),stat=ierr) - call ArrayStart('grot-matrix',ierr,1,kind(f_t),localrootsize) + call co_create_distr_array(grot(k1,k2)%me, ncontr) ! - call co_read_matrix_distr(grot(k1,k2)%me, ncontr, co_startdim, co_enddim, fileh_slice) - ! - call divided_slice_close_mpi(islice,fileh_slice,'g_rot') + call divided_slice_open(islice,sliceHandler,'g_rot',job%matelem_suffix) + call sliceHandler%read(grot(k1,k2)%me, ncontr) + call divided_slice_close(islice,sliceHandler,'g_rot') ! enddo enddo endif - !deallocate(mat_) - !call ArrayStop('PThamiltonian_contract: mat_') ! if (job%verbose>=4) write(out,"(' ...done!')") ! @@ -8810,16 +8636,12 @@ subroutine PTrestore_rot_kinetic_matrix_elements_mpi(jrot, task, fileh, dimen, & ! if (mpi_rank .eq. 0 .and. job%verbose>=4) write(out,"(' Read and process Coriolis part...')") ! - !allocate(mat_(maxcontr,maxcontr),stat=ierr) - !call ArrayStart('PThamiltonian_contract: mat_',ierr,1,kind(mat_),rootsize2_) - ! if (.not.job%IOmatelem_split) then ! - call MPI_File_read_all(fileh, readbuf, 5, mpi_character, mpi_status_ignore, ierr) + call ioHandler%read(readbuf(1:5)) if (readbuf(1:5)/='g_cor') then if (mpi_rank .eq. 0) write (out,"(' Vib. kinetic checkpoint file ',a,': g_cor is missing ',a)") job%kinetmat_file,readbuf(1:5) - call mpi_barrier(mpi_comm_world, ierr) - stop 'MPI_PTrestore_rot_kinetic_matrix_elements - in file - g_cor missing' + stop 'PTrestore_rot_kinetic_matrix_elements - in file - g_cor missing' end if ! allocate(gcor(3),stat=ierr) @@ -8829,11 +8651,9 @@ subroutine PTrestore_rot_kinetic_matrix_elements_mpi(jrot, task, fileh, dimen, & do k1 = 1,3 ! islice = islice + 1 - ! - allocate(gcor(k1)%me(localmatrix_y,co_startdim:co_enddim),stat=ierr) - call ArrayStart('gcor-matrix',ierr,1,kind(f_t),localrootsize) - ! - call co_read_matrix_distr(gcor(k1)%me, ncontr, co_startdim, co_enddim, fileh) + + call co_create_distr_array(gcor(k1)%me, ncontr) + call ioHandler%read(gcor(k1)%me, ncontr) ! enddo ! @@ -8847,22 +8667,15 @@ subroutine PTrestore_rot_kinetic_matrix_elements_mpi(jrot, task, fileh, dimen, & ! islice = islice + 1 ! - allocate(gcor(k1)%me(localmatrix_y,co_startdim:co_enddim),stat=ierr) - call ArrayStart('gcor-matrix',ierr,1,kind(f_t),localrootsize) - ! - call divided_slice_open_mpi(islice,fileh_slice,'g_cor',job%matelem_suffix) - ! - call co_read_matrix_distr(gcor(k1)%me, ncontr, co_startdim, co_enddim, fileh_slice) - ! - call divided_slice_close_mpi(islice,fileh_slice,'g_cor') + call co_create_distr_array(gcor(k1)%me, ncontr) ! + call divided_slice_open(islice,sliceHandler,'g_cor',job%matelem_suffix) + call sliceHandler%read(gcor(k1)%me, ncontr) + call divided_slice_close(islice,sliceHandler,'g_cor') enddo ! endif ! - !deallocate(mat_) - !call ArrayStop('PThamiltonian_contract: mat_') - ! if (mpi_rank .eq. 0 .and. job%verbose>=4) write(out,"(' ...done!')") ! case('vib') @@ -8871,652 +8684,56 @@ subroutine PTrestore_rot_kinetic_matrix_elements_mpi(jrot, task, fileh, dimen, & ! if (.not.job%IOmatelem_split.and.( (.not.FLrotation.or.jrot==0).and.trim(job%IOkinet_action)/='VIB_READ' ) ) then ! - !allocate(mat_t(maxcontr,maxcontr),stat=ierr) - !call ArrayStart('PThamiltonian_contract: mat_t',ierr,1,kind(mat_t),rootsize2_) - ! - call MPI_File_read_all(fileh, readbuf, 5, mpi_character, mpi_status_ignore, ierr) + call ioHandler%read(readbuf(1:5)) if (readbuf(1:5)/='g_rot') then if (mpi_rank .eq. 0) write (out,"(' Vib. kinetic checkpoint file ',a,': g_rot is missing ',a)") job%kinetmat_file,readbuf(1:5) - call mpi_barrier(mpi_comm_world, ierr) - stop 'MPI_PTrestore_rot_kinetic_matrix_elements - in file - g_rot missing' + stop 'PTrestore_rot_kinetic_matrix_elements - in file - g_rot missing' end if ! - file_offset = 9*int(ncontr,MPI_OFFSET_KIND)*ncontr*mpi_real_size - call MPI_File_seek(fileh, file_offset, MPI_SEEK_CUR) + file_offset = 9*int(ncontr,MPI_OFFSET_KIND)*ncontr*sizeof(0.0_rk) + call ioHandler%seek(file_offset) ! - call MPI_File_read_all(fileh, readbuf, 5, mpi_character, mpi_status_ignore, ierr) + call ioHandler%read(readbuf(1:5)) if (readbuf(1:5)/='g_cor') then if (mpi_rank .eq. 0) write (out,"(' Vib. kinetic checkpoint file ',a,': g_cor is missing ',a)") job%kinetmat_file,readbuf(1:5) - call mpi_barrier(mpi_comm_world, ierr) - stop 'MPI_PTrestore_rot_kinetic_matrix_elements - in file - g_cor missing' + stop 'PTrestore_rot_kinetic_matrix_elements - in file - g_cor missing' end if ! - file_offset = 3*int(ncontr,MPI_OFFSET_KIND)*ncontr*mpi_real_size - call MPI_File_seek(fileh, file_offset, MPI_SEEK_CUR) - ! - !deallocate(mat_t) - !call ArrayStop('mat_t') + file_offset = 3*int(ncontr,MPI_OFFSET_KIND)*ncontr*sizeof(0.0_rk) + call ioHandler%seek(file_offset) ! endif ! - if (mpi_rank .eq. 0 .and. job%verbose>=6) write(out,"(' rootsize_,rootsize = ',2i9)") rootsize_,rootsize - ! - call MPI_File_read_all(fileh, readbuf, 4, mpi_character, mpi_status_ignore, ierr) + call ioHandler%read(readbuf(1:4)) if (readbuf(1:4)/='hvib') then if (mpi_rank .eq. 0) write (out,"(' Vib. kinetic checkpoint file ',a,': hvib is missing ',a)") job%kinetmat_file,readbuf(1:4) if (mpi_rank .eq. 0 .and. readbuf(1:4)=='g_ro') then write (out,"(' Most likely non-divided chk-points used with MATELEM READ SPLIT')") write (out,"(' Re-do MATELEM SAVE SPLIT or use MATELEM SPLIT READ')") endif - call mpi_barrier(mpi_comm_world, ierr) - stop 'MPI_PTrestore_rot_kinetic_matrix_elements - in file - hvib or End missing' + stop 'PTrestore_rot_kinetic_matrix_elements - in file - hvib or End missing' end if + + call co_create_distr_array(hvib%me, ncontr) ! - !allocate(mat_(maxcontr,maxcontr),stat=ierr) - !call ArrayStart('PThamiltonian_contract: mat_',ierr,1,kind(mat_),rootsize2_) - ! - allocate(hvib%me(localmatrix_y,co_startdim:co_enddim),stat=ierr) - call ArrayStart('hvib-matrix',ierr,1,kind(f_t),localrootsize) - ! - call co_read_matrix_distr(hvib%me, ncontr, co_startdim, co_enddim, fileh) - !file_offset = int(ncontr,MPI_OFFSET_KIND)*ncontr*mpi_real_size - !call MPI_File_seek(fileh, file_offset, MPI_SEEK_CUR) - !call MPI_File_read_all(fileh, mat_(1:ncontr,1:ncontr), ncontr*ncontr, mpi_double_precision, mpi_status_ignore, ierr) - ! - !hvib%me(1:ncontr,1:ncontr) = mat_(1:ncontr,1:ncontr) - ! - !deallocate(mat_) - !call ArrayStop('PThamiltonian_contract: mat_') - ! - call MPI_File_read_all(fileh, readbuf, 16, mpi_character, mpi_status_ignore, ierr) + call ioHandler%read(hvib%me, ncontr) + + call ioHandler%read(readbuf(1:16)) if (readbuf(1:16)/='End Kinetic part') then if (mpi_rank .eq. 0) write (out,"(' Vib. kinetic checkpoint file ',a,' has bogus footer: ',a)") job%kinetmat_file,readbuf(1:16) - call mpi_barrier(mpi_comm_world, ierr) - stop 'MPI_PTrestore_rot_kinetic_matrix_elements - bogus file format' + stop 'PTrestore_rot_kinetic_matrix_elements - bogus file format' end if ! - call close_chkptfile_mpi(fileh) + deallocate(ioHandler) ! if (mpi_rank .eq. 0 .and. job%verbose>=4) write(out,"(' ...done!')") - ! - case('top-icontr') - write(*,*) "CASE TOP-ICONTR" - !! ! - !! nprocs = 1 - !! tid = 0 - !! !$omp parallel private(tid) - !! tid = omp_get_thread_num() - !! if (tid==0) then - !! nprocs = omp_get_num_threads() - !! endif - !! !$omp end parallel - !! ! - !! if (FLrotation.and.jrot/=0) then - !! ! - !! allocate(grot(3,3),stat=ierr) - !! ! - !! rootsize2_ = int(maxcontr,hik)*int(PT%max_deg_size,hik)*9_hik*int(nprocs,hik) - !! call ArrayStart('PThamiltonian_contract: grot',ierr,1,rk,rootsize2_) - !! ! - !! allocate(gcor(3),stat=ierr) - !! rootsize2_ = int(maxcontr,hik)*int(PT%max_deg_size,hik)*int(PT%Nmodes,hik)*3_hik*int(nprocs,hik) - !! call ArrayStart('PThamiltonian_contract: grot',ierr,1,rk,rootsize2_) - !! ! - !! if (job%vib_rot_contr) then - !! ! - !! do islice = 1,9+3*PT%Nmodes - !! call divided_slice_open_vib_rot(islice,job%matelem_suffix) - !! enddo - !! ! - !! endif - !! ! - !! endif - !! ! - case('rot-icontr') ! rotational part for the vib-rot contraction scheme - write(*,*) "CASE ROT-ICONTR" - !! ! - !! ! Read the rotational part only - !! ! - !! if (.not.job%IOmatelem_split ) then - !! ! - !! write (out,"('PTrestore_rot_kinetic_matrix_elements: vib-rot can be used with MATELEM SAVE SPLIT only ')") - !! stop 'PTrestore_rot_kinetic_matrix_elements: vib-rot can be used with SPLIT only' - !! ! - !! endif - !! ! - !! if (.not.job%vib_rot_contr) return - !! ! - !! icontr1 = PT%Ncontr02icase0(icontr,1) - !! icontr2 = PT%Ncontr02icase0(icontr,2) - !! ! - !! islice = 0 - !! ! - !! do k1 = 1,3 - !! do k2 = 1,3 - !! ! - !! islice = islice + 1 - !! ! - !! nullify(grot(k1,k2)%me) - !! ! - !! if (associated(grot(k1,k2)%me)) deallocate(grot(k1,k2)%me) - !! ! - !! allocate(grot(k1,k2)%me(maxcontr,icontr1:icontr2),stat=ierr) - !! ! - !! write(job_is,"('single swap_matrix #',i8)") islice - !! call IOStart(trim(job_is),chkptIO_) - !! ! - !! read(chkptIO_) grot(k1,k2)%me - !! ! - !! enddo - !! enddo - !! ! - case('cor-icontr') ! corriolis part for the vib-rot contraction scheme - write(*,*) "CASE COR-ICONTR" - !! ! - !! ! Read the Corriolis part only - !! ! - !! ! - !! if (.not.job%IOmatelem_split ) then - !! ! - !! write (out,"('PTrestore_rot_kinetic_matrix_elements: vib-rot can be used with MATELEM SAVE SPLIT only ')") - !! stop 'PTrestore_rot_kinetic_matrix_elements: vib-rot can be used with SPLIT only' - !! ! - !! endif - !! ! - !! if (.not.job%vib_rot_contr) then - !! ! - !! write (out,"('PTrestore_rot_kinetic_matrix_elements: vib-rot can be used with cor-icontr only ')") - !! stop 'PTrestore_rot_kinetic_matrix_elements: vib-rot can be used with cor-icontr only' - !! ! - !! endif - !! ! - !! icontr1 = PT%Ncontr02icase0(icontr,1) - !! icontr2 = PT%Ncontr02icase0(icontr,2) - !! ! - !! islice = 9 - !! ! - !! do k1 = 1,3 - !! ! - !! islice = islice + 1 - !! ! - !! nullify(gcor(k1)%me) - !! ! - !! if (associated(gcor(k1)%me)) deallocate(gcor(k1)%me) - !! ! - !! allocate(gcor(k1)%me(maxcontr,icontr1:icontr2),stat=ierr) - !! ! - !! write(job_is,"('single swap_matrix #',i8)") islice - !! call IOStart(trim(job_is),chkptIO_) - !! ! - !! read(chkptIO_) gcor(k1)%me - !! ! - !! enddo - !! ! - case('vib-icontr') ! vibrational part for the vib-rot contraction scheme - write(*,*) "CASE VIB-ICONTR" - !! ! - !! ! - !! !if (job%verbose>=4.and.irow==0) write(out,"(' Read and process vibrational part...')") - !! ! - !! if (.not.job%IOmatelem_split.and.( (.not.FLrotation.or.jrot==0).and.trim(job%IOkinet_action)/='VIB_READ' ) ) then - !! ! - !! if (job%vib_rot_contr) then - !! write (out,"('PTrestore_rot_kinetic_matrix_elements: vib-rot can be used with MATELEM SAVE SPLIT only ')") - !! stop 'PTrestore_rot_kinetic_matrix_elements: vib-rot can be used with SPLIT only' - !! endif - !! ! - !! endif - !! ! - !! if (icontr==1) then - !! ! - !! read(chkptIO) readbuf(1:4) - !! if (readbuf(1:4)/='hvib') then - !! write (out,"(' Vib. kinetic checkpoint file ',a,': hvib is missing ',a)") job%kinetmat_file,readbuf(1:4) - !! if (readbuf(1:4)=='g_ro') write (out,"(' Most likely non-divided chk-points used with MATELEM READ DIVIDED')") - !! write (out,"(' Re-do MATELEM SAVE DIVIDE or use MATELEM READ!')") - !! stop 'PTrestore_rot_kinetic_matrix_elements - in file - hvib or End missing' - !! end if - !! ! - !! endif - !! ! - !! if (associated(hvib%me)) deallocate(hvib%me) - !! ! - !! call ArrayStart('hvib-matrix',0,1,4) - !! call ArrayStart('grot-matrix',0,1,4) - !! call ArrayStart('gcor-matrix',0,1,4) - !! call ArrayStop('hvib-matrix') - !! call ArrayStop('grot-matrix') - !! call ArrayStop('gcor-matrix') - !! ! - !! icontr1 = PT%Ncontr02icase0(icontr,1) - !! icontr2 = PT%Ncontr02icase0(icontr,2) - !! ! - !! if (job%verbose>=6) write(out,"('allocate hvib for ',i9,' x ',i8,' -> ',i8)") maxcontr,icontr1,icontr2 - !! ! - !! allocate(hvib%me(maxcontr,icontr1:icontr2),stat=ierr) - !! call ArrayStart('hvib-matrix',ierr,1,kind(f_t),rootsize2_) - !! ! - !! read(chkptIO) hvib%me - !! ! - !! maxcontr0 = size(PT%Ncontr02icase0,dim=1) - !! ! - !! if (icontr==maxcontr0) then - !! ! - !! read(chkptIO) readbuf(1:4) - !! if (readbuf(1:4)/='hvib') then - !! write (out,"(' Vib. kinetic checkpoint file ',a,' has bogus footer: ',a)") job%kinetmat_file,readbuf(1:16) - !! stop 'PTrestore_rot_kinetic_matrix_elements - bogus file format' - !! end if - !! ! - !! close(chkptIO,status='keep') - !! ! - !! if (job%verbose>=4) write(out,"(' ...done!')") - !! ! - !! endif - !! ! - end select -#endif - end subroutine - - subroutine divided_slice_open_mpi(islice,fileh,chkpt_type,suffix) -#ifdef TROVE_USE_MPI_ - use mpi_f08 -#endif - use mpi_aux - implicit none - integer(ik),intent(in) :: islice - type(MPI_File),intent(inout) :: fileh - character(len=*),intent(in) :: chkpt_type,suffix - -#ifdef TROVE_USE_MPI_ - integer(ik) :: ilen - character(len=cl) :: readbuf,filename,jchar - integer :: ierr - - if (.not.job%IOmatelem_split) return - ! - !write(job_is,"('single swap_matrix')") - ! - !call IOStart(trim(job_is),chkptIO) - ! - write(jchar, '(i4)') islice - ! - filename = trim(suffix)//trim(adjustl(jchar))//'.chk' - ! - call open_chkptfile_mpi(fileh, filename, 'read') - ! - ilen = LEN_TRIM(chkpt_type) - ! - call MPI_File_read_all(fileh, readbuf, ilen, mpi_character, mpi_status_ignore, ierr) - if ( trim(readbuf(1:ilen))/=trim(chkpt_type) ) then - if (mpi_rank .eq. 0) write (out,"(' [MPI] kinetic checkpoint slice ',a20,': header is missing or wrong',a)") filename,readbuf(1:ilen) - stop 'PTrestore_rot_kinetic_matrix_elements_mpi - in slice - header missing or wrong' - end if -#endif - end subroutine divided_slice_open_mpi - - subroutine divided_slice_close_mpi(islice,fileh,chkpt_type) -#ifdef TROVE_USE_MPI_ - use mpi_f08 -#endif - - use mpi_aux - integer(ik),intent(in) :: islice - type(MPI_File),intent(inout) :: fileh - character(len=*),intent(in) :: chkpt_type - -#ifdef TROVE_USE_MPI_ - integer(ik) :: ilen - character(len=cl) :: readbuf - integer :: ierr - - if (.not.job%IOmatelem_split) return - ! - ilen = LEN_TRIM(chkpt_type) - ! - call MPI_File_read_all(fileh, readbuf, ilen, mpi_character, mpi_status_ignore, ierr) - if ( trim(readbuf(1:ilen))/=trim(chkpt_type) ) then - if(mpi_rank .eq. 0) write (out,"(' divided_slice_close_mpi, kinetic checkpoint slice: footer is missing or wrong',a)") readbuf(1:ilen) - stop 'divided_slice_close_mpi - in slice - footer missing or wrong' - end if - ! - call close_chkptfile_mpi(fileh) -#endif - end subroutine divided_slice_close_mpi - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!! POSIX IO !!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! - ! Here we restore the vibrational (J=0) matrix elements of the rotational kinetic part G_rot and G_cor - ! - subroutine PTrestore_rot_kinetic_matrix_elements(jrot,task,chkptIO,dimen,ncontr,maxcontr,icontr) - ! - integer(ik),intent(in) :: jrot - character(len=cl),intent(in) :: task - integer(ik),intent(inout) :: chkptIO - integer(ik),intent(in) :: dimen - integer(ik),intent(inout),optional :: ncontr - integer(ik),intent(inout),optional :: maxcontr - integer(ik),intent(in),optional :: icontr - ! - integer(ik) :: alloc - character(len=cl) :: job_is,filename - character(len=18) :: buf18 - integer(ik) :: k1,k2,j,ib0,islice,chkptIO_,ideg,icontr0,istart,iend,maxcontr0 - integer(ik),allocatable :: imat_t(:,:) - integer(hik) :: rootsize,rootsize_,rootsize2,rootsize2_,nprocs,tid,icontr1,icontr2 - integer :: OMP_GET_NUM_THREADS,omp_get_thread_num - real(rk),allocatable :: mat_t(:,:),mat_(:,:) - real(rk) :: f_t - ! - ! Return if there is nothing to read - ! - if ( trim(job%IOkinet_action)/='VIB_READ'.and.trim(job%IOkinet_action)/='READ'.and..not.job%contrci_me_fast) return - ! - if ( trim(job%IOswap_matelem)/='NONE') return - ! - rootsize = int(ncontr*(ncontr+1)/2,hik) - rootsize_ = int(maxcontr*(maxcontr+1)/2,hik) - ! - rootsize2 = int(ncontr*ncontr,hik) - rootsize2_ = int(maxcontr*maxcontr,hik) - ! - !dimen = max(min(int(PT%Maxcontracts*job%compress),PT%Maxcontracts),1) - if (job%verbose>=6.and.present(icontr)) write(out,"('icontr = ',i9)") icontr - ! - select case (trim(task)) - ! - case('top') - ! - job_is ='Vib. matrix elements of the rot. kinetic' - call IOStart(trim(job_is),chkptIO) - ! - !open(chkptIO,form='unformatted',action='read',position='rewind',status='old',file=job%kinetmat_file,recordtype='variable') - ! - open(chkptIO,form='unformatted',action='read',position='rewind',status='old',file=job%kinetmat_file) - ! - read(chkptIO) buf18 - if (buf18/='Start Kinetic part') then - write (out,"(' Vib. kinetic checkpoint file ',a,' has bogus header: ',a)") job%kinetmat_file,buf18 - stop 'PTrestore_rot_kinetic_matrix_elements - bogus file format' - end if - ! - read(chkptIO) ncontr - ! - if (jrot==0.and.PT%Maxcontracts/=ncontr.and.(.not.trove%triatom_sing_resolve ) ) then - write (out,"(' Vib. kinetic checkpoint file ',a)") job%kinetmat_file - write (out,"(' Actual and stored basis sizes at J=0 do not agree ',2i0)") PT%Maxcontracts,ncontr - stop 'PTrestore_rot_kinetic_matrix_elements - in file - illegal nroots ' - end if - ! - rootsize = int(ncontr*(ncontr+1)/2,hik) - rootsize2 = int(ncontr*ncontr,hik) - ! - if (job%verbose>=6) write(out,"(/'Restore_rot_kin...: Number of elements: ',i8)") PT%Maxcontracts - ! - ! Read the indexes of the J=0 contracted basis set. - ! - if (.not.FLrotation.or.jrot==0) then - ! - allocate (imat_t(0:PT%Nclasses,ncontr),stat=alloc) - call ArrayStart('mat_t',alloc,size(imat_t),kind(imat_t)) - ! - read(chkptIO) buf18(1:10) - if (buf18(1:10)/='icontr_cnu') then - write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_cnu is missing ',a)") job%kinetmat_file,buf18(1:10) - stop 'PTrestore_rot_kinetic_matrix_elements - in file - icontr_cnu missing' - end if - ! - read(chkptIO) imat_t(0:PT%Nclasses,1:ncontr) - ! - read(chkptIO) buf18(1:11) - if (buf18(1:11)/='icontr_ideg') then - write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_ideg is missing ',a)") job%kinetmat_file,buf18(1:11) - stop 'PTrestore_rot_kinetic_matrix_elements - in file - icontr_ideg missing' - end if - ! - read(chkptIO) imat_t(0:PT%Nclasses,1:ncontr) - ! - deallocate(imat_t) - ! - else - ! - allocate (PT%icontr_cnu(0:PT%Nclasses,ncontr),PT%icontr_ideg(0:PT%Nclasses,ncontr),stat=alloc) - call ArrayStart('PT%contractive_space',alloc,size(PT%icontr_cnu),kind(PT%icontr_cnu)) - call ArrayStart('PT%contractive_space',alloc,size(PT%icontr_ideg),kind(PT%icontr_ideg)) - ! - deallocate(PT%icontr2icase) - call Arraystop('PT%contractive_space') - ! - allocate (PT%icontr2icase(ncontr,2),stat=alloc) - call ArrayStart('PT%contractive_space',alloc,size(PT%icontr2icase),kind(PT%icontr2icase)) - ! - read(chkptIO) buf18(1:10) - if (buf18(1:10)/='icontr_cnu') then - write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_cnu is missing ',a)") job%kinetmat_file,buf18(1:10) - stop 'PTrestore_rot_kinetic_matrix_elements - in file - icontr_cnu missing' - end if - ! - read(chkptIO) PT%icontr_cnu(0:PT%Nclasses,1:ncontr) - ! - read(chkptIO) buf18(1:11) - if (buf18(1:11)/='icontr_ideg') then - write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_ideg is missing ',a)") job%kinetmat_file,buf18(1:11) - stop 'PTrestore_rot_kinetic_matrix_elements - in file - icontr_ideg missing' - end if - ! - read(chkptIO) PT%icontr_ideg(0:PT%Nclasses,1:ncontr) - ! - endif - ! - if (job%vib_rot_contr.and.PTvibrational_me_calc) then - ! - read(chkptIO) buf18(1:7) - if (buf18(1:7)/='vib-rot') then - write (out,"(' Vib. kinetic checkpoint file ',a,': label vib-rot is missing ',a)") job%kinetmat_file,buf18(1:7) - stop 'PTrestore_rot_kinetic_matrix_elements - in file - vib-rot missing' - end if - ! - close(chkptIO) - ! - write(job_is,"('single swap_matrix')") - ! - call IOStart(trim(job_is),chkptIO) - ! - filename = trim(job%matelem_suffix)//"0"//'.chk' - ! - open(chkptIO,form='unformatted',action='read',status='old',file=filename) - ! - endif - ! - ! reconstruct the correlation between the vib. indices for J=0 and current J - ! - call Find_groundstate_icontr(maxcontr) - ! - if (job%verbose>=4) write(out,"(' ...done!')") - ! - if (job%verbose>=4.and.maxcontr/=ncontr) then - write (out,"(' The contracted basis set is reduced: ',i8,' -> ',i8)") ncontr,maxcontr - end if - ! - if (maxcontr>ncontr) then - write (out,"(' Actual and stored basis sizes at J=0 do not agree (maxcontr,ncontr) ',2i8)") maxcontr,ncontr - stop 'PTrestore_rot_kinetic_matrix_elements - in file - illegal ncontr ' - end if - ! - case('rot') - ! - if (job%verbose>=4) write(out,"(' Read and process rotational part...')") - ! - ! Read the rotational part only - ! - allocate(mat_(maxcontr,maxcontr),stat=alloc) - call ArrayStart('PThamiltonian_contract: mat_',alloc,1,kind(mat_),rootsize2_) - ! - if (.not.job%IOmatelem_split) then - ! - read(chkptIO) buf18(1:5) - if (buf18(1:5)/='g_rot') then - write (out,"(' Vib. kinetic checkpoint file ',a,': g_rot is missing ',a)") job%kinetmat_file,buf18(1:5) - stop 'PTrestore_rot_kinetic_matrix_elements - in file - g_rot missing' - end if - ! - endif - ! - allocate(grot(3,3),stat=alloc) - ! - islice = 0 - chkptIO_ = chkptIO - ! - do k1 = 1,3 - ! - do k2 = 1,3 - ! - islice = islice + 1 - ! - call divided_slice_open(islice,chkptIO_,'g_rot',job%matelem_suffix) - ! - allocate(grot(k1,k2)%me(maxcontr,maxcontr),stat=alloc) - call ArrayStart('grot-matrix',alloc,1,kind(f_t),rootsize2_) - ! - read(chkptIO_) mat_ - ! - grot(k1,k2)%me(1:ncontr,1:ncontr) = mat_(1:ncontr,1:ncontr) - ! - call divided_slice_close(islice,chkptIO_,'g_rot') - ! - enddo - enddo - ! - deallocate(mat_) - call ArrayStop('PThamiltonian_contract: mat_') - ! - if (job%verbose>=4) write(out,"(' ...done!')") - ! - case('cor') - ! - if (job%verbose>=4) write(out,"(' Read and process Coriolis part...')") - ! - allocate(mat_(maxcontr,maxcontr),stat=alloc) - call ArrayStart('PThamiltonian_contract: mat_',alloc,1,kind(mat_),rootsize2_) - ! - if (.not.job%IOmatelem_split) then - ! - read(chkptIO) buf18(1:5) - if (buf18(1:5)/='g_cor') then - write (out,"(' Vib. kinetic checkpoint file ',a,': g_cor is missing ',a)") job%kinetmat_file,buf18(1:5) - stop 'PTrestore_rot_kinetic_matrix_elements - in file - g_cor missing' - end if - ! - endif - ! - allocate(gcor(3),stat=alloc) - ! - islice = 9 - chkptIO_ = chkptIO - ! - do k1 = 1,3 - ! - islice = islice + 1 - ! - allocate(gcor(k1)%me(maxcontr,maxcontr),stat=alloc) - call ArrayStart('gcor-matrix',alloc,1,kind(f_t),rootsize2_) - ! - call divided_slice_open(islice,chkptIO_,'g_cor',job%matelem_suffix) - ! - read(chkptIO_) mat_ - ! - gcor(k1)%me(1:ncontr,1:ncontr) = mat_(1:ncontr,1:ncontr) - ! - call divided_slice_close(islice,chkptIO_,'g_cor') - ! - enddo - ! - deallocate(mat_) - call ArrayStop('PThamiltonian_contract: mat_') - ! - if (job%verbose>=4) write(out,"(' ...done!')") - ! - case('vib') - ! - if (job%verbose>=4) write(out,"(' Read and process vibrational part...')") - ! - if (.not.job%IOmatelem_split.and.( (.not.FLrotation.or.jrot==0).and.trim(job%IOkinet_action)/='VIB_READ' ) ) then - ! - allocate(mat_t(maxcontr,maxcontr),stat=alloc) - call ArrayStart('PThamiltonian_contract: mat_t',alloc,1,kind(mat_t),rootsize2_) - ! - read(chkptIO) buf18(1:5) - if (buf18(1:5)/='g_rot') then - write (out,"(' Vib. kinetic checkpoint file ',a,': g_rot is missing ',a)") job%kinetmat_file,buf18(1:5) - stop 'PTrestore_rot_kinetic_matrix_elements - in file - g_rot missing' - end if - ! - do k1 = 1,3 - do k2 = 1,3 - ! - read(chkptIO) mat_t - ! - enddo - enddo - ! - read(chkptIO) buf18(1:5) - if (buf18(1:5)/='g_cor') then - write (out,"(' Vib. kinetic checkpoint file ',a,': g_cor is missing ',a)") job%kinetmat_file,buf18(1:5) - stop 'PTrestore_rot_kinetic_matrix_elements - in file - g_cor missing' - end if - ! - !do k1 = 1,PT%Nmodes - do k2 = 1,3 - ! - read(chkptIO) mat_t - ! - enddo - !enddo - ! - deallocate(mat_t) - call ArrayStop('mat_t') - ! - endif - ! - if (job%verbose>=6) write(out,"(' rootsize_,rootsize = ',2i0)") rootsize_,rootsize - ! - read(chkptIO) buf18(1:4) - if (buf18(1:4)/='hvib') then - write (out,"(' Vib. kinetic checkpoint file ',a,': hvib is missing ',a)") job%kinetmat_file,buf18(1:4) - if (buf18(1:4)=='g_ro') then - write (out,"(' Most likely non-divided chk-points used with MATELEM READ SPLIT')") - write (out,"(' Re-do MATELEM SAVE SPLIT or use MATELEM SPLIT READ')") - endif - stop 'PTrestore_rot_kinetic_matrix_elements - in file - hvib or End missing' - end if - ! - allocate(mat_(maxcontr,maxcontr),stat=alloc) - call ArrayStart('PThamiltonian_contract: mat_',alloc,1,kind(mat_),rootsize2_) - ! - allocate(hvib%me(maxcontr,maxcontr),stat=alloc) - call ArrayStart('hvib-matrix',alloc,1,kind(f_t),rootsize2_) - ! - read(chkptIO) mat_ - ! - hvib%me(1:ncontr,1:ncontr) = mat_(1:ncontr,1:ncontr) - ! - deallocate(mat_) - call ArrayStop('PThamiltonian_contract: mat_') - ! - read(chkptIO) buf18(1:16) - if (buf18(1:16)/='End Kinetic part') then - write (out,"(' Vib. kinetic checkpoint file ',a,' has bogus footer: ',a)") job%kinetmat_file,buf18(1:16) - stop 'PTrestore_rot_kinetic_matrix_elements - bogus file format' - end if - ! - close(chkptIO,status='keep') - ! - if (job%verbose>=4) write(out,"(' ...done!')") - ! case('top-icontr') - ! + if (blacs_size .gt. 1) then + if (mpi_rank .eq. 0) write (out,*) "task=top-icontr is not implemented for MPI." + stop 'PTrestore_rot_kinetic_matrix_elements - unimplemented MPI' + endif + nprocs = 1 tid = 0 !$omp parallel private(tid) @@ -9528,14 +8745,14 @@ subroutine PTrestore_rot_kinetic_matrix_elements(jrot,task,chkptIO,dimen,ncontr, ! if (FLrotation.and.jrot/=0) then ! - allocate(grot(3,3),stat=alloc) + allocate(grot(3,3),stat=ierr) ! rootsize2_ = int(maxcontr,hik)*int(PT%max_deg_size,hik)*9_hik*int(nprocs,hik) - call ArrayStart('PThamiltonian_contract: grot',alloc,1,rk,rootsize2_) + call ArrayStart('PThamiltonian_contract: grot',ierr,1,rk,rootsize2_) ! - allocate(gcor(3),stat=alloc) + allocate(gcor(3),stat=ierr) rootsize2_ = int(maxcontr,hik)*int(PT%max_deg_size,hik)*int(PT%Nmodes,hik)*3_hik*int(nprocs,hik) - call ArrayStart('PThamiltonian_contract: grot',alloc,1,rk,rootsize2_) + call ArrayStart('PThamiltonian_contract: grot',ierr,1,rk,rootsize2_) ! if (job%vib_rot_contr) then ! @@ -9546,11 +8763,8 @@ subroutine PTrestore_rot_kinetic_matrix_elements(jrot,task,chkptIO,dimen,ncontr, endif ! endif - ! + case('rot-icontr') ! rotational part for the vib-rot contraction scheme - ! - ! Read the rotational part only - ! if (.not.job%IOmatelem_split ) then ! write (out,"('PTrestore_rot_kinetic_matrix_elements: vib-rot can be used with MATELEM SAVE SPLIT only ')") @@ -9574,21 +8788,18 @@ subroutine PTrestore_rot_kinetic_matrix_elements(jrot,task,chkptIO,dimen,ncontr, ! if (associated(grot(k1,k2)%me)) deallocate(grot(k1,k2)%me) ! - allocate(grot(k1,k2)%me(maxcontr,icontr1:icontr2),stat=alloc) + allocate(grot(k1,k2)%me(maxcontr,icontr1:icontr2),stat=ierr) ! - write(job_is,"('single swap_matrix #',i8)") islice - call IOStart(trim(job_is),chkptIO_) + write(job_id,"('single swap_matrix #',i8)") islice + call IOStart(trim(job_id),chkptIO_) ! read(chkptIO_) grot(k1,k2)%me ! enddo enddo - ! + case('cor-icontr') ! corriolis part for the vib-rot contraction scheme - ! - ! Read the Corriolis part only - ! - ! + ! Read the Coriolis part only if (.not.job%IOmatelem_split ) then ! write (out,"('PTrestore_rot_kinetic_matrix_elements: vib-rot can be used with MATELEM SAVE SPLIT only ')") @@ -9616,17 +8827,16 @@ subroutine PTrestore_rot_kinetic_matrix_elements(jrot,task,chkptIO,dimen,ncontr, ! if (associated(gcor(k1)%me)) deallocate(gcor(k1)%me) ! - allocate(gcor(k1)%me(maxcontr,icontr1:icontr2),stat=alloc) + allocate(gcor(k1)%me(maxcontr,icontr1:icontr2),stat=ierr) ! - write(job_is,"('single swap_matrix #',i8)") islice - call IOStart(trim(job_is),chkptIO_) + write(job_id,"('single swap_matrix #',i8)") islice + call IOStart(trim(job_id),chkptIO_) ! read(chkptIO_) gcor(k1)%me ! enddo - ! + case('vib-icontr') ! vibrational part for the vib-rot contraction scheme - ! ! !if (job%verbose>=4.and.irow==0) write(out,"(' Read and process vibrational part...')") ! @@ -9641,10 +8851,10 @@ subroutine PTrestore_rot_kinetic_matrix_elements(jrot,task,chkptIO,dimen,ncontr, ! if (icontr==1) then ! - read(chkptIO) buf18(1:4) - if (buf18(1:4)/='hvib') then - write (out,"(' Vib. kinetic checkpoint file ',a,': hvib is missing ',a)") job%kinetmat_file,buf18(1:4) - if (buf18(1:4)=='g_ro') write (out,"(' Most likely non-divided chk-points used with MATELEM READ DIVIDED')") + read(chkptIO) readbuf(1:4) + if (readbuf(1:4)/='hvib') then + write (out,"(' Vib. kinetic checkpoint file ',a,': hvib is missing ',a)") job%kinetmat_file,readbuf(1:4) + if (readbuf(1:4)=='g_ro') write (out,"(' Most likely non-divided chk-points used with MATELEM READ DIVIDED')") write (out,"(' Re-do MATELEM SAVE DIVIDE or use MATELEM READ!')") stop 'PTrestore_rot_kinetic_matrix_elements - in file - hvib or End missing' end if @@ -9665,8 +8875,8 @@ subroutine PTrestore_rot_kinetic_matrix_elements(jrot,task,chkptIO,dimen,ncontr, ! if (job%verbose>=6) write(out,"('allocate hvib for ',i9,' x ',i8,' -> ',i8)") maxcontr,icontr1,icontr2 ! - allocate(hvib%me(maxcontr,icontr1:icontr2),stat=alloc) - call ArrayStart('hvib-matrix',alloc,1,kind(f_t),rootsize2_) + allocate(hvib%me(maxcontr,icontr1:icontr2),stat=ierr) + call ArrayStart('hvib-matrix',ierr,1,kind(f_t),rootsize2_) ! read(chkptIO) hvib%me ! @@ -9674,9 +8884,9 @@ subroutine PTrestore_rot_kinetic_matrix_elements(jrot,task,chkptIO,dimen,ncontr, ! if (icontr==maxcontr0) then ! - read(chkptIO) buf18(1:4) - if (buf18(1:4)/='hvib') then - write (out,"(' Vib. kinetic checkpoint file ',a,' has bogus footer: ',a)") job%kinetmat_file,buf18(1:4) + read(chkptIO) readbuf(1:4) + if (readbuf(1:4)/='hvib') then + write (out,"(' Vib. kinetic checkpoint file ',a,' has bogus footer: ',a)") job%kinetmat_file,readbuf(1:4) stop 'PTrestore_rot_kinetic_matrix_elements - bogus file format' end if ! @@ -9687,68 +8897,9 @@ subroutine PTrestore_rot_kinetic_matrix_elements(jrot,task,chkptIO,dimen,ncontr, endif ! end select - ! - contains - ! - subroutine divided_slice_open(islice,chkptIO,name,suffix) - ! - implicit none - integer(ik),intent(in) :: islice - integer(ik),intent(inout) :: chkptIO - character(len=*),intent(in) :: name,suffix - character(len=4) :: jchar - character(len=cl) :: buf,filename,job_is - integer(ik) :: ilen - logical :: ifopened - ! - if (.not.job%IOmatelem_split) return - ! - write(job_is,"('single swap_matrix')") - ! - call IOStart(trim(job_is),chkptIO) - ! - write(jchar, '(i4)') islice - ! - filename = trim(suffix)//trim(adjustl(jchar))//'.chk' - ! - open(chkptIO,form='unformatted',action='read',position='rewind',status='old',file=filename) - ! - ilen = LEN_TRIM(name) - ! - read(chkptIO) buf(1:ilen) - if ( trim(buf(1:ilen))/=trim(name) ) then - write (out,"(' kinetic checkpoint slice ',a20,': header is missing or wrong',a)") filename,buf(1:ilen) - stop 'PTrestore_rot_kinetic_matrix_elements - in slice - header missing or wrong' - end if - ! - end subroutine divided_slice_open - ! - subroutine divided_slice_close(islice,chkptIO,name) - ! - integer(ik),intent(in) :: islice - integer(ik),intent(inout) :: chkptIO - character(len=*),intent(in) :: name - character(len=4) :: jchar - character(len=cl) :: buf,filename - integer(ik) :: ilen - logical :: ifopened - ! - if (.not.job%IOmatelem_split) return - ! - ilen = LEN_TRIM(name) - ! - read(chkptIO) buf(1:ilen) - if ( trim(buf(1:ilen))/=trim(name) ) then - write (out,"(' divided_slice_close, kinetic checkpoint slice ',a,': footer is missing or wrong',a)") & - filename,buf(1:ilen) - stop 'divided_slice_close - in slice - footer missing or wrong' - end if - ! - close(chkptIO) - ! - end subroutine divided_slice_close - ! + contains + subroutine divided_slice_open_vib_rot(islice,suffix) ! implicit none @@ -9810,9 +8961,65 @@ subroutine divided_slice_close_vib_rot(islice,name,chkptIO) ! end subroutine divided_slice_close_vib_rot - end subroutine PTrestore_rot_kinetic_matrix_elements - ! - ! + subroutine divided_slice_open(islice,sliceHandler,chkpt_type,suffix) + use mpi_aux + implicit none + integer(ik),intent(in) :: islice + class(ioHandlerBase), allocatable, intent(out) :: sliceHandler + character(len=*),intent(in) :: chkpt_type,suffix + + integer(ik) :: ilen + character(len=cl) :: readbuf,filename,jchar + type(ErrorType) :: err + integer :: ierr + + if (.not.job%IOmatelem_split) return + ! + !write(job_is,"('single swap_matrix')") + ! + !call IOStart(trim(job_is),chkptIO) + ! + write(jchar, '(i4)') islice + ! + filename = trim(suffix)//trim(adjustl(jchar))//'.chk' + ! + call openFile(sliceHandler, filename, err, action='read', & + position='rewind', status='old', form='unformatted') + HANDLE_ERROR(err) + ! + ilen = LEN_TRIM(chkpt_type) + ! + call sliceHandler%read(readbuf(1:ilen)) + if ( trim(readbuf(1:ilen))/=trim(chkpt_type) ) then + if (mpi_rank .eq. 0) write (out,"('divided_slice_open, kinetic checkpoint slice ',a20,': header is missing or wrong',a)") filename,readbuf(1:ilen) + stop 'PTrestore_rot_kinetic_matrix_elements- in slice - header missing or wrong' + end if + end subroutine divided_slice_open + + subroutine divided_slice_close(islice,sliceHandler,chkpt_type) + use mpi_aux + integer(ik),intent(in) :: islice + class(ioHandlerBase), allocatable, intent(inout) :: sliceHandler + character(len=*),intent(in) :: chkpt_type + + integer(ik) :: ilen + character(len=cl) :: readbuf + type(ErrorType) :: err + integer :: ierr + + if (.not.job%IOmatelem_split) return + ! + ilen = LEN_TRIM(chkpt_type) + ! + call sliceHandler%read(readbuf(1:ilen)) + if ( trim(readbuf(1:ilen))/=trim(chkpt_type) ) then + if(mpi_rank .eq. 0) write (out,"('divided_slice_close, kinetic checkpoint slice: footer is missing or wrong',a)") readbuf(1:ilen) + stop 'divided_slice_close - in slice - footer missing or wrong' + end if + deallocate(sliceHandler) + end subroutine divided_slice_close + + end subroutine ! ! We construct the Hamiltonian matrix in symm. adapted representaion ! for the K-factorized rotational basis @@ -9847,8 +9054,11 @@ subroutine symm_mat_element_vector(jrot,irow,ijterm,func,mat_t,no_diagonalizatio hcontr = 0.0 ! do jrow = 1,irow - if ( present(no_diagonalization).and.no_diagonalization.and.jrow/=irow ) cycle - + ! + if (present(no_diagonalization)) then + if (no_diagonalization.and.jrow/=irow ) cycle + endif + ! cnu_j(:) = PT%contractive_space(:,jrow) jsize = PT%Index_deg(jrow)%size1 @@ -9977,153 +9187,6 @@ subroutine symm_mat_element_vector(jrot,irow,ijterm,func,mat_t,no_diagonalizatio !call TimerStop('Symmetrized Hamiltonian - one column') ! end subroutine symm_mat_element_vector - - ! - ! In this version of the routine we construct the Hamiltonian matrix in symm. adapted representaion - ! for rotational basis which is not factorized with "K". - ! - !subroutine symm_mat_element_vector(jrot,irow,ijterm,func,mat_t,no_diagonalization) - - !integer(ik),intent(in) :: jrot,irow,ijterm(:,:) - !real(rk),external :: func - !real(rk),intent(out) :: mat_t(:,:,:) - !logical,optional,intent(in) :: no_diagonalization - !! - !integer(ik) :: cnu_i(0:PT%Nclasses),cnu_j(0:PT%Nclasses) - !integer(ik) :: deg_i(0:PT%Nclasses),deg_j(0:PT%Nclasses) - !real(rk) :: mat_elem - !integer(ik) :: isize,jsize,ielem,jelem - !integer(ik) :: jrow,ideg,jdeg,isym,jsym,iL,iR,iterm,jterm,icontr,jcontr - !real(rk) :: vec_i(PT%max_deg_size),vec_j(PT%max_deg_size) - - !real(rk), dimension(:,:,:), allocatable :: hcontr - - !allocate(hcontr(PT%max_deg_size,PT%max_deg_size,irow)) - !! - !mat_t = 0 - !! - !cnu_i(:) = PT%contractive_space(:,irow) - !! - !isize = PT%Index_deg(irow)%size1 - !! - !do jrow = 1,irow - !! - !if ( present(no_diagonalization).and.no_diagonalization.and.jrow/=irow ) cycle - !! - !cnu_j(:) = PT%contractive_space(:,jrow) - !! - !jsize = PT%Index_deg(jrow)%size1 - !! - !do ideg = 1,isize - !! - !deg_i(:) = PT%Index_deg(irow)%icoeffs(:,ideg) - !! - !do jdeg = 1,jsize - !! - !deg_j(:) = PT%Index_deg(jrow)%icoeffs(:,jdeg) - !! - !icontr = PT%icase2icontr(irow,ideg) - !jcontr = PT%icase2icontr(jrow,jdeg) - !! - !! Matrix elements - !! - !hcontr(ideg,jdeg,jrow) = func(icontr,jcontr,jrot,cnu_i(0),cnu_j(0),deg_i(0),deg_j(0)) - !! - !enddo - !! - !enddo - !! - !end do - - !do jrow=1,irow - !jsize = PT%Index_deg(jrow)%size1 - !do isym = 1,sym%Nrepresen - !! - !iterm = ijterm(irow,isym) - !! - !do jsym = 1,isym - !! - !jterm = ijterm(jrow,jsym) - !! - !do ielem = 1,PT%irr(isym)%N(irow) - !! - !vec_i(1:isize) = PT%irr(isym)%repres(iterm+ielem,1,1:isize) - !! - !do jelem = 1,PT%irr(jsym)%N(jrow) - !! - !vec_j(1:jsize) = PT%irr(jsym)%repres(jterm+jelem,1,1:jsize) - !! - !vec_j(1:isize) = matmul(hcontr(1:isize,1:jsize, jrow),vec_j(1:jsize)) - !! - !mat_elem = dot_product(vec_i(1:isize),vec_j(1:isize)) - !! - !if (isym==jsym) then - !! - !iL = ielem - !iR = jterm+jelem - !! - !if (iterm+ielem(10.0_rk)**(-(rk-3))) then - !! - !! We print out non-zero mat. elements between different symmetries, which have to be zero. - !! - !if (job%verbose>=6) & - !write(out,"('<',a4,2i6,'|H|',a4,2i6,'> = ',g18.10)") & - !sym%label(isym),irow,iterm+ielem,sym%label(jsym),jrow,jterm+jelem,mat_elem - !! - !! if this error is too big - we stop - !! - !if(abs(mat_elem)>1.0_rk) then - !!write(out,"(/'A non-diagonal mat. element between different symmetries:')") - !write(out,"(/'<',a4,3i6,'|H|',a4,3i6,'> = ',g18.10,a)") & - !sym%label(isym),irow,ielem,iterm+ielem,sym%label(jsym),jrow,jelem,jterm+jelem,mat_elem,& - !' Non-diagonal element (euler) between different symmetries is too large!' - !stop 'non-zero element between two symmetries - symm_mat_element_vector' - !endif - !! - !endif - !endif - !! - !enddo - !enddo - !enddo - !! - !enddo - !! - !enddo - !! - !! - !!call TimerStop('Symmetrized Hamiltonian - one column') - !! - !! - !end subroutine symm_mat_element_vector - ! ! In this version of the routine we construct the Hamiltonian matrix in symm. adapted representaion ! for rotational basis which is not factorized with "K" and using the icontr-based sorting @@ -11350,7 +10413,7 @@ subroutine diagonalization_contract(jrot,gamma,dimen_s,mat,zpe,rlevel,total_root ! spur = spur*exp(-beta*mat0) ! - write(out, '(/1x, a, 1x, es16.8)'), 'qpart = ', spur + write(out, '(/1x, a, 1x, es16.8)') 'qpart = ', spur ! !mat = mat / (-planck * vellgt) * (boltz * intensity%temperature) !do ielem = 1, dimen_s @@ -11363,7 +10426,7 @@ subroutine diagonalization_contract(jrot,gamma,dimen_s,mat,zpe,rlevel,total_root ! if (gamma==sym%Nrepresen) then ! - write(out, '(/1x, a, 1x, es16.8)'), 'partition function value is', job%partfunc%value + write(out, '(/1x, a, 1x, es16.8)') 'partition function value is', job%partfunc%value ! endif ! @@ -16097,7 +15160,9 @@ subroutine PTcontracted_matelem_class(jrot) integer :: startdim, enddim, blocksize_, ierr, b, req_count, offset type(MPI_Request),allocatable :: reqs(:) type(MPI_Status) :: reqstat - type(MPI_File) :: chkptMPIIO, chkptMPIIO_ + class(ioHandlerBase), allocatable :: ioHandler + class(ioHandlerBase), allocatable :: sliceHandler + class(ioHandlerBase), allocatable :: extFmatHandler integer(kind=MPI_Offset_kind) :: mpioffset integer :: mpisz ! @@ -16117,7 +15182,6 @@ subroutine PTcontracted_matelem_class(jrot) type(PTcoeffsT) :: tmat(PT%Nclasses),mat_tt(PT%Nclasses) type(PTcoeffT),pointer :: fl - class(ioHandlerBase), allocatable :: ioHandler type(ErrorType) :: err ! ! @@ -16222,17 +15286,11 @@ subroutine PTcontracted_matelem_class(jrot) ! ! Prepare the checkpoint file ! - job_is ='Vib. matrix elements of the rot. kinetic part' - call IOStart(trim(job_is),chkptIO) + !job_is ='Vib. matrix elements of the rot. kinetic part' + !call IOStart(trim(job_is),chkptIO) ! -#ifdef TROVE_USE_MPI_ - allocate(ioHandlerMPI::ioHandler) -#else - allocate(ioHandlerFTN::ioHandler) -#endif - call ioHandler%open(& - job%kinetmat_file, err, & - action='write', position='rewind', status='replace', form='unformatted') + call openFile(ioHandler, job%kinetmat_file, err, action='write', & + position='rewind', status='replace', form='unformatted') HANDLE_ERROR(err) call ioHandler%write('Start Kinetic part') @@ -16469,11 +15527,7 @@ subroutine PTcontracted_matelem_class(jrot) if (trim(job%IOkinet_action)=='SAVE') then if (job%IOmatelem_split) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call write_divided_slice_mpi(islice,'g_rot',job%matelem_suffix,mdimen,grot_t) - else - call write_divided_slice(islice,'g_rot',job%matelem_suffix,mdimen,grot_t) - endif + call write_divided_slice(islice,'g_rot',job%matelem_suffix,mdimen,grot_t) ! else ! @@ -16554,11 +15608,7 @@ subroutine PTcontracted_matelem_class(jrot) ! if (job%IOmatelem_divide) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call write_divided_slice_mpi(islice,'g_cor',job%matelem_suffix,mdimen,grot_t) - else - call write_divided_slice(islice,'g_cor',job%matelem_suffix,mdimen,grot_t) - endif + call write_divided_slice(islice,'g_cor',job%matelem_suffix,mdimen,grot_t) ! else ! @@ -16576,11 +15626,7 @@ subroutine PTcontracted_matelem_class(jrot) ! if (job%IOmatelem_split) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call write_divided_slice_mpi(islice,'g_cor',job%matelem_suffix,mdimen,gcor_t) - else - call write_divided_slice(islice,'g_cor',job%matelem_suffix,mdimen,gcor_t) - endif + call write_divided_slice(islice,'g_cor',job%matelem_suffix,mdimen,gcor_t) ! else ! @@ -16652,39 +15698,35 @@ subroutine PTcontracted_matelem_class(jrot) if (job%verbose>=2) write(out,"('...end!')") ! if (treat_rotation.and.trim(job%IOkinet_action)=='SAVE') then - if (trim(job%kinetmat_format).eq.'MPIIO') then - write(out,*) "TODO implement MPI-IO version !POSIXIO!" - stop "Not yet implemented" - else - ! - ! store the rotational matrix elements - ! - call ioHandler%write('g_rot') - ! - do k1 = 1,3 - do k2 = 1,3 - ! - call ioHandler%write(grot_(k1,k2,:,:), mdimen_) - ! - enddo + write(out,*) "WARNING: MPI verion of treat_rotation is not verified" + ! + ! store the rotational matrix elements + ! + call ioHandler%write('g_rot') + ! + do k1 = 1,3 + do k2 = 1,3 + ! + call ioHandler%write(grot_(k1,k2,:,:), mdimen_) + ! enddo - ! - call ioHandler%write('g_cor') - ! - ! store the Coriolis matrix elements - ! - do k1 = 1,PT%Nmodes - do k2 = 1,3 - ! - call ioHandler%write(gcor_(k1,k2,:,:), mdimen_) - ! - enddo + enddo + ! + call ioHandler%write('g_cor') + ! + ! store the Coriolis matrix elements + ! + do k1 = 1,PT%Nmodes + do k2 = 1,3 + ! + call ioHandler%write(gcor_(k1,k2,:,:), mdimen_) + ! enddo - ! - deallocate(grot_,gcor_) - call ArrayStop('grot-gcor-fields') - ! - endif + enddo + ! + deallocate(grot_,gcor_) + call ArrayStop('grot-gcor-fields') + ! endif ! else ! if (.not.job%IOmatelem_split.or.job%iswap(1)==0 ) then @@ -16785,11 +15827,7 @@ subroutine PTcontracted_matelem_class(jrot) enddo !$omp end parallel do ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call write_divided_slice_mpi(islice,'g_vib',job%matelem_suffix,mdimen,gvib_t) - else - call write_divided_slice(islice,'g_vib',job%matelem_suffix,mdimen,gvib_t) - endif + call write_divided_slice(islice,'g_vib',job%matelem_suffix,mdimen,gvib_t) ! else ! @@ -16872,11 +15910,7 @@ subroutine PTcontracted_matelem_class(jrot) ! islice = (PT%Nmodes+3)*3+PT%Nmodes**2+1 ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call write_divided_slice_mpi(islice,'g_vib',job%matelem_suffix,mdimen,gvib_t) - else - call write_divided_slice(islice,'g_vib',job%matelem_suffix,mdimen,gvib_t) - endif + call write_divided_slice(islice,'g_vib',job%matelem_suffix,mdimen,gvib_t) ! if (job%IOmatelem_split.and.job%iswap(1)==1) job%iswap(1)=0 ! @@ -16894,53 +15928,32 @@ subroutine PTcontracted_matelem_class(jrot) ! f_t = -0.5_rk ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - do islice = iterm1,iterm2 - ! - if (islice==(PT%Nmodes+3)*3+PT%Nmodes**2+1) f_t = 1.0_rk - ! - !! TODO !! - write(*,*) "TODO: NEEDS VERIFICATION" - call divided_slice_open_mpi(islice,chkptMPIIO_,'g_vib',job%matelem_suffix) - ! - call co_read_matrix_distr_ordered(gvib_t, mdimen, startdim, enddim, chkptMPIIO_) - ! - do b=1,comm_size - if (send_or_recv(b).ge.0) then - !$omp parallel do private(icoeff,jcoeff) shared(b,grot_t) schedule(static) - do icoeff=startdim,enddim - do jcoeff=((b-1)*mdimen_p)+1,b*mdimen_p - if (jcoeff .gt. PT%Maxcontracts) cycle - hvib_t(jcoeff,icoeff) = hvib_t(jcoeff,icoeff)+f_t*gvib_t(jcoeff,icoeff) - enddo + do islice = iterm1,iterm2 + ! + if (islice==(PT%Nmodes+3)*3+PT%Nmodes**2+1) f_t = 1.0_rk + ! +#ifdef TROVE_USE_MPI_ + write(out,*) "WARNING: slice handling of g_vib using MPI is unverified" +#endif + call divided_slice_open(islice,sliceHandler,'g_vib',job%matelem_suffix) + ! + call sliceHandler%read(gvib_t, mdimen) + ! + do b=1,comm_size + if (send_or_recv(b).ge.0) then + !$omp parallel do private(icoeff,jcoeff) shared(b,grot_t) schedule(static) + do icoeff=startdim,enddim + do jcoeff=((b-1)*mdimen_p)+1,b*mdimen_p + if (jcoeff .gt. PT%Maxcontracts) cycle + hvib_t(jcoeff,icoeff) = hvib_t(jcoeff,icoeff)+f_t*gvib_t(jcoeff,icoeff) enddo - endif - enddo - ! - call divided_slice_close_mpi(islice,chkptMPIIO_,'g_vib') - ! - enddo - else - do islice = iterm1,iterm2 - ! - if (islice==(PT%Nmodes+3)*3+PT%Nmodes**2+1) f_t = 1.0_rk - ! - call divided_slice_open(islice,chkptIO_,'g_vib',job%matelem_suffix) - ! - read(chkptIO_) gvib_t - ! - !$omp parallel do private(icoeff,jcoeff) shared(hvib_t) schedule(dynamic) - do icoeff=1,mdimen - do jcoeff=1,icoeff - hvib_t(jcoeff,icoeff) = hvib_t(jcoeff,icoeff)+f_t*gvib_t(jcoeff,icoeff) - enddo - enddo - !$omp end parallel do - ! - call divided_slice_close(islice,chkptIO_,'g_vib') - ! + enddo + endif enddo - endif + ! + call divided_slice_close(islice,sliceHandler,'g_vib') + ! + enddo ! gvib_t = 0 ! @@ -17028,43 +16041,12 @@ subroutine PTcontracted_matelem_class(jrot) ! Prepare the checkpoint file ! job_is ='external field contracted matrix elements for J=0' - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - call MPI_File_open(mpi_comm_world, job%extFmat_file, mpi_mode_wronly+mpi_mode_create, mpi_info_null, chkptMPIIO, ierr) - mpioffset=0 - call MPI_File_set_size(chkptMPIIO, mpioffset, ierr) - if (mpi_rank.eq.0) then !AT - call TimerStart('mpiiosingle') !AT - - call MPI_File_write_shared(chkptMPIIO,'[MPIIO]Start external field',27,mpi_character,mpi_status_ignore,ierr) - !call MPI_File_write_shared(chkptMPIIO,'Start external field',20,mpi_character,mpi_status_ignore,ierr) - call MPI_File_write_shared(chkptMPIIO, PT%Maxcontracts, 1, mpi_integer, mpi_status_ignore, ierr) - ! - ! store the bookkeeping information about the contr. basis set - ! - !call PTstoreMPI_icontr_cnu(PT%Maxcontracts,chkptMPIIO,job%IOkinet_action) + call openFile(extFmatHandler, job%extFmat_file, err, action='write', & + form='unformatted',position='rewind',status='replace') + HANDLE_ERROR(err) - call TimerStop('mpiiosingle') !AT - ! - !else - ! call TimerStart('mpiiosingle') !AT - ! ! - ! call MPI_File_write_shared(chkptMPIIO,'',0,mpi_character,mpi_status_ignore,ierr) - ! call MPI_File_write_shared(chkptMPIIO,'',0,mpi_character,mpi_status_ignore,ierr) - ! ! - ! call TimerStop('mpiiosingle') !AT - endif -#endif - else - call IOStart(trim(job_is),chkptIO) - ! - open(chkptIO,form='unformatted',action='write',position='rewind',status='replace',file=job%extFmat_file) - write(chkptIO) 'Start external field' - ! - ! store the matrix elements - ! - write(chkptIO) PT%Maxcontracts - endif + call extFmatHandler%write('Start external field') + call extFmatHandler%write(PT%Maxcontracts) ! endif ! @@ -17110,30 +16092,14 @@ subroutine PTcontracted_matelem_class(jrot) ! if (job%IOextF_divide) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call write_divided_slice_mpi(imu,'extF',job%extmat_suffix,mdimen,extF_t) - else - call write_divided_slice(imu,'extF',job%extmat_suffix,mdimen,extF_t) - endif + call write_divided_slice(imu,'extF',job%extmat_suffix,mdimen,extF_t) ! else ! ! always store the matrix elements of the extF moment ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - if(mpi_rank.eq.0) then - call MPI_File_write_shared(chkptMPIIO,imu,1,mpi_integer,mpi_status_ignore,ierr) - endif - !call mpi_barrier(mpi_comm_world,ierr) - ! - call co_write_matrix_distr(extF_t,mdimen, startdim, enddim,chkptMPIIO) -#endif - else - write(chkptIO) imu - ! - write(chkptIO) extF_t - endif + call extFmatHandler%write(imu) + call extFmatHandler%write(extF_t, mdimen) ! endif ! @@ -17176,24 +16142,7 @@ subroutine PTcontracted_matelem_class(jrot) ! endif ! -#ifdef TROVE_USE_MPI_ - call mpi_barrier(mpi_comm_world,ierr) -#endif - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - if (mpi_rank.eq.0) then !AT - if(.not.job%IOextF_divide) then - !call MPI_File_seek(chkptMPIIO, mpioffset, MPI_SEEK_END) - call MPI_File_write_shared(chkptMPIIO,'End external field',18,mpi_character,mpi_status_ignore,ierr) - !else - ! call MPI_File_write_shared(chkptMPIIO,'',0,mpi_character,mpi_status_ignore,ierr) - endif - endif - call MPI_File_close(chkptMPIIO, ierr) -#endif - else - if (.not.job%IOextF_divide) write(chkptIO) 'End external field' - endif + if (.not.job%IOextF_divide) call extFmatHandler%write('End external field') ! endif ! @@ -17281,213 +16230,103 @@ subroutine PTcontracted_matelem_class(jrot) ! stop 'cannot proceede for extF = divide-join with matelem /= read' !endif ! + contains subroutine write_divided_slice(islice,name,suffix,N,field) - ! - integer(ik),intent(in) :: islice - character(len=*),intent(in) :: name,suffix - integer(ik),intent(in) :: N - real(rk),intent(in) :: field(N,N) - character(len=4) :: jchar - integer(ik) :: chkptIO - character(len=cl) :: filename - character(len=cl) :: job_is - ! - write(job_is,"('single swap_matrix')") - ! - call IOStart(trim(job_is),chkptIO) - ! - write(jchar, '(i4)') islice - ! - filename = trim(suffix)//trim(adjustl(jchar))//'.chk' - ! - open(chkptIO,form='unformatted',action='write',position='rewind',status='replace',file=filename) - ! - write(chkptIO) trim(name) - ! - write(chkptIO) field - ! - write(chkptIO) trim(name) - ! - close(chkptIO) - ! - end subroutine write_divided_slice - - subroutine write_divided_slice_mpi(islice,name,suffix,N,field) ! integer(ik),intent(in) :: islice character(len=*),intent(in) :: name,suffix integer(ik),intent(in) :: N real(rk),dimension(:,:),intent(in) :: field -#ifdef TROVE_USE_MPI_ character(len=4) :: jchar character(len=cl) :: filename - character(len=cl) :: job_is - type(MPI_File) :: chkptMPIIO - integer(kind=MPI_OFFSET_KIND) :: offset + character(len=cl) :: job_id + class(ioHandlerBase), allocatable :: ioHandler integer :: ierr ! - write(job_is,"('single swap_matrix')") + write(job_id,"('single swap_matrix')") ! - !!call IOStart(trim(job_is),chkptIO) + !!call IOStart(trim(job_id),chkptIO) ! write(jchar, '(i4)') islice ! filename = trim(suffix)//trim(adjustl(jchar))//'.chk' ! - offset = 0 - call MPI_File_open(mpi_comm_world, filename, mpi_mode_wronly+mpi_mode_create, mpi_info_null, chkptMPIIO, ierr) - call MPI_File_set_size(chkptMPIIO, offset, ierr) - call mpi_barrier(mpi_comm_world, ierr) + call openFile(ioHandler, filename, err, action='write', & + form='unformatted',position='rewind',status='replace') + HANDLE_ERROR(err) + + call ioHandler%write(name) + call ioHandler%write(field, N) + call ioHandler%write(name) + + deallocate(ioHandler) + end subroutine write_divided_slice + + subroutine divided_slice_open(islice,ioHandler,name,suffix) ! - if(mpi_rank .eq. 0) call MPI_File_write_shared(chkptMPIIO,name,len(trim(name)),mpi_character,mpi_status_ignore,ierr) - call mpi_barrier(mpi_comm_world, ierr) + implicit none + integer(ik),intent(in) :: islice + class(ioHandlerBase),intent(inout), allocatable :: ioHandler + character(len=*),intent(in) :: name,suffix + + character(len=4) :: jchar + character(len=cl) :: buf,filename,job_id + integer(ik) :: ilen + integer :: ierr ! - call co_write_matrix_distr(field, N, co_startdim, co_enddim,chkptMPIIO) + if (.not.job%IOmatelem_split) return ! - if(mpi_rank .eq. 0) call MPI_File_write_shared(chkptMPIIO,name,len(trim(name)),mpi_character,mpi_status_ignore,ierr) + write(job_id,"('single swap_matrix')") ! - call MPI_File_close(chkptMPIIO, ierr) + !!call IOStart(trim(job_is),chkptIO) ! -#endif - end subroutine write_divided_slice_mpi - - - subroutine divided_slice_open(islice,chkptIO,name,suffix) - ! - implicit none - integer(ik),intent(in) :: islice - integer(ik),intent(inout) :: chkptIO - character(len=*),intent(in) :: name,suffix - character(len=4) :: jchar - character(len=cl) :: buf,filename,job_is - integer(ik) :: ilen - logical :: ifopened - ! - if (.not.job%IOmatelem_split) return - ! - write(job_is,"('single swap_matrix')") - ! - call IOStart(trim(job_is),chkptIO) - ! - write(jchar, '(i4)') islice - ! - filename = trim(suffix)//trim(adjustl(jchar))//'.chk' - ! - open(chkptIO,form='unformatted',action='read',position='rewind',status='old',file=filename,err=10) - ! - ilen = LEN_TRIM(name) - ! - read(chkptIO) buf(1:ilen) - if ( trim(buf(1:ilen))/=trim(name) ) then - write (out,"(' kinetic checkpoint slice ',a20,': header is missing or wrong',a)") filename,buf(1:ilen) - stop 'PTrestore_rot_kinetic_matrix_elements - in slice - header missing or wrong' - end if - ! - return - ! - 10 write(out,"('divided_slice_open-error: The split-file ',a,' does not exist')") trim(filename) - stop 'divided_slice_open-error: The split-file does not exist' - ! - end subroutine divided_slice_open + write(jchar, '(i4)') islice + ! + filename = trim(suffix)//trim(adjustl(jchar))//'.chk' - subroutine divided_slice_open_mpi(islice,chkptIO,name,suffix) - ! - implicit none - integer(ik),intent(in) :: islice - type(MPI_File),intent(inout) :: chkptIO - character(len=*),intent(in) :: name,suffix + call openFile(ioHandler, filename, err, action='read', & + form='unformatted',position='rewind',status='old') + HANDLE_ERROR(err) + ! + ilen = LEN_TRIM(name) -#ifdef TROVE_USE_MPI_ - character(len=4) :: jchar - character(len=cl) :: buf,filename,job_is - integer(ik) :: ilen - integer :: ierr - ! - if (.not.job%IOmatelem_split) return - ! - write(job_is,"('single swap_matrix')") - ! - !!call IOStart(trim(job_is),chkptIO) - ! - write(jchar, '(i4)') islice - ! - filename = trim(suffix)//trim(adjustl(jchar))//'.chk' - ! - call MPI_File_open(mpi_comm_world, filename, mpi_mode_rdonly, mpi_info_null, chkptMPIIO, ierr) - if (ierr.ne.0) then - if(mpi_rank.eq.0) write(out,"('divided_slice_open-error: The split-file ',a,' does not exist')") trim(filename) - stop 'divided_slice_open-error: The split-file does not exist' - endif - ! - ilen = LEN_TRIM(name) - ! - if (mpi_rank.eq.0) then - call MPI_File_read_shared(chkptIO, buf, ilen, mpi_character, mpi_status_ignore, ierr) + call ioHandler%read(buf(1:ilen)) if ( trim(buf(1:ilen))/=trim(name) ) then write (out,"(' kinetic checkpoint slice ',a20,': header is missing or wrong',a)") filename,buf(1:ilen) +#ifdef TROVE_USE_MPI_ call MPI_Abort(mpi_comm_world, 1) - !stop 'PTrestore_rot_kinetic_matrix_elements - in slice - header missing or wrong' - endif - endif #endif - end subroutine divided_slice_open_mpi - ! - subroutine divided_slice_close(islice,chkptIO,name) - ! - integer(ik),intent(in) :: islice - integer(ik),intent(inout) :: chkptIO - character(len=*),intent(in) :: name - character(len=4) :: jchar - character(len=cl) :: buf,filename - integer(ik) :: ilen - logical :: ifopened - ! - if (.not.job%IOmatelem_split) return - ! - ilen = LEN_TRIM(name) - ! - read(chkptIO) buf(1:ilen) - if ( trim(buf(1:ilen))/=trim(name) ) then - write (out,"(' divided_slice_close, kinetic checkpoint slice ',a,': footer is missing or wrong',a)") filename,buf(1:ilen) - stop 'divided_slice_close - in slice - footer missing or wrong' - end if - ! - close(chkptIO) - ! - end subroutine divided_slice_close + stop 'PTrestore_rot_kinetic_matrix_elements - in slice - header missing or wrong' + endif + end subroutine divided_slice_open - subroutine divided_slice_close_mpi(islice,chkptIO,name) - ! - integer(ik),intent(in) :: islice - type(MPI_File),intent(inout) :: chkptIO - character(len=*),intent(in) :: name - character(len=cl) :: buf - integer(ik) :: ilen - integer :: ierr - ! - if (.not.job%IOmatelem_split) return - ! - ilen = LEN_TRIM(name) - ! - if(mpi_rank .eq. 0) then -#ifdef TROVE_USE_MPI_ - call MPI_File_read_shared(chkptIO, buf, ilen, mpi_character, mpi_status_ignore, ierr) + subroutine divided_slice_close(islice,ioHandler,name) + + integer(ik),intent(in) :: islice + class(ioHandlerBase),intent(inout), allocatable :: ioHandler + character(len=*),intent(in) :: name + character(len=cl) :: buf + integer(ik) :: ilen + integer :: ierr + ! + if (.not.job%IOmatelem_split) return + ! + ilen = LEN_TRIM(name) + ! + call ioHandler%read(buf(1:ilen)) if ( trim(buf(1:ilen))/=trim(name) ) then - write (out,"(' divided_slice_close, kinetic checkpoint slice ',a,': footer is missing or wrong',a)") trim(name),buf(1:ilen) - call MPI_Abort(mpi_comm_world, 1) - !stop 'divided_slice_close - in slice - footer missing or wrong' - endif -#endif - endif - ! + write (out,"(' divided_slice_close, kinetic checkpoint slice ',a,': footer is missing or wrong',a)") trim(name),buf(1:ilen) #ifdef TROVE_USE_MPI_ - call MPI_File_close(chkptIO, ierr) + call MPI_Abort(mpi_comm_world, 1) #endif - ! - end subroutine divided_slice_close_mpi + stop 'PTrestore_rot_kinetic_matrix_elements - in slice - footer missing or wrong' + endif + + deallocate(ioHandler) + end subroutine divided_slice_close ! ! This procedure is thought to make the calculations of the contracted mat. elements @@ -17860,11 +16699,10 @@ subroutine PTcontracted_matelem_class_fast_II(jrot) job_is ='Vib. matrix elements of the rot. kinetic part' call IOStart(trim(job_is),chkptIO) - allocate(ioHandlerFTN::ioHandler) - call ioHandler%open( & - job%kinetmat_file, err, & - action='write', position='rewind', status='replace', form='unformatted') + call openFile(ioHandler, job%kinetmat_file, err, action='write', & + position='rewind', status='replace', form='unformatted') HANDLE_ERROR(err) + call ioHandler%write('Start Kinetic part') ! ! store the bookkeeping information about the contr. basis set @@ -29434,6 +28272,8 @@ subroutine PTTest_eigensolution(jmin,jmax) integer(ik) :: iunit integer(ik) :: ncontr,maxcontr character(len=cl) :: task + ! + class(ioHandlerBase), allocatable :: ioHandler if (job%verbose>=2) write (out,"(//'Test the read/write procedure of the eigenvectors and all auxiliary information.')") @@ -29812,7 +28652,7 @@ subroutine PTTest_eigensolution(jmin,jmax) !call restore_vib_matrix_elements ! task = 'top' - call PTrestore_rot_kinetic_matrix_elements(j,task,iunit,ncontr,maxcontr) + call PTrestore_rot_kinetic_matrix_elements(j,task,ioHandler,ncontr,maxcontr) ! ! obtain the rotational matrix elements ! @@ -34428,7 +33268,7 @@ subroutine partfunc_matexp_taylor(dimen,m,norm_thresh,max_deg,spur_thresh,max_or ! ! perform squaring ! - write(out, '(/1x, a/1x, a, 1x, a)'), 'perform squaring', 'deg of 2', 'norm' + write(out, '(/1x, a/1x, a, 1x, a)') 'perform squaring', 'deg of 2', 'norm' ! if (job%verbose>=2) call TimerStart('Partition function my mat-exp') ! @@ -34443,13 +33283,13 @@ subroutine partfunc_matexp_taylor(dimen,m,norm_thresh,max_deg,spur_thresh,max_or ! do if (deg > max_deg) then - write(out, '(/1x, a, 1x, i3, 1x, a)'), 'max degree of 2', max_deg, 'is reached' + write(out, '(/1x, a, 1x, i3, 1x, a)') 'max degree of 2', max_deg, 'is reached' exit end if ! norm = norm / real(2**deg, kind = rk) ! - write(out, '(1x, i3, 1x, es16.8)'), deg, norm + write(out, '(1x, i3, 1x, es16.8)') deg, norm ! if (abs(norm) <= norm_thresh) exit deg = deg + 1 @@ -34507,14 +33347,14 @@ subroutine partfunc_matexp_taylor(dimen,m,norm_thresh,max_deg,spur_thresh,max_or spur = real(dimen, kind = rk) spur0 = spur ! - write(out, '(/1x, a/1x, a, 13x, a)'), 'compute exponential', 'ord', 'spur' + write(out, '(/1x, a/1x, a, 13x, a)') 'compute exponential', 'ord', 'spur' ! ! loop over Taylor series ! do iorder = iorder + 1 if (iorder > max_order) then - write(out, '(/1x, a, 1x, i3, 1x, a)'), 'max exp degree', max_order, 'is reached' + write(out, '(/1x, a, 1x, i3, 1x, a)') 'max exp degree', max_order, 'is reached' exit end if ! @@ -34751,7 +33591,7 @@ subroutine partfunc_matexp_taylor(dimen,m,norm_thresh,max_deg,spur_thresh,max_or end do !$omp end parallel do ! - write(out, '(1x, i3, 1x, es16.8)'), iorder, spur + write(out, '(1x, i3, 1x, es16.8)') iorder, spur ! if (abs(spur - spur0) <= spur_thresh) exit spur0 = spur @@ -34859,7 +33699,7 @@ subroutine PTstore_icontr_cnu(Maxcontracts,ioHandler,dir) if (Maxcontracts/=ncontr) then write (out,"(' Vib. kinetic checkpoint file ',a)") job%kinetmat_file write (out,"(' Actual and stored basis sizes at J=0 do not agree ',2i8)") PT%Maxcontracts,ncontr - stop 'PTrestore_rot_kinetic_matrix_elements - in file - illegal nroots ' + stop 'PTstore_icontr_cnu - in file - illegal nroots ' end if ! allocate (imat_t(0:PT%Nclasses,ncontr),stat=alloc) @@ -34868,7 +33708,7 @@ subroutine PTstore_icontr_cnu(Maxcontracts,ioHandler,dir) call ioHandler%read(buf18(1:10)) if (buf18(1:10)/='icontr_cnu') then write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_cnu is missing ',a)") job%kinetmat_file,buf18(1:10) - stop 'PTrestore_rot_kinetic_matrix_elements - in file - icontr_cnu missing' + stop 'PTstore_icontr_cnu - in file - icontr_cnu missing' end if ! call ioHandler%read(imat_t(0:PT%Nclasses,1:ncontr)) @@ -34876,7 +33716,7 @@ subroutine PTstore_icontr_cnu(Maxcontracts,ioHandler,dir) call ioHandler%read(buf18(1:11)) if (buf18(1:11)/='icontr_ideg') then write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_ideg is missing ',a)") job%kinetmat_file,buf18(1:11) - stop 'PTrestore_rot_kinetic_matrix_elements - in file - icontr_ideg missing' + stop 'PTstore_icontr_cnu - in file - icontr_ideg missing' end if ! call ioHandler%read(imat_t(0:PT%Nclasses,1:ncontr)) @@ -38280,10 +37120,8 @@ subroutine PTcontracted_matelem_class_fast(jrot) job_is ='Vib. matrix elements of the rot. kinetic part' call IOStart(trim(job_is),chkptIO) ! - allocate(ioHandlerFTN::ioHandler) - call ioHandler%open(& - job%kinetmat_file, err, & - action='write', position='rewind', status='replace', form='unformatted') + call openFile(ioHandler, job%kinetmat_file, err, action='write', & + position='rewind', status='replace', form='unformatted') HANDLE_ERROR(err) call ioHandler%write('Start Kinetic part') diff --git a/test/regression/scripts/H2CO/compare_results.sh b/test/regression/scripts/H2CO/compare_results.sh index 1c1ffdc..e4cd39a 100755 --- a/test/regression/scripts/H2CO/compare_results.sh +++ b/test/regression/scripts/H2CO/compare_results.sh @@ -11,12 +11,13 @@ quantum_files="eigen_descr0_1.chk eigen_descr0_2.chk eigen_descr0_3.chk eigen_de python compare_results.py --kind quantum --folder1 "$folder1" --folder2 "$folder2" $quantum_files -if [[ ${USE_MPI} -ne 1 ]]; then - python compare_results.py --kind intensity --precision 5e-6 --folder1 "$folder1" --folder2 "$folder2" file_intensity.out -fi - python compare_results.py --kind column --column 3 --precision 5e-3 --folder1 "$folder1" --folder2 "$folder2" external.chk python compare_results.py --kind column --column 2 --precision 1e-8 --folder1 "$folder1" --folder2 "$folder2" potential.chk python compare_results.py --kind column --column 4 --precision 1e-8 --folder1 "$folder1" --folder2 "$folder2" kinetic.chk + +if [[ ${USE_MPI} -ne 1 ]]; then + python compare_results.py --kind intensity --precision 5e-6 --folder1 "$folder1" --folder2 "$folder2" file_intensity.out +fi + diff --git a/test/regression/scripts/H2CO/run_benchmark.sh b/test/regression/scripts/H2CO/run_benchmark.sh index ad2b899..31e0704 100755 --- a/test/regression/scripts/H2CO/run_benchmark.sh +++ b/test/regression/scripts/H2CO/run_benchmark.sh @@ -9,18 +9,19 @@ exe=$2 ## SYSTEM OPTIONS -export OMP_NUM_THREADS=$nproc - # Ensure stacksize unlimited (for fortran) ulimit -d unlimited -if [[ ${USE_MPI} -eq 1 ]]; then - LAUNCH="time mpirun -ppn 1 -np $nproc" +if [[ ${USE_MPI} == 1 ]]; then + echo "MPI enabled" + LAUNCH="time mpirun -np $nproc" ./set_io_format.sh enable - echo "Will run with MPI" + export OMP_NUM_THREADS=1 else + echo "MPI disabled" LAUNCH="time" - echo "Will run without MPI" + ./set_io_format.sh disable + export OMP_NUM_THREADS=$nproc fi echo "Time: `date`" @@ -28,7 +29,7 @@ echo "Current directory: `pwd`" echo "Using ${nproc} process(es)" files_to_check=(file{1..12}) -if [[ ${USE_MPI} -ne 1 ]]; then +if [[ ${USE_MPI} == 0 ]]; then # The intensity file does not work with MPI at the moment files_to_check+=(file_intensity) fi diff --git a/test/unit/test_mpi_io.pf b/test/unit/test_mpi_io.pf index 537c347..7830625 100644 --- a/test/unit/test_mpi_io.pf +++ b/test/unit/test_mpi_io.pf @@ -21,7 +21,7 @@ module test_mpi_io integer,external :: INDXL2G logical :: is_mpi_initialised = .false. - integer, parameter :: totalTestCount = 2 + integer, parameter :: totalTestCount = 7 ! CHANGE ME TO NUMBER OF TESTS integer :: currentTestCount = 0 contains @@ -43,18 +43,11 @@ module test_mpi_io end subroutine tearDown @test - subroutine testMPIWriting(this) + subroutine testMPIWritingScalars(this) class(TestMPI), intent(inout) :: this type(ioHandlerMPI) :: ioHandler - integer, parameter :: array2DNRow = 4 - integer, parameter :: array2DNCol = 3 - real(rk), allocatable :: array2D(:,:) - real(rk) :: in_array2D(array2DNRow,array2DNCol) - integer :: array2D_descr(9) = 0 - type(MPI_Datatype) :: array2D_block_type - real :: true_real = 4.0, in_real complex :: true_complex = (5.0, 1.0), in_complex integer :: true_integer = 5, in_integer @@ -67,15 +60,15 @@ module test_mpi_io character(len=*), parameter :: form = "unformatted" character(len=*), parameter :: access = "sequential" - integer i, j - integer :: gi = 0, gj = 0, MB, NB, RSRC, CSRC type(ErrorType) :: err integer :: ierr, rank, allocinfo = 0 + ! Set up MPI call MPI_Comm_rank(MPI_COMM_WORLD, rank, ierr) if(ierr.ne.0) print *, "Error: could not get rank" + ! Write call ioHandler%open(fname, err, action='write', & form=form, access=access, status=status, position=position) HANDLE_ERROR(err) @@ -86,7 +79,65 @@ module test_mpi_io call ioHandler%write(true_complex) ! complex call ioHandler%write(true_str) ! string - ! Test writing an array + call ioHandler%close() + + ! Only test result on main process + if(rank == 0) then + open(newunit=iounit, iostat=stat, action='read', file=fname, & + form=form, access=access, status=status, position=position) + + read(iounit) in_integer + @assertTrue(in_integer == true_integer) + + read(iounit) in_real + @assertTrue(in_real == true_real) + + read(iounit) in_complex + @assertTrue(in_complex == true_complex) + + read(iounit) in_str + @assertTrue(in_str == true_str) + + ! Remove test file + if (stat == 0) close(iounit, status='delete') + + endif + + end subroutine + + @test + subroutine testMPIWritingBlacsDistArray(this) + class(TestMPI), intent(inout) :: this + + type(ioHandlerMPI) :: ioHandler + + integer, parameter :: array2DNRow = 4 + integer, parameter :: array2DNCol = 3 + real(rk), allocatable :: array2D(:,:) + real(rk) :: in_array2D(array2DNRow,array2DNCol) + integer :: array2D_descr(9) = 0 + type(MPI_Datatype) :: array2D_block_type + + integer :: true_integer = 5, in_integer + + integer :: iounit, stat + character(len=*), parameter :: fname = "test.dat" + character(len=*), parameter :: position = "asis" + character(len=*), parameter :: status = "unknown" + character(len=*), parameter :: form = "unformatted" + character(len=*), parameter :: access = "sequential" + + integer i, j + integer :: gi = 0, gj = 0, MB, NB, RSRC, CSRC + type(ErrorType) :: err + + integer :: ierr, rank, allocinfo = 0 + + ! Set up MPI + call MPI_Comm_rank(MPI_COMM_WORLD, rank, ierr) + if(ierr.ne.0) print *, "Error: could not get rank" + + ! Set up 2D Array call co_block_type_init(array2D, array2DNCol, array2DNRow, array2D_descr, allocinfo, array2D_block_type) if(allocinfo.ne.0) print *, "ERROR: couldn't allocate array" @@ -104,17 +155,20 @@ module test_mpi_io end do end do + ! Write + call ioHandler%open(fname, err, action='write', & + form=form, access=access, status=status, position=position) + HANDLE_ERROR(err) + + ! Test writing an array call ioHandler%write(array2D, array2D_descr, array2D_block_type) ! Test writing something after array - call ioHandler%write(true_integer) ! int + call ioHandler%write(true_integer) ! Test writing another array call ioHandler%write(array2D, array2D_descr, array2D_block_type) - ! Test writing something after 2nd array - call ioHandler%write(true_integer) ! int - call ioHandler%close() ! Only test result on main process @@ -122,19 +176,6 @@ module test_mpi_io open(newunit=iounit, iostat=stat, action='read', file=fname, & form=form, access=access, status=status, position=position) - read(iounit) in_integer - @assertTrue(in_integer == true_integer) - - read(iounit) in_real - @assertTrue(in_real == true_real) - - read(iounit) in_complex - @assertTrue(in_complex == true_complex) - - read(iounit) in_str - print *, sizeof(in_str) - @assertTrue(in_str == true_str) - read(iounit) in_array2D do i=1,array2DNRow do j=1,array2DNCol @@ -152,24 +193,21 @@ module test_mpi_io end do end do - read(iounit) in_integer - @assertTrue(in_integer == true_integer) - ! Remove test file if (stat == 0) close(iounit, status='delete') endif - end subroutine testMPIWriting + end subroutine @test subroutine testMPIWritingColumnDistArray(this) class(TestMPI), intent(inout) :: this type(ioHandlerMPI) :: ioHandler - integer :: dimen = 12 - integer :: startdim, enddim, blocksize, mdimen_p, mdimen_b, mdimen + integer, parameter :: dimen = 57 integer :: b, icoeff, jcoeff + integer :: startdim_, enddim_, blocksize_ integer :: iounit, stat character(len=*), parameter :: fname = "test.dat" @@ -178,9 +216,9 @@ module test_mpi_io character(len=*), parameter :: form = "unformatted" character(len=*), parameter :: access = "sequential" - real(rk), allocatable :: grot_t(:,:) - real(rk), allocatable :: grot_full(:,:) - real(rk),allocatable :: recvbuf(:,:,:) + real(rk) :: grot_full(dimen, dimen) + real(rk), pointer :: grot_t(:,:) + real(rk), allocatable :: recvbuf(:,:,:) integer :: true_integer = 5, in_integer @@ -194,23 +232,17 @@ module test_mpi_io call MPI_Comm_rank(MPI_COMM_WORLD, rank, ierr) if(ierr.ne.0) print *, "Error: could not get rank" - call co_init_distr(dimen, startdim, enddim, blocksize) - - mdimen = dimen - mdimen_p = int(1+real(dimen/comm_size)) - mdimen_b = comm_size*mdimen_p + call co_init_distr(dimen, startdim_, enddim_, blocksize_) - allocate(recvbuf(mdimen_p,mdimen_p,comm_size)) - allocate(grot_t(mdimen_b,startdim:startdim+mdimen_p-1)) - allocate(grot_full(dimen, dimen)) + call co_create_distr_array(grot_t, dimen) grot_t = 0 ! Fill local chunk of symmetric matrix do b=1,comm_size if (send_or_recv(b).ge.0) then - do icoeff=startdim,enddim - do jcoeff=((b-1)*mdimen_p)+1,b*mdimen_p + do icoeff=co_startdim,co_enddim + do jcoeff=((b-1)*co_localsize)+1,b*co_localsize grot_t(jcoeff,icoeff) = jcoeff*icoeff enddo enddo @@ -218,7 +250,8 @@ module test_mpi_io enddo ! Distribute between processes and save - call co_distr_data(grot_t, recvbuf, mdimen_p, startdim, enddim) + allocate(recvbuf(co_localsize,co_localsize,comm_size)) + call co_distr_data(grot_t, recvbuf, co_localsize, co_startdim, co_enddim) call ioHandler%open(fname, err, action='write', & form=form, access=access, status=status, position=position) @@ -228,13 +261,13 @@ module test_mpi_io call ioHandler%write(true_integer) ! Test writing an array - call ioHandler%write(grot_t, mdimen) + call ioHandler%write(grot_t, dimen) ! Test writing something after array call ioHandler%write(true_integer) ! Test writing another array - call ioHandler%write(grot_t, mdimen) + call ioHandler%write(grot_t, dimen) ! Test writing something after second array call ioHandler%write(true_integer) @@ -273,4 +306,347 @@ module test_mpi_io endif end subroutine testMPIWritingColumnDistArray + + @test + subroutine testMPIReadScalars(this) + class(TestMPI), intent(inout) :: this + type(ioHandlerMPI) :: ioHandler + + real :: true_real = 4.0, in_real + complex :: true_complex = (5.0, 1.0), in_complex + integer :: true_integer = 5, in_integer + character (len=11) :: true_str = "test string", in_str + + integer :: iounit, stat + character(len=*), parameter :: fname = "test.dat" + character(len=*), parameter :: position = "asis" + character(len=*), parameter :: status = "unknown" + character(len=*), parameter :: form = "unformatted" + character(len=*), parameter :: access = "sequential" + + integer :: ierr, rank, allocinfo = 0 + type(ErrorType) :: err + + ! Set up MPI + call MPI_Comm_rank(MPI_COMM_WORLD, rank, ierr) + if(ierr.ne.0) print *, "Error: could not get rank" + + ! Write test file + if(rank == 0) then + open(newunit=iounit, iostat=stat, action='write', file=fname, & + form=form, access=access, status=status, position=position) + + write(iounit) true_integer ! int + write(iounit) true_real ! double + write(iounit) true_complex ! complex + write(iounit) true_str ! string + + if (stat == 0) close(iounit) + endif + + call MPI_Barrier(MPI_COMM_WORLD) + + ! Read test file + call ioHandler%open(fname, err, \ + action='read', form=form, access=access, status=status, position=position) + HANDLE_ERROR(err) + + call ioHandler%read(in_integer) + @assertTrue(in_integer == true_integer) + + call ioHandler%read(in_real) + @assertTrue(in_real == true_real) + + call ioHandler%read(in_complex) + @assertTrue(in_complex == true_complex) + + call ioHandler%read(in_str) + @assertTrue(in_str == true_str) + + call ioHandler%close() + + if(rank == 0) then + ! Cleanup test file + open(newunit=iounit, iostat=stat, action='read', file=fname) + if (stat == 0) close(iounit, status='delete') + endif + + end subroutine + + @test + subroutine testMPIReadColumnDistArray(this) + class(TestMPI), intent(inout) :: this + type(ioHandlerMPI) :: ioHandler + + integer, parameter :: dimen = 57 + integer :: b, icoeff, jcoeff + integer :: startdim_, enddim_, blocksize_ + + real(rk) :: true_array(dimen,dimen) + real(rk), pointer :: in_array(:,:) + + integer :: true_integer = 5, in_integer + + integer :: iounit, stat + character(len=*), parameter :: fname = "test.dat" + character(len=*), parameter :: position = "asis" + character(len=*), parameter :: status = "unknown" + character(len=*), parameter :: form = "unformatted" + character(len=*), parameter :: access = "sequential" + + integer i, j + integer :: ierr, rank, allocinfo = 0 + type(ErrorType) :: err + + ! Set up array + do i=1,dimen + do j=1,dimen + true_array(i,j) = i*j + end do + end do + + ! Set up MPI + call MPI_Comm_rank(MPI_COMM_WORLD, rank, ierr) + if(ierr.ne.0) print *, "Error: could not get rank" + + ! Set up 2D array block type + call co_init_distr(dimen, startdim_, enddim_, blocksize_) + call co_create_distr_array(in_array, dimen) + + ! Write test file + if(rank == 0) then + open(newunit=iounit, iostat=stat, action='write', file=fname, & + form=form, access=access, status=status, position=position) + + write(iounit) true_integer ! int + write(iounit) true_array ! 2D array + write(iounit) true_integer ! int + write(iounit) true_array ! 2D array + write(iounit) true_integer ! int + + if (stat == 0) close(iounit) + endif + + call MPI_Barrier(MPI_COMM_WORLD) + + ! Read test file + call ioHandler%open(fname, err, \ + action='read', form=form, access=access, status=status, position=position) + HANDLE_ERROR(err) + + call ioHandler%read(in_integer) + @assertTrue(in_integer == true_integer) + + call ioHandler%read(in_array, dimen) + do b=1,comm_size + if (send_or_recv(b).ge.0) then + do icoeff=co_startdim,co_enddim + do jcoeff=co_startdim,co_enddim + @assertTrue(in_array(jcoeff,icoeff) == true_array(jcoeff,icoeff)) + enddo + enddo + endif + enddo + + call ioHandler%read(true_integer) + @assertTrue(in_integer == true_integer) + + call ioHandler%read(in_array, dimen) + do b=1,comm_size + if (send_or_recv(b).ge.0) then + do icoeff=co_startdim,co_enddim + do jcoeff=co_startdim,co_enddim + @assertTrue(in_array(jcoeff,icoeff) == true_array(jcoeff,icoeff)) + enddo + enddo + endif + enddo + + call ioHandler%read(true_integer) + @assertTrue(in_integer == true_integer) + + call ioHandler%close() + + if(rank == 0) then + ! Cleanup test file + open(newunit=iounit, iostat=stat, action='read', file=fname) + if (stat == 0) close(iounit, status='delete') + endif + + end subroutine + + @test + subroutine testMPIReadBlacsDistArray(this) + class(TestMPI), intent(inout) :: this + type(ioHandlerMPI) :: ioHandler + + integer, parameter :: array2DNRow = 13 + integer, parameter :: array2DNCol = 13 + real(rk) :: array2D(array2DNRow,array2DNCol) + real(rk), allocatable :: in_array2D(:,:) + integer :: array2D_descr(9) = 0 + type(MPI_Datatype) :: array2D_block_type + + integer :: true_integer = 5, in_integer + + integer :: iounit, stat + character(len=*), parameter :: fname = "test.dat" + character(len=*), parameter :: position = "asis" + character(len=*), parameter :: status = "unknown" + character(len=*), parameter :: form = "unformatted" + character(len=*), parameter :: access = "sequential" + + integer i, j + integer :: gi = 0, gj = 0, MB, NB, RSRC, CSRC + integer :: ierr, rank, allocinfo = 0 + type(ErrorType) :: err + + ! Set up 2D array + do i=1,array2DNRow + do j=1,array2DNCol + array2D(i,j) = array2DNCol*(i-1) + j + end do + end do + + ! Set up MPI + call MPI_Comm_rank(MPI_COMM_WORLD, rank, ierr) + if(ierr.ne.0) print *, "Error: could not get rank" + + ! Set up 2D array block type + call co_block_type_init(in_array2D, array2DNCol, array2DNRow, array2D_descr, allocinfo, array2D_block_type) + if(allocinfo.ne.0) print *, "ERROR: couldn't allocate array" + + MB = array2D_descr(5) + NB = array2D_descr(6) + + RSRC = array2D_descr(7) + CSRC = array2D_descr(8) + + ! Write test file + if(rank == 0) then + open(newunit=iounit, iostat=stat, action='write', file=fname, & + form=form, access=access, status=status, position=position) + + write(iounit) true_integer ! int + write(iounit) array2D ! 2D array + write(iounit) true_integer ! int + write(iounit) array2D ! 2D array + write(iounit) true_integer ! int + + if (stat == 0) close(iounit) + endif + + call MPI_Barrier(MPI_COMM_WORLD) + + ! Read test file + call ioHandler%open(fname, err, \ + action='read', form=form, access=access, status=status, position=position) + HANDLE_ERROR(err) + + call ioHandler%read(in_integer) + @assertTrue(in_integer == true_integer) + + call ioHandler%read(in_array2D, array2D_descr, array2D_block_type) + do i=1,size(in_array2D,1) + do j=1,size(in_array2D,2) + gi = INDXL2G (i, MB, myprow, RSRC, nprow) + gj = INDXL2G (j, NB, mypcol, CSRC, npcol) + @assertTrue(in_array2D(i,j) == array2DNCol*(gi-1) + (gj)) + end do + end do + + call ioHandler%read(true_integer) + @assertTrue(in_integer == true_integer) + + call ioHandler%read(in_array2D, array2D_descr, array2D_block_type) + do i=1,size(in_array2D,1) + do j=1,size(in_array2D,2) + gi = INDXL2G (i, MB, myprow, RSRC, nprow) + gj = INDXL2G (j, NB, mypcol, CSRC, npcol) + @assertTrue(in_array2D(i,j) == array2DNCol*(gi-1) + (gj)) + end do + end do + + call ioHandler%read(true_integer) + @assertTrue(in_integer == true_integer) + + call ioHandler%close() + + if(rank == 0) then + ! Cleanup test file + open(newunit=iounit, iostat=stat, action='read', file=fname) + if (stat == 0) close(iounit, status='delete') + endif + + end subroutine + + @test + subroutine testSeekWhileReading(this) + class(TestMPI), intent(inout) :: this + type(ioHandlerMPI) :: ioHandler + + real :: true_real = 4.0, in_real + integer :: dummy_arr(1:5) + integer :: true_integer = 5, in_integer + + integer :: iounit, stat + character(len=*), parameter :: fname = "test.dat" + character(len=*), parameter :: position = "asis" + character(len=*), parameter :: status = "unknown" + character(len=*), parameter :: form = "unformatted" + character(len=*), parameter :: access = "sequential" + + integer i + + integer :: ierr, rank, allocinfo = 0 + type(ErrorType) :: err + + integer(kind=4) :: offset + + ! Set up MPI + call MPI_Comm_rank(MPI_COMM_WORLD, rank, ierr) + if(ierr.ne.0) print *, "Error: could not get rank" + + ! Setup dummy array + do i=1,5 + dummy_arr(i) = 10 + enddo + + ! Write test file + if(rank == 0) then + open(newunit=iounit, iostat=stat, action='write', file=fname, & + form=form, access=access, status=status, position=position) + + write(iounit) true_real ! double + write(iounit) dummy_arr + write(iounit) true_integer ! int + + if (stat == 0) close(iounit) + endif + + call MPI_Barrier(MPI_COMM_WORLD) + + ! Read test file + call ioHandler%open(fname, err, \ + action='read', form=form, access=access, status=status, position=position) + HANDLE_ERROR(err) + + call ioHandler%read(in_real) + @assertTrue(in_real == true_real) + + offset = size(dummy_arr)*sizeof(0) + call ioHandler%seek(offset) + + call ioHandler%read(in_integer) + @assertTrue(in_integer == true_integer) + + call ioHandler%close() + + if(rank == 0) then + ! Cleanup test file + open(newunit=iounit, iostat=stat, action='read', file=fname) + if (stat == 0) close(iounit, status='delete') + endif + end subroutine + end module diff --git a/tran.f90 b/tran.f90 index db13c45..d61f630 100644 --- a/tran.f90 +++ b/tran.f90 @@ -22,6 +22,7 @@ module tran #ifdef TROVE_USE_MPI_ use io_handler_mpi #endif + use io_factory use errors private @@ -1323,10 +1324,11 @@ subroutine TRconvert_matel_j0_eigen(jrot) integer(ik),allocatable :: ijterm(:,:) double precision,parameter :: alpha = 1.0d0,beta=0.0d0 character(len=cl) :: jchar,filename - type(MPI_File) :: fileh, fileh_w - integer(kind=MPI_OFFSET_KIND) :: mpioffset,read_offset,write_offset integer :: ierr - class(ioHandlerBase), allocatable :: ioHandler + class(ioHandlerBase), allocatable :: kineteigenHandler + class(ioHandlerBase), allocatable :: kinetmatHandler + class(ioHandlerBase), allocatable :: extFmatHandler + class(ioHandlerBase), allocatable :: exteigenHandler type(ErrorType) :: err type(MPI_Datatype) :: gmat_block_type, psi_block_type, mat_t_block_type, mat_s_block_type, extF_block_type @@ -1391,10 +1393,10 @@ subroutine TRconvert_matel_j0_eigen(jrot) else call co_block_type_init(psi, Neigenroots, dimen, desc_psi, info) call ArrayStart('psi',info,1,kind(psi),int(size(psi),hik)) - + !shape(psi_t) == shape(psi^T) == shape(mat_t) call co_block_type_init(psi_t, dimen, Neigenroots, desc_mat_t, info) - call ArrayStart('psi_t',info,1,kind(mat_t),int(size(mat_t),hik)) + call ArrayStart('psi_t',info,1,kind(psi_t),int(size(psi_t),hik)) call co_block_type_init(mat_t, dimen, Neigenroots, desc_mat_t, info) call ArrayStart('mat_t',info,1,kind(mat_t),int(size(mat_t),hik)) endif @@ -1612,22 +1614,16 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! job_is ='Eigen-vib. matrix elements of the rot. kinetic part' ! - call IOStart(trim(job_is),chkptIO) -#ifdef TROVE_USE_MPI_ - allocate(ioHandlerMPI::ioHandler) -#else - allocate(ioHandlerFTN::ioHandler) -#endif - call ioHandler%open(& - job%kineteigen_file, err, & - action='write', position='rewind', status='replace', form='unformatted') + call openFile(kineteigenHandler, job%kineteigen_file, err, action='write', & + position='rewind', status='replace', form='unformatted') HANDLE_ERROR(err) - call ioHandler%write('Start Kinetic part') + + call kineteigenHandler%write('Start Kinetic part') treat_vibration = .false. - call PTstore_icontr_cnu(Neigenroots,ioHandler,job%IOj0matel_action) + call PTstore_icontr_cnu(Neigenroots,kineteigenHandler,job%IOj0matel_action) if (job%vib_rot_contr) then - call ioHandler%write('vib-rot') + call kineteigenHandler%write('vib-rot') endif ! endif @@ -1661,34 +1657,23 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! endif ! - -#ifdef TROVE_USE_MPI_ - call MPI_File_open(mpi_comm_world, job%kinetmat_file, mpi_mode_rdonly, mpi_info_null, fileh, ierr) -#endif ! The eigen-vibrational (J=0) matrix elements of the rotational and coriolis ! kinetic parts are being computed here. ! if (job%verbose>=3) write(out,"(/' Transform grot to J0-repres...')") ! - if (.not.job%IOmatelem_split) then - ! - task = 'rot' - ! - call ioHandler%write('g_rot') - ! - call restore_rot_kinetic_matrix_elements(jrot,treat_vibration,task,iunit) - ! - else - ! + if (job%IOmatelem_split) then task = 'top' - ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call restore_rot_kinetic_matrix_elements_mpi(jrot,treat_vibration,task,fileh) - else - call restore_rot_kinetic_matrix_elements(jrot,treat_vibration,task,iunit) - endif - ! + else + task = 'rot' + call kineteigenHandler%write('g_rot') endif + + call openFile(kinetmatHandler, job%kinetmat_file, err, action='read', & + position='rewind', status='old', form='unformatted') + HANDLE_ERROR(err) + + call restore_rot_kinetic_matrix_elements(jrot,treat_vibration,task,kinetmatHandler) ! if (job%verbose>=5) call TimerStart('J0-convertion for g_rot') ! @@ -1696,13 +1681,6 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! islice = 0 ! - if ((.not.job%IOmatelem_split) .and. blacs_size.gt.1) then -#ifdef TROVE_USE_MPI_ - call MPI_File_get_position(fileh, read_offset, ierr) - call MPI_File_set_view(fileh, read_offset, mpi_byte, gmat_block_type, "native", MPI_INFO_NULL, ierr) -#endif - endif - ! do k1 = 1,3 ! do k2 = 1,3 @@ -1713,25 +1691,15 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! if (job%IOmatelem_split.and..not.job%vib_rot_contr) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call divided_slice_read_mpi(islice,'g_rot',job%matelem_suffix,dimen,gmat,gmat_block_type,ierror) - else - call divided_slice_read(islice,'g_rot',job%matelem_suffix,dimen,gmat,ierror) - endif + call divided_slice_read(islice,'g_rot',job%matelem_suffix,gmat,desc_gmat,gmat_block_type,ierror) ! elseif (job%IOmatelem_split.and.job%vib_rot_contr) then ! - call divided_slice_read_vibrot(islice,job%matelem_suffix,dimen,gmat) + call divided_slice_read_vibrot(islice,job%matelem_suffix,dimen,gmat,desc_gmat,gmat_block_type) ! else ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - call MPI_File_read_all(fileh, gmat, size(gmat), mpi_double_precision, mpi_status_ignore, ierr) -#endif - else - read(iunit) gmat - endif + call kinetmatHandler%read(gmat, desc_gmat, gmat_block_type) ! endif ! @@ -1752,35 +1720,21 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! if (job%IOmatelem_split.and..not.job%vib_rot_contr) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call divided_slice_write_mpi(islice,'g_rot',job%j0matelem_suffix,Neigenroots,mat_s,mat_s_block_type) - else - if(mpi_rank.eq.0) call divided_slice_write(islice,'g_rot',job%j0matelem_suffix,Neigenroots,mat_s) - endif + call divided_slice_write(islice,'g_rot',job%j0matelem_suffix,mat_s, desc_mat_s,mat_s_block_type) ! elseif (job%IOmatelem_split.and.job%vib_rot_contr) then ! - if(mpi_rank.eq.0) call divided_slice_write_vibrot(islice,job%j0matelem_suffix,Neigenroots,mat_s) + call divided_slice_write_vibrot(islice,job%j0matelem_suffix,Neigenroots,mat_s,desc_mat_s,mat_s_block_type) ! else ! - call ioHandler%write(mat_s, desc_mat_s, mat_s_block_type) + call kineteigenHandler%write(mat_s, desc_mat_s, mat_s_block_type) ! endif ! enddo ! enddo - ! - ! Reset view to flat file - if ((.not.job%IOmatelem_split) .and. blacs_size.gt.1) then -#ifdef TROVE_USE_MPI_ - read_offset = read_offset + 9*int(dimen,MPI_OFFSET_KIND)*dimen*mpi_real_size - call MPI_File_set_view(fileh, int(0,MPI_OFFSET_KIND), mpi_byte, mpi_byte, "native", MPI_INFO_NULL, ierr) - call MPI_File_seek(fileh, read_offset, MPI_SEEK_SET) -#endif - endif - ! if (job%verbose>=5) call TimerStop('J0-convertion for g_rot') ! @@ -1790,17 +1744,8 @@ subroutine TRconvert_matel_j0_eigen(jrot) if (.not.job%IOmatelem_split) then ! task = 'cor' - if (trim(job%kinetmat_format).eq.'MPIIO') then - ! -#ifdef TROVE_USE_MPI_ - call restore_rot_kinetic_matrix_elements_mpi(jrot,treat_vibration,task,fileh) - ! - call MPI_Barrier(MPI_COMM_WORLD, ierr) ! May no longer be needed? -#endif - else - call restore_rot_kinetic_matrix_elements(jrot,treat_vibration,task,iunit) - endif - call ioHandler%write('g_cor') + call restore_rot_kinetic_matrix_elements(jrot,treat_vibration,task,kinetmatHandler) + call kineteigenHandler%write('g_cor') ! endif ! @@ -1810,13 +1755,6 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! islice = 9 ! - if ((.not.job%IOmatelem_split) .and. blacs_size.gt.1) then -#ifdef TROVE_USE_MPI_ - call MPI_File_get_position(fileh, read_offset, ierr) - call MPI_File_set_view(fileh, read_offset, mpi_byte, gmat_block_type, "native", MPI_INFO_NULL, ierr) -#endif - endif - ! !do k1 = 1,FLNmodes ! !if (job%contrci_me_fast.and.k1>1) cycle @@ -1829,25 +1767,15 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! if (job%IOmatelem_split.and..not.job%vib_rot_contr) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call divided_slice_read_mpi(islice,'g_cor',job%matelem_suffix,dimen,gmat,gmat_block_type,ierror) - else - call divided_slice_read(islice,'g_cor',job%matelem_suffix,dimen,gmat,ierror) - endif + call divided_slice_read(islice,'g_cor',job%matelem_suffix,gmat,desc_gmat,gmat_block_type,ierror) ! elseif (job%IOmatelem_split.and.job%vib_rot_contr) then ! - call divided_slice_read_vibrot(islice,job%matelem_suffix,dimen,gmat) + call divided_slice_read_vibrot(islice,job%matelem_suffix,dimen,gmat,desc_gmat,gmat_block_type) ! else ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - call MPI_File_read_all(fileh, gmat, size(gmat), mpi_double_precision, mpi_status_ignore, ierr) -#endif - else - read(iunit) gmat - endif + call kinetmatHandler%read(gmat,desc_gmat,gmat_block_type) ! endif ! @@ -1866,72 +1794,39 @@ subroutine TRconvert_matel_j0_eigen(jrot) endif ! ! - if (job%IOmatelem_split.and..not.job%vib_rot_contr) then - ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call divided_slice_write_mpi(islice,'g_cor',job%j0matelem_suffix,Neigenroots,mat_s,mat_s_block_type) + if (job%IOmatelem_split) then + if(job%vib_rot_contr) then + call divided_slice_write_vibrot(islice,job%j0matelem_suffix,Neigenroots,mat_s,desc_mat_s,mat_s_block_type) else - if(mpi_rank.eq.0) call divided_slice_write(islice,'g_cor',job%j0matelem_suffix,Neigenroots,mat_s) + call divided_slice_write(islice,'g_cor',job%j0matelem_suffix,mat_s,desc_mat_s,mat_s_block_type) endif - ! - elseif (job%IOmatelem_split.and.job%vib_rot_contr) then - ! - if(mpi_rank.eq.0) call divided_slice_write_vibrot(islice,job%j0matelem_suffix,Neigenroots,mat_s) - ! else - ! - call ioHandler%write(mat_s, desc_mat_s, mat_s_block_type) - ! + call kineteigenHandler%write(mat_s, desc_mat_s, mat_s_block_type) endif ! enddo ! !enddo ! - ! Reset view to flat file - if ((.not.job%IOmatelem_split) .and. blacs_size.gt.1) then -#ifdef TROVE_USE_MPI_ - read_offset = read_offset + 3*int(dimen,MPI_OFFSET_KIND)*dimen*mpi_real_size - call MPI_File_set_view(fileh, int(0,MPI_OFFSET_KIND), mpi_byte, mpi_byte, "native", MPI_INFO_NULL, ierr) - call MPI_File_seek(fileh, read_offset, MPI_SEEK_SET) -#endif - endif - ! if (job%verbose>=5) call TimerStop('J0-convertion for g_cor') ! - if ((.not.job%IOmatelem_split.or.job%iswap(1)==1).and.(mpi_rank.eq.0)) then - call ioHandler%write('End Kinetic part') - endif - ! - if (.not.job%vib_rot_contr) then - deallocate(ioHandler) + if ((.not.job%IOmatelem_split.or.job%iswap(1)==1)) then + call kineteigenHandler%write('End Kinetic part') endif ! task = 'end' ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call restore_rot_kinetic_matrix_elements_mpi(jrot,treat_vibration,task,fileh) - else - call restore_rot_kinetic_matrix_elements(jrot,treat_vibration,task,iunit) - endif + call restore_rot_kinetic_matrix_elements(jrot,treat_vibration,task,kinetmatHandler) ! if (allocated(gmat)) deallocate(gmat) + if (allocated(kinetmatHandler)) deallocate(kinetmatHandler) call ArrayStop('gmat-fields') ! if (job%verbose>=3) write(out,"(' ...done!')") ! endif ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - call MPI_File_close(fileh, ierr) -#endif - else - ! Should this be close(iunit) instead of close(chkptIO)? - ! In which case, we don't want to deallocate the ioHandler (formerly chkptIO)... - !close(chkptIO,status='keep') - if (allocated(ioHandler)) deallocate(ioHandler) - endif + if (allocated(kineteigenHandler)) deallocate(kineteigenHandler) ! ! External field part ! @@ -1960,49 +1855,22 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! filename = job%extFmat_file ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - ! -#ifdef TROVE_USE_MPI_ - call MPI_File_open(mpi_comm_world, filename, mpi_mode_rdonly, mpi_info_null, fileh, ierr) - ! - call MPI_File_read_all(fileh, buf20, 7, mpi_character, mpi_status_ignore, ierr) - if (buf20(1:7)/='[MPIIO]') then - write (out,"(' Vib. kinetic checkpoint file ',a,' is not an MPIIO file: ',a)") filename, buf20 - stop 'restore_vib_matric_elements - Not an MPIIO file' - end if - ! - call MPI_File_read_all(fileh, buf20, 20, mpi_character, mpi_status_ignore, ierr) - if (buf20/='Start external field') then - write (out,"(' restore_vib_matrix_elements ',a,' has bogus header: ',a)") filename,buf20 - stop 'restore_vib_matrix_elements - bogus file format' - end if - ! - call MPI_File_read_all(fileh, ncontr_t, 1, mpi_integer, mpi_status_ignore, ierr) - ! - if (bset_contr(1)%Maxcontracts/=ncontr_t) then - write (out,"(' Dipole moment checkpoint file ',a)") filename - write (out,"(' Actual and stored basis sizes at J=0 do not agree ',2i8)") bset_contr(1)%Maxcontracts,ncontr_t - stop 'restore_Extvib_matrix_elements - in file - illegal ncontracts ' - end if - ! -#endif - else - open(iunit,form='unformatted',action='read',position='rewind',status='old',file=filename) - ! - read(iunit) buf20 - if (buf20/='Start external field') then - write (out,"(' restore_vib_matrix_elements ',a,' has bogus header: ',a)") filename,buf20 - stop 'restore_vib_matrix_elements - bogus file format' - end if - ! - read(iunit) ncontr_t - if (bset_contr(1)%Maxcontracts/=ncontr_t) then - write (out,"(' Dipole moment checkpoint file ',a)") filename - write (out,"(' Actual and stored basis sizes at J=0 do not agree ',2i8)") bset_contr(1)%Maxcontracts,ncontr_t - stop 'restore_Extvib_matrix_elements - in file - illegal ncontracts ' - end if - ! - endif + call openFile(extFmatHandler, filename, err, action='read', & + form='unformatted',position='rewind',status='old') + HANDLE_ERROR(err) + + call extFmatHandler%read(buf20(1:20)) + if (buf20/='Start external field') then + write (out,"(' restore_vib_matrix_elements ',a,' has bogus header: ',a)") filename,buf20 + stop 'restore_vib_matrix_elements - bogus file format' + end if + + call extFmatHandler%read(ncontr_t) + if (bset_contr(1)%Maxcontracts/=ncontr_t) then + write (out,"(' Dipole moment checkpoint file ',a)") filename + write (out,"(' Actual and stored basis sizes at J=0 do not agree ',2i8)") bset_contr(1)%Maxcontracts,ncontr_t + stop 'restore_Extvib_matrix_elements - in file - illegal ncontracts ' + end if ! endif ! @@ -2013,40 +1881,15 @@ subroutine TRconvert_matel_j0_eigen(jrot) job_is ='external field contracted matrix elements for J=0' call IOStart(trim(job_is),chkptIO) ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - ! - call mpi_file_open(mpi_comm_world, job%exteigen_file, mpi_mode_wronly+mpi_mode_create, mpi_info_null, fileh_w, ierr) - call mpi_file_set_errhandler(fileh_w, mpi_errors_are_fatal) - ! - mpioffset = 0 - call mpi_file_set_size(fileh_w, mpioffset, ierr) - ! - if(mpi_rank.eq.0) then - call mpi_file_write_shared(fileh_w, '[MPIIO]', 7, mpi_character, mpi_status_ignore, ierr) - call mpi_file_write_shared(fileh_w, 'Start external field', 20, mpi_character, mpi_status_ignore, ierr) - endif - ! - ! store the matrix elements - ! - if(mpi_rank.eq.0) call mpi_file_write(fileh_w, neigenroots, 1, mpi_integer, mpi_status_ignore, ierr) - call mpi_barrier(mpi_comm_world, ierr) - call mpi_file_seek(fileh_w, int(0,mpi_offset_kind), mpi_seek_end) - ! -#endif - else - ! - open(chkptio,form='unformatted',action='write',position='rewind',status='replace',file=job%exteigen_file) - write(chkptio) 'Start external field' - ! - ! store the matrix elements - ! - write(chkptio) neigenroots - endif + call openFile(exteigenHandler, job%exteigen_file, err, action='write', & + form='unformatted',position='rewind',status='replace') + HANDLE_ERROR(err) + call exteigenHandler%write('Start external field') + call exteigenHandler%write(neigenroots) ! endif ! - if (trim(job%kinetmat_format).ne.'MPIIO'.and.job%IOextF_divide) close(iunit) + if (job%IOextF_divide) deallocate(exteigenHandler) ! rootsize = int(ncontr_t*(ncontr_t+1)/2,hik) rootsize2= int(ncontr_t*ncontr_t,hik) @@ -2061,27 +1904,13 @@ subroutine TRconvert_matel_j0_eigen(jrot) call ArrayStart('extF_me',info,1,kind(extF_me),int(size(extF_me),hik)) endif ! - if ((.not.job%IOextF_divide) .and. blacs_size.gt.1) then -#ifdef TROVE_USE_MPI_ - call MPI_File_get_position(fileh, read_offset, ierr) - call MPI_File_set_view(fileh, read_offset, mpi_byte, extF_block_type, "native", MPI_INFO_NULL, ierr) - - call MPI_File_get_position_shared(fileh_w, write_offset, ierr) - call MPI_File_set_view(fileh_w, write_offset, mpi_byte, mat_s_block_type, "native", MPI_INFO_NULL, ierr) -#endif - endif - ! do imu = fitting%iparam(1),fitting%iparam(2) ! if (job%verbose>=4) write(out,"(' imu = ',i8)",advance='NO') imu ! if (job%IOextF_divide) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call divided_slice_read_mpi(imu,'extF',job%extmat_suffix,dimen,extF_me,extF_block_type,ierror) - else - call divided_slice_read(imu,'extF',job%extmat_suffix,dimen,extF_me,ierror) - endif + call divided_slice_read(imu,'extF',job%extmat_suffix,extF_me,desc_extF,extF_block_type,ierror) ! if (ierror==1) cycle ! @@ -2089,17 +1918,8 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! else ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - call MPI_File_read_all(fileh, imu_t, 1, mpi_integer, mpi_status_ignore, ierr) - ! - call MPI_File_read_all(fileh, extF_me, size(extF_me), mpi_double_precision, mpi_status_ignore, ierr) -#endif - else - read(iunit) imu_t - ! - read(iunit) extF_me - endif + call extFmatHandler%read(imu_t) + call extFmatHandler%read(extF_me, desc_extF, extF_block_type) ! endif ! @@ -2160,87 +1980,37 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! if (.not.job%IOextF_divide.or.job%IOextF_stitch) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - if(mpi_rank.eq.0) call MPI_File_write(fileh_w, imu, 1, mpi_integer, mpi_status_ignore, ierr) - call MPI_Barrier(mpi_comm_world, ierr) - call MPI_File_seek_shared(fileh_w, int(0,MPI_OFFSET_KIND), MPI_SEEK_END) - call MPI_File_write_all(fileh_w, mat_s, size(mat_s), mpi_double_precision, mpi_status_ignore, ierr) -#endif - else - write(chkptIO) imu - write(chkptIO) mat_s - endif + call exteigenHandler%write(imu) + call exteigenHandler%write(mat_s, desc_mat_s, mat_s_block_type) ! else ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call divided_slice_write_mpi(imu,'extF',job%j0extmat_suffix,Neigenroots,mat_s,mat_s_block_type) - else - if(mpi_rank.eq.0) call divided_slice_write(imu,'extF',job%j0extmat_suffix,Neigenroots,mat_s) - endif + call divided_slice_write(imu,'extF',job%j0extmat_suffix,mat_s,desc_mat_s,mat_s_block_type) ! endif ! enddo - ! Reset view to flat file - if ((.not.job%IOextF_divide) .and. blacs_size.gt.1) then -#ifdef TROVE_USE_MPI_ - read_offset = read_offset + (fitting%iparam(2)-fitting%iparam(1)+1)*int(ncontr_t,MPI_OFFSET_KIND)*ncontr_t*mpi_real_size & - + (fitting%iparam(2)-fitting%iparam(1)+1)*mpi_int_size - call MPI_File_set_view(fileh, int(0,MPI_OFFSET_KIND), mpi_byte, mpi_byte, "native", MPI_INFO_NULL, ierr) - call MPI_File_seek(fileh, read_offset, MPI_SEEK_SET) - - write_offset = write_offset + 3*int(Neigenroots,MPI_OFFSET_KIND)*Neigenroots*mpi_real_size - call MPI_File_set_view(fileh_w, int(0,MPI_OFFSET_KIND), mpi_byte, mpi_byte, "native", MPI_INFO_NULL, ierr) - call MPI_File_seek_shared(fileh_w, write_offset, MPI_SEEK_SET) - !write_offset = 0 - !call MPI_File_seek_shared(fileh_w, write_offset, MPI_SEEK_END) -#endif - endif ! if (allocated(extF_me)) deallocate(extF_me) call ArrayStop('extF_me') ! if (.not.job%IOextF_divide) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - call MPI_File_read_all(fileh, buf20, 18, mpi_character, mpi_status_ignore, ierr) -#endif - else - read(iunit) buf20(1:18) - endif + call extFmatHandler%read(buf20(1:18)) ! if (buf20(1:18)/='End external field') then write (out,"(' restore_Extvib_matrix_elements ',a,' has bogus footer: ',a)") filename,buf20(1:18) stop 'restore_Extvib_matrix_elements - bogus file format' end if ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - call MPI_File_close(fileh, ierr) -#endif - else - close(iunit,status='keep') - endif - ! - !job_is ='external field contracted matrix elements for J=0' - !call IOStart(trim(job_is),iunit) + deallocate(extFmatHandler) ! endif ! if (.not.job%IOextF_divide.or.job%IOextF_stitch) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - if(mpi_rank.eq.0) call MPI_File_write(fileh_w, 'End external field', 18, mpi_character, mpi_status_ignore, ierr) - call MPI_File_close(fileh_w, ierr) -#endif - else - write(chkptIO) 'End external field' - close(chkptIO,status='keep') - endif + call exteigenHandler%write('End external field') + deallocate(exteigenHandler) ! endif ! @@ -2270,97 +2040,52 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! contains ! - subroutine divided_slice_write(islice,name,suffix,N,field) + subroutine divided_slice_write(islice,name,suffix,field, field_desc, block_type) ! implicit none ! integer(ik),intent(in) :: islice character(len=*),intent(in) :: name,suffix - integer(ik),intent(in) :: N - real(rk),intent(in) :: field(N,N) + real(rk),intent(in) :: field(:,:) + integer, intent(in) :: field_desc(9) + type(MPI_Datatype),intent(in) :: block_type character(len=4) :: jchar - integer(ik) :: chkptIO + class(ioHandlerBase), allocatable :: ioHandler character(len=cl) :: buf,filename,job_is ! write(job_is,"('single swap_matrix')") ! - call IOStart(trim(job_is),chkptIO) - ! write(jchar, '(i4)') islice ! filename = trim(suffix)//trim(adjustl(jchar))//'.chk' ! - open(chkptIO,form='unformatted',action='write',position='rewind',status='replace',file=filename) - ! - write(chkptIO) trim(name) + call openFile(ioHandler, filename, err, action='write', & + position='rewind', status='replace', form='unformatted') + HANDLE_ERROR(err) + + call ioHandler%write(name) ! - write(chkptIO) field + call ioHandler%write(field, field_desc, block_type) ! - write(chkptIO) trim(name) + call ioHandler%write(name) ! - close(chkptIO) + deallocate(ioHandler) ! end subroutine divided_slice_write ! ! - subroutine divided_slice_write_mpi(islice,name,suffix,N,field,block_type) - ! - implicit none - ! - integer(ik),intent(in) :: islice - character(len=*),intent(in) :: name,suffix - integer(ik),intent(in) :: N - real(rk),dimension(:,:),intent(in) :: field - type(MPI_Datatype),intent(in) :: block_type - -#ifdef TROVE_USE_MPI_ - character(len=4) :: jchar - character(len=cl) :: filename - type(MPI_File) :: chkptMPIIO - integer(kind=MPI_OFFSET_KIND) :: offset - integer :: ierr - ! - write(job_is,"('single swap_matrix MPI')") - ! - ! - write(jchar, '(i4)') islice - ! - filename = trim(suffix)//trim(adjustl(jchar))//'.chk' - ! - offset = 0 - call MPI_File_open(mpi_comm_world, filename, mpi_mode_wronly+mpi_mode_create, mpi_info_null, chkptMPIIO, ierr) - call MPI_File_set_size(chkptMPIIO, offset, ierr) - ! - if(mpi_rank .eq. 0) call MPI_File_write(chkptMPIIO,name,len(trim(name)),mpi_character,mpi_status_ignore,ierr) - ! - !call MPI_File_get_position(chkptMPIIO, offset, ierr) - offset = int(len(trim(name)),MPI_OFFSET_KIND) - call MPI_File_set_view(chkptMPIIO, offset, mpi_byte, block_type, "native", MPI_INFO_NULL, ierr) - ! - call MPI_File_write_all(chkptMPIIO, field, size(field), mpi_double_precision, mpi_status_ignore, ierr) - ! - offset = offset + N*int(N,MPI_OFFSET_KIND)*mpi_real_size - call MPI_File_set_view(chkptMPIIO, offset, mpi_byte, mpi_byte, "native", MPI_INFO_NULL, ierr) - ! - if(mpi_rank .eq. 0) call MPI_File_write(chkptMPIIO,name,len(trim(name)),mpi_character,mpi_status_ignore,ierr) - ! - call MPI_File_close(chkptMPIIO, ierr) - ! -#endif - end subroutine divided_slice_write_mpi - ! - ! - subroutine divided_slice_read(islice,name,suffix,N,field,ierror) + subroutine divided_slice_read(islice,name,suffix,field,field_desc,block_type,ierror) ! implicit none ! integer(ik),intent(in) :: islice character(len=*),intent(in) :: name,suffix - integer(ik),intent(in) :: N - real(rk),intent(out) :: field(N,N) + real(rk),intent(out) :: field(:,:) + integer, intent(in) :: field_desc(9) + type(MPI_Datatype),intent(in) :: block_type integer(ik),intent(out) :: ierror character(len=4) :: jchar - integer(ik) :: chkptIO + class(ioHandlerBase), allocatable :: ioHandler character(len=cl) :: buf,filename,job_is integer(ik) :: ilen logical :: ifopened @@ -2369,124 +2094,35 @@ subroutine divided_slice_read(islice,name,suffix,N,field,ierror) ! write(job_is,"('single swap_matrix')") ! - call IOStart(trim(job_is),chkptIO) - ! write(jchar, '(i4)') islice ! filename = trim(suffix)//trim(adjustl(jchar))//'.chk' ! - open(chkptIO,form='unformatted',action='read',position='rewind',status='old',file=filename,err=15) + call openFile(ioHandler, filename, err, action='read', & + position='rewind', status='old', form='unformatted') + HANDLE_ERROR(err) ! ilen = LEN_TRIM(name) ! - read(chkptIO) buf(1:ilen) + call ioHandler%read(buf(1:ilen)) if ( trim(buf(1:ilen))/=trim(name) ) then write (out,"(' kinetic divided_slice_read in slice ',a20,': header is missing or wrong',a)") filename,buf(1:ilen) stop 'divided_slice_read - in slice - header missing or wrong' end if ! - read(chkptIO) field + call ioHandler%read(field, field_desc, block_type) ! - read(chkptIO) buf(1:ilen) + call ioHandler%read(buf(1:ilen)) if ( trim(buf(1:ilen))/=trim(name) ) then - write (out,"(' kinetic divided_slice_read in slice ',a20,': header is missing or wrong',a)") filename,buf(1:ilen) - stop 'divided_slice_read - in slice - header missing or wrong' + write (out,"(' kinetic divided_slice_read in slice ',a20,': footer is missing or wrong',a)") filename,buf(1:ilen) + stop 'divided_slice_read - in slice - footer missing or wrong' end if ! - close(chkptIO) - ! - return - ! - ! This error code will allow simply skipping the corresponding record/file without crashing the program - ! - 15 ierror = 1 - ! - ! we allow to skip opening the file only for the external matrix elements - ! - if (trim(name)/="extF") then - write (out,"(' kinetic divided_slice_read in slice ',a20,': file does not exist')") filename - stop 'divided_slice_read - in slice - file does not exist' - endif - ! - if (job%verbose>=4) write (out,"(' (skipped).')",advance='YES') - ! + deallocate(ioHandler) end subroutine divided_slice_read ! ! - subroutine divided_slice_read_mpi(islice,name,suffix,N,field,block_type,ierr) - ! - implicit none - ! - integer(ik),intent(in) :: islice - character(len=*),intent(in) :: name,suffix - integer(ik),intent(in) :: N - real(rk),dimension(:,:),intent(out) :: field - type(MPI_Datatype),intent(in) :: block_type - integer(ik),intent(out) :: ierr -#ifdef TROVE_USE_MPI_ - character(len=4) :: jchar - type(MPI_File) :: chkptMPIIO - character(len=cl) :: buf,filename,job_is - integer(ik) :: ilen - integer(MPI_OFFSET_KIND) :: offset - ! - ierror = 0 - ! - write(job_is,"('single swap_matrix MPI')") - ! - write(jchar, '(i4)') islice - ! - filename = trim(suffix)//trim(adjustl(jchar))//'.chk' - ilen = LEN_TRIM(name) - ! - call MPI_File_open(mpi_comm_world, filename, mpi_mode_rdonly, mpi_info_null, chkptMPIIO, ierr) - ! - if (ierr .ne. 0) then - ! - ! This error code will allow simply skipping the corresponding record/file without crashing the program - ! - ierr = 1 - ! - ! we allow to skip opening the file only for the external matrix elements - ! - if (trim(name)/="extF") then - write (out,"(' kinetic divided_slice_read_mpi in slice ',a20,': file does not exist')") filename - stop 'divided_slice_read_mpi - in slice - file does not exist' - endif - ! - if (job%verbose>=4) write (out,"(' (skipped).')",advance='YES') - ! - return - endif - ! - call MPI_File_read_all(chkptMPIIO,buf,ilen,mpi_character,mpi_status_ignore,ierr) - ! - if ( trim(buf(1:ilen))/=trim(name) ) then - write (out,"(' kinetic divided_slice_read_mpi in slice ',a20,': header is missing or wrong: ',a)") filename,buf(1:ilen) - stop 'divided_slice_read_mpi - in slice - header missing or wrong' - end if - ! - offset = ilen - call MPI_File_set_view(chkptMPIIO, int(ilen,MPI_OFFSET_KIND), mpi_byte, block_type, "native", MPI_INFO_NULL, ierr) - call MPI_File_read_all(chkptMPIIO, field, size(field), mpi_double_precision, mpi_status_ignore, ierr) - ! - offset = offset + N*int(N,MPI_OFFSET_KIND)*mpi_real_size - call MPI_File_set_view(chkptMPIIO, offset, mpi_byte, mpi_byte, "native", MPI_INFO_NULL, ierr) - ! - call MPI_File_read_all(chkptMPIIO,buf,ilen,mpi_character,mpi_status_ignore,ierr) - ! - if ( trim(buf(1:ilen))/=trim(name) ) then - write (out,"(' kinetic divided_slice_read_mpi in slice ',a20,': header is missing or wrong',a)") filename,buf(1:ilen) - stop 'divided_slice_read_mpi - in slice - header missing or wrong' - end if - ! - call MPI_File_close(chkptMPIIO, ierr) - ! -#endif - end subroutine divided_slice_read_mpi - ! - ! - subroutine divided_slice_read_vibrot(islice,suffix,N,field) + subroutine divided_slice_read_vibrot(islice,suffix,N,field,field_desc,block_type) ! implicit none ! @@ -2494,12 +2130,20 @@ subroutine divided_slice_read_vibrot(islice,suffix,N,field) character(len=*),intent(in) :: suffix integer(ik),intent(in) :: N real(rk),intent(out) :: field(N,N) + integer, intent(in) :: field_desc(9) + type(MPI_Datatype),intent(in) :: block_type character(len=4) :: jchar integer(ik) :: chkptIO character(len=cl) :: buf,filename,job_is integer(ik) :: rec_len,icontr,icontr1,icontr2 logical :: ifopened real(rk) :: f_t + +#ifdef TROVE_USE_MPI_ + if(blacs_size > 1) then + stop "MPI output not implemented in tran.f90:divided_slice_read_vibrot" + endif +#endif ! if (.not.job%IOmatelem_split) return ! @@ -2533,7 +2177,7 @@ subroutine divided_slice_read_vibrot(islice,suffix,N,field) end subroutine divided_slice_read_vibrot ! ! - subroutine divided_slice_write_vibrot(islice,suffix,N,field) + subroutine divided_slice_write_vibrot(islice,suffix,N,field,field_desc,block_type) ! implicit none ! @@ -2541,6 +2185,8 @@ subroutine divided_slice_write_vibrot(islice,suffix,N,field) character(len=*),intent(in) :: suffix integer(ik),intent(in) :: N real(rk),intent(in) :: field(N,N) + integer, intent(in) :: field_desc(9) + type(MPI_Datatype),intent(in) :: block_type character(len=4) :: jchar integer(ik) :: chkptIO character(len=cl) :: buf,filename,job_is @@ -2585,40 +2231,38 @@ end subroutine TRconvert_matel_j0_eigen ! Here we restore the vibrational (J=0) matrix elements of the rotational kinetic part G_rot and G_cor ! - subroutine restore_rot_kinetic_matrix_elements(jrot,treat_vibration,kinetic_part,chkptIO) - ! - implicit none - ! - integer(ik),intent(in) :: jrot - logical,intent(in) :: treat_vibration - character(len=3),intent(in) :: kinetic_part - integer(ik),intent(inout) :: chkptIO - ! - integer(ik) :: iclasses,alloc,ilevel,ib,max_deg_size,nclasses,islice - character(len=cl) :: job_is - character(len=18) :: buf18 - integer(ik) :: ncontr_t,rootsize - integer(ik),allocatable :: imat_t(:,:) - real(rk),allocatable :: mat_t(:) - ! - nclasses = bset_contr(1)%nclasses - ! - select case (kinetic_part) - ! - case('rot','top') + subroutine restore_rot_kinetic_matrix_elements(jrot,treat_vibration,kinetic_part,fileHandler) + ! + implicit none + ! + integer(ik),intent(in) :: jrot + logical,intent(in) :: treat_vibration + character(len=3),intent(in) :: kinetic_part + class(ioHandlerBase), intent(in) :: fileHandler + type(ErrorType) :: err + ! + integer(ik) :: iclasses,alloc,ilevel,ib,max_deg_size,nclasses,islice + character(len=cl) :: job_is + character(len=18) :: buf18 + integer(ik) :: ncontr_t,rootsize + integer(ik),allocatable :: imat_t(:,:) + real(rk),allocatable :: mat_t(:) + ! + nclasses = bset_contr(1)%nclasses + ! + select case (kinetic_part) + ! + case('rot','top') ! job_is ='Vib. matrix elements of the rot. kinetic' - call IOStart(trim(job_is),chkptIO) ! - open(chkptIO,form='unformatted',action='read',position='rewind',status='old',file=job%kinetmat_file) - ! - read(chkptIO) buf18 + call fileHandler%read(buf18) if (buf18/='Start Kinetic part') then write (out,"(' Vib. kinetic checkpoint file ',a,' has bogus header: ',a)") job%kinetmat_file,buf18 stop 'PTcontracted_matelem_class - bogus file format' end if ! - read(chkptIO) ncontr_t + call fileHandler%read(ncontr_t) ! if (bset_contr(1)%Maxcontracts/=ncontr_t) then write (out,"(' Vib. kinetic checkpoint file ',a)") job%kinetmat_file @@ -2635,13 +2279,13 @@ subroutine restore_rot_kinetic_matrix_elements(jrot,treat_vibration,kinetic_part allocate (imat_t(0:nclasses,ncontr_t),stat=alloc) call ArrayStart('contractive_space',alloc,size(imat_t),kind(imat_t)) ! - read(chkptIO) buf18(1:10) + call fileHandler%read(buf18(1:10)) if (buf18(1:10)/='icontr_cnu') then write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_cnu is missing ',a)") job%kinetmat_file,buf18(1:10) stop 'restore_rot_kinetic_matrix_elements - in file - icontr_cnu missing' end if ! - read(chkptIO) imat_t(0:nclasses,1:ncontr_t) + call fileHandler%read(imat_t(0:nclasses,1:ncontr_t)) ! if (job%vib_rot_contr) then ! @@ -2651,13 +2295,13 @@ subroutine restore_rot_kinetic_matrix_elements(jrot,treat_vibration,kinetic_part ! endif ! - read(chkptIO) buf18(1:11) + call fileHandler%read(buf18(1:11)) if (buf18(1:11)/='icontr_ideg') then write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_ideg is missing ',a)") job%kinetmat_file,buf18(1:11) stop 'restore_rot_kinetic_matrix_elements - in file - icontr_ideg missing' end if ! - read(chkptIO) imat_t(0:nclasses,1:ncontr_t) + call fileHandler%read(imat_t(0:nclasses,1:ncontr_t)) ! deallocate(imat_t) ! @@ -2667,8 +2311,8 @@ subroutine restore_rot_kinetic_matrix_elements(jrot,treat_vibration,kinetic_part ! if (trim(kinetic_part)=='rot') then ! - read(chkptIO) buf18(1:4) - if (buf18(1:4)/='g_ro') then + call fileHandler%read(buf18(1:5)) + if (buf18(1:5)/='g_rot') then write (out,"(' Vib. kinetic checkpoint file ',a,': g_rot is missing ',a)") trim(job%kinetmat_file),buf18(1:5) ! if (buf18(1:4)=='hvib'.or.buf18(1:3)=='End') & @@ -2681,13 +2325,6 @@ subroutine restore_rot_kinetic_matrix_elements(jrot,treat_vibration,kinetic_part endif ! if (job%vib_rot_contr) then - ! - !inquire(iolength=rec_len) f_t - !rec_len = rec_len*ncontr*PT%max_deg_size - !! - !do islice = 1,9+3*PT%Nmodes - ! call divided_slice_open_vib_rot(islice,rec_len,job%matelem_suffix) - !enddo ! call find_groundstate_icontr(bset_contr(1)%nclasses) ! @@ -2695,7 +2332,7 @@ subroutine restore_rot_kinetic_matrix_elements(jrot,treat_vibration,kinetic_part ! case('cor') ! - read(chkptIO) buf18(1:5) + call fileHandler%read(buf18(1:5)) if (buf18(1:5)/='g_cor') then write (out,"(' Vib. kinetic checkpoint file ',a,': g_cor is missing ',a)") job%kinetmat_file,buf18(1:5) stop 'restore_rot_kinetic_matrix_elements_posix - in file - g_cor missing' @@ -2711,14 +2348,12 @@ subroutine restore_rot_kinetic_matrix_elements(jrot,treat_vibration,kinetic_part ! endif ! - read(chkptIO) buf18(1:4) + call fileHandler%read(buf18(1:4)) if (buf18(1:4)/='hvib'.and.buf18(1:3)/='End'.and.buf18(1:4)/='vib-') then write (out,"(' Vib. kinetic checkpoint file ',a,': hvib, vib-rot, End is missing ',a)") job%kinetmat_file,buf18(1:4) stop 'restore_rot_kinetic_matrix_elements - in file - hvib or End missing' end if ! - close(chkptIO,status='keep') - ! end select ! contains @@ -2831,252 +2466,6 @@ end subroutine find_groundstate_icontr end subroutine restore_rot_kinetic_matrix_elements ! ! - ! [MPIIO] Here we restore the vibrational (J=0) matrix elements of the rotational kinetic part G_rot and G_cor - ! - subroutine restore_rot_kinetic_matrix_elements_mpi(jrot,treat_vibration,kinetic_part,fileh) - use mpi_aux - ! - implicit none - ! - integer(ik),intent(in) :: jrot - logical,intent(in) :: treat_vibration - character(len=3),intent(in) :: kinetic_part - type(MPI_File),intent(in) :: fileh -#ifdef TROVE_USE_MPI_ - ! - integer(ik) :: iclasses,alloc,ilevel,ib,max_deg_size,nclasses,islice - character(len=cl) :: job_is - character(len=18) :: buf18 - integer(ik) :: ncontr_t,rootsize - integer(ik),allocatable :: imat_t(:,:) - real(rk),allocatable :: mat_t(:) - integer :: ierr - ! - nclasses = bset_contr(1)%nclasses - ! - select case (kinetic_part) - ! - case('rot','top') - ! - job_is ='Vib. matrix elements of the rot. kinetic [mpi]' - ! - call MPI_File_read_all(fileh, buf18, 7, mpi_character, mpi_status_ignore, ierr) - if (buf18(1:7)/='[MPIIO]') then - write (out,"(' Vib. kinetic checkpoint file ',a,' is not an MPIIO file: ',a)") job%kinetmat_file,buf18 - stop 'restore_rot_kinetic_matrix_elements_mpi - Not an MPIIO file' - end if - ! - call MPI_File_read_all(fileh, buf18, 18, mpi_character, mpi_status_ignore, ierr) - if (buf18/='Start Kinetic part') then - write (out,"(' Vib. kinetic checkpoint file ',a,' has bogus header: ',a)") job%kinetmat_file,buf18 - stop 'restore_rot_kinetic_matrix_elements_mpi - bogus file format' - end if - ! - call MPI_File_read_all(fileh, ncontr_t, 1, mpi_integer, mpi_status_ignore, ierr) - ! - if (bset_contr(1)%Maxcontracts/=ncontr_t) then - write (out,"(' Vib. kinetic checkpoint file ',a)") job%kinetmat_file - write (out,"(' Actual and stored basis sizes at J=0 do not agree ',2i0)") bset_contr(1)%Maxcontracts,ncontr_t - stop 'restore_rot_kinetic_matrix_elements_mpi - in file - illegal nroots ' - end if - ! - rootsize = bset_contr(1)%Maxcontracts*(bset_contr(1)%Maxcontracts+1)/2 - ! - if (job%verbose>=6) write(out,"(/'Restore_rot_kin...: Number of elements: ',i8)") bset_contr(1)%Maxcontracts - ! - ! Read the indexes of the J=0 contracted basis set. - ! - allocate (imat_t(0:nclasses,ncontr_t),stat=alloc) - call ArrayStart('contractive_space',alloc,size(imat_t),kind(imat_t)) - ! - call MPI_File_read_all(fileh, buf18, 10, mpi_character, mpi_status_ignore, ierr) - if (buf18(1:10)/='icontr_cnu') then - write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_cnu is missing ',a)") job%kinetmat_file,buf18(1:10) - stop 'restore_rot_kinetic_matrix_elements_mpi - in file - icontr_cnu missing' - end if - ! - call MPI_File_read_all(fileh, imat_t, (nclasses+1)*ncontr_t, mpi_integer, mpi_status_ignore, ierr) - ! - if (job%vib_rot_contr) then - ! - allocate (TRicontr_cnu(0:nclasses,ncontr_t),stat=alloc) - call ArrayStart('TRicontr_cnu',alloc,size(TRicontr_cnu),kind(TRicontr_cnu)) - TRicontr_cnu = imat_t - ! - endif - ! - call MPI_File_read_all(fileh, buf18, 11, mpi_character, mpi_status_ignore, ierr) - if (buf18(1:11)/='icontr_ideg') then - write (out,"(' Vib. kinetic checkpoint file ',a,': icontr_ideg is missing ',a)") trim(job%kinetmat_file),buf18(1:11) - stop 'restore_rot_kinetic_matrix_elements_mpi - in file - icontr_ideg missing' - end if - ! - call MPI_File_read_all(fileh, imat_t, (nclasses+1)*ncontr_t, mpi_integer, mpi_status_ignore, ierr) - ! - deallocate(imat_t) - ! - call arraystop('contractive_space') - ! - ! Read the rotational part only if its really needed. - ! - if (trim(kinetic_part)=='rot') then - ! - call MPI_File_read_all(fileh, buf18, 5, mpi_character, mpi_status_ignore, ierr) - if (buf18(1:5)/='g_rot') then - write (out,"(' Vib. kinetic checkpoint file ',a,': g_rot is missing ',a)") trim(job%kinetmat_file),buf18(1:5) - ! - if (buf18(1:4)=='hvib'.or.buf18(1:3)=='End') & - write (out,"(a,a)") & - ' Most likely the split chk-points are supposed to be used.',& - 'Re-do MATELEM SAVE or use SPLIT in MATELEM !' - stop 'restore_rot_kinetic_matrix_elements_mpi - in file - g_rot missing' - end if - ! - endif - ! - if (job%vib_rot_contr) then - ! - call find_groundstate_icontr(bset_contr(1)%nclasses) - ! - endif - ! - case('cor') - ! - call MPI_File_read_all(fileh, buf18, 5, mpi_character, mpi_status_ignore, ierr) - if (buf18(1:5)/='g_cor') then - write (out,"(' Vib. kinetic checkpoint file ',a,': g_cor is missing ',a)") job%kinetmat_file,buf18(1:5) - stop 'restore_rot_kinetic_matrix_elements_mpi - in file - g_cor missing' - end if - ! - case('end') - ! - if (job%vib_rot_contr) then - ! - do islice = 1,9+3*FLNmodes - call divided_slice_close_vib_rot(islice) - enddo - ! - endif - ! - call MPI_File_read_all(fileh, buf18, 4, mpi_character, mpi_status_ignore, ierr) - if (buf18(1:4)/='hvib'.and.buf18(1:3)/='End'.and.buf18(1:4)/='vib-') then - write (out,"(' Vib. kinetic checkpoint file ',a,': hvib, vib-rot, End is missing ',a)") job%kinetmat_file,buf18(1:4) - stop 'restore_rot_kinetic_matrix_elements_mpi - in file - hvib or End missing' - end if - ! - end select - ! - contains - ! - subroutine divided_slice_open_vib_rot(islice,rec_len,suffix) - ! - implicit none - integer(ik),intent(in) :: islice,rec_len - character(len=*),intent(in) :: suffix - integer(ik) :: chkptIO - character(len=4) :: jchar - character(len=cl) :: buf,filename,job_is - integer(ik) :: ilen - logical :: ifopened - ! - if (.not.job%IOmatelem_split) return - ! - write(jchar, '(i4)') islice - ! - filename = trim(suffix)//trim(adjustl(jchar))//'.chk' - ! - write(job_is,"('single swap_matrix #',i8)") islice - ! - call IOStart(trim(job_is),chkptIO) - ! - ! use direct access format for vib-rot - ! - open(chkptIO,access='direct',action='read',status='old',recl=rec_len,file=filename) - ! - end subroutine divided_slice_open_vib_rot - ! - subroutine divided_slice_close_vib_rot(islice) - ! - implicit none - integer(ik),intent(in) :: islice - integer(ik) :: chkptIO - character(len=4) :: jchar - character(len=cl) :: filename,job_is - integer(ik) :: ilen - ! - if (.not.job%IOmatelem_split) return - ! - write(jchar, '(i4)') islice - ! - write(job_is,"('single swap_matrix #',i8)") islice - ! - call IOStart(trim(job_is),chkptIO) - ! - ! use direct access format for vib-rot - ! - close(chkptIO,status='keep') - ! - end subroutine divided_slice_close_vib_rot - ! - ! find correspondence between contracted quantum numbers: current and for J=0 - ! - subroutine find_groundstate_icontr(Nclasses) - ! - integer(ik), intent(in) :: Nclasses - integer(ik) :: maxcontr - ! - integer(ik) :: icontr,iterm,cnu_(1:bset_contr(1)%nclasses),icontr0,alloc,maxcontr0 - ! - if (job%verbose>=4) write(out,"(' Find correlation between J=0 and J/=0 contr. basis functions...')") - ! - maxcontr = size(TRicontr_cnu,dim=2) - ! - ! count icontr0-s - icontr0 = 0 - cnu_ = 0 - ! - do icontr = 1,maxcontr - ! - if (any(TRicontr_cnu(1:Nclasses,icontr)/=cnu_(1:Nclasses))) then - ! - icontr0 = icontr0 + 1 - cnu_ = TRicontr_cnu(1:Nclasses,icontr) - ! - endif - ! - enddo - ! - maxcontr0 = icontr0 - ! - allocate(Ncontr02icase0(maxcontr0,2),stat=alloc) - call ArrayStart('Ncontr02icase0',alloc,size(Ncontr02icase0),kind(Ncontr02icase0)) - ! - ! count ideg0 - icontr0 = 0 - cnu_ = 0 - ! - Ncontr02icase0(1,1) = 1 - Ncontr02icase0(maxcontr0,2) = maxcontr - ! - do icontr = 1,maxcontr - ! - if (any(TRicontr_cnu(1:Nclasses,icontr)/=cnu_(1:Nclasses))) then - ! - icontr0 = icontr0 + 1 - cnu_ = TRicontr_cnu(1:Nclasses,icontr) - Ncontr02icase0(icontr0,1) = icontr - if (icontr0>1) Ncontr02icase0(icontr0-1,2) = icontr-1 - ! - endif - ! - enddo - ! - end subroutine find_groundstate_icontr - ! -#endif - end subroutine restore_rot_kinetic_matrix_elements_mpi - ! - ! subroutine eigen_vib_matelem_vector(iparity,ilevelI,irootI,nlevels,nroots,cdimenmax,icoeff,fcoeff,cdimen,field,mat) ! implicit none