diff --git a/f90/3D_MT/CMakeLists.txt b/f90/3D_MT/CMakeLists.txt index 0911f4a6..64db53c6 100644 --- a/f90/3D_MT/CMakeLists.txt +++ b/f90/3D_MT/CMakeLists.txt @@ -24,6 +24,7 @@ target_sources(${MODEM_EXE} #PRIVATE DataIO_HDF5.f90 PRIVATE DataIO.f90 PRIVATE EMfieldInterp.f90 + PRIVATE ESolnManager.f90 PRIVATE ForwardSolver.f90 PRIVATE GridCalc.f90 #PRIVATE GridCalcS.f90 diff --git a/f90/3D_MT/ESolnManager.f90 b/f90/3D_MT/ESolnManager.f90 new file mode 100644 index 00000000..9e38f962 --- /dev/null +++ b/f90/3D_MT/ESolnManager.f90 @@ -0,0 +1,667 @@ +! EsolnManager Module +! ==================== +! +! The EsolnManager module is a module that helps with transferring electric +! and adjoint fields to and from worker tasks. It does so in two ways: +! +! 1. MPI Transferring - Original ModEM MPI Behavior +! +! If specified, main tasks and worker tasks will communicate their +! electric fields/adjoint solutions to each other using MPI send and +! receives. +! +! This will cause the main task to allocate space for eAll. +! +! 2. Saving Electric fields/Adjoin Solutions to File +! +! If specified, the worker tasks will read and write out electric fields +! when they are told to access them by the main task rather than communicating +! the electric fields back to the main task. +! +! +! How to use this Module: +! +! To use this module, you will want to first call EsMgr_init with either: +! +! save_in_file = .true. or .false. depending on how you wish electric fields to +! be saved. +! +! Regardless of value of save_in_file, you will need to call EsMgr_get and EsMgr_save +! in order to complete a 'transaction'. +! +! With save_in_file = .false., EsMgr_get will call a MPI_Recv and EsMgr_save will call +! a corresponding MPI_Send +! +! With save_in_file = .ture., the worker tasks will write out their e fields to disk via +! EsMgr_save and will read e fields using EsMgr_get. The main task will still need to call +! EsMgr_get/EsMgr_save as the worker task will communicate to the main task when it is done +! reading/writing. + +module ESolnManager + + use mpi + use Declaration_MPI + use GridDef + use SolnSpace + use utilities + + implicit None + + private + + type (ModEM_mpi_context_t), pointer :: EsMgr_ctx + + character(len=*), parameter :: FTYPE_ASCII = "ascii" + character(len=*), parameter :: FTYPE_BINARY = "binary" + + character(len=*), parameter :: E_FIELD_TYPE_FWD = "FWD" + character(len=*), parameter :: E_FIELD_TYPE_JMULTT = "JmultT" + + type (grid_t), pointer :: EsMgr_grid => null() + character(len=25) :: EsMgr_ftype + character(len=256) :: EsMgr_prefix + logical :: EsMgr_save_in_file + + character, pointer, dimension(:) :: esmgr_holder => null() + integer :: esmgr_holder_bytes + logical :: esmgr_holder_allocated + + public :: FTYPE_ASCII, FTYPE_BINARY + public :: E_FIELD_TYPE_FWD, E_FIELD_TYPE_JMULTT + public :: EsMgr_init + public :: EsMgr_create_solnVectorMTX, EsMgr_create_e + public :: EsMgr_get + public :: EsMgr_save + public :: EsMgr_save_in_file + +contains + + ! EsMgr_init - Initializes the ESolnManager sets several module variables + ! and determines if the Esoln files should be saved or not. The main task and all + ! and worker tasks must call this before calling any other subroutines in this + ! module. + ! + ! Arguments: + ! * grid - grid_t target + ! The grid + ! * context - ModEM_mpi_context_t target + ! An initalized ModEM_mpi_context_t which the ESolnManager should use + ! when doing MPI communication routines. + ! * save_in_file - logical - optional - Default: .false. + ! If save_in_file is present and .true., the electric fields will be saved in a file + ! rather than passed from the main task and worker tasks via MPI. + ! + ! If .false. or not present, the electric fields will be passed from the main and + ! worker task via MPI. + ! * prefix - character - optional - Default: 'esoln' + ! A character string prefix which to label the electric field/adjoint solution + ! filenames if they are written out to disk. + ! + ! If save_in_file is .false. or not present, prefix has no effect. + ! * ftype - character - optional - Default: 'FTYPE_BINARY' + ! Determines the type of file the cvectors should be written as if `save_in_file` + ! is true. + ! + ! Options are above in this module and are: + ! * FTYPE_ASCII - ASCII representation + ! * FTYPE_BINARY - Fortran unformatted binary + subroutine EsMgr_init(grid, context, save_in_file, prefix, ftype) + + implicit none + + type (grid_t), target, intent(in) :: grid + type (ModEM_mpi_context_t), target, intent(in) :: context + logical, intent(in), optional :: save_in_file + character(len=*), intent(in), optional :: prefix + character(len=*), intent(in), optional :: ftype + + logical :: save_in_file_lcl + character(len=256) :: prefix_lcl + character(len=256) :: ftype_lcl + + EsMgr_ctx => context + + if (present(save_in_file)) then + save_in_file_lcl = save_in_file + else + save_in_file_lcl = .false. + end if + + if ( .not. (save_in_file_lcl .and. present(prefix)) ) then + if (EsMgr_ctx % rank_world == 0) then + write(0,*) "Warning: Argument 'prefix' was passed, but 'save_in_file' was not present" + write(0,*) "Warning: 'prefix' will not have an effect. Set 'save_in_file' to true to save" + write(0,*) "Warning: esolns in files" + end if + end if + + if ( .not. (save_in_file_lcl .and. present(ftype)) ) then + if (EsMgr_ctx % rank_world == 0) then + write(0,*) "Warning: Argument 'ftype' was passed, but 'save_in_file' was not present" + write(0,*) "Warning: 'ftype' will not have an effect. Set 'save_in_file' to true when calling" + write(0,*) "WARNING: EsMgr_init to save esolns in files" + end if + end if + + if (present(prefix)) then + prefix_lcl = prefix + else + prefix_lcl = "esoln" + end if + + if (present(ftype)) then + ftype_lcl = ftype + else + ftype_lcl = FTYPE_ASCII + end if + + select case (ftype_lcl) + case (FTYPE_ASCII : FTYPE_BINARY) + case DEFAULT + write(0,*) "ERROR: ", trim(ftype_lcl), " is not a valid file type for Esoln" + write(0,*) "ERROR: Valid options are: [", trim(FTYPE_ASCII), " | ", trim(FTYPE_BINARY), "]" + call ModEM_abort() + end select + + EsMgr_grid => grid + EsMgr_save_in_file = save_in_file_lcl + EsMgr_prefix = prefix_lcl + EsMgr_ftype = ftype_lcl + + call MPI_Barrier(MPI_COMM_WORLD, ierr) + + end subroutine EsMgr_init + + ! EsMgr_create_eAll + ! + ! Convenience function to initialize an solnVectorMTX_t with solnVectors. + ! + ! If EsMgr_init was called with save_in_file=.true., then the solnVectors_t + ! created in this subroutine will be created as 'place_holders' and their + ! pol cvectors variables will also be created as 'place_holders'. These + ! cvectors will not have their x, y, z arrays will not be allocated. + ! + ! However, if save_in_file=.false. then these arrays will be allocated. + ! + ! Note: EsMgr_init should be called before calling this function. + ! + ! Arguments: + ! eAll - solnVectorMTX_T + ! The solnVectorMTX_t to hold eAll. See description above. + ! nTx - integer + ! The number of transmitters to create for the SolnVectorMTX_t + ! grid - grid_t target + ! The grid that should be associated with this SolnVectorMTX_t + ! place_holder - logical - Defaults to save_in_file passed in during init + ! Overrides the default action described above. + ! + subroutine EsMgr_create_solnVectorMTX(eAll, nTx, grid, place_holder) + + implicit none + + type (solnVectorMTX_t), intent(inout) :: eAll + integer, intent(in) :: nTx + type(grid_t), intent(in), target :: grid + logical, intent(in), optional :: place_holder + + integer :: iTx + type(solnVector_t) :: e + logical :: place_holder_lcl + + if (present(place_holder)) then + place_holder_lcl = place_holder + else + place_holder_lcl = EsMgr_save_in_file + end if + + call create_solnVectorMTX(nTx, eAll) + + do iTx = 1, nTx + call create_solnVector(grid, iTx, eall % solns(iTx), place_holder=place_holder_lcl) + end do + + end subroutine EsMgr_create_solnVectorMTX + + ! EsMgr_create_e + ! + ! Convenience function to create a single solnVector_t. + ! + ! If EsMgr_init was called with save_in_file=.true., then the solnVector_t + ! created in this subroutine by the main task will be created as 'place_holders' + ! and their pol cvectors variables will also be created as 'place_holders'. + ! These cvectors will not have their x, y, z arrays will not be allocated. + ! + ! If save_in_file=.true. and this routine is called by a worker task, then + ! the this will routine will completly allocate the cvectors associated with + ! the solnVector_t. + ! + ! However, if save_in_file=.false. then these arrays will be allocated. + ! + ! Note: EsMgr_init should be called before calling this function. + ! + ! Arguments: + ! e - solnVector_t + ! The solnVector_t that will be created + ! iTx - integer + ! The transmitter number to associate with this solnVector_t + ! place_holder - logical - Defaults to save_in_file passed in during init + ! Overides the default action described above. + ! + subroutine EsMgr_create_e(e, iTx, place_holder) + + implicit none + + type (solnVector_t), intent(inout) :: e + integer, intent(in) :: iTx + logical, intent(in), optional :: place_holder + + logical :: place_holder_lcl + + if (present(place_holder)) then + place_holder_lcl = place_holder + else + if (EsMgr_ctx % rank_current == 0) then + place_holder_lcl = EsMgr_save_in_file + else + ! This is a worker task - Don't make this a place holder + place_holder_lcl = .false. + end if + end if + + call create_solnVector(EsMgr_grid, iTx, e, place_holder=place_holder_lcl) + + end subroutine EsMgr_create_e + + ! create_prefix - Generate a prefix for a solnVector file name + ! + ! This routine will create a prefix that can be used to name the + ! solnVector file. + ! + ! The prefix will always start with EsMgr_prefix, which was + ! passed in during the init call, the rest will be in one of + ! depending on which optional arguments are present or not: + ! + ! - $EsMgr_prefix.$E_field_type.$label + ! - $EsMgr_prefix.$E_field_type + ! - $EsMgr_prefix.$label + ! - $EsMgr_prefix + ! + subroutine create_prefix(prefix, E_field_type, label) + + implicit none + + character(len=*), intent(out) :: prefix + character(len=*), intent(in), optional :: E_field_type + character(len=*), intent(in), optional :: label + + character(len=512) :: PREFIX_FNAME_FORMAT + logical :: use_label + + use_label = .false. + + if (present(label)) then + use_label = trim(label) /= "NULL" .and. len_trim(label) /= 0 + end if + + if (present(E_field_type) .and. use_label) then + PREFIX_FNAME_FORMAT = '(A, A, A, A, A)' + write(prefix, PREFIX_FNAME_FORMAT) trim(EsMgr_prefix), ".", trim(E_field_type), ".", trim(label) + else if (present(E_field_type) .and. .not. use_label) then + PREFIX_FNAME_FORMAT = '(A, A, A)' + write(prefix, PREFIX_FNAME_FORMAT) trim(EsMgr_prefix), ".", trim(E_field_type) + else if (.not. present(E_field_type) .and. use_label) then + PREFIX_FNAME_FORMAT = '(A, A, A)' + write(prefix, PREFIX_FNAME_FORMAT) trim(EsMgr_prefix), ".", trim(label) + else + PREFIX_FNAME_FORMAT = '(A, A)' + write(prefix, PREFIX_FNAME_FORMAT) trim(EsMgr_prefix) + end if + + end subroutine create_prefix + + ! EsMgr_get + ! + ! Populate e with the transmitter iTx and polarization pol_index + ! via MPI or by reading a file. + ! + ! EsMgr_save_in_file == .false. + ! ====================== + ! + ! If `EsMgr_save_in_file == .false.`, then a task will receive the + ! solnVector via MPI from the MPI tasks specified in the `from`, + ! argument. + ! + ! Note: The task in the `from` argument must call EsMgr_save + ! with a corosponding task number in the `to` argument. + ! + ! EsMgr_save_in_file == .true. + ! ====================== + ! + ! If the EsolnManager was initialized with `EsMgr_save_in_file == .true.`, + ! then worker_tasks that call this function will read the corresponding + ! file that matches the requested arguments. + ! + ! After reading, the worker task will communicate with the main task + ! to indicate that it is done reading. + ! + ! Main Tasks: + ! + ! If EsMgr_save_in_file is .true., then the main task will wait for + ! the task speciried in from to finish writing it's solnVector_t + ! that it made in ESoln_save. + ! + ! Note: EsMgr_init should be called before calling this function. + ! + ! Arguments: + ! e - solnVector_t + ! The solnVector_t to read the polarization pol_index of iTx into + ! iTx - integer + ! The transmitter number to read + ! pol_index - integer - optional + ! The polarization to read + ! from - integer - optional + ! The MPI task to receive the solnVector_t from if EsMgr_save_in_file is .false. + ! E_field_type - character - optional + ! The type of the electric field that the SolnVector was saved with, + ! normally this is either E_FIELD_TYPE_FWD or E_FIELD_TYPE_JMULTT. + ! label - character - optional + ! The additonal label that was used to save this SolnVector. Normally + ! this is used to label the specific line search calulation (e.g. + ! trial, fwd01, fwd02 etc.). + subroutine EsMgr_get(e, iTx, pol_index, from, E_field_type, label) + + implicit none + + type (solnVector_t), intent(inout) :: e + integer, intent(in) :: iTx + integer, intent(in), optional :: pol_index + integer, intent(in), optional :: from + character(len=*), intent(in), optional :: E_field_type + character(len=*), intent(in), optional :: label + character(len=512) :: prefix + + e % tx = iTx + + ! If we are reading and writing files, do nothing + if (EsMgr_save_in_file .and. EsMgr_ctx % rank_world == 0) then + call wait_on_task(from) + return + end if + + if (EsMgr_save_in_file .and. EsMgr_ctx % rank_world /= 0) then + call create_prefix(prefix, E_field_type, label) + call read_esoln_from_file(e, iTx, pol_index, prefix=prefix) + return + end if + + call EsMgr_recv_e(e, from) + + end subroutine EsMgr_get + + ! EsMgr_save + ! + ! Save `e` by either writing it to disk or sending it to the main task via MPI. + ! + ! EsMgr_save_in_file = .false. + ! ============================= + ! + ! If `EsMgr_save_in_file == .false.`, then a calling task will send `e` to the task + ! specified in the `to` argument. + ! + ! To receive this `e`, a task will need to call EsMgr_get with the correspond from. + ! + ! EsMgr_save_in_file = .true. + ! ============================= + ! + ! If EsMgr_save_in_file is .true., then when a worker task that calls this + ! subroutine will write out e to a corresponding file. The worker task will then + ! communicate to the main task to declare that it is done writing. + ! + ! If the main task calls this routine with EsMgr_save_in_file == .true., then + ! it will return immediately as no work is needed. + ! + ! Arguments + ! e - solnVector_t + ! The solnVector_t to save + ! to - integer - optional + ! The task to send e to + ! E_field_type - character - optional + ! The type of the electric field, normally this is either E_FIELD_TYPE_FWD + ! or E_FIELD_TYPE_JMULTT, but could be anything. However, the same E_field_type + ! will need to be used in EsMgr_get to read the same solnVector. + ! label - character - optional + ! An additional label to use to label the solnVector file (if + ! EsMgr_save_in_file == .true.). Normally, this is used to label the specific + ! line search calculations (e.g. trial, fwd01, fwd02, etc.). Again, the same + ! label will need to be specified when calling the coorosponding EsMgr_get. + ! + subroutine EsMgr_save(e, to, E_field_type, label) + + implicit none + + type (solnVector_t), intent(inout) :: e + integer, intent(in), optional :: to + character(len=*), intent(in), optional :: E_field_type + character(len=*), intent(in), optional :: label + + character(len=512) :: prefix + + if (EsMgr_save_in_file .and. EsMgr_ctx % rank_world == 0) then + return + end if + + if (EsMgr_save_in_file .and. .not. EsMgr_ctx % rank_world == 0) then + call create_prefix(prefix, E_field_type, label) + call EsMgr_write_to_file(e, prefix) + call communicate_file_done_writing() + return + end if + + call EsMgr_send_e(e, to) + + end subroutine EsMgr_save + + subroutine communicate_file_done_writing() + + implicit none + + integer :: dummy_buffer + integer :: tag + integer :: count + integer :: to + + count = 0 + tag = 0 + to = 0 + + call MPI_Send(dummy_buffer, count, MPI_INTEGER, to, tag, EsMgr_ctx % comm_current, ierr) + + end subroutine communicate_file_done_writing + + subroutine wait_on_task(task) + + implicit none + + integer, intent(in) :: task + + integer :: dummy_buff + integer :: count + integer :: tag + + count = 0 + tag = 0 + + call MPI_Recv(dummy_buff, count, MPI_INTEGER, task, tag, EsMgr_ctx % comm_current, MPI_STATUS_IGNORE, ierr) + + end subroutine wait_on_task + + subroutine EsMgr_write_to_file(e, prefix, iPol, ftype) + + implicit none + + type (solnVector_t), intent(in) :: e + character(len=*), intent(in), optional :: prefix + integer, intent(in), optional :: iPol + character(len=*), optional, intent(in) :: ftype + + integer :: iPol_lcl + character(len=256) :: prefix_lcl + character(len=256) :: ftype_lcl + + if (present(iPol)) then + iPol_lcl = iPol + else + iPol_lcl = 1 + end if + + if (present(prefix)) then + prefix_lcl = prefix + else + prefix_lcl = "esoln" + end if + + if (present(ftype)) then + ftype_lcl = ftype + else + ftype_lcl = EsMgr_ftype + endif + + call write_solnVector(e, trim(prefix_lcl), ftype=ftype_lcl, pol_index=iPol_lcl) + + end subroutine EsMgr_write_to_file + + subroutine read_esoln_from_file(e, iTx, iPol, prefix) + + implicit none + + type (solnVector_t), intent(inout) :: e + integer, intent(in) :: iTx + integer, intent(in), optional :: iPol + character(len=*), intent(in), optional :: prefix + + integer :: iPol_lcl + character(len=256) :: prefix_lcl + + if (present(iPol)) then + iPol_lcl = iPol + else + iPol_lcl = 1 + end if + + if (present(prefix)) then + prefix_lcl = prefix + else + prefix_lcl = "esoln" + end if + + e % tx = iTx + call read_solnVector(e, trim(prefix_lcl), ftype=EsMgr_ftype, pol_index=iPol_lcl) + + end subroutine read_esoln_from_file + + subroutine EsMgr_send_e(e, to) + + implicit none + + type (solnVector_t), intent(inout) :: e + integer, intent(in) :: to + integer :: tag + character, pointer, dimension(:) :: e_packed => null() + + if (to == 0) then + tag = FROM_WORKER + else + tag = FROM_MASTER + end if + + call int_create_e_param_place_holder(e) + call int_Pack_packed_vec(e) + call MPI_Send(esmgr_holder, esmgr_holder_bytes, MPI_PACKED, to, tag, EsMgr_ctx % comm_current, ierr) + + end subroutine EsMgr_send_e + + subroutine EsMgr_recv_e(e, from) + + implicit none + + type (solnVector_t), intent(inout) :: e + integer, intent(in) :: from + integer :: tag + + if (from == 0) then + tag = FROM_MASTER + else + tag = FROM_WORKER + end if + + call int_create_e_param_place_holder(e) + call MPI_Recv(esmgr_holder, esmgr_holder_bytes, MPI_PACKED, from, tag, EsMgr_ctx % comm_current, STATUS, ierr) + call int_Unpack_packed_vec(e) + + end subroutine EsMgr_recv_e + + subroutine int_create_e_param_place_holder(e) + + implicit none + type(solnVector_t), intent(in) :: e + integer :: Ex_size,Ey_size,Ez_size,Nbytes1,Nbytes2,Nbytes3,Nbytes4 + + if (esmgr_holder_allocated) then + return + end if + + Ex_size=size(e%pol(1)%x) + CALL MPI_PACK_SIZE(Ex_size, MPI_DOUBLE_COMPLEX, MPI_COMM_WORLD, Nbytes1, ierr) + Ey_size=size(e%pol(1)%y) + CALL MPI_PACK_SIZE(Ey_size, MPI_DOUBLE_COMPLEX, MPI_COMM_WORLD, Nbytes2, ierr) + Ez_size=size(e%pol(1)%z) + CALL MPI_PACK_SIZE(Ez_size, MPI_DOUBLE_COMPLEX, MPI_COMM_WORLD, Nbytes3, ierr) + CALL MPI_PACK_SIZE(1, MPI_INTEGER, MPI_COMM_WORLD, Nbytes4, ierr) + esmgr_holder_bytes=((Nbytes1+Nbytes2+Nbytes3+Nbytes4))+1 + + allocate(esmgr_holder(esmgr_holder_bytes)) + esmgr_holder_allocated = .true. + + end subroutine int_create_e_param_place_holder + + subroutine int_pack_packed_vec(e) + + implicit none + + type(solnVector_t), intent(in) :: e + integer index,Ex_size,Ey_size,Ez_size + + Ex_size=size(e%pol(1)%x) + Ey_size=size(e%pol(1)%y) + Ez_size=size(e%pol(1)%z) + index=1 + + call MPI_Pack(e%pol(which_pol)%x(1,1,1),Ex_size, MPI_DOUBLE_COMPLEX, esmgr_holder, esmgr_holder_bytes, index, MPI_COMM_WORLD, ierr) + call MPI_Pack(e%pol(which_pol)%y(1,1,1),Ey_size, MPI_DOUBLE_COMPLEX, esmgr_holder, esmgr_holder_bytes, index, MPI_COMM_WORLD, ierr) + call MPI_Pack(e%pol(which_pol)%z(1,1,1),Ez_size, MPI_DOUBLE_COMPLEX, esmgr_holder, esmgr_holder_bytes, index, MPI_COMM_WORLD, ierr) + call MPI_Pack(e%tx, 1, MPI_INTEGER, esmgr_holder, esmgr_holder_bytes, index, MPI_COMM_WORLD, ierr) + + end subroutine int_pack_packed_vec + + + subroutine int_unpack_packed_vec(e) + + implicit none + + type(solnVector_t), intent(inout) :: e + + integer index,Ex_size,Ey_size,Ez_size + + Ex_size=size(e%pol(1)%x) + Ey_size=size(e%pol(1)%y) + Ez_size=size(e%pol(1)%z) + index=1 + + call MPI_Unpack(esmgr_holder, esmgr_holder_bytes, index, e%pol(which_pol)%x(1,1,1),Ex_size, MPI_DOUBLE_COMPLEX,MPI_COMM_WORLD, ierr) + call MPI_Unpack(esmgr_holder, esmgr_holder_bytes, index, e%pol(which_pol)%y(1,1,1),Ey_size, MPI_DOUBLE_COMPLEX,MPI_COMM_WORLD, ierr) + call MPI_Unpack(esmgr_holder, esmgr_holder_bytes, index, e%pol(which_pol)%z(1,1,1),Ez_size, MPI_DOUBLE_COMPLEX,MPI_COMM_WORLD, ierr) + call MPI_Unpack(esmgr_holder, esmgr_holder_bytes, index, e%tx,1, MPI_INTEGER,MPI_COMM_WORLD, ierr) + end subroutine int_unpack_packed_vec + + +end module ESolnManager diff --git a/f90/3D_MT/SolnSpace.f90 b/f90/3D_MT/SolnSpace.f90 index 63a41115..397fc83b 100755 --- a/f90/3D_MT/SolnSpace.f90 +++ b/f90/3D_MT/SolnSpace.f90 @@ -45,7 +45,6 @@ module SolnSpace MODULE PROCEDURE dotProd_sparseVsolnV end interface - type :: solnVector_t !! Generic solution type, same name must be used to allow !! use of higher level inversion modules on different problems. @@ -81,6 +80,8 @@ module SolnSpace !! avoid memory leaks: set this to true for function outputs only logical :: temporary = .false. + logical :: place_holder = .false. + end type solnVector_t type :: solnVectorMTX_t @@ -165,7 +166,7 @@ module SolnSpace ! Basic solnVector methods !********************************************************************** - subroutine create_solnVector(grid,iTx,e) + subroutine create_solnVector(grid, iTx, e, place_holder) ! generic routine for creating the solnVector type for 3D problems: ! number of polarization obtained from the transmitter dictionary @@ -174,9 +175,17 @@ subroutine create_solnVector(grid,iTx,e) type(grid_t), intent(in), target :: grid integer, intent(in) :: iTx type (solnVector_t), intent(inout) :: e + logical, intent(in), optional :: place_holder ! local variables integer :: k,istat,iPol + logical :: place_holder_lcl + + if (present(place_holder)) then + place_holder_lcl = place_holder + else + place_holder_lcl = .false. + end if if (e%allocated) then if (associated(e%grid, target=grid) .and. (e%tx == iTx)) then @@ -190,11 +199,11 @@ subroutine create_solnVector(grid,iTx,e) e%nPol = txDict(iTx)%nPol allocate(e%Pol_index(e%nPol), STAT=istat) allocate(e%Pol_name(e%nPol), STAT=istat) - + do iPol=1,e%nPol e%Pol_index(iPol)=iPol end do - + ! set up the mode names based on transmitter type; ! for now, only set up for MT. Used for I/O. if (trim(txDict(iTx)%tx_type) .eq. 'MT') then @@ -207,9 +216,12 @@ subroutine create_solnVector(grid,iTx,e) end if allocate(e%pol(e%nPol), STAT=istat) + do k = 1,e%nPol - call create_cvector(grid,e%pol(k),EDGE) + call create_cvector(grid,e%pol(k),EDGE, place_holder=place_holder_lcl) enddo + + e % place_holder = place_holder_lcl e%tx = iTx e%grid => grid @@ -343,6 +355,175 @@ function dotProd_solnVector(FV1,FV2,Conj_Case) result(c) end function dotProd_solnVector + !********************************************************************** + subroutine write_solnVector(e, prefix, pol_index, ftype) + + implicit none + + type (solnVector_t), intent(in) :: e + character(len=*), intent(in) :: prefix + integer, optional, intent(in) :: pol_index + character(len=*), optional, intent(in) :: ftype + + character(len=512) :: fname, ftype_lcl, form + character(len=512) :: iomsg + integer :: iostat + integer :: fid + integer :: pol_index_lcl + + if (.not. e % allocated) then + write(0,*) "ERROR: solnVector_t (argument e) msut be allocated before calling write_solnVector" + write(0,*) "ERROR: Allocate it by calling `create_solnVector` first" + call ModEM_abort() + end if + + if (present(pol_index)) then + pol_index_lcl = pol_index + else + pol_index_lcl = 1 + end if + + if (present(ftype)) then + ftype_lcl = ftype + else + ftype_lcl = 'binary' + end if + + if (ftype_lcl == 'binary') then + form = 'unformatted' + else if (ftype_lcl == 'ascii') then + form = 'formatted' + else + ftype_lcl = 'binary' + form = 'unformatted' + end if + + fname = construct_esoln_fname(prefix, e % tx, trim(e % pol_name(e % pol_index(pol_index_lcl)))) + open(newunit=fid, file=trim(fname), action='write', form=form, status='replace', & + iostat=iostat, iomsg=iomsg) + if (iostat /= 0) then + write(0,'(A, A, A)') "ERROR: There was an issue when opening: '", trim(fname), "' for writing:" + write(0,'(A, A)') "ERROR: Reason: ", trim(iomsg) + call ModEM_abort() + end if + + write(6, '(A, i4.4, A, A, A, A)') "Saving electric solution for Tx: ", e % tx, " pol: '", & + trim(e % pol_name(pol_index_lcl)), "' to file: ", trim(fname) + + call write_cvector(fid, e % pol(pol_index_lcl), ftype_lcl) + close(fid) + + end subroutine write_solnVector + + !********************************************************************** + subroutine read_solnVector(e, prefix, pol_index, ftype) + + implicit none + + type (solnVector_t), intent(inout) :: e + character(len=*), intent(in) :: prefix + integer, optional, intent(in) :: pol_index + character(len=*), optional, intent(in) :: ftype + + character(len=512) :: fname, ftype_lcl, form + character(len=512) :: iomsg + integer :: fid + integer :: pol_index_lcl + integer :: iostat + logical :: file_exists + + if (.not. e % allocated) then + write(0,*) "ERROR: solnVector_t (argument e) msut be allocated before calling read_solnVector" + write(0,*) "ERROR: Allocate it by calling `create_solnVector` first" + call ModEM_abort() + end if + + if (present(pol_index)) then + pol_index_lcl = pol_index + else + pol_index_lcl = 1 + end if + + if (present(ftype)) then + ftype_lcl = ftype + else + ftype_lcl = 'binary' + end if + + if (ftype_lcl == 'binary') then + form = 'unformatted' + else if (ftype_lcl == 'ascii') then + form = 'formatted' + else + ftype_lcl = 'binary' + form = 'unformatted' + end if + + fname = construct_esoln_fname(prefix, e % tx, trim(e % pol_name(pol_index_lcl))) + + if (.not. does_esoln_file_exist(e, prefix)) then + write(0,*) "ERROR: The file for this solnVector_t (argument e) does not exist" + write(0,*) "ERROR: Could not find: ", trim(fname) + call ModEM_abort() + end if + + open(newunit=fid, file=trim(fname), action='read', form=form, status='old', & + iostat=iostat, iomsg=iomsg) + if (iostat /= 0) then + write(0,'(A, A, A)') "ERROR: There was an issue when opening: '", trim(fname), "' for reading:" + write(0,'(A, A)') "ERROR: Reason: ", trim(iomsg) + call ModEM_abort() + end if + + call read_cvector(fid, e % pol(pol_index_lcl), ftype_lcl) + close(fid) + + end subroutine read_solnVector + + !********************************************************************** + function does_esoln_file_exist(e, prefix, pol_index) result(file_exists) + + implicit none + + type (solnVector_t), intent(in) :: e + character(len=*), intent(in) :: prefix + integer, optional, intent(in) :: pol_index + character(len=512) :: fname + logical :: file_exists + + integer :: pol_index_lcl + + if (present(pol_index)) then + pol_index_lcl = pol_index + else + pol_index_lcl = 1 + end if + + if (.not. e % allocated) then + write(0,*) "ERROR: solnVector_t (argument e) msut be allocated before calling does_esoln_file_eixst" + write(0,*) "ERROR: Allocate it by calling `create_solnVector` first" + call ModEM_abort() + end if + + fname = construct_esoln_fname(prefix, e % tx, trim(e % pol_name(pol_index_lcl))) + inquire(file=trim(fname), exist=file_exists) + + end function does_esoln_file_exist + + !********************************************************************** + function construct_esoln_fname(prefix, iTx, pol_name) result(fname) + + implicit none + + character(len=*), intent(in) :: prefix + integer, intent(in) :: iTx + character(len=*), intent(in) :: pol_name + character(len=512) :: fname + + write(fname, '(A, A, I4.4, A, A, A)') trim(prefix), '.iTx.', iTx, '.', trim(pol_name), '.cvec' + + end function construct_esoln_fname + !********************************************************************** ! Basic solnVectorMTX methods !********************************************************************** diff --git a/f90/3D_MT/Sub_MPI.f90 b/f90/3D_MT/Sub_MPI.f90 index facdd40c..ba52c6fe 100644 --- a/f90/3D_MT/Sub_MPI.f90 +++ b/f90/3D_MT/Sub_MPI.f90 @@ -71,6 +71,18 @@ subroutine get_nPol_MPI(emsoln) nPol_MPI= emsoln%nPol end subroutine get_nPol_MPI + +function get_nPol(iPer) result(nPol) + + implicit none + + integer, intent(in) :: iPer + integer :: nPol + nPol = 0 + + nPol = txDict(iPer) % nPol + +end function get_nPol !**************************************************************************** subroutine count_number_of_messages_to_RECV(eAll1) diff --git a/f90/FIELDS/FiniteDiff3D/sg_vector.f90 b/f90/FIELDS/FiniteDiff3D/sg_vector.f90 index 6eb8e7c0..60294ea3 100755 --- a/f90/FIELDS/FiniteDiff3D/sg_vector.f90 +++ b/f90/FIELDS/FiniteDiff3D/sg_vector.f90 @@ -200,6 +200,8 @@ module sg_vector ! ISO/IEC 15581 - the "allocatable array extension") logical :: temporary = .false. + logical :: place_holder = .false. + ! pointer to parent grid type (grid_t), pointer :: grid @@ -415,16 +417,24 @@ end subroutine create_rvector ! create_rvector ! using grid definition in structure "grid" ; ! allocates memory in x,y,z component arrays ! gridType is a character string to describe intended usage - subroutine create_cvector(igrid, E, gridType) + subroutine create_cvector(igrid, E, gridType, place_holder) implicit none type(grid_t), target, intent(in) :: igrid ! the grid for which an edge/ face node field is being initialized type (cvector), intent(inout) :: E + logical, intent(in), optional :: place_holder integer :: status,nx,ny,nz character (len=80), intent(in) :: gridType + logical :: place_holder_lcl + + if (present(place_holder)) then + place_holder_lcl = place_holder + else + place_holder_lcl = .false. + end if ! First deallocate anything, that's allocated call deall_cvector(E) @@ -443,6 +453,17 @@ subroutine create_cvector(igrid, E, gridType) ! gridType E%gridType = gridType + E%place_holder = place_holder_lcl + + ! Initialize everything to null + E % x => null() + E % y => null() + E % z => null() + + if (E % place_holder .eqv. .true.) then + ! Don't allocate anything for place holders + return + end if ! allocate memory for x,y,z ; ! E%allocated will be true if all allocations succeed @@ -654,6 +675,7 @@ subroutine write_cvector(fid, E, ftype) complex (kind(E%x)), allocatable, dimension(:,:,:) :: x, y, z logical :: ok, hasname, binary character(80) :: fname, isbinary + character(512) :: iomsg if(.not. E%allocated) then write(0, *) 'cvector must be allocated before call to write_cvector' @@ -679,10 +701,30 @@ subroutine write_cvector(fid, E, ftype) ! write binary to unformatted files if (binary) then - write(fid) E%nx,E%ny,E%nz,E%gridType - write(fid) E%x - write(fid) E%y - write(fid) E%z + write(fid, iostat=istat, iomsg=iomsg) E%nx,E%ny,E%nz,E%gridType + if (istat /= 0) then + write(0,*) "ERROR: There was an writing file '", trim(fname), "' header" + write(0,*) "ERROR: Reason: ", trim(iomsg) + call ModEM_abort() + end if + write(fid, iostat=istat, iomsg=iomsg) E%x + if (istat /= 0) then + write(0,*) "ERROR: There was an writing file '", trim(fname), "' E % x" + write(0,*) "ERROR: Reason: ", trim(iomsg) + call ModEM_abort() + end if + write(fid, iostat=istat, iomsg=iomsg) E%y + if (istat /= 0) then + write(0,*) "ERROR: There was an writing file '", trim(fname), "' E % y" + write(0,*) "ERROR: Reason: ", trim(iomsg) + call ModEM_abort() + end if + write(fid, iostat=istat, iomsg=iomsg) E%z + if (istat /= 0) then + write(0,*) "ERROR: There was an writing file '", trim(fname), "' E % z" + write(0,*) "ERROR: Reason: ", trim(iomsg) + call ModEM_abort() + end if return end if @@ -848,6 +890,7 @@ subroutine read_cvector(fid, E, ftype) complex (kind(E%x)), allocatable, dimension(:,:,:) :: x, y, z logical :: ok, hasname, binary character(80) :: fname, isbinary + character(512) :: iomsg if (.not. present(ftype)) then binary = .false. @@ -887,9 +930,38 @@ subroutine read_cvector(fid, E, ftype) if (binary) then ! read binary from unformatted files - read(fid) E%x - read(fid) E%y - read(fid) E%z + read(fid, iostat=istat, iomsg=iomsg) E%x + if (istat < 0) then + write(1, '(A, A, A)') "ERROR: EOF reached when reading E % x in file: ", trim(fname), ' aborting!' + write(0,*) 'ERROR: Error read iostat code:', istat + call ModEM_abort() + else if (istat > 0) then + write(0, '(A, A, A)') "ERROR: Erorr when reading E % x: ", trim(fname), ' aborting!' + write(0, '(A, i4.4, A, A)') 'ERROR: Error read iostat code:', istat, ' Message: ', trim(iomsg) + call ModEM_abort() + end if + + read(fid, iostat=istat, iomsg=iomsg) E%y + if (istat < 0) then + write(0, '(A, A, A)') "ERROR: EOF reached when reading E % y in file: ", trim(fname), ' aborting!' + write(0, '(A, i4.4, A, A)') 'ERROR: Error read iostat code:', istat, ' Message: ', trim(iomsg) + call ModEM_abort() + else if (istat > 0) then + write(0, '(A, A, A)') "ERROR: Erorr when reading E % y: ", trim(fname), ' aborting!' + write(0, '(A, i4.4, A, A)') 'ERROR: Error read iostat code:', istat, ' Message: ', trim(iomsg) + call ModEM_abort() + end if + + read(fid, iostat=istat, iomsg=iomsg) E%z + if (istat < 0) then + write(0, '(A, A, A)') "ERROR: EOF reached when reading E % z in file: ", trim(fname), ' aborting!' + write(0, '(A, i4.4, A, A)') 'ERROR: Error read iostat code:', istat, ' Message: ', trim(iomsg) + call ModEM_abort() + else if (istat > 0) then + write(0, '(A, A, A)') "ERROR: Erorr when reading E % z: ", trim(fname), ' aborting!' + write(0, '(A, i4.4, A, A)') 'ERROR: Error read iostat code:', istat, ' Message: ', trim(iomsg) + call ModEM_abort() + end if return end if @@ -1003,7 +1075,7 @@ subroutine copy_cvector(E2, E1) integer :: status ! check to see if RHS (E1) is active (allocated) - if(.not.E1%allocated) then + if(.not.E1%allocated .and. .not. E1 % place_holder) then write(0,*) 'RHS not allocated yet for copy_cvector' else @@ -1012,11 +1084,15 @@ subroutine copy_cvector(E2, E1) if (E1%gridType == E2%gridType) then ! just copy components - E2%x = E1%x - E2%y = E1%y - E2%z = E1%z + if (.not. E1 % place_holder) then + E2%x = E1%x + E2%y = E1%y + E2%z = E1%z + end if + E2%gridType = E1%gridType E2%grid => E1%grid + E2%place_holder = E1%place_holder else write (0, *) 'not compatible usage for copy_cvector' @@ -1030,13 +1106,18 @@ subroutine copy_cvector(E2, E1) end if ! then allocate E2 as correct size ... - Call create_cvector(E1%grid, E2, E1%gridType) + Call create_cvector(E1%grid, E2, E1%gridType, place_holder=E1 % place_holder) ! .... and copy E1 - E2%x = E1%x - E2%y = E1%y - E2%z = E1%z + + if (.not. E1 % place_holder) then + E2%x = E1%x + E2%y = E1%y + E2%z = E1%z + end if + E2%gridType = E1%gridType E2%grid => E1%grid + E2%place_holder = E1%place_holder end if diff --git a/f90/INV/INVcore.f90 b/f90/INV/INVcore.f90 index 40269338..78a0b429 100755 --- a/f90/INV/INVcore.f90 +++ b/f90/INV/INVcore.f90 @@ -62,7 +62,7 @@ end subroutine printf !********************************************************************** - subroutine func(lambda,d,m0,mHat,F,mNorm,dHat,eAll,RMS) + subroutine func(lambda,d,m0,mHat,F,mNorm,dHat,eAll,RMS,label) ! Compute the full penalty functional F ! Also output the predicted data and the EM solution @@ -76,12 +76,20 @@ subroutine func(lambda,d,m0,mHat,F,mNorm,dHat,eAll,RMS) type(dataVectorMTX_t), optional, intent(inout) :: dHat type(solnVectorMTX_t), optional, intent(inout) :: eAll real(kind=prec), optional, intent(out) :: RMS + character(len=*), optional, intent(in) :: label ! local variables type(dataVectorMTX_t) :: res,Nres type(modelParam_t) :: m,JTd real(kind=prec) :: SS integer :: Ndata, Nmodel + character(len=80) :: label_lcl + + if (present(label)) then + label_lcl = label + else + label_lcl = "NULL" + end if ! compute the smoothed model parameter vector call CmSqrtMult(mHat,m) @@ -96,7 +104,7 @@ subroutine func(lambda,d,m0,mHat,F,mNorm,dHat,eAll,RMS) ! also sets up forward solutions for all transmitters in eAll ! (which is created on the fly if it doesn't exist) #ifdef MPI - call Master_Job_fwdPred(m,dHat,eAll) + call Master_Job_fwdPred(m,dHat,eAll,label=label_lcl) #else call fwdPred(m,dHat,eAll) #endif @@ -138,7 +146,7 @@ subroutine func(lambda,d,m0,mHat,F,mNorm,dHat,eAll,RMS) end subroutine func !********************************************************************** - subroutine gradient(lambda,d,m0,mHat,grad,dHat,eAll) + subroutine gradient(lambda,d,m0,mHat,grad,dHat,eAll,label) ! Computes the gradient of the penalty functional, ! using EM solution (eAll) and the predicted data (dHat) @@ -157,13 +165,19 @@ subroutine gradient(lambda,d,m0,mHat,grad,dHat,eAll) type(modelParam_t), intent(inout) :: grad type(dataVectorMTX_t), intent(inout) :: dHat type(solnVectorMTX_t), intent(inout) :: eAll + character(len=*), intent(in), optional :: label ! local variables real(kind=prec) :: Ndata,Nmodel type(dataVectorMTX_t) :: res type(modelParam_t) :: m,JTd,CmJTd + character(len=80) :: label_lcl - ! integer :: j, Ny, NzEarth + if (present(label)) then + label_lcl = label + else + label_lcl = "NULL" + end if ! compute the smoothed model parameter vector call CmSqrtMult(mHat,m) @@ -180,7 +194,7 @@ subroutine gradient(lambda,d,m0,mHat,grad,dHat,eAll) ! multiply by J^T call CdInvMult(res) #ifdef MPI - call Master_job_JmultT(m,res,JTd,eAll) + call Master_job_JmultT(m,res,JTd,eAll,label=label) #else call JmultT(m,res,JTd,eAll) #endif diff --git a/f90/INV/LBFGS.f90 b/f90/INV/LBFGS.f90 index f6f7b11b..c9f005d9 100644 --- a/f90/INV/LBFGS.f90 +++ b/f90/INV/LBFGS.f90 @@ -64,6 +64,10 @@ module LBFGS type(LBFGSiterControl_t), private, save :: iterControl type(LBFGSiterCache_t), private, save :: iterCache + character(len=80), parameter :: STARTUP_LABEL = 'STARTUP' + character(len=80), parameter :: TRIAL_FWD_F_1_LABEL = 'F_1' + character(len=80), parameter :: FWD_F_LABEL = 'F' + Contains !********************************************************************** @@ -604,6 +608,7 @@ subroutine lineSearchQuadratic(lambda,d,m0,h,alpha,mHat,f,grad, & type(dataVectorMTX_t) :: dHat_1 type(solnVectorMTX_t) :: eAll_1 character(100) :: logFile + character(80) :: label ! parameters c = iterControl%c @@ -642,7 +647,7 @@ subroutine lineSearchQuadratic(lambda,d,m0,h,alpha,mHat,f,grad, & call linComb(ONE,mHat_0,alpha_1,h,mHat_1) ! compute the penalty functional and predicted data at mHat_1 - call func(lambda,d,m0,mHat_1,f_1,mNorm_1,dHat_1,eAll_1,rms_1) + call func(lambda,d,m0,mHat_1,f_1,mNorm_1,dHat_1,eAll_1,rms_1,label=TRIAL_FWD_F_1_LABEL) call printf('STARTLS',lambda,alpha,f_1,mNorm_1,rms_1) call printf('STARTLS',lambda,alpha,f_1,mNorm_1,rms_1,logFile) niter = niter + 1 @@ -671,7 +676,7 @@ subroutine lineSearchQuadratic(lambda,d,m0,h,alpha,mHat,f,grad, & ! alpha = alpha_i/TWO ! reset alpha to ensure progress ! end if call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('QUADLS',lambda,alpha,f,mNorm,rms) call printf('QUADLS',lambda,alpha,f,mNorm,rms,logFile) niter = niter + 1 @@ -701,22 +706,25 @@ subroutine lineSearchQuadratic(lambda,d,m0,h,alpha,mHat,f,grad, & end if if (starting_guess) then - alpha = alpha_1 - dHat = dHat_1 - eAll = eAll_1 - mHat = mHat_1 - rms = rms_1 - f = f_1 + label=TRIAL_FWD_F_1_LABEL + alpha = alpha_1 + dHat = dHat_1 + eAll = eAll_1 + mHat = mHat_1 + rms = rms_1 + f = f_1 + else + label=FWD_F_LABEL end if ! compute gradient of the full penalty functional and exit if (relaxation) then - call linComb(ONE,mHat_0,gamma*alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) - call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms) - call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms,logFile) - end if - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call linComb(ONE,mHat_0,gamma*alpha,h,mHat) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms) + call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms,logFile) + end if + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=label) write(*,'(a39)') 'Gradient computed, line search finished' write(ioLog,'(a39)') 'Gradient computed, line search finished' @@ -798,6 +806,7 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & type(dataVectorMTX_t) :: dHat_1 type(solnVectorMTX_t) :: eAll_1 character(100) :: logFile + character(80) :: label ! parameters c = iterControl%c @@ -832,7 +841,7 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & ! compute the trial mHat, f, dHat, eAll, rms mHat_1 = mHat_0 call linComb(ONE,mHat_0,alpha_1,h,mHat_1) - call func(lambda,d,m0,mHat_1,f_1,mNorm_1,dHat_1,eAll_1,rms_1) + call func(lambda,d,m0,mHat_1,f_1,mNorm_1,dHat_1,eAll_1,rms_1,label=TRIAL_FWD_F_1_LABEL) call printf('STARTLS',lambda,alpha,f_1,mNorm_1,rms_1) call printf('STARTLS',lambda,alpha,f_1,mNorm_1,rms_1,logFile) niter = niter + 1 @@ -860,7 +869,7 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms) call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms,logFile) end if - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=TRIAL_FWD_F_1_LABEL) write(*,'(a45)') 'Quadratic has no minimum, exiting line search' write(ioLog,'(a45)') 'Quadratic has no minimum, exiting line search' call deall_dataVectorMTX(dHat_1) @@ -873,21 +882,24 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & ! otherwise compute the functional at the minimizer of the quadratic alpha = - b/(TWO*a) call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('QUADLS',lambda,alpha,f,mNorm,rms) call printf('QUADLS',lambda,alpha,f,mNorm,rms,logFile) niter = niter + 1 ! check whether the solution satisfies the sufficient decrease condition if (f < f_0 + c * alpha * g_0) then ! if the initial guess was better than what we found, take it - if (f_1 < f) then - starting_guess = .true. - alpha = alpha_1 - dHat = dHat_1 - eAll = eAll_1 - mHat = mHat_1 - rms = rms_1 - f = f_1 + if (f_1 < f) then + label = TRIAL_FWD_F_1_LABEL + starting_guess = .true. + alpha = alpha_1 + dHat = dHat_1 + eAll = eAll_1 + mHat = mHat_1 + rms = rms_1 + f = f_1 + else + label = FWD_F_LABEL end if ! compute the gradient and exit if (relaxation) then @@ -897,7 +909,7 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms,logFile) end if - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=label) write(*,'(a60)') 'Sufficient decrease condition satisfied, exiting line search' write(ioLog,'(a60)') 'Sufficient decrease condition satisfied, exiting line search' call deall_dataVectorMTX(dHat_1) @@ -939,7 +951,7 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & ! end if ! compute the penalty functional call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('CUBICLS',lambda,alpha,f,mNorm,rms) call printf('CUBICLS',lambda,alpha,f,mNorm,rms,logFile) niter = niter + 1 @@ -967,12 +979,15 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & ! if the initial guess was better than what we found, take it if (starting_guess) then - alpha = alpha_1 - dHat = dHat_1 - eAll = eAll_1 - mHat = mHat_1 - rms = rms_1 - f = f_1 + label = TRIAL_FWD_F_1_LABEL + alpha = alpha_1 + dHat = dHat_1 + eAll = eAll_1 + mHat = mHat_1 + rms = rms_1 + f = f_1 + else + label = FWD_F_LABEL end if ! compute gradient of the full penalty functional and exit @@ -982,7 +997,7 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms) call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms,logFile) end if - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=label) write(*,'(a39)') 'Gradient computed, line search finished' write(ioLog,'(a39)') 'Gradient computed, line search finished' @@ -1132,7 +1147,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & mHat_1 = mHat_0 ! mHat_1 = mHat_0 + dir*step call linComb(ONE,mHat_0,alpha_1,h,mHat_1) - call func(lambda,d,m0,mHat_1,f_1,mNorm_1,dHat_1,eAll_1,rms_1) + call func(lambda,d,m0,mHat_1,f_1,mNorm_1,dHat_1,eAll_1,rms_1,label=TRIAL_FWD_F_1_LABEL) call printf('STARTLS',lambda,alpha_1,f_1,mNorm_1,rms_1) call printf('STARTLS',lambda,alpha_1,f_1,mNorm_1,rms_1,logFile) niter = niter + 1 @@ -1161,7 +1176,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & rms = rms_1 f = f_1 ! compute the gradient and exit - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=TRIAL_FWD_F_1_LABEL) write(*,'(a45)') 'Quadratic has no minimum, exiting line search' write(ioLog,'(a45)') 'Quadratic has no minimum, exiting line search' if (present(flag)) then @@ -1175,7 +1190,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & end if ! otherwise compute the functional at the minimizer of the quadratic call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('QUADLS',lambda,alpha,f,mNorm,rms) call printf('QUADLS',lambda,alpha,f,mNorm,rms,logFile) niter = niter + 1 @@ -1183,7 +1198,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & ! Strong Wolfe's condition needs the gradient ! well, we are going to calculate it anyway - so why don't we do it now? if (f <= f_1) then - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=FWD_F_LABEL) g_1 = dotProd(grad, h) write(*,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0 write(*,'(a4,es12.5)') ' g1=',g_1 @@ -1216,7 +1231,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & else if (f_1 < f_0) then! is the initial making any progress? ! Test if the initial guess is good for Strong Wolfe condition - call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1) + call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1,label=TRIAL_FWD_F_1_LABEL) g_1 = dotProd(grad, h) write(*,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0 write(*,'(a4,es12.5)') ' g1=',g_1 @@ -1317,12 +1332,12 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & ! end if ! compute the penalty functional call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('CUBICLS',lambda,alpha,f,mNorm,rms) call printf('CUBICLS',lambda,alpha,f,mNorm,rms,logFile) niter = niter + 1 ibracket = ibracket + 1 - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=FWD_F_LABEL) g_1 = dotProd(grad, h) write(*,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0 write(*,'(a4,es12.5)') ' g1=',g_1 @@ -1367,7 +1382,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & starting_guess = .false. else starting_guess = .true. - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=TRIAL_FWD_F_1_LABEL) endif exit endif @@ -1382,7 +1397,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & starting_guess = .false. else starting_guess = .true. - call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1) + call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1,label=TRIAL_FWD_F_1_LABEL) endif exit endif @@ -1538,7 +1553,7 @@ subroutine lineSearchWolfe2(lambda,d,m0,h,alpha,mHat,f,grad, & ! mHat_1 = mHat_0 + dir*step call linComb(ONE,mHat_0,alpha_1,h,mHat_1) ! compute the trial m, f, dHat, eAll, rms - call func(lambda,d,m0,mHat_1,f_1,mNorm,dHat_1,eAll_1,rms_1) + call func(lambda,d,m0,mHat_1,f_1,mNorm,dHat_1,eAll_1,rms_1,label=trim(TRIAL_FWD_F_1_LABEL)) call printf('STARTLS',lambda,alpha_1,f_1,mNorm,rms_1) call printf('STARTLS',lambda,alpha_1,f_1,mNorm,rms_1,logFile) niter = niter + 1 @@ -1552,7 +1567,7 @@ subroutine lineSearchWolfe2(lambda,d,m0,h,alpha,mHat,f,grad, & call ModEM_abort() end if ! calculate the gradient at the first guess - call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1) + call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1,label=trim(TRIAL_FWD_F_1_LABEL)) g_1 = dotProd(grad, h) grad_1 = grad write(*,'(a29,es12.5)',advance='no') & @@ -1678,12 +1693,12 @@ subroutine lineSearchWolfe2(lambda,d,m0,h,alpha,mHat,f,grad, & gPrev = g ! evaluate the function and derivetive at the new alpha call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=trim(FWD_F_LABEL)) call printf('CUBICLS',lambda,alpha,f,mNorm,rms) call printf('CUBICLS',lambda,alpha,f,mNorm,rms,logFile) niter = niter + 1 nbracket = nbracket +1 - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=trim(FWD_F_LABEL)) g = dotProd(grad, h) if (f.lt.f_1) then alpha_1 = alpha @@ -1744,7 +1759,7 @@ subroutine lineSearchWolfe2(lambda,d,m0,h,alpha,mHat,f,grad, & alpha = alphaNext ! compute the penalty functional call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('QUADLS',lambda,alpha,f,mNorm,rms) call printf('QUADLS',lambda,alpha,f,mNorm,rms,logFile) else if ((istrapped .eq. 2).and.(nbracket.eq.0)) then ! try quadratic @@ -1754,7 +1769,7 @@ subroutine lineSearchWolfe2(lambda,d,m0,h,alpha,mHat,f,grad, & alpha = alphaNext ! compute the penalty functional call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('QUADLS',lambda,alpha,f,mNorm,rms) call printf('QUADLS',lambda,alpha,f,mNorm,rms,logFile) else if ((istrapped .eq. -1).and.(nbracket.eq.0)) then ! jump @@ -1764,7 +1779,7 @@ subroutine lineSearchWolfe2(lambda,d,m0,h,alpha,mHat,f,grad, & alpha = alphaNext ! compute the penalty functional call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('QUADLS',lambda,alpha,f,mNorm,rms) call printf('QUADLS',lambda,alpha,f,mNorm,rms,logFile) else ! cubic @@ -1774,7 +1789,7 @@ subroutine lineSearchWolfe2(lambda,d,m0,h,alpha,mHat,f,grad, & alpha = alphaNext ! compute the penalty functional call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('CUBICLS',lambda,alpha,f,mNorm,rms) call printf('CUBICLS',lambda,alpha,f,mNorm,rms,logFile) endif @@ -1784,7 +1799,7 @@ subroutine lineSearchWolfe2(lambda,d,m0,h,alpha,mHat,f,grad, & fPrev = f gPrev = g !calculatie gradient to test the Wolfe condition - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=FWD_F_LABEL) g = dotProd(grad, h) if (f.lt.f_1) then alpha_1 = alpha diff --git a/f90/INV/NLCG.f90 b/f90/INV/NLCG.f90 index bb367757..563f60ee 100755 --- a/f90/INV/NLCG.f90 +++ b/f90/INV/NLCG.f90 @@ -53,6 +53,10 @@ module NLCG type(NLCGiterControl_t), private, save :: iterControl + character(len=*), parameter :: STARTUP_LABEL = 'STARTUP' + character(len=*), parameter :: TRIAL_FWD_F_1_LABEL = 'F_1' + character(len=*), parameter :: FWD_F_LABEL = 'F' + Contains !********************************************************************** @@ -284,7 +288,6 @@ subroutine NLCGsolver(d,lambda,m0,m,fname) write(ioLog,'(a41,es8.1)') 'The initial damping parameter lambda is ',lambda write(ioLog,'(a55,f12.6)') 'The initial line search step size (in model units) is ',startdm - ! starting from the prior hardcoded by setting mHat = 0 and m = m0 ! m = m0 ! mHat = m0 @@ -559,14 +562,15 @@ subroutine lineSearchQuadratic(lambda,d,m0,h,alpha,mHat,f,grad, & ! local variables real(kind=prec) :: alpha_1,alpha_i,mNorm - logical :: starting_guess - logical :: relaxation + logical :: starting_guess + logical :: relaxation real(kind=prec) :: eps,k,c,a,b real(kind=prec) :: g_0,f_0,f_1,f_i,rms_1,mNorm_1 - type(modelParam_t) :: mHat_0,mHat_1 - type(dataVectorMTX_t) :: dHat_1 - type(solnVectorMTX_t) :: eAll_1 - character(100) :: logFile + type(modelParam_t) :: mHat_0,mHat_1 + type(dataVectorMTX_t) :: dHat_1 + type(solnVectorMTX_t) :: eAll_1 + character(100) :: logFile + character(len=80) :: label ! parameters c = iterControl%c @@ -601,11 +605,11 @@ subroutine lineSearchQuadratic(lambda,d,m0,h,alpha,mHat,f,grad, & ! initialize mHat_1 = mHat_0 - ! compute the trial parameter mHat_1 + ! compute the parameter mHat_1 call linComb(ONE,mHat_0,alpha_1,h,mHat_1) ! compute the penalty functional and predicted data at mHat_1 - call func(lambda,d,m0,mHat_1,f_1,mNorm_1,dHat_1,eAll_1,rms_1) + call func(lambda,d,m0,mHat_1,f_1,mNorm_1,dHat_1,eAll_1,rms_1,label=TRIAL_FWD_F_1_LABEL) call printf('STARTLS',lambda,alpha,f_1,mNorm_1,rms_1) call printf('STARTLS',lambda,alpha,f_1,mNorm_1,rms_1,logFile) niter = niter + 1 @@ -634,7 +638,7 @@ subroutine lineSearchQuadratic(lambda,d,m0,h,alpha,mHat,f,grad, & ! alpha = alpha_i/TWO ! reset alpha to ensure progress ! end if call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('QUADLS',lambda,alpha,f,mNorm,rms) call printf('QUADLS',lambda,alpha,f,mNorm,rms,logFile) niter = niter + 1 @@ -664,12 +668,15 @@ subroutine lineSearchQuadratic(lambda,d,m0,h,alpha,mHat,f,grad, & end if if (starting_guess) then - alpha = alpha_1 - dHat = dHat_1 - eAll = eAll_1 - mHat = mHat_1 - rms = rms_1 - f = f_1 + alpha = alpha_1 + dHat = dHat_1 + eAll = eAll_1 + mHat = mHat_1 + label = TRIAL_FWD_F_1_LABEL + rms = rms_1 + f = f_1 + else + label=FWD_F_LABEL end if ! compute gradient of the full penalty functional and exit @@ -679,7 +686,7 @@ subroutine lineSearchQuadratic(lambda,d,m0,h,alpha,mHat,f,grad, & call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms) call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms,logFile) end if - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll, label=label) write(*,'(a39)') 'Gradient computed, line search finished' write(ioLog,'(a39)') 'Gradient computed, line search finished' @@ -752,14 +759,15 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & ! local variables real(kind=prec) :: alpha_1,alpha_i,alpha_j,mNorm - logical :: starting_guess - logical :: relaxation + logical :: starting_guess + logical :: relaxation real(kind=prec) :: eps,k,c,a,b,q1,q2,q3 real(kind=prec) :: g_0,f_0,f_1,f_i,f_j,rms_1,mNorm_1 - type(modelParam_t) :: mHat_0,mHat_1 - type(dataVectorMTX_t) :: dHat_1 - type(solnVectorMTX_t) :: eAll_1 - character(100) :: logFile + type(modelParam_t) :: mHat_0,mHat_1 + type(dataVectorMTX_t) :: dHat_1 + type(solnVectorMTX_t) :: eAll_1 + character(100) :: logFile + character(len=80) :: label ! parameters c = iterControl%c @@ -791,10 +799,10 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & relaxation = .false. end if - ! compute the trial mHat, f, dHat, eAll, rms + ! compute the mHat, f, dHat, eAll, rms mHat_1 = mHat_0 call linComb(ONE,mHat_0,alpha_1,h,mHat_1) - call func(lambda,d,m0,mHat_1,f_1,mNorm_1,dHat_1,eAll_1,rms_1) + call func(lambda,d,m0,mHat_1,f_1,mNorm_1,dHat_1,eAll_1,rms_1, label=TRIAL_FWD_F_1_LABEL) call printf('STARTLS',lambda,alpha,f_1,mNorm_1,rms_1) call printf('STARTLS',lambda,alpha,f_1,mNorm_1,rms_1,logFile) niter = niter + 1 @@ -807,22 +815,26 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & a = (f_1 - f_0 - g_0*alpha_1)/(alpha_1**2) b = g_0 ! if the curvature is -ve, there is no minimum; take the initial guess + label = FWD_F_LABEL if (a < 0) then - starting_guess = .true. - alpha = alpha_1 - dHat = dHat_1 + write(0,*) " USE STARTING UGESS" + stop + starting_guess = .true. + label = TRIAL_FWD_F_1_LABEL + alpha = alpha_1 + dHat = dHat_1 eAll = eAll_1 - mHat = mHat_1 - rms = rms_1 - f = f_1 + mHat = mHat_1 + rms = rms_1 + f = f_1 ! compute the gradient and exit if (relaxation) then call linComb(ONE,mHat_0,gamma*alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms) call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms,logFile) end if - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=label) write(*,'(a45)') 'Quadratic has no minimum, exiting line search' write(ioLog,'(a45)') 'Quadratic has no minimum, exiting line search' call deall_dataVectorMTX(dHat_1) @@ -835,22 +847,26 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & ! otherwise compute the functional at the minimizer of the quadratic alpha = - b/(TWO*a) call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('QUADLS',lambda,alpha,f,mNorm,rms) call printf('QUADLS',lambda,alpha,f,mNorm,rms,logFile) niter = niter + 1 ! check whether the solution satisfies the sufficient decrease condition if (f < f_0 + c * alpha * g_0) then ! if the initial guess was better than what we found, take it - if (f_1 < f) then - starting_guess = .true. - alpha = alpha_1 - dHat = dHat_1 - eAll = eAll_1 - mHat = mHat_1 - rms = rms_1 - f = f_1 + if (f_1 < f) then + starting_guess = .true. + label = TRIAL_FWD_F_1_LABEL + alpha = alpha_1 + dHat = dHat_1 + eAll = eAll_1 + mHat = mHat_1 + rms = rms_1 + f = f_1 + else + label=FWD_F_LABEL end if + ! compute the gradient and exit if (relaxation) then call linComb(ONE,mHat_0,gamma*alpha,h,mHat) @@ -859,7 +875,7 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms,logFile) end if - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=label) write(*,'(a60)') 'Sufficient decrease condition satisfied, exiting line search' write(ioLog,'(a60)') 'Sufficient decrease condition satisfied, exiting line search' call deall_dataVectorMTX(dHat_1) @@ -901,7 +917,7 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & ! end if ! compute the penalty functional call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('CUBICLS',lambda,alpha,f,mNorm,rms) call printf('CUBICLS',lambda,alpha,f,mNorm,rms,logFile) niter = niter + 1 @@ -929,12 +945,15 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & ! if the initial guess was better than what we found, take it if (starting_guess) then - alpha = alpha_1 - dHat = dHat_1 - eAll = eAll_1 - mHat = mHat_1 - rms = rms_1 - f = f_1 + alpha = alpha_1 + dHat = dHat_1 + eAll = eAll_1 + mHat = mHat_1 + rms = rms_1 + f = f_1 + label=TRIAL_FWD_F_1_LABEL + else + label=FWD_F_LABEL end if ! compute gradient of the full penalty functional and exit @@ -944,7 +963,7 @@ subroutine lineSearchCubic(lambda,d,m0,h,alpha,mHat,f,grad, & call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms) call printf('RELAX',lambda,gamma*alpha,f,mNorm,rms,logFile) end if - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=label) write(*,'(a39)') 'Gradient computed, line search finished' write(ioLog,'(a39)') 'Gradient computed, line search finished' @@ -1092,11 +1111,11 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & ! alpha_1 is the initial step size, which is set in LBFGS alpha_1 = alpha - ! compute the trial mHat, f, dHat, eAll, rms + ! compute the mHat, f, dHat, eAll, rms mHat_1 = mHat_0 ! mHat_1 = mHat_0 + dir*step call linComb(ONE,mHat_0,alpha_1,h,mHat_1) - call func(lambda,d,m0,mHat_1,f_1,mNorm_1,dHat_1,eAll_1,rms_1) + call func(lambda,d,m0,mHat_1,f_1,mNorm_1,dHat_1,eAll_1,rms_1,label=TRIAL_FWD_F_1_LABEL) call printf('STARTLS',lambda,alpha_1,f_1,mNorm_1,rms_1) call printf('STARTLS',lambda,alpha_1,f_1,mNorm_1,rms_1,logFile) niter = niter + 1 @@ -1136,7 +1155,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & end if ! otherwise compute the functional at the minimizer of the quadratic call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('QUADLS',lambda,alpha,f,mNorm,rms) call printf('QUADLS',lambda,alpha,f,mNorm,rms,logFile) niter = niter + 1 @@ -1144,7 +1163,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & ! Strong Wolfe's condition needs the gradient ! well, we are going to calculate it anyway - so why don't we do it now? if (f <= f_1) then - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=FWD_F_LABEL) g_1 = dotProd(grad, h) write(*,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0 write(*,'(a4,es12.5)') ' g1=',g_1 @@ -1174,7 +1193,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & else if (f_1 < f_0) then! is the initial making any progress? ! Test if the initial guess is good for Strong Wolfe condition - call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1) + call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1,label=TRIAL_FWD_F_1_LABEL) g_1 = dotProd(grad, h) write(*,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0 write(*,'(a4,es12.5)') ' g1=',g_1 @@ -1262,12 +1281,12 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & endif ! compute the penalty functional call linComb(ONE,mHat_0,alpha,h,mHat) - call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms) + call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms,label=FWD_F_LABEL) call printf('CUBICLS',lambda,alpha,f,mNorm,rms) call printf('CUBICLS',lambda,alpha,f,mNorm,rms,logFile) niter = niter + 1 ibracket = ibracket + 1 - call gradient(lambda,d,m0,mHat,grad,dHat,eAll) + call gradient(lambda,d,m0,mHat,grad,dHat,eAll,label=FWD_F_LABEL) g_1 = dotProd(grad, h) write(*,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0 write(*,'(a4,es12.5)') ' g1=',g_1 @@ -1306,7 +1325,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & starting_guess = .false. else starting_guess = .true. - call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1) + call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1,label=TRIAL_FWD_F_1_LABEL) endif exit endif @@ -1318,7 +1337,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, & starting_guess = .false. else starting_guess = .true. - call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1) + call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1,label=TRIAL_FWD_F_1_LABEL) endif exit endif diff --git a/f90/MPI/Declaration_MPI.f90 b/f90/MPI/Declaration_MPI.f90 index 529fa823..027df5f6 100644 --- a/f90/MPI/Declaration_MPI.f90 +++ b/f90/MPI/Declaration_MPI.f90 @@ -121,6 +121,7 @@ Module Declaration_MPI type :: define_worker_job SEQUENCE character*80 :: what_to_do='NOTHING' + character*80 :: label='NOTHING' Integer :: per_index,Stn_index,pol_index,data_type_index,data_type Integer :: taskid logical :: keep_E_soln=.false. @@ -132,25 +133,37 @@ Module Declaration_MPI type(define_worker_job), save :: worker_job_task !******************************************************************** +type :: ModEM_mpi_context_t + integer :: comm_current + integer :: comm_world + integer :: comm_local + integer :: rank_world + integer :: rank_local + integer :: rank_current + integer :: size_local + integer :: number_of_workers +end type ModEM_mpi_context_t + +type (ModEM_mpi_context_t), target, save :: modem_ctx + Contains !########################################################################## subroutine create_worker_job_task_place_holder - implicit none - integer index,Nbytes1,Nbytes2,Nbytes3 + implicit none + integer index,Nbytes1,Nbytes2,Nbytes3,NBytes4 - CALL MPI_PACK_SIZE(80, MPI_CHARACTER, MPI_COMM_WORLD, Nbytes1, ierr) - ! 2019.05.10, Liu Zhongyin, replace 6 with 7, for the iSite - ! CALL MPI_PACK_SIZE(6, MPI_INTEGER, MPI_COMM_WORLD, Nbytes2, ierr) - CALL MPI_PACK_SIZE(7, MPI_INTEGER, MPI_COMM_WORLD, Nbytes2, ierr) - CALL MPI_PACK_SIZE(3, MPI_LOGICAL, MPI_COMM_WORLD, Nbytes3, ierr) + CALL MPI_PACK_SIZE(80, MPI_CHARACTER, MPI_COMM_WORLD, Nbytes1, ierr) + CALL MPI_PACK_SIZE(80, MPI_CHARACTER, MPI_COMM_WORLD, Nbytes2, ierr) + CALL MPI_PACK_SIZE(7, MPI_INTEGER, MPI_COMM_WORLD, Nbytes3, ierr) + CALL MPI_PACK_SIZE(3, MPI_LOGICAL, MPI_COMM_WORLD, Nbytes4, ierr) - Nbytes=(Nbytes1+Nbytes2+Nbytes3)+1 + Nbytes=(Nbytes1+Nbytes2+Nbytes3+NBytes4)+1 - if(.not. associated(worker_job_package)) then - allocate(worker_job_package(Nbytes)) - end if + if(.not. associated(worker_job_package)) then + allocate(worker_job_package(Nbytes)) + end if end subroutine create_worker_job_task_place_holder !******************************************************************************* @@ -162,6 +175,7 @@ subroutine Pack_worker_job_task index=1 call MPI_Pack(worker_job_task%what_to_do,80, MPI_CHARACTER, worker_job_package, Nbytes, index, MPI_COMM_WORLD, ierr) + call MPI_Pack(worker_job_task%label,80, MPI_CHARACTER, worker_job_package, Nbytes, index, MPI_COMM_WORLD, ierr) call MPI_Pack(worker_job_task%per_index ,1 , MPI_INTEGER , worker_job_package, Nbytes, index, MPI_COMM_WORLD, ierr) call MPI_Pack(worker_job_task%Stn_index ,1 , MPI_INTEGER , worker_job_package, Nbytes, index, MPI_COMM_WORLD, ierr) call MPI_Pack(worker_job_task%pol_index ,1 , MPI_INTEGER , worker_job_package, Nbytes, index, MPI_COMM_WORLD, ierr) @@ -183,6 +197,7 @@ subroutine Unpack_worker_job_task integer index index=1 call MPI_Unpack(worker_job_package, Nbytes, index, worker_job_task%what_to_do,80, MPI_CHARACTER,MPI_COMM_WORLD, ierr) + call MPI_Unpack(worker_job_package, Nbytes, index, worker_job_task%label,80, MPI_CHARACTER,MPI_COMM_WORLD, ierr) call MPI_Unpack(worker_job_package, Nbytes, index, worker_job_task%per_index ,1 , MPI_INTEGER,MPI_COMM_WORLD, ierr) call MPI_Unpack(worker_job_package, Nbytes, index, worker_job_task%Stn_index ,1 , MPI_INTEGER,MPI_COMM_WORLD, ierr) diff --git a/f90/MPI/Main_MPI.f90 b/f90/MPI/Main_MPI.f90 index e9adf8d6..100a1376 100644 --- a/f90/MPI/Main_MPI.f90 +++ b/f90/MPI/Main_MPI.f90 @@ -18,6 +18,9 @@ module Main_MPI use Declaration_MPI use Sub_MPI + + use EsolnManager + use ModEM_timers ! use ioascii implicit none @@ -28,6 +31,9 @@ module Main_MPI type(rhsVector_t) , save, private :: b0,comb type(grid_t), target, save, private :: grid + private Master_job_DataResp_main_only + private Master_job_DataResp_IO + Contains @@ -79,6 +85,14 @@ Subroutine Constructor_MPI size_gpu = 0 ! tic-toc previous_time = MPI_Wtime() + + modem_ctx % comm_world = comm_world + modem_ctx % comm_local = comm_local + modem_ctx % rank_world = rank_world + modem_ctx % rank_local = rank_local + modem_ctx % size_local = size_local + modem_ctx % number_of_workers = number_of_workers + End Subroutine Constructor_MPI !---------------------------------------------------------------------------- @@ -606,13 +620,14 @@ end subroutine Master_job_Regroup !---------------------------------------------------------------------------- !########################## Master_job_fwdPred ############################ -Subroutine Master_job_fwdPred(sigma,d1,eAll,comm) +Subroutine Master_job_fwdPred(sigma,d1,eAll,comm,label) implicit none - type(modelParam_t), intent(in) :: sigma - type(dataVectorMTX_t), intent(inout) :: d1 - type(solnVectorMTX_t), intent(inout) :: eAll - integer, intent(inout),optional :: comm + type(modelParam_t), intent(in) :: sigma + type(dataVectorMTX_t), intent(inout) :: d1 + type(solnVectorMTX_t), intent(inout) :: eAll + integer, intent(inout), optional :: comm + character(len=*), intent(in), optional :: label integer nTx @@ -620,6 +635,13 @@ Subroutine Master_job_fwdPred(sigma,d1,eAll,comm) Integer :: iper, comm_current Integer :: per_index,pol_index,stn_index,iTx,i,iDt,j character(80) :: job_name + character(len=80) :: label_lcl + + if (present(label)) then + label_lcl = label + else + label_lcl = 'NULL' + end if ! nTX is number of transmitters; nTx = d1%nTx @@ -639,43 +661,22 @@ Subroutine Master_job_fwdPred(sigma,d1,eAll,comm) end if end if + modem_ctx % comm_current = comm_current + ! First, distribute the current model to all workers call Master_job_Distribute_Model(sigma) + ! call Master_job_Distribute_Data(d1) if(.not. eAll%allocated) then - ! call deall(eAll) - ! end if - call create_solnVectorMTX(d1%nTx,eAll) - do iTx=1,nTx - call create_solnVector(grid,iTx,e0) - call copy_solnVector(eAll%solns(iTx),e0) - end do + call EsMgr_create_solnVectorMTX(eAll, d1 % nTx, grid=grid) end if - job_name= 'FORWARD' - call Master_job_Distribute_Taskes(job_name,nTx,sigma,eAll,comm_current) - ! Initialize only those grid elements on the master that are used in - ! EMfieldInterp - ! (obviously a quick patch, needs to be fixed in a major way) - ! A.Kelbert 2018-01-28 + job_name= 'FORWARD' + call Master_job_Distribute_Taskes(job_name,nTx,sigma,eAll,comm_current, label=label_lcl) - Call EdgeLength(grid, l_E) - Call FaceArea(grid, S_F) + ! Calculate Data Response + call Master_job_DataResp(nTx, sigma, eAll, d1, label=label_lcl) - ! Compute the model Responces - do iTx=1,nTx - do i = 1,d1%d(iTx)%nDt - d1%d(iTx)%data(i)%errorBar = .false. - iDt = d1%d(iTx)%data(i)%dataType - do j = 1,d1%d(iTx)%data(i)%nSite - call dataResp(eAll%solns(iTx),sigma,iDt,d1%d(iTx)%data(i)%rx(j),d1%d(iTx)%data(i)%value(:,j), & - d1%d(iTx)%data(i)%orient(j)) - end do - end do - end do - ! clean up the grid elements stored in GridCalc on the master node - call deall_rvector(l_E) - call deall_rvector(S_F) write(ioMPI,*)'FWD: Finished calculating for (', nTx , ') Transmitters ' endtime=MPI_Wtime() time_used = endtime-starttime @@ -684,6 +685,180 @@ Subroutine Master_job_fwdPred(sigma,d1,eAll,comm) end subroutine Master_job_fwdPred +!---------------------------------------------------------------------------- +!########################## Master_job_DataResp ############################ +! Compute the data response by either using Master_job_DataResp_main_only +! or Master_job_DataResp_IO: +! +! * Master_job_DataResp_main_only - Runs if EsMgr_save_in_file == .false. +! * Master_job_DataResp_IO - Runs if EsMgr_save_in_file == .true. +! +subroutine Master_job_DataResp(nTx, sigma, eAll, d1, label) + + implicit none + + integer, intent(in) :: nTx + type(modelParam_t), intent(in) :: sigma + type(solnVectorMTX_t), intent(in) :: eAll + type(dataVectorMTX_t), intent(inout) :: d1 + character(len=*), intent(in), optional :: label + + integer :: iTx, i, j + integer :: iDt + + character(len=80) :: label_lcl + + if (present(label)) then + label_lcl = label + else + label_lcl = "NULL" + end if + + call ModEM_timers_create('DataResp') + call ModEM_timers_start('DataResp') + + if (EsMgr_save_in_file) then + ! Use the worker jobs to calculate the data resposne by reading + ! the electric field files + call Master_job_DataResp_IO(nTx, sigma, d1, label=trim(label_lcl)) + else + ! The main task only computes the data response + call Master_job_DataResp_main_only(nTx, sigma, d1, eAll) + end if + + call ModEM_timers_stop('DataResp') + +end subroutine Master_job_DataResp + +! Master_job_DataResp_main_only - Calculate the data response on +! the main task only. The main task will need to have eAll completly +! present. +! +! Run this function if EsMgr_save_in_file == .false. +subroutine Master_job_DataResp_main_only(nTx, sigma, d1, eAll) + + implicit none + + integer, intent(in) :: nTx + type (modelParam_t), intent(in) :: sigma + type (dataVectorMTX_t), intent(inout) :: d1 + type(solnVectorMTX_t), intent(in) :: eAll + + integer :: iTx, i, j + integer :: iDt + logical :: trial_lcl + + ! Initialize only those grid elements on the master that are used in + ! EMfieldInterp + ! (obviously a quick patch, needs to be fixed in a major way) + ! A.Kelbert 2018-01-28 + call EdgeLength(grid, l_E) + call FaceArea(grid, S_F) + + do iTx=1,nTx + do i = 1,d1%d(iTx)%nDt + d1%d(iTx)%data(i)%errorBar = .false. + iDt = d1%d(iTx)%data(i)%dataType + do j = 1,d1%d(iTx)%data(i)%nSite + call dataResp(eAll%solns(iTx),sigma,iDt,d1%d(iTx)%data(i)%rx(j),d1%d(iTx)%data(i)%value(:,j), & + d1%d(iTx)%data(i)%orient(j)) + end do + end do + end do + + ! clean up the grid elements stored in GridCalc on the master node + call deall_rvector(l_E) + call deall_rvector(S_F) + +end subroutine Master_job_DataResp_main_only + +! Master_job_DataResp_IO +! +! Use the worker jobs to calculate the data response by +! reading in the electric fields from files created by the ESolnManager. +! +! This function should only be used if EsMgr_init was initalized with +! `save_in_file = .true.`/EsMgr_save_in_file = .true.. +subroutine Master_job_DataResp_IO(nTx, sigma, d, label) + + implicit none + + integer, intent(in) :: nTx + type (modelParam_t), intent(in) :: sigma + type (dataVectorMTX_t), intent(inout) :: d + character(len=*), intent(in), optional :: label + + character (len=*), parameter :: JOB_NAME = "DATARESP" + + integer :: dest, nTasks, remainder, iTx + integer :: iTx_min, iTx_max, i, j, k + integer :: comm_current + integer :: size_current + character(len=80) :: label_lcl + + if (present(label)) then + label_lcl = label + else + label_lcl = trim('NULL') + endif + + if (para_method.eq.0) then + comm_current = comm_world + else + comm_current = comm_leader + end if + + modem_ctx % comm_current = comm_current + call MPI_COMM_SIZE( comm_current, size_current, ierr ) + + call create_worker_job_task_place_holder + + nTasks = nTx / (size_leader - 1) + remainder = modulo(nTx, size_leader - 1) + iTx_max = 0 + + do dest = 1, size_current - 1 + iTx_min = iTx_max + 1 + iTx_max = iTx_min + nTasks - 1 + + if (remainder > 0) then + iTx_max = iTx_max + 1 + remainder = remainder - 1 + end if + + if ( iTx_max >= iTx_min ) then + worker_job_task % what_to_do = trim(job_name) + worker_job_task % per_index = iTx_min + worker_job_task % pol_index = iTx_max + worker_job_task % label = trim(label_lcl) + + call Pack_worker_job_task + call MPI_Send(worker_job_package, Nbytes, MPI_PACKED, dest, FROM_MASTER, comm_current, ierr) + write(ioMPI, '(a10,a16,i5,a8,i5,a11,i5)') trim(job_name), ': Send Per from ', iTx_min, ' to', iTx_max, ' to ', dest + end if + end do + + remainder = modulo(nTx, size_current - 1) + iTx_max = 0 + + do dest = 1, size_current - 1 + iTx_min = iTx_max + 1 + iTx_max = iTx_min + nTasks - 1 + + if (remainder > 0) then + iTx_max = iTx_max + 1 + remainder = remainder - 1 + end if + + if (iTx_max >= iTx_min) then + call create_data_vec_place_holder(d, start_iTx=iTx_min, end_iTx=iTx_max) + call MPI_Recv(data_para_vec, Nbytes, MPI_PACKED, dest, FROM_WORKER, comm_current, MPI_STATUS_IGNORE, ierr) + call UnPack_data_para_vec(d, start_iTx=iTx_min, end_iTx=iTx_max) + end if + end do + +end subroutine Master_job_DataResp_IO + !######################### Master_job_Compute_J ########################## @@ -934,7 +1109,7 @@ end subroutine Master_job_calcJ !######################### Master_job_JmultT ############################ -Subroutine Master_job_JmultT(sigma,d,dsigma,eAll,s_hat,comm) +Subroutine Master_job_JmultT(sigma,d,dsigma,eAll,s_hat,comm,label) implicit none type(modelParam_t), intent(in) :: sigma @@ -942,7 +1117,8 @@ Subroutine Master_job_JmultT(sigma,d,dsigma,eAll,s_hat,comm) type(modelParam_t), intent(Out) :: dsigma type(solnVectorMTX_t), intent(in), optional :: eAll type(modelParam_t),intent(inout),pointer,dimension(:), optional :: s_hat - integer, intent(in),optional :: comm + integer, intent(in), optional :: comm + character(len=*), intent(in), optional :: label ! Local type(modelParam_t) :: dsigma_temp @@ -956,6 +1132,14 @@ Subroutine Master_job_JmultT(sigma,d,dsigma,eAll,s_hat,comm) Integer :: per_index,pol_index,stn_index character(80) :: job_name,file_name Integer :: comm_current + character(80) :: label_lcl + + + if (present(label)) then + label_lcl = label + else + label_lcl = 'NULL' + end if savedSolns = present(eAll) returne_m_vectors= present(s_hat) @@ -973,6 +1157,9 @@ Subroutine Master_job_JmultT(sigma,d,dsigma,eAll,s_hat,comm) comm_current = comm_leader endif end if + + modem_ctx % comm_current = comm_current + ! nTX is number of transmitters; nTx = d%nTx if(.not. eAll_temp%allocated) then @@ -993,7 +1180,6 @@ Subroutine Master_job_JmultT(sigma,d,dsigma,eAll,s_hat,comm) end do end if - if (returne_m_vectors) then if (.not. associated(s_hat)) then allocate(s_hat(nTx)) @@ -1016,25 +1202,15 @@ Subroutine Master_job_JmultT(sigma,d,dsigma,eAll,s_hat,comm) call zero(Qcomb) job_name= 'JmultT' call Master_job_Distribute_Taskes(job_name,nTx,sigma,eAll_out, & - & comm_current, eAll_temp) + & comm_current, eAll_temp, label=label_lcl) file_name='e0.soln' !call write_solnVectorMTX(10,file_name,eAll_temp) file_name='e.soln' !call write_solnVectorMTX(20,file_name,eAll_out) - do iper=1,nTx - !e0=eAll%solns(iper) - !e =eAll_out%solns(iper) - call PmultT(eAll_temp%solns(iper),sigma,eAll_out%solns(iper) & - & ,dsigma_temp) - call QmultT(eAll_temp%solns(iper),sigma,d%d(iper),Qcomb) - call scMultAdd(ONE,Qcomb,dsigma_temp) - if (returne_m_vectors) then - s_hat(iper)=dsigma_temp - end if - call linComb_modelParam(ONE,dsigma,ONE,dsigma_temp,dsigma) - end do + ! Calculate PQMult + call Master_job_PQMult(nTx, d, eAll_temp, eAll_out, sigma, dsigma, s_hat, label=label_lcl) endtime=MPI_Wtime() time_used = endtime-starttime @@ -1055,6 +1231,211 @@ Subroutine Master_job_JmultT(sigma,d,dsigma,eAll,s_hat,comm) end Subroutine Master_job_JmultT +!########################## Master_job_PQMult ############################ +! Compute the model response by either using Master_job_PQMult_main_only +! or Master_job_PQMult_IO: +! +! * Master_job_PQMult_main_only - Runs if EsMgr_save_in_file == .false. +! * Master_job_PQMult_IO - Runs if EsMgr_save_in_file == .true. +! +subroutine Master_job_PQMult(nTx, d, eAll_temp, eAll_out, sigma, dsigma, s_hat, label) + + implicit none + + integer, intent(in) :: nTx + type(dataVectorMTX_t), intent(in) :: d + type(solnVectorMTX_t), intent(in) :: eAll_out + type(solnVectorMTX_t), intent(in) :: eAll_temp + type (modelParam_t), intent(in) :: sigma + type (modelParam_t), intent(inout) :: dsigma + type(modelParam_t),intent(inout),pointer,dimension(:), optional :: s_hat + character(len=*), intent(in), optional :: label + character(len=80) :: label_lcl + + if (present(label)) then + label_lcl = label + else + label_lcl = 'NULL' + end if + + call ModEM_timers_create('PQMult') + call ModEM_timers_start('PQMult') + + if (EsMgr_save_in_file) then + call Master_job_PQMult_IO(nTx, sigma, dsigma, label=label_lcl) + else + call Master_job_PQMult_main_only(nTx, d, eAll_temp, eAll_out, sigma, dsigma, s_hat) + end if + + call ModEM_timers_stop('PQMult') + +end subroutine Master_job_PQMult + + +! Master_job_PQMult_main_only - Calculate the model response on +! the main task only. The main task will need to have eAll_temp set to eAll. +! +! Run this function if ESmgr_save_in_file == f.alse. +subroutine Master_job_PQMult_main_only(nTx, d, eAll_temp, eAll_out, sigma, dsigma, s_hat) + + implicit none + + integer, intent(in) :: nTx + type(dataVectorMTX_t), intent(in) :: d + type(solnVectorMTX_t), intent(in) :: eAll_out + type(solnVectorMTX_t), intent(in) :: eAll_temp + type (modelParam_t), intent(in) :: sigma + type (modelParam_t), intent(inout) :: dsigma + type(modelParam_t),intent(inout),pointer,dimension(:), optional :: s_hat + + type(modelParam_t) :: dsigma_temp + type(modelParam_t) :: Qcomb + type(dataVectorMTX_t) :: d_temp + logical :: returne_m_vectors + integer :: iper + + returne_m_vectors = present(s_hat) + + dsigma_temp = sigma + dsigma = sigma + Qcomb = sigma + + call zero(dsigma_temp) + call zero(dsigma) + call zero(Qcomb) + + do iper=1,nTx + call PmultT(eAll_temp%solns(iper),sigma,eAll_out%solns(iper), dsigma_temp) + call QmultT(eAll_temp%solns(iper),sigma,d%d(iper),Qcomb) + call scMultAdd(ONE,Qcomb,dsigma_temp) + if (returne_m_vectors) then + s_hat(iper)=dsigma_temp + end if + call linComb_modelParam(ONE,dsigma,ONE,dsigma_temp,dsigma) + end do + + call deall_modelParam(dsigma_temp) + call deall_modelParam(Qcomb) + call deall_dataVectorMTX(d_temp) + +end subroutine Master_job_PQMult_main_only + +! Master_job_PQMult_IO +! +! Use the worker jobs to calculate the model response by +! reading in the electric fields from files created by the ESolnManager and +! the adjoint files. +! +! This function should only be used if EsMgr_init was initalized with +! `save_in_file = .true.`/EsMgr_save_in_file = .true.. +subroutine Master_job_PQMult_IO(nTx, sigma, dsigma, label) + + implicit none + + integer, intent(in) :: nTx + type (modelParam_t), intent(in) :: sigma + type (modelParam_t), intent(inout) :: dsigma + character(len=*), intent(in), optional :: label + + character(len=*), parameter :: job_name = "PQMULT" + type(modelParam_t), dimension(nTx) :: dsigma_recv + character(len=80) :: label_lcl + + integer :: dest, nTasks, remainder, iTx + integer :: iTx_min, iTx_max, i, j, k + logical :: flag + integer :: comm_current, size_current + + logical, dimension(number_of_workers) :: task_is_working + logical, dimension(nTx) :: transmitters_processing, transmitters_done + logical :: sending + + if (present(label)) then + label_lcl = label + else + label_lcl = 'NULL' + end if + + if (para_method.eq.0) then + comm_current = comm_world + else + comm_current = comm_leader + end if + + modem_ctx % comm_current = comm_current + call MPI_COMM_SIZE( comm_current, size_current, ierr ) + + sending = .true. + + task_is_working(:) = .false. + transmitters_processing(:) = .false. + transmitters_done(:) = .false. + + call zero(dsigma) + + do while(sending) + ! Send out jobs + do iTx = 1, nTx + + if (transmitters_processing(iTx) .or. transmitters_done(iTx)) then + cycle + end if + + do dest = 1, size_current - 1 + if (task_is_working(dest)) then + cycle + end if + + ! Send the task the job + worker_job_task % what_to_do=trim(job_name) + worker_job_task % per_index = iTx + worker_job_task % pol_index = -1 + worker_job_task % label = trim(label_lcl) + + call create_worker_job_task_place_holder + call Pack_worker_job_task + call MPI_Send(worker_job_package, Nbytes, MPI_PACKED, dest, FROM_MASTER, comm_current, ierr) + transmitters_processing(iTx) = .true. + task_is_working(dest) = .true. + exit + end do + end do + + ! Recv any jobs + ! See if anyone is sending us a message... + call MPI_Iprobe(MPI_ANY_SOURCE, FROM_WORKER, comm_current, flag, MPI_STATUS_IGNORE, ierr) + if (flag) then ! Someone is sending us a message + call create_worker_job_task_place_holder + call MPI_Recv(worker_job_package, Nbytes, MPI_PACKED, MPI_ANY_SOURCE, FROM_WORKER, comm_current, STATUS, ierr) + call Unpack_worker_job_task + + dest = worker_job_task % taskid + iTx = worker_job_task % per_index + + dsigma_recv(iTx) = sigma + call zero(dsigma_recv(iTx)) + + call create_model_param_place_holder(dsigma_recv(iTx)) + call MPI_Recv(sigma_para_vec, Nbytes, MPI_PACKED, dest, FROM_WORKER, comm_current, MPI_STATUS_IGNORE, ierr) + call Unpack_model_para_values(dsigma_recv(iTx)) + + transmitters_done(iTx) = .true. + task_is_working(dest) = .false. + end if + + if (all(transmitters_done(:))) then + sending = .false. + exit + end if + end do + + do iTx = 1, nTx + call linComb_modelParam(ONE, dsigma, ONE, dsigma_recv(iTx), dsigma) + call deall_modelParam(dsigma_recv(iTx)) + end do + +end subroutine Master_job_PQMult_IO + !######################## Master_job_Jmult ############################## Subroutine Master_job_Jmult(mHat,m,d,eAll,comm) @@ -1535,7 +1916,7 @@ end Subroutine Master_job_Stop_MESSAGE !********************** Master Distribute Tasks ***************************** subroutine Master_job_Distribute_Taskes(job_name,nTx,sigma,eAll_out, & - & comm, eAll_in) + & comm, eAll_in, label) implicit none character(80) , intent(in) :: job_name @@ -1544,13 +1925,15 @@ subroutine Master_job_Distribute_Taskes(job_name,nTx,sigma,eAll_out, & type(solnVectorMTX_t), intent(in), optional :: eAll_in integer, intent(in), optional :: comm type(solnVectorMTX_t), intent(inout), optional :: eAll_out + character(len=*), intent(in), optional :: label !Local Integer :: iper,ipol,ipol1,ijob,total_jobs Integer :: bcounter, tcounter, robin Integer :: per_index,pol_index,des_index - logical :: keep_soln,savedSolns,ascend + logical :: keep_soln,savedSolns,ascend,trial_lcl integer :: comm_current, size_current double precision :: now, just, time_passed + character(len=80) :: label_lcl savedSolns = present(eAll_in) ! over-ride the default communicator, if needed @@ -1563,6 +1946,14 @@ subroutine Master_job_Distribute_Taskes(job_name,nTx,sigma,eAll_out, & else comm_current = comm_world end if + + if (present(label)) then + label_lcl = label + else + label_lcl = 'NULL' + end if + + call get_nPol_MPI(eAll_out%solns(1)) if (rank_local.eq.-1) then ! first run! ! run initial regroup -- note this requires the comm to be @@ -1576,6 +1967,7 @@ subroutine Master_job_Distribute_Taskes(job_name,nTx,sigma,eAll_out, & endif call MPI_COMM_SIZE( comm_current, size_current, ierr ) who = 0 + worker_job_task%what_to_do=trim(job_name) call count_number_of_messages_to_RECV(eAll_out) total_jobs = answers_to_receive @@ -1609,6 +2001,8 @@ subroutine Master_job_Distribute_Taskes(job_name,nTx,sigma,eAll_out, & end if worker_job_task%per_index = per_index worker_job_task%pol_index = pol_index + worker_job_task%label=trim(label_lcl) + call create_worker_job_task_place_holder call Pack_worker_job_task call MPI_SEND(worker_job_package,Nbytes, MPI_PACKED, who, & @@ -1623,10 +2017,7 @@ subroutine Master_job_Distribute_Taskes(job_name,nTx,sigma,eAll_out, & which_per=per_index do ipol1=1,nPol_MPI which_pol=ipol1 - call create_e_param_place_holder(eAll_in%solns(which_per)) - call Pack_e_para_vec(eAll_in%solns(which_per)) - call MPI_SEND(e_para_vec, Nbytes, MPI_PACKED, who, & - & FROM_MASTER, comm_current, ierr) + call EsMgr_save(eAll_in % solns(which_per), to=who) end do end if write(ioMPI,'(a10,a16,i5,a8,i5,a11,i5)')trim(job_name), & @@ -1652,14 +2043,9 @@ subroutine Master_job_Distribute_Taskes(job_name,nTx,sigma,eAll_out, & who=worker_job_task%taskid which_per=worker_job_task%per_index which_pol=worker_job_task%pol_index - - call create_e_param_place_holder(eAll_out%solns(which_per)) - call MPI_RECV(e_para_vec, Nbytes, MPI_PACKED, who,FROM_WORKER, & - & comm_current, STATUS, ierr) - ! call get_nPol_MPI(eAll_out%solns(which_per)) - ! if (nPol_MPI==1) which_pol=1 - call Unpack_e_para_vec(eAll_out%solns(which_per)) + modem_ctx % comm_current = comm_current + call EsMgr_get(eAll_out % solns(which_per), which_pol, 1, from=who) write(ioMPI,'(a10,a16,i5,a8,i5,a11,i5)')trim(job_name) , & & ': Receive Per # ',which_per ,' and Pol # ', which_pol ,' from ',& @@ -1703,7 +2089,8 @@ subroutine Master_job_Distribute_Taskes(job_name,nTx,sigma,eAll_out, & worker_job_task%per_index= per_index worker_job_task%pol_index= pol_index - worker_job_task%what_to_do= trim(job_name) + worker_job_task%what_to_do= trim(job_name) + worker_job_task%label = trim(label_lcl) call create_worker_job_task_place_holder call Pack_worker_job_task @@ -1716,11 +2103,7 @@ subroutine Master_job_Distribute_Taskes(job_name,nTx,sigma,eAll_out, & call get_nPol_MPI(eAll_out%solns(per_index)) do ipol1=1,nPol_MPI which_pol=ipol1 - call create_e_param_place_holder(eAll_in%solns( & - & which_per)) - call Pack_e_para_vec(eAll_in%solns(which_per)) - call MPI_SEND(e_para_vec, Nbytes, MPI_PACKED, & - & who, FROM_MASTER, comm_current, ierr) + call EsMgr_save(eAll_in % solns(which_per), who) end do end if write(ioMPI,'(a10,a16,i5,a8,i5,a11,i5)')trim(job_name), & @@ -1748,6 +2131,7 @@ subroutine Master_job_Distribute_Taskes(job_name,nTx,sigma,eAll_out, & end subroutine Master_job_Distribute_Taskes + !################## find next job -- from back or front ################## Subroutine find_next_job(nTx,total_jobs,counter,fromhead, eAll_out, & & per_index,pol_index) @@ -1822,8 +2206,9 @@ Subroutine Worker_job(sigma,d) ! Local type(modelParam_t) :: delSigma type(userdef_control) :: ctrl - Integer :: nTx,m_dimension,ndata - Integer :: itx, ndt, dt, dt_index + Integer :: nTx,m_dimension,ndata, start_iTx, end_iTx + Integer :: itx, ndt, dt, dt_index, iDt + Integer :: j Integer :: iper,ipol,i,des_index Integer :: per_index,per_index_pre @@ -1847,7 +2232,9 @@ Subroutine Worker_job(sigma,d) ! 2019.05.08, Liu Zhongyin, add isite for rx in dataBlock_t integer :: isite - + ! PQMult + type (modelParam_t) :: dsigma_temp, dsigma_send, Qcomb + ! time DOUBLE PRECISION :: time_passed, now, just DOUBLE PRECISION, pointer,dimension(:) :: time_buff @@ -1905,6 +2292,10 @@ Subroutine Worker_job(sigma,d) write(6,'(a12,a12,a30,a16,i5)') node_info,' MPI TASK [', & & trim(worker_job_task%what_to_do),'] received from ', & & STATUS(MPI_SOURCE) + + modem_ctx % comm_current = comm_current + modem_ctx % rank_current = rank_current + ! for debug ! write(6,*) 'source = ', MPI_SOURCE ! write(6,*) 'tag = ', MPI_TAG @@ -1914,9 +2305,14 @@ Subroutine Worker_job(sigma,d) ! '; several TX = ',worker_job_task%several_Tx,']' if (trim(worker_job_task%what_to_do) .eq. 'FORWARD') then + ! forward modelling per_index=worker_job_task%per_index pol_index=worker_job_task%pol_index + + modem_ctx % comm_current = comm_current + modem_ctx % rank_current = rank_current + if ((size_local.gt.1).and.(para_method.gt.0).and. & & (rank_local.eq.0)) then ! group leader passing the command to workers @@ -1964,13 +2360,10 @@ Subroutine Worker_job(sigma,d) call Pack_worker_job_task call MPI_SEND(worker_job_package,Nbytes, MPI_PACKED,0, & & FROM_WORKER, comm_current, ierr) - ! Create e0_temp package (one Period and one Polarization) - ! and send it to the master + ! Use EsMgr_save to send e0 back to main task or save it to disk which_pol=1 - call create_e_param_place_holder(e0) - call Pack_e_para_vec(e0) - call MPI_SEND(e_para_vec, Nbytes, MPI_PACKED, 0, & - & FROM_WORKER, comm_current, ierr) + + call EsMgr_save(e0, to=0, E_field_type='fwd', label=worker_job_task % label) end if ! so long! call reset_e_soln(e0) @@ -1978,6 +2371,93 @@ Subroutine Worker_job(sigma,d) time_passed = now - previous_time previous_time = now + elseif (trim(worker_job_task%what_to_do) .eq. 'DATARESP') then + + start_iTx = worker_job_task % per_index + end_iTx = worker_job_task % pol_index + worker_job_task % taskid = taskid + + call zero_solnvector(e0) + call zero_solnvector(e) + + call EsMgr_create_e(e, per_index) + + do per_index = start_iTx, end_iTx + call zero_solnvector(e0) + e0 % tx = per_index + e % tx = per_index + + do pol_index = 1, get_nPol(per_index) + call EsMgr_get(e0, e0 % tx, pol_index=pol_index, E_field_type='fwd', label=worker_job_task % label) + end do + + do i = 1, d % d(per_index) % nDt + d % d(per_index) % data(i) % errorBar = .false. + iDt = d % d(per_index) % data(i) % dataType + + do j = 1, d % d(per_index) % data(i) % nSite + call dataResp(e0, sigma, iDt, d % d(per_index) % data(i) % rx(j), & + d % d(per_index) % data(i) % value(:,j), & + d % d(per_index) % data(i) % orient(j)) + end do + end do + end do + + call create_data_vec_place_holder(d, start_iTx=start_iTx, end_iTx=end_iTx) + call Pack_data_para_vec(d, start_iTx=start_iTx, end_iTx=end_iTx) + call MPI_Send(data_para_vec, NBytes, MPI_PACKED, 0, FROM_WORKER, comm_current, ierr) + deallocate(data_para_vec) + data_para_vec => null() + + elseif (trim(worker_job_task%what_to_do) .eq. 'PQMULT') then + + per_index = worker_job_task % per_index + worker_job_task % taskid = taskid + + + call create_solnVector(grid, 1, e0) + call create_solnVector(grid, 1, e) + + call zero_solnvector(e0) + call zero_solnvector(e) + + dsigma_temp = sigma + Qcomb = sigma + dsigma_send = sigma + + call zero(dsigma_temp) + call zero(dsigma_send) + call zero(QComb) + + e0 % tx = per_index + e % tx = per_index + + do pol_index = 1, get_nPol(per_index) + call EsMgr_get(e0, e0 % tx, pol_index=pol_index, E_field_type='fwd', label=trim(worker_job_task % label)) + call EsMgr_get(e, e % tx, pol_index=pol_index, E_field_type='JmultT', label=trim(worker_job_task % label)) + end do + + call PmultT(e0, sigma, e, dsigma_temp) + call QmultT(e0, sigma, d%d(per_index), Qcomb) + call scMultAdd(ONE, Qcomb, dsigma_temp) + + ! Send to Main that we are done + call create_worker_job_task_place_holder + worker_job_task % what_to_do = 'PQMULT-RESP' + worker_job_task % taskid = rank_current + worker_job_task % per_index = per_index + call Pack_worker_job_task + call MPI_Send(worker_job_package, Nbytes, MPI_PACKED, 0, & + FROM_WORKER, comm_current, ierr) + + call create_model_param_place_holder(dsigma_temp) + call pack_model_para_values(dsigma_temp) + call MPI_Send(sigma_para_vec, NBytes, MPI_PACKED, 0, FROM_WORKER, comm_current, ierr) + call deall(QComb) + call deall(dsigma_temp) + call deall(dsigma_send) + + elseif (trim(worker_job_task%what_to_do) .eq. 'COMPUTE_J') then ! compute (explicit) J if ((size_local.gt.1).and.(para_method.gt.0).and. & @@ -2151,6 +2631,7 @@ Subroutine Worker_job(sigma,d) ! calculate JmultT (a.k.a. adjoint) per_index=worker_job_task%per_index pol_index=worker_job_task%pol_index + if ((size_local.gt.1).and.(para_method.gt.0).and. & & (rank_local.eq.0)) then ! group leader passing the command to workers @@ -2170,10 +2651,9 @@ Subroutine Worker_job(sigma,d) & ' Start Receiving ' , orginal_nPol, ' from Master' do ipol=1,nPol_MPI which_pol=ipol - call create_e_param_place_holder(e0) - call MPI_RECV(e_para_vec, Nbytes, MPI_PACKED, 0, & - & FROM_MASTER,comm_current, STATUS, ierr) - call Unpack_e_para_vec(e0) + call EsMgr_get(e0, per_index, pol_index=ipol, from=0, & + E_field_type='fwd', & + label=worker_job_task % label) end do call initSolverWithOutE0(per_index,sigma,grid,size_local,& e,comb) @@ -2217,11 +2697,7 @@ Subroutine Worker_job(sigma,d) call MPI_SEND(worker_job_package,Nbytes, MPI_PACKED,0, & & FROM_WORKER, comm_current, ierr) which_pol=1 - call create_e_param_place_holder(e) - call Pack_e_para_vec(e) - call MPI_SEND(e_para_vec, Nbytes, MPI_PACKED, 0, & - & FROM_WORKER, comm_current, ierr) - !deallocate(e_para_vec,worker_job_package) + call EsMgr_save(e, to=0, E_field_type="JmultT", label=worker_job_task % label) end if ! hasta la vista! now = MPI_Wtime() @@ -2453,6 +2929,7 @@ Subroutine Worker_job(sigma,d) call Pack_worker_job_task call MPI_SEND(worker_job_package,Nbytes, MPI_PACKED,0, & & FROM_WORKER, comm_current, ierr) + elseif (trim(worker_job_task%what_to_do) .eq. 'REGROUP') then ! calculate the time between two regroup events if ((size_local.gt.1).and.(para_method.gt.0).and. & @@ -2626,14 +3103,39 @@ End Subroutine Worker_job !****************************************************************************** -subroutine create_data_vec_place_holder(d) +subroutine create_data_vec_place_holder(d, start_iTx, end_iTx) implicit none integer Nbytes1,Nbytes2,ndata,iper,ndt,sum1,sum2 type(dataVectorMTX_t), intent(in) :: d + integer, optional, intent(in) :: start_iTx + integer, optional, intent(in) :: end_iTx + integer :: start_iTx_lcl + integer :: end_iTx_lcl + + ! .neqv. here is the exclusive Or + if (present(start_iTx) .neqv. present(end_iTx)) then + write(0,*) "ERROR: In create_data_vec_place_holder - Either 'start_iTx' or 'end_iTx' was specified, but not the other" + write(0,*) "ERROR: Please either: Specify BOTH, or, specify NONE" + call ModEM_Abort() + end if + + if (present(start_iTx)) then + start_iTx_lcl = start_iTx + else + start_iTx_lcl = 1 + end if + + if (present(end_iTx)) then + end_iTx_lcl = end_iTx + else + end_iTx_lcl = d % nTx + end if + sum1=0 sum2=0 - do iper=1,d%nTx + + do iper = start_iTx_lcl, end_iTX_lcl do ndt=1,d%d(iper)%ndt ndata=size(d%d(iper)%data(ndt)%value) CALL MPI_PACK_SIZE(ndata, MPI_DOUBLE_PRECISION, & @@ -2651,14 +3153,30 @@ subroutine create_data_vec_place_holder(d) end subroutine create_data_vec_place_holder !****************************************************************************** - subroutine Pack_data_para_vec(d) +subroutine Pack_data_para_vec(d, start_iTx, end_iTx) implicit none type(dataVectorMTX_t), intent(in) :: d + integer, optional, intent(in) :: start_iTx + integer, optional, intent(in) :: end_iTx + integer :: start_iTx_lcl, end_iTx_lcl integer index integer ndata,iper,ndt + + if (present(start_iTx)) then + start_iTx_lcl = start_iTx + else + start_iTx_lcl = 1 + end if + + if (present(end_iTx)) then + end_iTx_lcl = end_iTx + else + end_iTx_lcl = d % nTx + end if + index=1 - do iper=1,d%nTx + do iper=start_iTx_lcl, end_iTx_lcl do ndt=1,d%d(iper)%ndt ndata=size(d%d(iper)%data(ndt)%value) call MPI_Pack(d%d(iper)%data(ndt)%value(1,1),ndata, & @@ -2676,15 +3194,37 @@ subroutine Pack_data_para_vec(d) end subroutine Pack_data_para_vec !****************************************************************************** -subroutine UnPack_data_para_vec(d) +subroutine UnPack_data_para_vec(d, start_iTx, end_iTx) implicit none type(dataVectorMTX_t), intent(inout) :: d + integer, optional, intent(in) :: start_iTx + integer, optional, intent(in) :: end_iTx ! Local - integer index - integer ndata,iper,ndt + integer :: start_iTx_lcl, end_iTx_lcl + integer :: index + integer :: ndata,iper,ndt + + ! .neqv. here is the exclusive Or + if (present(start_iTx) .neqv. present(end_iTx)) then + write(0,*) "ERROR: In UnPack_Data_para_vec - Either 'start_iTx' or 'end_iTx' was specified, but not the other" + write(0,*) "ERROR: Please either: Specify BOTH, or, specify NONE" + call ModEM_Abort() + end if + + if (present(start_iTx)) then + start_iTx_lcl = start_iTx + else + start_iTx_lcl = 1 + end if + + if (present(end_iTx)) then + end_iTx_lcl = end_iTx + else + end_iTx_lcl = d % nTx + end if index=1 - do iper=1,d%nTx + do iper=start_iTx_lcl, end_iTx_lcl do ndt=1,d%d(iper)%ndt ndata=size(d%d(iper)%data(ndt)%value) call MPI_Unpack(data_para_vec, Nbytes, index, & @@ -2703,7 +3243,8 @@ subroutine UnPack_data_para_vec(d) end do if (associated(data_para_vec)) then - deallocate (data_para_vec) + deallocate(data_para_vec) + data_para_vec => null() end if end subroutine UnPack_data_para_vec diff --git a/f90/Mod3DMT.f90 b/f90/Mod3DMT.f90 index 1c46b911..fa6af649 100755 --- a/f90/Mod3DMT.f90 +++ b/f90/Mod3DMT.f90 @@ -62,7 +62,13 @@ program Mod3DMT call initGlobalData(cUserDef) ! set the grid for the numerical computations #ifdef MPI - call setGrid_MPI(grid) + call setGrid_MPI(grid) + + call EsMgr_init(grid, context=modem_ctx, & + save_in_file=cUserDef % storeSolnsInFile, & + prefix=cUserDef % prefix, & + ftype=FTYPE_BINARY) + ! Check if a large grid file with E field is defined: ! NOTE: right now both grids share the same transmitters. ! This why, reading and setting the large grid and its E solution comes after setting the trasnmitters Dictionary. diff --git a/f90/UserCtrl.f90 b/f90/UserCtrl.f90 index 7eeeff67..699b70a0 100644 --- a/f90/UserCtrl.f90 +++ b/f90/UserCtrl.f90 @@ -60,8 +60,8 @@ module UserCtrl ! Choose the sort of test / procedure variant you wish to perform character(80) :: option - ! Out-of-core file prefix for storing working E-field solutions (NCI) - character(80) :: prefix + ! Out-of-core file prefix for storing working E-field solutions (NCI) + character(80) :: prefix ! Specify damping parameter for the inversion real(8) :: lambda @@ -72,14 +72,14 @@ module UserCtrl ! Specify the magnitude for random perturbations real(8) :: delta - ! Specify the Covariance Type used in 3D (reserved for future use) - integer :: CovType + ! Specify the Covariance Type used in 3D (reserved for future use) + integer :: CovType ! Indicate how much output you want integer :: output_level - ! Reduce master memory usage by storing E-fields in files (NCI) - logical :: storeSolnsInFile + ! Reduce master memory usage by storing E-fields in files (NCI) + logical :: storeSolnsInFile end type userdef_control @@ -144,7 +144,8 @@ end subroutine initUserCtrl subroutine parseArgs(program,ctrl) character(1) :: job - integer (kind=4) :: iargc,narg,k + integer (kind=4) :: iargc, narg, k + integer (kind=4) :: nOptArgs, optArgSize integer :: istat logical :: exists character*80 gridType, header,arg, verbose, paramtype @@ -169,20 +170,32 @@ subroutine parseArgs(program,ctrl) ! by the verbose level; ignore the others verbose = 'regular' if(narg .gt. 1) then - k=1 - search_arg: & - do - k=k+1 - if (k .eq. narg) exit search_arg - call getarg ( k, arg ) - if(arg(1:1).eq.'-') then - if (arg(2:2).eq.'v') then - call getarg ( k+1, verbose ) - end if - narg=k-1 - exit search_arg + k=1 + nOptArgs = 0 + search_arg: & + do + k=k+1 + if (k .eq. narg) exit search_arg + call getarg(k, arg) + + if(arg(1:1).eq.'-') then + ! Set `optArgSize` to the number of argument parts: e.g. `-P [file_prefix]` + ! has an optArgSize == 2, one part for the `-P` and one part for `[file_prefix]`. + if (arg(2:2).eq.'v') then + call getarg(k+1, verbose) + optArgSize = 2 end if - end do search_arg + + if (arg(2:2) .eq. 'P') then + call getarg(k+1, ctrl % prefix) + ctrl % storeSolnsInFile = .true. + optArgSize = 2 + end if + + nOptArgs = nOptArgs + optArgSize + end if + end do search_arg + narg = narg - nOptArgs end if ! set the level of output based on the user input @@ -524,6 +537,15 @@ subroutine parseArgs(program,ctrl) write(0,*) ' rFile_Cov = the model covariance configuration file' write(0,*) ' rFile_dModel = the starting model parameter perturbation' write(0,*) + write(0,*) 'Optional final arguments:' + write(0,*) ' -v [debug|full|regular|compact|result|none]' + write(0,*) ' indicates the desired level of output to screen and to files' + write(0,*) + write(0,*) ' -P [file_prefix]' + write(0,*) ' indicates that the partial solutions should be stored as ' + write(0,*) ' temporary files on disk rather in memory with file prefix' + write(0,*) ' [file_prefix]' + write(0,*) call ModEM_abort() else ctrl%search = temp(1)