diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 4c6f29e..f132e5c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -11,6 +11,7 @@ jobs: env: OMP_NUM_THREADS: 1 USE_MPI: ${{ matrix.mpi && 1 || 0 }} + COMPILER: ${{ matrix.compiler }} runs-on: ${{ matrix.os }} strategy: fail-fast: false @@ -67,6 +68,17 @@ jobs: if: ${{ !matrix.mpi }} run: sudo apt-get install -y intel-mkl + - name: Build pFUnit + run: | + if [[ -d "/opt/intel/oneapi" ]]; then + source /opt/intel/oneapi/compiler/latest/env/vars.sh + source /opt/intel/oneapi/setvars.sh + if [ ${USE_MPI} = 1 ]; then + source /opt/intel/oneapi/mkl/latest/env/vars.sh + fi + fi + make COMPILER=${COMPILER} install-pfunit + - name: Build (gfortran) if: ${{ matrix.compiler == 'gfortran' }} run: | @@ -98,4 +110,4 @@ jobs: source /opt/intel/oneapi/mkl/latest/env/vars.sh fi fi - make test + make COMPILER=${COMPILER} USE_MPI=${USE_MPI} test diff --git a/.gitignore b/.gitignore index 2bade0a..5ec92f5 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,15 @@ wigxjpf-1.5/lib/ wigxjpf-1.5/bin/ wigxjpf-1.5/mod/ +# test executables +test/unit/test_io +test/unit/test_mpi_io + +# regression test files +test/regression/benchmarks +test/regression/benchmarks.tar.gz +test/regression/outputs + # Windows image file caches Thumbs.db ehthumbs.db diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..43ec5aa --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "lib/pFUnit"] + path = lib/pFUnit + url = https://github.com/Goddard-Fortran-Ecosystem/pFUnit.git diff --git a/errors.f90 b/errors.f90 new file mode 100644 index 0000000..27389e0 --- /dev/null +++ b/errors.f90 @@ -0,0 +1,12 @@ +module Errors + implicit none + + integer, parameter :: ERR_None = 0, & + ERR_Default = 1, & + ERR_FileNotFound = 2 + + type:: ErrorType + integer :: code = ERR_None + character(len=256) :: message = "" + end type +end module Errors diff --git a/errors.fpp b/errors.fpp new file mode 100644 index 0000000..0eafedb --- /dev/null +++ b/errors.fpp @@ -0,0 +1,16 @@ +! Error.fpp + +! From http://www.luckingtechnotes.com/fortran-error-handling-techniques on 23/06/2021 + +! Macros for error handling. +! Enables user to store errors and exit the subroutine in single statement. +! Fortran preprocessor must be enabled: -fpp. + +! Raise Error +! Store the error code and info (only if the current code is zero). +! Return from the subroutine. +#define RAISE_ERROR(msg, err) if (err%Code == ERR_None) then; err = ErrorType(Code=ERR_Default, Message=msg); end if; return; + +! Pass Error +! Returns if there's an error. +#define HANDLE_ERROR(err) if (err%Code /= ERR_None) then; print *, err%message; return; end if; diff --git a/io_handler_base.f90 b/io_handler_base.f90 new file mode 100644 index 0000000..eb4e9ff --- /dev/null +++ b/io_handler_base.f90 @@ -0,0 +1,83 @@ +module io_handler_base + use mpi_aux + use errors + + implicit none + + type, abstract :: ioHandlerBase + contains + procedure(open), deferred :: open + 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 + procedure(readScalar), deferred :: readScalar + procedure(read1DArray), deferred :: read1DArray + procedure(read2DArray), deferred :: read2DArray + end type ioHandlerBase + + abstract interface + subroutine open(this, fname, err, action, position, status, form, access) + import ioHandlerBase + import ErrorType + class(ioHandlerBase) :: this + type(ErrorType), intent(inout) :: err + character (len = *), intent(in) :: fname + character (len = *), intent(in) :: action + character (len = *), intent(in), optional :: position, status, form, access + end subroutine + subroutine writeScalar(this, object) + import ioHandlerBase + class(ioHandlerBase) :: this + class(*), intent(in) :: object + end subroutine + subroutine write1DArray(this, object) + import ioHandlerBase + class(ioHandlerBase) :: this + class(*), dimension(:), intent(in) :: object + end subroutine + subroutine write2DArray(this, object) + import ioHandlerBase + class(ioHandlerBase) :: this + class(*), dimension(:,:), intent(in) :: object + end subroutine + subroutine write2DArrayDistBlacs(this, object, descr, block_type) + import ioHandlerBase + import MPI_Datatype + class(ioHandlerBase) :: this + class(*), dimension(:,:), intent(in) :: object + integer, intent(in) :: descr(9) + type(MPI_Datatype), intent(in) :: block_type + end subroutine + subroutine write2DArrayDistColumn(this, object, mdimen) + import ioHandlerBase + import MPI_Datatype + class(ioHandlerBase) :: this + class(*), dimension(:,:), intent(in) :: object + integer, intent(in) :: mdimen + end subroutine + subroutine readScalar(this, object) + import ioHandlerBase + class(ioHandlerBase) :: this + class(*), intent(out) :: object + end subroutine + subroutine read1DArray(this, object) + import ioHandlerBase + class(ioHandlerBase) :: this + class(*), dimension(:), intent(out) :: object + end subroutine + subroutine read2DArray(this, object) + import ioHandlerBase + class(ioHandlerBase) :: this + class(*), dimension(:,:), intent(out) :: object + end subroutine + end interface + + private + + public :: ioHandlerBase + +end module diff --git a/io_handler_ftn.f90 b/io_handler_ftn.f90 new file mode 100644 index 0000000..fb48e26 --- /dev/null +++ b/io_handler_ftn.f90 @@ -0,0 +1,239 @@ +#include "errors.fpp" + +module io_handler_ftn + use mpi_aux + use io_handler_base + use errors + use iso_fortran_env + + implicit none + + type, extends(ioHandlerBase) :: ioHandlerFTN + integer :: iounit = 0 + integer :: stat = 0 + logical :: isOpen = .false. + contains + procedure :: writeScalar => writeScalarFTN + procedure :: write1DArray => write1DArrayFTN + procedure :: write2DArray => write2DArrayFTN + procedure :: write2DArrayDistBlacs => write2DArrayDistBlacsFTN + procedure :: write2DArrayDistColumn => write2DArrayDistColumnFTN + procedure :: readScalar => readScalarFTN + procedure :: read1DArray => read1DArrayFTN + procedure :: read2DArray => read2DArrayFTN + procedure :: open + procedure :: close + 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() + end subroutine + + subroutine open(this, fname, err, action, position, status, form, access) + class(ioHandlerFTN) :: this + type(ErrorType), intent(inout) :: err + character (len = *), intent(in) :: fname + character (len = *), intent(in) :: action + character (len = *), intent(in), optional :: position, status, form, access + character (len = 20) :: positionVal, statusVal, formVal, accessVal + + if (this%isOpen) then + RAISE_ERROR("ERROR: Tried to open second file", err) + endif + + if (present(position)) then + positionVal = position + else + positionVal = 'asis' + end if + + if (present(status)) then + statusVal = status + else + statusVal = 'unknown' + end if + + if (present(form)) then + formVal = form + else + formVal = 'unformatted' + end if + + if (present(access)) then + accessVal = access + else + accessVal = 'sequential' + end if + + print *, "FTN: Opening ", trim(fname), " with ", \ + trim(positionVal), " ", trim(statusVal), " ", trim(formVal), " ", trim(accessVal) + + open(newunit=this%iounit, action=action,\ + form=formVal, position=positionVal, status=statusVal, file=fname,\ + iostat=this%stat) + + if (this%stat == 0) then + this%isOpen = .true. + else + RAISE_ERROR("ERROR: Could not open file", err) + endif + end subroutine + + subroutine close(this) + class(ioHandlerFTN) :: this + if (this%isOpen) then + close(this%iounit) + this%isOpen = .false. + endif + 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 + type is (real) + write(this%iounit) object + type is (complex) + write(this%iounit) object + type is (character(len=*)) + write(this%iounit) object + class default + print *, "ERROR: 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 + type is (real) + write(this%iounit) object + type is (complex) + write(this%iounit) object + class default + print *, "ERROR: 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 + type is (integer(int64)) + write(this%iounit) object + type is (real(real32)) + write(this%iounit) object + type is (real(real64)) + write(this%iounit) object + type is (complex(kind=4)) + write(this%iounit) object + type is (complex(kind=8)) + write(this%iounit) object + class default + print *, "ERROR: Tried to write unsupported type" + end select + end subroutine + + subroutine write2DArrayDistBlacsFTN(this, object, descr, block_type) + ! Write arrays distributed using blacs + + class(ioHandlerFTN) :: this + class(*), dimension(:,:), intent(in) :: object + integer, intent(in) :: descr(9) + type(MPI_Datatype), intent(in) :: block_type + + ! Using the fortran io_handler means array isn't distributed, just write normally + call this%write2DArray(object) + end subroutine + + subroutine write2DArrayDistColumnFTN(this, object, mdimen) + ! Write arrays distributed as columns using co_distr_data + + class(ioHandlerFTN) :: this + class(*), intent(in) :: object(:,:) + integer, intent(in) :: mdimen ! Dimension of entire distributed array + + ! Using the fortran io_handler means array isn't distributed, just write normally + call this%write2DArray(object) + end subroutine + + subroutine readScalarFTN(this, object) + class(ioHandlerFTN) :: this + class(*), intent(out) :: object + print *, "reading object with FTN IO" + select type(object) + type is (integer) + read(this%iounit) object + type is (real) + read(this%iounit) object + type is (complex) + read(this%iounit) object + class default + print *, "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) + read(this%iounit) object + type is (real) + read(this%iounit) object + type is (complex) + read(this%iounit) object + class default + print *, "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) + read(this%iounit) object + type is (real) + read(this%iounit) object + type is (complex) + read(this%iounit) object + class default + print *, "Unsupported type!" + end select + end subroutine +end module diff --git a/io_handler_mpi.f90 b/io_handler_mpi.f90 new file mode 100644 index 0000000..4d12362 --- /dev/null +++ b/io_handler_mpi.f90 @@ -0,0 +1,343 @@ +#include "errors.fpp" + +#define MPI_WRAPPER(function, handle, obj, size, mytype, status, err) \ +select type(obj); \ + type is (integer(ik)); \ + call function(handle, obj, size, mytype, status, err); \ + type is (integer(hik)); \ + call function(handle, obj, size, mytype, status, err); \ + type is (real); \ + call function(handle, obj, size, mytype, status, err); \ + type is (real(rk)); \ + call function(handle, obj, size, mytype, status, err); \ + type is (real(ark)); \ + call function(handle, obj, size, mytype, status, err); \ + type is (complex); \ + call function(handle, obj, size, mytype, status, err); \ + type is (character(len=*)); \ + call function(handle, obj, size, mytype, status, err); \ + class default; \ + write(*,*) 'Not covered'; \ +end select + + +module io_handler_mpi + use mpi_f08 + use mpi_aux + use io_handler_base + use errors + use accuracy, only : rk, ark, ik, hik + + implicit none + + type, extends(ioHandlerBase) :: ioHandlerMPI + integer (kind=MPI_Offset_kind) :: offset = 0 + integer :: rank = 0 + type(MPI_File) :: fileh + logical :: isOpen = .false. + contains + procedure :: writeScalar => writeScalarMPI + procedure :: write1DArray => write1DArrayMPI + procedure :: write2DArray => write2DArrayMPI + procedure :: write2DArrayDistBlacs => write2DArrayDistBlacsMPI + procedure :: write2DArrayDistColumn => write2DArrayDistColumnMPI + procedure :: readScalar => readScalarMPI + procedure :: read1DArray => read1DArrayMPI + procedure :: read2DArray => read2DArrayMPI + procedure :: open + procedure :: close + 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() + end subroutine + + subroutine open(this, fname, err, action, position, status, form, access) + ! writer MPI constructor + class(ioHandlerMPI) :: this + character (len = *), intent(in) :: fname + 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 + integer :: ierr + + if (this%isOpen) then + RAISE_ERROR("ERROR: Tried to open second file", err) + endif + + if (present(position)) then + positionVal = position + else + positionVal = 'append' + end if + + if (present(status)) then + statusVal = status + else + statusVal = 'unknown' + end if + + if (present(form)) then + formVal = form + else + formVal = 'formatted' + end if + + if (present(access)) then + accessVal = access + else + accessVal = 'sequential' + end if + + print *, "MPI: Opening ", trim(fname), " with ", \ + trim(positionVal), " ", trim(statusVal), " ", trim(formVal), " ", trim(accessVal) + + ! FIXME use above flags to change open behaviour + + call MPI_Comm_rank(MPI_COMM_WORLD, this%rank, ierr) + + ! 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 + call MPI_File_open(MPI_COMM_WORLD, fname, MPI_MODE_RDONLY, MPI_INFO_NULL, this%fileh, ierr) + endif + + call MPI_File_set_errhandler(this%fileh, MPI_ERRORS_ARE_FATAL) + ! FIXME handle error + end subroutine open + + subroutine close(this) + class(ioHandlerMPI) :: this + integer :: ierr + if (this%isOpen) then + call MPI_File_close(this%fileh, ierr) + this%isOpen = .false. + endif + ! FIXME handle error + end subroutine close + + subroutine getMPIVarInfo(object, byteSize, mpiType) + class(*), intent(in) :: object + integer, intent(out) :: byteSize + type(MPI_Datatype), intent(out) :: mpiType + + select type(object) + type is (integer(kind=4)) + byteSize = sizeof(object) + mpiType = MPI_INTEGER + type is (integer(kind=8)) + byteSize = sizeof(object) + mpiType = MPI_LONG + type is (real(kind=4)) + byteSize = sizeof(object) + mpiType = MPI_FLOAT + type is (real(kind=8)) + byteSize = sizeof(object) + mpiType = MPI_DOUBLE + type is (complex(kind=4)) + byteSize = sizeof(object) + mpiType = MPI_COMPLEX + type is (complex(kind=8)) + byteSize = sizeof(object) + mpiType = MPI_DOUBLE_COMPLEX + type is (character(len=*)) + byteSize = len(object) * sizeof('a') + mpiType = MPI_CHARACTER + class default + print *, "ERROR: Unknown type" + end select + end subroutine + + subroutine writeScalarMPI(this, object) + class(ioHandlerMPI) :: this + class(*), intent(in) :: object + + integer :: byteSize, ierr, length + type(MPI_Datatype) :: mpiType + + call getMPIVarInfo(object, byteSize, mpiType) + this%offset = this%offset + 4+byteSize+4 + + select type(object) + type is (character(len=*)) + length = len(object) + class default + length = 1 + end select + + if (this%rank /= 0) then + return + end if + + 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) + end subroutine + + subroutine write1DArrayMPI(this, object) + class(ioHandlerMPI) :: this + class(*), intent(in) :: object(:) + print *, "ERROR: 1D array saving not currently supported" + end subroutine + + subroutine write2DArrayMPI(this, object) + class(ioHandlerMPI) :: this + class(*), intent(in) :: object(:,:) + + type(MPI_Datatype) :: mpiType + integer :: byteSize, globalSize, ierr, length + + integer(kind = MPI_OFFSET_KIND) :: arrSizeBytes + + globalSize = size(object) + + call getMPIVarInfo(object(1,1), byteSize, mpiType) + arrSizeBytes = globalSize*byteSize + + this%offset = this%offset + 4 + arrSizeBytes + 4 + + if (this%rank /= 0) then + return + 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) + end subroutine + + subroutine write2DArrayDistBlacsMPI(this, object, descr, block_type) + class(ioHandlerMPI) :: this + class(*), intent(in) :: 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, globalSize, ierr + integer(kind = MPI_OFFSET_KIND) :: arrSizeBytes + + integer :: dims(2) + + dims(:) = descr(3:4) + globalSize = dims(1)*dims(2) + + 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, & + '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, & + 'native', MPI_INFO_NULL, ierr) + end subroutine + + subroutine write2DArrayDistColumnMPI(this, object, mdimen) + class(ioHandlerMPI) :: this + class(*), intent(in) :: object(:,:) + integer, intent(in) :: mdimen ! Dimension of entire distributed array + + type(MPI_Datatype) :: mpiType + integer :: byteSize, globalSize, ierr, writestat + integer(kind = MPI_OFFSET_KIND) :: arrSizeBytes + + 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) + ! 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) + end subroutine + + subroutine readScalarMPI(this, object) + class(ioHandlerMPI) :: this + class(*), intent(out) :: object + print *, "reading object to MPI IO" + ! Example object handling + select type(object) + type is (integer) + print *, object + type is (real) + print *, object + type is (complex) + print *, object + class default + print *, "Unsupported type!" + end select + end subroutine + + subroutine read1DArrayMPI(this, object) + class(ioHandlerMPI) :: this + class(*), dimension(:), intent(out) :: object + print *, "reading 1D array to MPI IO" + end subroutine + + subroutine read2DArrayMPI(this, object) + class(ioHandlerMPI) :: this + class(*), dimension(:,:), intent(out) :: object + print *, "reading 2D array to MPI IO" + end subroutine + +end module diff --git a/lib/pFUnit b/lib/pFUnit new file mode 160000 index 0000000..bc4de1f --- /dev/null +++ b/lib/pFUnit @@ -0,0 +1 @@ +Subproject commit bc4de1fa39f1eb8843cb60e2f418787301dde06e diff --git a/makefile b/makefile index a6d8e43..bff92c5 100644 --- a/makefile +++ b/makefile @@ -17,7 +17,7 @@ USE_MPI ?= 0 # Intel ####### ifeq ($(strip $(COMPILER)),intel) - FOR = ifort + FC = ifort FFLAGS = -cpp -ip -align -ansi-alias -mcmodel=medium -parallel -nostandard-realloc-lhs -qopenmp -module $(OBJDIR) ifeq ($(strip $(MODE)),debug) @@ -36,7 +36,7 @@ ifeq ($(strip $(COMPILER)),intel) # gfortran ########## else ifeq ($(strip $(COMPILER)),gfortran) - FOR = gfortran + FC = gfortran FFLAGS = -cpp -std=gnu -fopenmp -march=native -ffree-line-length-512 -fcray-pointer -I$(OBJDIR) -J$(OBJDIR) GCC_VERSION_GT_10 := $(shell expr `gcc -dumpversion | cut -f1 -d.` \>= 10) @@ -65,14 +65,18 @@ endif CPPFLAGS = -D_EXTFIELD_DEBUG_ ifneq ($(strip $(USE_MPI)),0) - FOR = mpif90 + FC = mpif90 FFLAGS += -DTROVE_USE_MPI_ endif +export FC +export USE_MPI + ################################################################################ ## LIBRARIES ################################################################################ +PFUNIT_DIR = lib/pFUnit WIGXJPF_DIR = wigxjpf-1.5 WIGXJPF_LIB = $(WIGXJPF_DIR)/lib/libwigxjpf.a LIB = $(LAPACK) $(LIBS) $(WIGXJPF_LIB) $(ARPACK) @@ -85,6 +89,12 @@ BINDIR=. SRCDIR=. OBJDIR=. user_pot_dir=. +TARGET=$(BINDIR)/$(EXE) + +MPI_SRCS = +ifneq ($(strip $(USE_MPI)),0) + MPI_SRCS += io_handler_mpi.f90 +endif SRCS := timer.f90 accuracy.f90 diag.f90 dipole.f90 extfield.f90 fields.f90 fwigxjpf.f90 input.f90 kin_xy2.f90 lapack.f90 \ me_bnd.f90 me_numer.f90 me_rot.f90 me_str.f90 \ @@ -94,8 +104,11 @@ 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 \ - refinement.f90 richmol_data.f90 rotme_cart_tens.f90 symmetry.f90 tran.f90 trove.f90 $(pot_user).f90 + io_handler_base.f90 io_handler_ftn.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} +MPI_OBJS := ${MPI_SRCS:.f90=.o} VPATH = $(SRCDIR):$(user_pot_dir):$(OBJDIR) @@ -103,15 +116,15 @@ VPATH = $(SRCDIR):$(user_pot_dir):$(OBJDIR) ## TARGETS ################################################################################ -.PHONY: all, clean, cleanall, tarball, checkin, test +.PHONY: all, clean, cleanall, tarball, checkin, test, install-pfunit -all: $(BINDIR) $(OBJDIR) $(BINDIR)/$(EXE) +all: $(OBJDIR) $(TARGET) %.o : %.f90 - $(FOR) -c $(FFLAGS) $(CPPFLAGS) -o $(OBJDIR)/$@ $< + $(FC) -c $(FFLAGS) $(CPPFLAGS) -o $(OBJDIR)/$@ $< -$(BINDIR)/$(EXE): $(OBJS) $(WIGXJPF_LIB) - $(FOR) $(FFLAGS) -o $@ $(addprefix $(OBJDIR)/,$(OBJS)) $(LIB) +$(BINDIR)/$(EXE): $(BINDIR) $(OBJS) $(WIGXJPF_LIB) + $(FC) $(FFLAGS) -o $@ $(addprefix $(OBJDIR)/,$(OBJS)) $(LIB) $(WIGXJPF_LIB): $(MAKE) -C $(WIGXJPF_DIR) @@ -126,8 +139,16 @@ ifneq ($(BINDIR),.) mkdir -p $(BINDIR) endif +install-pfunit: + git submodule update --init # Make sure we have pfunit + mkdir $(PFUNIT_DIR)/build + cd $(PFUNIT_DIR)/build; cmake .. + $(MAKE) -C $(PFUNIT_DIR)/build + $(MAKE) -C $(PFUNIT_DIR)/build install + clean: - rm -rf $(BINDIR)/$(EXE) $(OBJDIR)/*.mod $(OBJDIR)/*.o + rm -rf $(TARGET) $(OBJDIR)/*.mod $(OBJDIR)/*.o + $(MAKE) -C test/unit clean cleanall: clean $(MAKE) -C $(WIGXJPF_DIR) clean @@ -138,8 +159,26 @@ tarball: checkin: ci -l Makefile *.f90 -test: $(BINDIR)/$(EXE) - cd test; ./run_regression_tests.sh +test: regression-tests unit-tests-nompi unit-tests-mpi + +regression-tests: $(TARGET) + echo "Running regression tests" + cd test/regression; ./run_regression_tests.sh + +unit-tests-nompi: $(TARGET) + $(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) + $(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 +else +unit-tests-mpi: $(TARGET) + echo "Skipping unit tests with MPI (USE_MPI not set)" +endif ################################################################################ ## DEPENDENCIES @@ -179,7 +218,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 +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) 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 @@ -202,3 +241,6 @@ 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 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 diff --git a/perturbation.f90 b/perturbation.f90 index 3d98e8c..eb78c07 100644 --- a/perturbation.f90 +++ b/perturbation.f90 @@ -2,6 +2,8 @@ ! This unit is not the perturbation theory anymore, it solves the Schroedinger equation ! variationally. ! +#include "errors.fpp" + module perturbation use accuracy use molecules, only : MOrepres_arkT,MLsymmetry_transform_func,polintark,MLrotsymmetry_func,MLrotsymmetry_generate @@ -13,6 +15,12 @@ module perturbation use symmetry , only : SymmetryInitialize,sym use me_numer use diag + use io_handler_base + use io_handler_ftn +#ifdef TROVE_USE_MPI_ + use io_handler_mpi +#endif + use errors ! use omp_lib @@ -23,7 +31,7 @@ module perturbation public PTcontracted_matelem_class,PTeigenfunction_orthogonality public PThamiltonian_contract,PTget_primitive_matelements,PTDVR_initialize public PTcheck_point_contracted_space,PT_conctracted_rotational_bset - public PTTest_eigensolution,PTanalysis_density,PTstore_icontr_cnu,PTstorempi_icontr_cnu + public PTTest_eigensolution,PTanalysis_density,PTstore_icontr_cnu public PTintcoeffsT,PTrotquantaT,PTNclasses,PTdefine_contr_from_eigenvect,PTeigenT,PTrepresT public PTstore_contr_matelem,PTcontracted_matelem_class_fast,PTstore_contr_matelem_II,PTcontracted_matelem_class_fast_II @@ -16108,6 +16116,9 @@ 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 ! ! call TimerStart('Contracted matelements-class') @@ -16212,44 +16223,23 @@ subroutine PTcontracted_matelem_class(jrot) ! Prepare the checkpoint file ! job_is ='Vib. matrix elements of the rot. kinetic part' - if (trim(job%kinetmat_format).eq.'MPIIO') then + call IOStart(trim(job_is),chkptIO) + ! #ifdef TROVE_USE_MPI_ - call MPI_File_open(mpi_comm_world, job%kinetmat_file, mpi_mode_wronly+mpi_mode_create, mpi_info_null, chkptMPIIO, ierr) - - call MPI_File_set_errhandler(chkptMPIIO, MPI_ERRORS_ARE_FATAL) - - mpioffset=0 - - call MPI_File_set_size(chkptMPIIO, mpioffset, ierr) - - call MPI_File_seek(chkptMPIIO, mpioffset, MPI_SEEK_CUR, ierr) - - - if (mpi_rank.eq.0) then !AT - call MPI_File_write(chkptMPIIO,'[MPIIO]',7,mpi_character,mpi_status_ignore,ierr) - call MPI_File_write(chkptMPIIO,'Start Kinetic part',18,mpi_character,mpi_status_ignore,ierr) - call TimerStart('mpiiosingle') !AT - ! - ! store the bookkeeping information about the contr. basis set - ! - call PTstoreMPI_icontr_cnu(PT%Maxcontracts,chkptMPIIO,job%IOkinet_action) - - call TimerStop('mpiiosingle') !AT - endif - call MPI_Barrier(mpi_comm_world, ierr) - call MPI_File_seek_shared(chkptMPIIO, mpioffset, MPI_SEEK_END, ierr) + allocate(ioHandlerMPI::ioHandler) +#else + allocate(ioHandlerFTN::ioHandler) #endif - else - call IOStart(trim(job_is),chkptIO) - ! - open(chkptIO,form='unformatted',action='write',position='rewind',status='replace',file=job%kinetmat_file) - write(chkptIO) 'Start Kinetic part' - ! - ! store the bookkeeping information about the contr. basis set - ! - call PTstore_icontr_cnu(PT%Maxcontracts,chkptIO,job%IOkinet_action) - ! - endif + call ioHandler%open(& + 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 + ! + call PTstore_icontr_cnu(PT%Maxcontracts,ioHandler,job%IOkinet_action) endif ! ! maximal size of the primitive matrix in all classes @@ -16426,18 +16416,7 @@ subroutine PTcontracted_matelem_class(jrot) ! if (trim(job%IOkinet_action)=='SAVE'.and..not.job%IOmatelem_split) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - if(mpi_rank.eq.0) then - !call MPI_File_seek(chkptMPIIO, mpioffset, MPI_SEEK_END) - call MPI_File_write_shared(chkptMPIIO,'g_rot',5,mpi_character,mpi_status_ignore,ierr) - !else - ! call MPI_File_write_shared(chkptMPIIO,'',0,mpi_character,mpi_status_ignore,ierr) - endif -#endif - else - write(chkptIO) 'g_rot' - endif + call ioHandler%write('g_rot') ! endif ! @@ -16500,11 +16479,7 @@ subroutine PTcontracted_matelem_class(jrot) ! ! store the matrix elements ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call co_write_matrix_distr(grot_t,mdimen, startdim, enddim,chkptMPIIO) - else - write(chkptIO) grot_t - endif + call ioHandler%write(grot_t, mdimen) ! endif endif @@ -16515,20 +16490,7 @@ subroutine PTcontracted_matelem_class(jrot) if (job%verbose>=2) write(out,"(/'Coriolis part of the Kinetic energy operator...')") ! if (trim(job%IOkinet_action)=='SAVE'.and..not.job%IOmatelem_split) then - ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - ! call MPI_File_seek(chkptMPIIO, mpioffset, MPI_SEEK_END) - if(mpi_rank.eq.0) then - call MPI_File_write_shared(chkptMPIIO,'g_cor',5,mpi_character,mpi_status_ignore,ierr) - !else - ! call MPI_File_write_shared(chkptMPIIO,'',0,mpi_character,mpi_status_ignore,ierr) - endif -#endif - else - write(chkptIO) 'g_cor' - endif - ! + call ioHandler%write('g_cor') endif ! ! Run the loop over all term of the expansion of the Hamiltonian @@ -16624,11 +16586,7 @@ subroutine PTcontracted_matelem_class(jrot) ! ! store the matrix elements ! - if (trim(job%kinetmat_format).eq.'MPIIO') then - call co_write_matrix_distr(gcor_t,mdimen, startdim, enddim,chkptMPIIO) - else - write(chkptIO) gcor_t - endif + call ioHandler%write(gcor_t, mdimen) ! endif endif @@ -16701,24 +16659,24 @@ subroutine PTcontracted_matelem_class(jrot) ! ! store the rotational matrix elements ! - write(chkptIO) 'g_rot' + call ioHandler%write('g_rot') ! do k1 = 1,3 do k2 = 1,3 ! - write(chkptIO) grot_(k1,k2,:,:) + call ioHandler%write(grot_(k1,k2,:,:), mdimen_) ! enddo enddo ! - write(chkptIO) 'g_cor' + call ioHandler%write('g_cor') ! ! store the Coriolis matrix elements ! do k1 = 1,PT%Nmodes do k2 = 1,3 ! - write(chkptIO) gcor_(k1,k2,:,:) + call ioHandler%write(gcor_(k1,k2,:,:), mdimen_) ! enddo enddo @@ -17030,20 +16988,8 @@ subroutine PTcontracted_matelem_class(jrot) (.not.job%IOmatelem_divide.or.job%iswap(1)==0) .and. & (.not.job%IOmatelem_split.or.job%iswap(1)==0)) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - if(mpi_rank.eq.0) then - !call MPI_File_seek(chkptMPIIO, mpioffset, MPI_SEEK_END) - call MPI_File_write_shared(chkptMPIIO,'hvib',4,mpi_character,mpi_status_ignore,ierr) - !else - ! call MPI_File_write_shared(chkptMPIIO,'',0,mpi_character,mpi_status_ignore,ierr) - endif - call co_write_matrix_distr(hvib%me,mdimen, startdim, enddim,chkptMPIIO) -#endif - else - write(chkptIO) 'hvib' - write(chkptIO) hvib%me - endif + call ioHandler%write('hvib') + call ioHandler%write(hvib%me, mdimen) ! endif ! @@ -17055,19 +17001,8 @@ subroutine PTcontracted_matelem_class(jrot) (.not.job%IOmatelem_divide.or.job%iswap(1)==0 ).and. & (.not.job%IOmatelem_split.or.job%iswap(1)==0 ) ) then ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - call mpi_barrier(mpi_comm_world, ierr) - call MPI_File_seek_shared(chkptMPIIO, mpioffset, MPI_SEEK_END) - if(mpi_rank.eq.0) then - call MPI_File_write_shared(chkptMPIIO,'End Kinetic part',16,mpi_character,mpi_status_ignore,ierr) - endif - call MPI_File_close(chkptMPIIO, ierr) -#endif - else - write(chkptIO) 'End Kinetic part' - close(chkptIO,status='keep') - endif + call ioHandler%write('End Kinetic part') + deallocate(ioHandler) ! endif ! @@ -17814,6 +17749,9 @@ subroutine PTcontracted_matelem_class_fast_II(jrot) integer(ik) :: Nsym,isym,Nsymi,Nsymj,jsymcoeff,Ntot integer(ik), allocatable :: isymcoeff_vs_isym(:,:) double precision,parameter :: a0 = -0.5d0 + + class(ioHandlerBase), allocatable :: ioHandler + type(ErrorType) :: err ! call TimerStart('Contracted matelements-class-fast') ! @@ -17922,15 +17860,19 @@ subroutine PTcontracted_matelem_class_fast_II(jrot) job_is ='Vib. matrix elements of the rot. kinetic part' call IOStart(trim(job_is),chkptIO) - open(chkptIO,form='unformatted',action='write',position='rewind',status='replace',file=job%kinetmat_file) - write(chkptIO) 'Start Kinetic part' + allocate(ioHandlerFTN::ioHandler) + call ioHandler%open( & + 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 ! - call PTstore_icontr_cnu(PT%Maxcontracts,chkptIO,job%IOkinet_action) + call PTstore_icontr_cnu(PT%Maxcontracts,ioHandler,job%IOkinet_action) ! if (job%vib_rot_contr) then - write(chkptIO) 'vib-rot' + call ioHandler%write('vib-rot') endif ! endif @@ -18042,7 +17984,7 @@ subroutine PTcontracted_matelem_class_fast_II(jrot) ! if (trim(job%IOkinet_action)=='SAVE'.and..not.job%IOmatelem_split) then ! - write(chkptIO) 'g_rot' + call ioHandler%write('g_rot') ! endif ! @@ -18122,7 +18064,7 @@ subroutine PTcontracted_matelem_class_fast_II(jrot) if (job%IOmatelem_split) then write(chkptIO_) grot_t else - write(chkptIO) grot_t + call ioHandler%write(grot_t, mdimen) endif ! endif @@ -18146,7 +18088,7 @@ subroutine PTcontracted_matelem_class_fast_II(jrot) ! if (trim(job%IOkinet_action)=='SAVE'.and..not.job%IOmatelem_split) then ! - write(chkptIO) 'g_cor' + call ioHandler%write('g_cor') ! endif ! @@ -18297,7 +18239,7 @@ subroutine PTcontracted_matelem_class_fast_II(jrot) if (job%IOmatelem_split) then write(chkptIO_) grot_t else - write(chkptIO) grot_t + call ioHandler%write(grot_t, mdimen) endif ! endif @@ -18342,7 +18284,7 @@ subroutine PTcontracted_matelem_class_fast_II(jrot) if ((trim(job%IOkinet_action)=='SAVE'.or.trim(job%IOkinet_action)=='VIB_SAVE').and. & (.not.job%IOmatelem_split.or.job%iswap(1)==0.or.job%iswap(1)==(PT%Nmodes+3)*3 ) ) then ! - write(chkptIO) 'hvib' + call ioHandler%write('hvib') ! endif ! @@ -18552,10 +18494,8 @@ subroutine PTcontracted_matelem_class_fast_II(jrot) enddo !$omp end parallel do ! - !write(chkptIO) hvib_t - ! if (trim(job%IOkinet_action)=='SAVE') then - write(chkptIO) hvib_t + call ioHandler%write(hvib_t) endif ! endif @@ -18592,8 +18532,8 @@ subroutine PTcontracted_matelem_class_fast_II(jrot) if ((trim(job%IOkinet_action)=='SAVE'.or.trim(job%IOkinet_action)=='VIB_SAVE').and.& (.not.job%IOmatelem_split.or.job%iswap(1)==0 )) then ! - write(chkptIO) 'End Kinetic part' - close(chkptIO,status='keep') + call ioHandler%write('End Kinetic part') + deallocate(ioHandler) ! endif ! @@ -34889,10 +34829,10 @@ subroutine matvec_sym(n,bterm,h,z,w) return end subroutine matvec_sym + subroutine PTstore_icontr_cnu(Maxcontracts,ioHandler,dir) - subroutine PTstore_icontr_cnu(Maxcontracts,iunit,dir) - - integer(ik),intent(in) :: Maxcontracts,iunit + integer(ik),intent(in) :: Maxcontracts + class(ioHandlerBase), intent(in) :: ioHandler character(len=18),intent(in) :: dir integer(ik) :: alloc character(len=18) :: buf18 @@ -34902,20 +34842,19 @@ subroutine PTstore_icontr_cnu(Maxcontracts,iunit,dir) selectcase (trim(dir)) ! case ('SAVE') - ! - write(iunit) Maxcontracts - ! - write(iunit) 'icontr_cnu' - ! - write(iunit) PT%icontr_cnu(0:PT%Nclasses,1:Maxcontracts) - ! - write(iunit) 'icontr_ideg' - ! - write(iunit) PT%icontr_ideg(0:PT%Nclasses,1:Maxcontracts) - ! + + call ioHandler%write(Maxcontracts) + call ioHandler%write('icontr_cnu') + call ioHandler%write(PT%icontr_cnu(0:PT%Nclasses,1:Maxcontracts)) + call ioHandler%write('icontr_ideg') + call ioHandler%write(PT%icontr_ideg(0:PT%Nclasses,1:Maxcontracts)) + case ('APPEND') ! - read(iunit) ncontr +#ifdef TROVE_USE_MPI_ + stop "APPEND in PTstore_icontr_cnu currently unsupported using MPI" +#endif + call ioHandler%read(ncontr) ! if (Maxcontracts/=ncontr) then write (out,"(' Vib. kinetic checkpoint file ',a)") job%kinetmat_file @@ -34926,21 +34865,21 @@ subroutine PTstore_icontr_cnu(Maxcontracts,iunit,dir) allocate (imat_t(0:PT%Nclasses,ncontr),stat=alloc) call ArrayStart('mat_t',alloc,size(imat_t),kind(imat_t)) ! - read(iunit) buf18(1:10) + 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' end if ! - read(iunit) imat_t(0:PT%Nclasses,1:ncontr) + call ioHandler%read(imat_t(0:PT%Nclasses,1:ncontr)) ! - read(iunit) buf18(1:11) + 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' end if ! - read(iunit) imat_t(0:PT%Nclasses,1:ncontr) + call ioHandler%read(imat_t(0:PT%Nclasses,1:ncontr)) ! deallocate(imat_t) ! @@ -34948,40 +34887,6 @@ subroutine PTstore_icontr_cnu(Maxcontracts,iunit,dir) ! end subroutine PTstore_icontr_cnu - subroutine PTstoreMPI_icontr_cnu(maxcontracts,iunit,dir) - use mpi_aux - - integer(ik),intent(in) :: maxcontracts - type(mpi_file),intent(in) :: iunit - character(len=18),intent(in) :: dir - integer(ik) :: alloc - character(len=18) :: buf18 - integer(ik) :: ncontr - integer(ik),allocatable :: imat_t(:,:) - integer::ierr - ! - select case(trim(dir)) - ! - case ('SAVE') - ! -#ifdef TROVE_USE_MPI_ - call mpi_file_write(iunit, maxcontracts, 1,mpi_integer,mpi_status_ignore,ierr) - ! - call mpi_file_write(iunit, 'icontr_cnu', 10,mpi_character,mpi_status_ignore,ierr) - ! - call mpi_file_write(iunit, pt%icontr_cnu(0:pt%nclasses,1:maxcontracts), (1+pt%nclasses)*maxcontracts, mpi_integer, & - mpi_status_ignore, ierr) - ! - call mpi_file_write(iunit, 'icontr_ideg', 11,mpi_character,mpi_status_ignore,ierr) - ! - call mpi_file_write(iunit, pt%icontr_ideg(0:pt%nclasses,1:maxcontracts), (1+pt%nclasses)*maxcontracts, mpi_integer, & - mpi_status_ignore, ierr) - -#endif - end select - - end subroutine PTstorempi_icontr_cnu - subroutine PTdefine_contr_from_eigenvect(nroots,Neigenlevels,eigen) @@ -38270,6 +38175,9 @@ subroutine PTcontracted_matelem_class_fast(jrot) type(me_type), allocatable :: vpot_me(:) type(me_type), allocatable :: pseudo_me(:) type(me_type), allocatable :: extF_me(:) + + class(ioHandlerBase), allocatable :: ioHandler + type(ErrorType) :: err ! call TimerStart('Contracted matelements-class-fast') ! @@ -38372,15 +38280,20 @@ subroutine PTcontracted_matelem_class_fast(jrot) job_is ='Vib. matrix elements of the rot. kinetic part' call IOStart(trim(job_is),chkptIO) ! - open(chkptIO,form='unformatted',action='write',position='rewind',status='replace',file=job%kinetmat_file) - write(chkptIO) 'Start Kinetic part' + allocate(ioHandlerFTN::ioHandler) + call ioHandler%open(& + 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 ! - call PTstore_icontr_cnu(PT%Maxcontracts,chkptIO,job%IOkinet_action) + call PTstore_icontr_cnu(PT%Maxcontracts,ioHandler,job%IOkinet_action) ! if (job%vib_rot_contr) then - write(chkptIO) 'vib-rot' + call ioHandler%write('vib-rot') endif ! endif @@ -38453,7 +38366,7 @@ subroutine PTcontracted_matelem_class_fast(jrot) ! if (trim(job%IOkinet_action)=='SAVE'.and..not.job%IOmatelem_split) then ! - write(chkptIO) 'g_rot' + call ioHandler%write('g_rot') ! endif ! @@ -38527,7 +38440,7 @@ subroutine PTcontracted_matelem_class_fast(jrot) if (job%IOmatelem_split) then write(chkptIO_) grot_t else - write(chkptIO) grot_t + call ioHandler%write(grot_t, mdimen) endif ! endif @@ -38551,7 +38464,7 @@ subroutine PTcontracted_matelem_class_fast(jrot) ! if (trim(job%IOkinet_action)=='SAVE'.and..not.job%IOmatelem_split) then ! - write(chkptIO) 'g_cor' + call ioHandler%write('g_cor') ! endif ! @@ -38668,7 +38581,7 @@ subroutine PTcontracted_matelem_class_fast(jrot) if (job%IOmatelem_split) then write(chkptIO_) grot_t else - write(chkptIO) grot_t + call ioHandler%write(grot_t, mdimen) endif ! endif @@ -38713,7 +38626,7 @@ subroutine PTcontracted_matelem_class_fast(jrot) if ((trim(job%IOkinet_action)=='SAVE'.or.trim(job%IOkinet_action)=='VIB_SAVE').and. & (.not.job%IOmatelem_split.or.job%iswap(1)==0.or.job%iswap(1)==(PT%Nmodes+3)*3 ) ) then ! - write(chkptIO) 'hvib' + call ioHandler%write('hivb') ! endif ! @@ -38848,10 +38761,8 @@ subroutine PTcontracted_matelem_class_fast(jrot) enddo !$omp end parallel do ! - !write(chkptIO) hvib_t - ! if (trim(job%IOkinet_action)=='SAVE') then - write(chkptIO) hvib_t + call ioHandler%write(hvib_t, mdimen) endif ! endif @@ -38888,8 +38799,8 @@ subroutine PTcontracted_matelem_class_fast(jrot) if ((trim(job%IOkinet_action)=='SAVE'.or.trim(job%IOkinet_action)=='VIB_SAVE').and.& (.not.job%IOmatelem_split.or.job%iswap(1)==0 )) then ! - write(chkptIO) 'End Kinetic part' - close(chkptIO,status='keep') + call ioHandler%write('End Kinetic part') + deallocate(ioHandler) ! endif ! diff --git a/test/compare_results.py b/test/regression/compare_results.py similarity index 100% rename from test/compare_results.py rename to test/regression/compare_results.py diff --git a/test/download_benchmarks.sh b/test/regression/download_benchmarks.sh similarity index 100% rename from test/download_benchmarks.sh rename to test/regression/download_benchmarks.sh diff --git a/test/run_regression_tests.sh b/test/regression/run_regression_tests.sh similarity index 98% rename from test/run_regression_tests.sh rename to test/regression/run_regression_tests.sh index 2fe003a..66e9852 100755 --- a/test/run_regression_tests.sh +++ b/test/regression/run_regression_tests.sh @@ -3,7 +3,7 @@ set -e exe_name=j-trove.x -exe=../$exe_name +exe=../../$exe_name # Use 1 process unless we have specified differently (e.g. in CI) nproc=${nproc:-1} diff --git a/test/scripts/H2CO/compare_results.sh b/test/regression/scripts/H2CO/compare_results.sh similarity index 100% rename from test/scripts/H2CO/compare_results.sh rename to test/regression/scripts/H2CO/compare_results.sh diff --git a/test/scripts/H2CO/run_benchmark.sh b/test/regression/scripts/H2CO/run_benchmark.sh similarity index 100% rename from test/scripts/H2CO/run_benchmark.sh rename to test/regression/scripts/H2CO/run_benchmark.sh diff --git a/test/scripts/set_io_format.sh b/test/regression/scripts/set_io_format.sh similarity index 100% rename from test/scripts/set_io_format.sh rename to test/regression/scripts/set_io_format.sh diff --git a/test/unit/.gitignore b/test/unit/.gitignore new file mode 100644 index 0000000..873ee29 --- /dev/null +++ b/test/unit/.gitignore @@ -0,0 +1,2 @@ +*.inc +*.F90 diff --git a/test/unit/makefile b/test/unit/makefile new file mode 100644 index 0000000..a7561ee --- /dev/null +++ b/test/unit/makefile @@ -0,0 +1,26 @@ +BASE_DIR := ../.. + +include $(BASE_DIR)/lib/pFUnit/build/PFUNIT.mk + +%.o : %.F90 + $(FC) -c $(FFLAGS) $< + +FFLAGS += $(PFUNIT_EXTRA_FFLAGS) +FFLAGS += -I../.. + +test_io_TESTS := test_io.pf +test_io_OTHER_SRCS := $(addprefix $(BASE_DIR)/, io_handler_ftn.f90 io_handler_base.f90 ) +$(eval $(call make_pfunit_test,test_io)) + +ifdef USE_MPI +FFLAGS += -DTROVE_USE_MPI_ +test_mpi_io_TESTS := test_mpi_io.pf +test_mpi_io_OTHER_LIBRARIES := ${LAPACK} # passed in from base makefile +test_mpi_io_OTHER_SRCS := $(addprefix $(BASE_DIR)/, io_handler_mpi.f90 io_handler_base.f90 mpi_aux.f90 timer.f90 accuracy.f90) +$(eval $(call make_pfunit_test,test_mpi_io)) +endif + +clean: + $(RM) *.o *.mod *.a *.inc + $(RM) test_io.F90 test_io + $(RM) test_mpi_io.F90 test_mpi_io diff --git a/test/unit/test_io.pf b/test/unit/test_io.pf new file mode 100644 index 0000000..749c34c --- /dev/null +++ b/test/unit/test_io.pf @@ -0,0 +1,164 @@ +#include "errors.fpp" + +module test_io + use funit + use io_handler_base + use io_handler_ftn + use errors + + implicit none + + contains + + @test + subroutine testFTNWriting() + type(ioHandlerFTN) :: ioHandler + + real, dimension(5) :: array1D + real, dimension(5, 5) :: array2D + real, dimension(5) :: in_array1D + real, dimension(5, 5) :: in_array2D + real :: true_real, in_real + complex :: true_complex, in_complex + integer :: true_integer, 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 + type(ErrorType) :: err + + call ioHandler%open(fname, err, \ + action='write', form=form, access=access, status=status, position=position) + HANDLE_ERROR(err) + + true_integer = 5 + call ioHandler%write(true_integer) ! int + + true_real = 4.0 + call ioHandler%write(true_real) ! double + + true_complex = (5.0, 1.0) + call ioHandler%write(true_complex) ! complex + + ! Array + array1D = (/2, 3, 4, 5, 6/) + call ioHandler%write(array1D) + + ! 2D array + do i=1,5 + do j=1,5 + array2D(i,j) = i+j + end do + end do + call ioHandler%write(array2D) + + call ioHandler%close() + + open(newunit=iounit, iostat=stat, action='read', file=fname, \ + form=form, access=access, status=status, position=position) + + read(iounit) in_integer + read(iounit) in_real + read(iounit) in_complex + read(iounit) in_array1D + read(iounit) in_array2D + + if (stat == 0) close(iounit, status='delete') + + @assertTrue(in_integer == true_integer) + @assertTrue(in_real == true_real) + @assertTrue(in_complex == true_complex) + do i=1,5 + @assertTrue(in_array1D(i) == array1D(i)) + end do + do i=1,5 + do j=1,5 + @assertTrue(in_array2D(i,j) == array2D(i,j)) + end do + end do + end subroutine + + @test + subroutine testFTNioHandler() + type(ioHandlerFTN) :: ioHandler + + real, dimension(5) :: array1D + real, dimension(5, 5) :: array2D + real, dimension(5) :: in_array1D + real, dimension(5, 5) :: in_array2D + real :: true_real, in_real + complex :: true_complex, in_complex + integer :: true_integer, 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 + type(ErrorType) :: err + + open(newunit=iounit, iostat=stat, action='write', file=fname, \ + form=form, access=access, status=status, position=position) + + true_integer = 5 + write(iounit) true_integer ! int + + true_real = 4.0 + write(iounit) true_real ! double + + true_complex = (5.0, 1.0) + write(iounit) true_complex ! complex + + ! Array + array1D = (/2, 3, 4, 5, 6/) + write(iounit) array1D ! 1D array + + ! 2D array + do i=1,5 + do j=1,5 + array2D(i,j) = i+j + end do + end do + write(iounit) array2D ! 2D array + + if (stat == 0) close(iounit) + + call ioHandler%open(fname, err, \ + action='read', form=form, access=access, status=status, position=position) + HANDLE_ERROR(err) + + call ioHandler%read(in_integer) + call ioHandler%read(in_real) + call ioHandler%read(in_complex) + call ioHandler%read(in_array1D) + call ioHandler%read(in_array2D) + + call ioHandler%close() + + @assertTrue(in_integer == true_integer) + @assertTrue(in_real == true_real) + @assertTrue(in_complex == true_complex) + do i=1,5 + @assertTrue(in_array1D(i) == array1D(i)) + end do + do i=1,5 + do j=1,5 + @assertTrue(in_array2D(i,j) == array2D(i,j)) + end do + end do + + ! Cleanup test file + open(newunit=iounit, iostat=stat, action='read', file=fname) + if (stat == 0) close(iounit, status='delete') + + end subroutine + +end module test_io diff --git a/test/unit/test_mpi_io.pf b/test/unit/test_mpi_io.pf new file mode 100644 index 0000000..537c347 --- /dev/null +++ b/test/unit/test_mpi_io.pf @@ -0,0 +1,276 @@ +#include "errors.fpp" + +module test_mpi_io + use mpi_f08 + use funit + use mpi_aux + use io_handler_base + use io_handler_mpi + use errors + use accuracy + + implicit none + + @TestCase + type, extends(TestCase) :: TestMPI + contains + procedure :: setUp ! overides generic + procedure :: tearDown ! overrides generic + end type TestMPI + + integer,external :: INDXL2G + + logical :: is_mpi_initialised = .false. + integer, parameter :: totalTestCount = 2 + integer :: currentTestCount = 0 + + contains + + subroutine setUp(this) + class(TestMPI), intent(inout) :: this + currentTestCount = currentTestCount + 1 + if (is_mpi_initialised .eqv. .false.) then + call co_init_comms() + is_mpi_initialised = .true. + endif + end subroutine setUp + + subroutine tearDown(this) + class(TestMPI), intent(inout) :: this + if (is_mpi_initialised .eqv. .true. .and. currentTestCount == totalTestCount) then + call co_finalize_comms() + endif + end subroutine tearDown + + @test + subroutine testMPIWriting(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 + 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 i, j + integer :: gi = 0, gj = 0, MB, NB, RSRC, CSRC + type(ErrorType) :: err + + integer :: ierr, rank, allocinfo = 0 + + call MPI_Comm_rank(MPI_COMM_WORLD, rank, ierr) + if(ierr.ne.0) print *, "Error: could not get rank" + + call ioHandler%open(fname, err, action='write', & + form=form, access=access, status=status, position=position) + HANDLE_ERROR(err) + + ! Test writing scalars + call ioHandler%write(true_integer) ! int + call ioHandler%write(true_real) ! double + call ioHandler%write(true_complex) ! complex + call ioHandler%write(true_str) ! string + + ! Test writing an 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" + + MB = array2D_descr(5) + NB = array2D_descr(6) + + RSRC = array2D_descr(7) + CSRC = array2D_descr(8) + + do i=1,size(array2D,1) + do j=1,size(array2D,2) + gi = INDXL2G (i, MB, myprow, RSRC, nprow) + gj = INDXL2G (j, NB, mypcol, CSRC, npcol) + array2D(i,j) = array2DNCol*(gi-1) + (gj) + end do + end do + + call ioHandler%write(array2D, array2D_descr, array2D_block_type) + + ! Test writing something after array + call ioHandler%write(true_integer) ! int + + ! 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 + 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 + print *, sizeof(in_str) + @assertTrue(in_str == true_str) + + read(iounit) in_array2D + do i=1,array2DNRow + do j=1,array2DNCol + @assertTrue(in_array2D(i,j) == array2DNCol*(i-1) + j) + end do + end do + + read(iounit) in_integer + @assertTrue(in_integer == true_integer) + + read(iounit) in_array2D + do i=1,array2DNRow + do j=1,array2DNCol + @assertTrue(in_array2D(i,j) == array2DNCol*(i-1) + j) + 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 + + @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 :: b, icoeff, jcoeff + + 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" + + real(rk), allocatable :: grot_t(:,:) + real(rk), allocatable :: grot_full(:,:) + real(rk),allocatable :: recvbuf(:,:,:) + + integer :: true_integer = 5, in_integer + + integer :: i,j + + type(MPI_File) :: chkptMPIIO + integer :: ierr, rank, allocinfo = 0 + + type(ErrorType) :: err + + 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 + + allocate(recvbuf(mdimen_p,mdimen_p,comm_size)) + allocate(grot_t(mdimen_b,startdim:startdim+mdimen_p-1)) + allocate(grot_full(dimen, 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 + grot_t(jcoeff,icoeff) = jcoeff*icoeff + enddo + enddo + endif + enddo + + ! Distribute between processes and save + call co_distr_data(grot_t, recvbuf, mdimen_p, startdim, enddim) + + call ioHandler%open(fname, err, action='write', & + form=form, access=access, status=status, position=position) + HANDLE_ERROR(err) + + ! Test writing something before array + call ioHandler%write(true_integer) + + ! Test writing an array + call ioHandler%write(grot_t, mdimen) + + ! Test writing something after array + call ioHandler%write(true_integer) + + ! Test writing another array + call ioHandler%write(grot_t, mdimen) + + ! Test writing something after second array + call ioHandler%write(true_integer) + + call ioHandler%close() + + ! Check output + 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) grot_full + do i=1,dimen + do j=1,dimen + @assertTrue(grot_full(i,j) == i*j) + end do + end do + + read(iounit) in_integer + @assertTrue(in_integer == true_integer) + + read(iounit) grot_full + do i=1,dimen + do j=1,dimen + @assertTrue(grot_full(i,j) == i*j) + end do + end do + + read(iounit) in_integer + @assertTrue(in_integer == true_integer) + + close(iounit, status='delete') + endif + + end subroutine testMPIWritingColumnDistArray +end module diff --git a/tran.f90 b/tran.f90 index ee92b9b..db13c45 100644 --- a/tran.f90 +++ b/tran.f90 @@ -2,6 +2,8 @@ ! defines different transformations of the eigenvectors given ! in terms of the contracted basis state representaion. ! +#include "errors.fpp" + module tran ! set tran_debug > 2 with small vibrational bases and small expansions only !#define tran_debug 1 @@ -14,7 +16,13 @@ module tran use moltype, only : manifold,intensity,extF use symmetry, only : sym - use perturbation, only : PTintcoeffsT,PTrotquantaT,PTNclasses,PTstore_icontr_cnu,PTeigenT,PTdefine_contr_from_eigenvect,PTrepresT,PTstorempi_icontr_cnu + use perturbation, only : PTintcoeffsT,PTrotquantaT,PTNclasses,PTstore_icontr_cnu,PTeigenT,PTdefine_contr_from_eigenvect,PTrepresT + use io_handler_base + use io_handler_ftn +#ifdef TROVE_USE_MPI_ + use io_handler_mpi +#endif + use errors private public read_contrind,read_eigenval, TReigenvec_unit, bset_contrT, & @@ -1318,6 +1326,8 @@ subroutine TRconvert_matel_j0_eigen(jrot) type(MPI_File) :: fileh, fileh_w integer(kind=MPI_OFFSET_KIND) :: mpioffset,read_offset,write_offset integer :: ierr + class(ioHandlerBase), allocatable :: ioHandler + type(ErrorType) :: err type(MPI_Datatype) :: gmat_block_type, psi_block_type, mat_t_block_type, mat_s_block_type, extF_block_type integer,dimension(9) :: desc_gmat, desc_mat_t, desc_mat_s, desc_psi, desc_extF @@ -1602,46 +1612,24 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! job_is ='Eigen-vib. matrix elements of the rot. kinetic part' ! - if (trim(job%kinetmat_format).eq.'MPIIO') then + call IOStart(trim(job_is),chkptIO) #ifdef TROVE_USE_MPI_ - call MPI_File_open(mpi_comm_world, job%kineteigen_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(fileh_w, '[MPIIO]', 7, mpi_character, mpi_status_ignore, ierr) - call MPI_File_write(fileh_w, 'Start Kinetic part', 18, mpi_character, mpi_status_ignore, ierr) - ! - treat_vibration = .false. - ! - call PTstorempi_icontr_cnu(Neigenroots,fileh_w,job%IOj0matel_action) - ! - if (job%vib_rot_contr) then - call MPI_File_write(fileh_w, 'vib-rot', 7, mpi_character, mpi_status_ignore, ierr) - endif - else - mpioffset = 0 - treat_vibration = .false. - endif + allocate(ioHandlerMPI::ioHandler) +#else + allocate(ioHandlerFTN::ioHandler) #endif - else - call IOStart(trim(job_is),chkptIO) - ! - open(chkptIO,form='unformatted',action='write',position='rewind',status='replace',file=job%kineteigen_file) - write(chkptIO) 'Start Kinetic part' - ! - treat_vibration = .false. - if(mpi_rank .eq. 0) then - ! - call PTstore_icontr_cnu(Neigenroots,chkptIO,job%IOj0matel_action) - ! - if (job%vib_rot_contr) then - write(chkptIO) 'vib-rot' - endif - endif + call ioHandler%open(& + job%kineteigen_file, err, & + action='write', position='rewind', status='replace', form='unformatted') + HANDLE_ERROR(err) + call ioHandler%write('Start Kinetic part') + + treat_vibration = .false. + call PTstore_icontr_cnu(Neigenroots,ioHandler,job%IOj0matel_action) + if (job%vib_rot_contr) then + call ioHandler%write('vib-rot') endif - ! + ! endif ! rootsize = int(bset_contr(1)%Maxcontracts*(bset_contr(1)%Maxcontracts+1)/2,hik) @@ -1686,19 +1674,9 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! task = 'rot' ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - call MPI_File_seek_shared(fileh_w, int(0,MPI_OFFSET_KIND), MPI_SEEK_END) - if(mpi_rank.eq.0) call MPI_File_write_shared(fileh_w, 'g_rot', 5, mpi_character, mpi_status_ignore, ierr) - call mpi_barrier(MPI_COMM_WORLD, ierr) + call ioHandler%write('g_rot') ! - call restore_rot_kinetic_matrix_elements_mpi(jrot,treat_vibration,task,fileh) -#endif - else - write(chkptIO) 'g_rot' - ! - call restore_rot_kinetic_matrix_elements(jrot,treat_vibration,task,iunit) - endif + call restore_rot_kinetic_matrix_elements(jrot,treat_vibration,task,iunit) ! else ! @@ -1722,10 +1700,6 @@ subroutine TRconvert_matel_j0_eigen(jrot) #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) - - !call MPI_File_seek_shared(fileh_w, int(0,MPI_OFFSET_KIND),MPI_SEEK_END) - 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 ! @@ -1790,13 +1764,7 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! else ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - call MPI_File_write_all(fileh_w, mat_s, size(mat_s), mpi_double_precision, mpi_status_ignore, ierr) -#endif - else - write (chkptIO) mat_s - endif + call ioHandler%write(mat_s, desc_mat_s, mat_s_block_type) ! endif ! @@ -1810,13 +1778,6 @@ subroutine TRconvert_matel_j0_eigen(jrot) 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) - - write_offset = write_offset + 9*int(Neigenroots,MPI_OFFSET_KIND)*Neigenroots*mpi_real_size - !call MPI_File_set_view(fileh_w, write_offset, mpi_byte, mpi_byte, "native", MPI_INFO_NULL, ierr) - 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 @@ -1834,15 +1795,12 @@ subroutine TRconvert_matel_j0_eigen(jrot) #ifdef TROVE_USE_MPI_ call restore_rot_kinetic_matrix_elements_mpi(jrot,treat_vibration,task,fileh) ! - if(mpi_rank.eq.0) call MPI_File_write_shared(fileh_w, 'g_cor', 5, mpi_character, 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_Barrier(MPI_COMM_WORLD, ierr) ! May no longer be needed? #endif else call restore_rot_kinetic_matrix_elements(jrot,treat_vibration,task,iunit) - ! - write(chkptIO) 'g_cor' endif + call ioHandler%write('g_cor') ! endif ! @@ -1856,9 +1814,6 @@ subroutine TRconvert_matel_j0_eigen(jrot) #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) - - 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 ! @@ -1925,13 +1880,7 @@ subroutine TRconvert_matel_j0_eigen(jrot) ! else ! - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - call MPI_File_write_all(fileh_w, mat_s, size(mat_s), mpi_double_precision, mpi_status_ignore, ierr) -#endif - else - write (chkptIO) mat_s - endif + call ioHandler%write(mat_s, desc_mat_s, mat_s_block_type) ! endif ! @@ -1945,35 +1894,17 @@ subroutine TRconvert_matel_j0_eigen(jrot) 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) - - write_offset = write_offset + 3*int(Neigenroots,MPI_OFFSET_KIND)*Neigenroots*mpi_real_size - !call MPI_File_set_view(fileh_w, write_offset, mpi_byte, mpi_byte, "native", MPI_INFO_NULL, ierr) - 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_END) - !write_offset = 0 #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 - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - call MPI_File_write_shared(fileh_w, 'End Kinetic part', 16, mpi_character, mpi_status_ignore, ierr) -#endif - else - write(chkptIO) 'End Kinetic part' - endif + call ioHandler%write('End Kinetic part') endif ! - if (.not.job%vib_rot_contr) then - if (trim(job%kinetmat_format).eq.'MPIIO') then -#ifdef TROVE_USE_MPI_ - call MPI_File_close(fileh_w, ierr) -#endif - else - close(chkptIO,status='keep') - endif + if (.not.job%vib_rot_contr) then + deallocate(ioHandler) endif ! task = 'end' @@ -1996,7 +1927,10 @@ subroutine TRconvert_matel_j0_eigen(jrot) call MPI_File_close(fileh, ierr) #endif else - close(chkptIO,status='keep') + ! 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 ! ! External field part