text:
! Copyright 2012-2018, University of Strathclyde
! Authors: Lawrence T. Campbell
! License: BSD-3-Clause

module RK4int

use ParallelInfoType
use TransformInfoType
!use FFTW_Constants

use Globals
use Derivative
use IO
use ParaField

implicit none

  REAL(KIND=WP), DIMENSION(:),ALLOCATABLE :: dadz_r0, dadz_i0
  REAL(KIND=WP), DIMENSION(:),ALLOCATABLE :: dadz_r1, dadz_i1
  REAL(KIND=WP), DIMENSION(:),ALLOCATABLE :: dadz_r2, dadz_i2

  REAL(KIND=WP), DIMENSION(:),ALLOCATABLE :: A_localtr0, A_localti0
  REAL(KIND=WP), DIMENSION(:),ALLOCATABLE :: A_localtr1, A_localti1
  REAL(KIND=WP), DIMENSION(:),ALLOCATABLE :: A_localtr2, A_localti2
  REAL(KIND=WP), DIMENSION(:),ALLOCATABLE :: A_localtr3, A_localti3

  REAL(KIND=WP), DIMENSION(:),ALLOCATABLE :: ac_rfield_in, ac_ifield_in

  REAL(KIND=WP), DIMENSION(:),ALLOCATABLE :: dxdx, dydx, dz2dx, dpxdx, dpydx, dpz2dx

  REAL(KIND=WP), DIMENSION(:), ALLOCATABLE :: dxm, dxt, xt    ! *t is 'temp', for use in next rhs call...
  REAL(KIND=WP), DIMENSION(:), ALLOCATABLE :: dym, dyt, yt
  REAL(KIND=WP), DIMENSION(:), ALLOCATABLE :: dpxm, dpxt, pxt
  REAL(KIND=WP), DIMENSION(:), ALLOCATABLE :: dpym, dpyt, pyt
  REAL(KIND=WP), DIMENSION(:), ALLOCATABLE :: dz2m, dz2t, z2t
  REAL(KIND=WP), DIMENSION(:), ALLOCATABLE :: dpz2m, dpz2t, pz2t

contains

subroutine rk4par(sZ,h,qD)

  implicit none
!
! Perform 4th order Runge-Kutta integration, tailored
! to Puffin and its method of parallelization:
! This is NOT a general, all-purpose RK4 routine, it
! is specific to Puffin. Includes MPI_gathers and
! scatters etc between calculation of derivatives for
! use with the parallel field derivative.
!
!                ARGUMENTS
!
! y       INPUT/OUTPUT   Electron values
! SA      INPUT/OUTPUT   Field values
! x       INPUT          Propagation distance zbar
! h       INPUT          Step size in zbar

!  REAL(KIND=WP),  DIMENSION(:), INTENT(INOUT) :: sA, A_local
  REAL(KIND=WP),  INTENT(IN)                  :: sZ
  REAL(KIND=WP),                INTENT(IN)  :: h
  LOGICAL, INTENT(INOUT) :: qD

!               LOCAL ARGS
!
! h6         Step size divided by 6
! hh         Half of the step size
! xh         x position incremented by half a step
! dym        Intermediate derivatives
! dyt        Intermediate derivatives
! yt         Incremental solution
! dAdx       Field derivative
! dydx       Electron derivatives

  INTEGER(KIND=IP) :: iy,idydx,iyout,i,p
  REAL(KIND=WP)    :: h6, hh, szh
  !REAL(KIND=WP), DIMENSION(size(y)) :: dym, dyt, yt

  REAL(KIND=WP), DIMENSION(:),ALLOCATABLE :: dAdx
  REAL(KIND=WP), DIMENSION(:),ALLOCATABLE :: A_localt
  INTEGER(KIND=IP) :: error, trans

!    Transverse nodes

  trans = (NX_G)*(NY_G)

!    Step sizes

  hh = h * 0.5_WP
  h6 = h / 6.0_WP
  szh = sz + hh

  dadz_r0 = 0_wp
  dadz_r1 = 0_wp
  dadz_r2 = 0_wp
  dadz_i0 = 0_wp
  dadz_i1 = 0_wp
  dadz_i2 = 0_wp

  A_localtr0 = 0_wp
  A_localtr1 = 0_wp
  A_localtr2 = 0_wp
  A_localtr3 = 0_wp
  A_localti0 = 0_wp
  A_localti1 = 0_wp
  A_localti2 = 0_wp
  A_localti3 = 0_wp

  xt = 0_wp
  yt = 0_wp
  z2t = 0_wp
  pxt = 0_wp
  pyt = 0_wp
  pz2t = 0_wp

  dxdx = 0.0_wp
  dydx = 0.0_wp
  dz2dx = 0.0_wp
  dpxdx = 0.0_wp
  dpydx = 0.0_wp
  dpz2dx = 0.0_wp

  dxm = 0.0_wp
  dxt = 0.0_wp
  dym = 0.0_wp
  dyt = 0.0_wp
  dpxm = 0.0_wp
  dpxt = 0.0_wp
  dpym = 0.0_wp
  dpyt = 0.0_wp
  dz2m = 0.0_wp
  dz2t = 0.0_wp
  dpz2m = 0.0_wp
  dpz2t = 0.0_wp

  A_localtr0 = ac_rfield_in
  A_localti0 = ac_ifield_in

!if (count(abs(ac_rfield) > 1.0E2) > 0) print*, 'HELP IM RUBBUSH AT START I habve ', &
!              (count(abs(ac_rfield) > 1.0E2) > 0), 'bigger than 100...'

!  allocate(DADx(2*local_rows))
!  allocate(A_localt(2*local_rows))

!    A_local from A_big

  if (qD) then

!    if (tTransInfo_G%qOneD) then
!       A_local(1:local_rows)=sA(fst_row:lst_row)
!       A_local(local_rows+1:2*local_rows)=&
!            sA(fst_row+iNumberNodes_G:lst_row+iNumberNodes_G)
!    ELSE
!       CALL getAlocalFS(sA,A_local)
!    END if
!
    qD = .false.
!
  end if

!    First step
!  iy = size(sElX_G)
!  idydx = size(dxdx)

!    Get derivatives

  call derivs(sZ, A_localtr0, A_localti0, &
              sElX_G, sElY_G, sElZ2_G, sElPX_G, sElPY_G, sElGam_G, &
              dxdx, dydx, dz2dx, dpxdx, dpydx, dpz2dx, &
              dadz_r0, dadz_i0)

!call mpi_finalize(error)
!stop

!  allocate(dAm(2*local_rows),dAt(2*local_rows))
  !print*, dpydx

!    Increment local electron and field values

  if (qPArrOK_G) then

!$OMP PARALLEL WORKSHARE
    xt = sElX_G      +  hh*dxdx
    yt = sElY_G      +  hh*dydx
    z2t = sElZ2_G    +  hh*dz2dx
    pxt = sElPX_G    +  hh*dpxdx
    pyt = sElPY_G    +  hh*dpydx
    pz2t = sElGam_G  +  hh*dpz2dx

    A_localtr1 = A_localtr0 + hh * dadz_r0
    A_localti1 = A_localti0 + hh * dadz_i0
!$OMP END PARALLEL WORKSHARE

!    Update large field array with new values
!  call local2globalA(A_localt,sA,recvs,displs,tTransInfo_G%qOneD)

    call upd8a(A_localtr1, A_localti1)

  end if

!    Second step
!    Get derivatives

  if (qPArrOK_G) &
    call derivs(szh, A_localtr1, A_localti1, &
       xt, yt, z2t, pxt, pyt, pz2t, &
       dxt, dyt, dz2t, dpxt, dpyt, dpz2t, &
       dadz_r1, dadz_i1)

!    Incrementing with newest derivative value...

  if (qPArrOK_G) then
!$OMP PARALLEL WORKSHARE
    xt = sElX_G      +  hh*dxt
    yt = sElY_G      +  hh*dyt
    z2t = sElZ2_G    +  hh*dz2t
    pxt = sElPX_G    +  hh*dpxt
    pyt = sElPY_G    +  hh*dpyt
    pz2t = sElGam_G  +  hh*dpz2t

    A_localtr2 = A_localtr0 + hh * dadz_r1
    A_localti2 = A_localti0 + hh * dadz_i1
!$OMP END PARALLEL WORKSHARE
!    Update full field array

!  call local2globalA(A_localt,sA,recvs,displs,tTransInfo_G%qOneD)

    call upd8a(A_localtr2, A_localti2)

  end if

!    Third step
!    Get derivatives

  if (qPArrOK_G) &
    call derivs(szh, A_localtr2, A_localti2, &
       xt, yt, z2t, pxt, pyt, pz2t, &
       dxm, dym, dz2m, dpxm, dpym, dpz2m, &
       dadz_r2, dadz_i2)

!    Incrementing

  if (qPArrOK_G) then
!$OMP PARALLEL WORKSHARE
    xt = sElX_G      +  h * dxm
    yt = sElY_G      +  h * dym
    z2t = sElZ2_G    +  h * dz2m
    pxt = sElPX_G    +  h * dpxm
    pyt = sElPY_G    +  h * dpym
    pz2t = sElGam_G  +  h * dpz2m

    A_localtr3 = A_localtr0 + h * dadz_r2
    A_localti3 = A_localti0 + h * dadz_i2
!$OMP END PARALLEL WORKSHARE
!  call local2globalA(A_localt, sA, recvs, displs, tTransInfo_G%qOneD)

    call upd8a(A_localtr3, A_localti3)

!$OMP PARALLEL WORKSHARE
    dxm = dxt + dxm
    dym = dyt + dym
    dz2m = dz2t + dz2m
    dpxm = dpxt + dpxm
    dpym = dpyt + dpym
    dpz2m = dpz2t + dpz2m

    dadz_r2 = dadz_r1 + dadz_r2
    dadz_i2 = dadz_i1 + dadz_i2

    dadz_r1 = 0_wp
    dadz_i1 = 0_wp
!$OMP END PARALLEL WORKSHARE
  end if

!    Fourth step

  szh = sz + h

!    Get derivatives

  if (qPArrOK_G) &
      call derivs(szh, A_localtr3, A_localti3, &
       xt, yt, z2t, pxt, pyt, pz2t, &
       dxt, dyt, dz2t, dpxt, dpyt, dpz2t, &
       dadz_r1, dadz_i1)

!    Accumulate increments with proper weights

  if (qPArrOK_G) then
!$OMP PARALLEL WORKSHARE
    sElX_G    = sElX_G   + h6 * ( dxdx   + dxt   + 2.0_WP * dxm  )
    sElY_G    = sElY_G   + h6 * ( dydx   + dyt   + 2.0_WP * dym  )
    sElZ2_G   = sElZ2_G  + h6 * ( dz2dx  + dz2t  + 2.0_WP * dz2m )
    sElPX_G   = sElPX_G  + h6 * ( dpxdx  + dpxt  + 2.0_WP * dpxm )
    sElPY_G   = sElPY_G  + h6 * ( dpydx  + dpyt  + 2.0_WP * dpym )
    sElGam_G  = sElGam_G + h6 * ( dpz2dx + dpz2t + 2.0_WP * dpz2m)

    ac_rfield_in = ac_rfield_in + h6 * (dadz_r0 + dadz_r1 + 2.0_WP * dadz_r2)
    ac_ifield_in = ac_ifield_in + h6 * (dadz_i0 + dadz_i1 + 2.0_WP * dadz_i2)
!$OMP END PARALLEL WORKSHARE
!  if (count(abs(dadz_r0) > 0.0_wp) <= 0) print*, 'HELP IM TOO RUBBUSH'

!  if (count(abs(ac_rfield) > 0.0_wp) <= 0) print*, 'HELP IM RUBBUSH'

!if (count(abs(ac_rfield) > 1.0E2) > 0) print*, 'HELP IM RUBBUSH for I habve ', &
!              count(abs(ac_rfield) > 1.0E2) , 'bigger than 100...', &
!              'and I am', tProcInfo_G%rank

    call upd8a(ac_rfield_in, ac_ifield_in)

  end if

!  call local2globalA(A_local,sA,recvs,displs,tTransInfo_G%qOneD)

!    Deallocating temp arrays

!  deallocate(dAm,dAt,A_localt)

!  deallocate(DADx)

!  deallocate(dadz_r0, dadz_i0)
!  deallocate(dadz_r1, dadz_i1)
!  deallocate(dadz_r2, dadz_i2)
!
!  deallocate(A_localtr0, A_localti0)
!  deallocate(A_localtr1, A_localti1)
!  deallocate(A_localtr2, A_localti2)
!  deallocate(A_localtr3, A_localti3)
!
!  deallocate(DxDx)
!  deallocate(DyDx)
!  deallocate(DpxDx)
!  deallocate(DpyDx)
!  deallocate(Dz2Dx)
!  deallocate(Dpz2Dx)

!   Set error flag and exit

  GOTO 2000

!   Error Handler - Error log Subroutine in CIO.f90 line 709

1000 CALL Error_log('Error in MathLib:rk4',tErrorLog_G)
  PRINT*,'Error in MathLib:rk4'
2000 CONTINUE

end subroutine rk4par

subroutine allact_rk4_arrs()

  integer(kind=ip) :: tllen43D

  tllen43D = tllen * ntrndsi_G

  allocate(DxDx(iNumberElectrons_G))
  allocate(DyDx(iNumberElectrons_G))
  allocate(DpxDx(iNumberElectrons_G))
  allocate(DpyDx(iNumberElectrons_G))
  allocate(Dz2Dx(iNumberElectrons_G))
  allocate(Dpz2Dx(iNumberElectrons_G))

  allocate(dadz_r0(tllen43D), dadz_i0(tllen43D))
  allocate(dadz_r1(tllen43D), dadz_i1(tllen43D))
  allocate(dadz_r2(tllen43D), dadz_i2(tllen43D))

  allocate(A_localtr0(tllen43D), A_localti0(tllen43D))
  allocate(A_localtr1(tllen43D), A_localti1(tllen43D))
  allocate(A_localtr2(tllen43D), A_localti2(tllen43D))
  allocate(A_localtr3(tllen43D), A_localti3(tllen43D))

  allocate(ac_rfield_in(tllen43D), ac_ifield_in(tllen43D))

  allocate(dxm(iNumberElectrons_G), &
    dxt(iNumberElectrons_G), xt(iNumberElectrons_G))
  allocate(dym(iNumberElectrons_G), &
    dyt(iNumberElectrons_G), yt(iNumberElectrons_G))
  allocate(dpxm(iNumberElectrons_G), &
    dpxt(iNumberElectrons_G), pxt(iNumberElectrons_G))
  allocate(dpym(iNumberElectrons_G), &
    dpyt(iNumberElectrons_G), pyt(iNumberElectrons_G))
  allocate(dz2m(iNumberElectrons_G), &
    dz2t(iNumberElectrons_G), z2t(iNumberElectrons_G))
  allocate(dpz2m(iNumberElectrons_G), &
    dpz2t(iNumberElectrons_G), pz2t(iNumberElectrons_G))

  allocate(dadz_w(iNumberElectrons_G))

  call outer2Inner(ac_rfield_in, ac_ifield_in)

  qInnerXYOK_G = .true.

end subroutine allact_rk4_arrs

subroutine deallact_rk4_arrs()

  call inner2Outer(ac_rfield_in, ac_ifield_in)

  deallocate(ac_rfield_in, ac_ifield_in)

  deallocate(dadz_r0, dadz_i0)
  deallocate(dadz_r1, dadz_i1)
  deallocate(dadz_r2, dadz_i2)

  deallocate(A_localtr0, A_localti0)
  deallocate(A_localtr1, A_localti1)
  deallocate(A_localtr2, A_localti2)
  deallocate(A_localtr3, A_localti3)

  deallocate(DxDx)
  deallocate(DyDx)
  deallocate(DpxDx)
  deallocate(DpyDx)
  deallocate(Dz2Dx)
  deallocate(Dpz2Dx)

  deallocate(dxm, &
    dxt, xt)
  deallocate(dym, &
    dyt, yt)
  deallocate(dpxm, &
    dpxt, pxt)
  deallocate(dpym, &
    dpyt, pyt)
  deallocate(dz2m, &
    dz2t, z2t)
  deallocate(dpz2m, &
    dpz2t, pz2t)

  deallocate(dadz_w)

end subroutine deallact_rk4_arrs

!subroutine RK4_inc

!end subroutine RK4_inc

! Note - the dxdz ad intermediates should all be global,
! if allocating outside of RK4 routine.
!
! They should be passed through to rhs / derivs, and
! be local in there, I think....
!
! All those vars being passed into rhs? GLOBAL.
! Only the arrays should be global.
! Same for the equations module...
!
! Scoop out preamble of rhs, defining temp vars.
! These can all be defined outside this routine to make
! it more readable
!
! label consistently - *_g for main global,
!                      *_rg for rhs global,
!                      *_dg for diffraction global
!
! Check 3D und eqns - are they general? i.e. can kx and ky be anything?
!
! Lj, field4elec (real and imag) and dp2f should be global?
!
! Interface for private var access??
! So have dp2f, field4elec and Lj private arrays in the equations module,
! and interface dp2f and field4elec to rhs.f90 to alter them...
! OR alter them through a subroutine....
! SO Lj is common to both field and e eqns...
! whereas field4elec and dp2f are e only.
! So only making Lj and dp2f private to eqns for now
! In fact, they are just globallay defined ATM until
! we get this working....
! Rename eqns to electron eqns or something...
!
! Make rhs / eqns vars global
! fix dp2f interface in rhs.f90 ... DONE
! only allocate / calc dp2f when it will be used!
!

end module rk4int

max_stars_repo_name: mightylorenzo/Puffin

max_issues_repo_name: mightylorenzo/Puffin

max_forks_repo_name: mightylorenzo/Puffin


!        Generated by TAPENADE     (INRIA, Ecuador team)
!  tapenade 3.x
!
!  Differentiation of set_optip in forward (tangent) mode:
!   variations   of useful results: *dbc_data *bcperiod *propunit
!   with respect to varying inputs: *dbc_data *bcperiod *propunit
!                pvalue
!   Plus diff mem management of: dbc_data:in bcperiod:in propunit:in
! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!>    @brief set parameter value with index "seedi"
!>    @param[in] seedi parameter component/index
!>    @param[in] pvalue value to set
!>    @param[in] ismpl local sample index
SUBROUTINE g_SET_OPTIP(seedi, pvalue, g_pvalue, ismpl)
  USE ARRAYS

  USE g_ARRAYS

  USE MOD_GENRL
  IMPLICIT NONE
  INTEGER :: ismpl
  INTEGER :: seedi, s_k, s_u
  DOUBLE PRECISION :: pvalue
  DOUBLE PRECISION :: g_pvalue
  IF (seedi .LT. 1 .OR. seedi .GT. mpara) THEN
    WRITE(*, *) 'error: in "set_optip", no valid index=', seedi, &
&   ' out of range'
    STOP
  ELSE
    s_k = seed_para(1, seedi)
    s_u = seed_para(2, seedi)
    IF (s_k .LE. lastidx .AND. s_k .GE. firstidx) THEN
!         parameter units
      g_propunit(s_u, s_k, ismpl) = g_pvalue
      propunit(s_u, s_k, ismpl) = pvalue
      CALL g_STAB_PARAM(propunit(s_u, s_k, ismpl), g_propunit(s_u, &
&                   s_k, ismpl), s_k, s_u)
    ELSE IF (s_k .LE. bc_lastidx .AND. s_k .GE. bc_firstidx) THEN
!        bc units
      g_propunit(s_u, s_k, ismpl) = g_pvalue
      propunit(s_u, s_k, ismpl) = pvalue
    ELSE IF (s_k .EQ. -1) THEN
!        tp units
!        warning, the second parameter need an addition of 1,
!          because of restricted range [1..2] instead of [1..3]
      g_bcperiod(opti_tp(1, s_u), opti_tp(2, s_u)+1, opti_tp(3, s_u), &
&     ismpl) = g_pvalue
      bcperiod(opti_tp(1, s_u), opti_tp(2, s_u)+1, opti_tp(3, s_u), &
&     ismpl) = pvalue
    ELSE IF (s_k .EQ. -2) THEN
!        single cell bc
      g_dbc_data(s_u, 1, ismpl) = g_pvalue
      dbc_data(s_u, 1, ismpl) = pvalue
    ELSE
      WRITE(*, *) 'error: "set_optip" has an incorrect index=', seedi, &
&     ' value !'
      STOP
    END IF
    RETURN
  END IF
END SUBROUTINE g_SET_OPTIP


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!
!!!!   MIT License
!!!!
!!!!   ParaMonte: plain powerful parallel Monte Carlo library.
!!!!
!!!!   Copyright (C) 2012-present, The Computational Data Science Lab
!!!!
!!!!   This file is part of the ParaMonte library.
!!!!
!!!!   Permission is hereby granted, free of charge, to any person obtaining a
!!!!   copy of this software and associated documentation files (the "Software"),
!!!!   to deal in the Software without restriction, including without limitation
!!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!!!!   and/or sell copies of the Software, and to permit persons to whom the
!!!!   Software is furnished to do so, subject to the following conditions:
!!!!
!!!!   The above copyright notice and this permission notice shall be
!!!!   included in all copies or substantial portions of the Software.
!!!!
!!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
!!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
!!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
!!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
!!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
!!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
!!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
!!!!
!!!!   ACKNOWLEDGMENT
!!!!
!!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
!!!!   As per the ParaMonte library license agreement terms, if you use any parts of
!!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
!!!!   work (education/research/industry/development/...) by citing the ParaMonte
!!!!   library as described on this page:
!!!!
!!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!>  \brief This module contains procedures for computing the parallel performance of the parallel algorithms.
!>  \author Amir Shahmoradi

module Parallelism_mod

    use Optimization_mod, only: PowellMinimum_type
    use Constants_mod, only: RK, IK
    use Err_mod, only: Err_type

    implicit none

    character(*), parameter :: MODULE_NAME = "@Parallelism_mod"

    !> The [Image_type](@ref image_type) class.
    type                            :: Image_type
        integer(IK)                 :: id                   !<  The ID of the runtime parallel process: 1, 2, 3, ...
        integer(IK)                 :: count                !<  The total number of runtime parallel processes available.
        logical                     :: isFirst = .false.    !<  A logical flag indicating whether the current process is ID #1.
        logical                     :: isNotFirst = .false. !<  A logical flag indicating whether the current process is NOT ID #1.
        logical                     :: isLeader = .false.   !<  A logical flag indicating whether the current process is a leader. This must be defined by the user.
        logical                     :: isRooter = .false.   !<  A logical flag indicating whether the current process is a follower. This must be defined by the user.
        character(:), allocatable   :: name                 !<  The name of the current process as a string with the format `"@process(id)"`.
    contains
        procedure, nopass           :: sync => syncImages
        procedure, pass             :: query => queryImage
        procedure, nopass           :: finalize => finalizeImages
    end type Image_type

    type, private :: Maximum_type
        real(RK)                    :: value            !<  The maximum speedup attained (or attainable).
        integer(IK)                 :: nproc            !<  The required number of processes for maximum speedup.
    end type Maximum_type

    type, private :: Speedup_type
        integer(IK)                 :: count            !<  The size of the `Scaling` vector.
        real(RK)                    :: current          !<  The speedup given the current available number of processes.
        type(Maximum_type)          :: Maximum          !<  An object of type [Maximum_type](@ref maximum_type), containing the predicted maximum speedup and process count.
        real(RK)    , allocatable   :: Scaling(:)       !<  A real vector of length `Speedup_type%count` containing the predicted speedup for a range process counts.
    end type Speedup_type

    type, private :: UniqueProcess_type
        integer(IK)                 :: count            !<  The sizes of the two vector components `Identity` and `Frequency` of `UniqueProcess_type`.
        integer(IK) , allocatable   :: Identity(:)      !<  A vector of size `UniqueProcess_type%count` containing the unique IDs of processes, i.e., the ranks of processes, starting from 1.
        integer(IK) , allocatable   :: Frequency(:)     !<  The frequency with which the process IDs have contributed to the simulation at hand.
    end type UniqueProcess_type

    type, private, extends(UniqueProcess_type) :: Contribution_type
        real(RK) , allocatable      :: LogFrequency(:)  !<  The natural logarithm of the `Frequency` vector component of the type [UniqueProcess_type](@ref uniqueprocess_type).
    end type Contribution_type

    type, private :: SuccessProb_type
        real(RK)                    :: current          !<  The probability of success in the Bernoulli trial. For example, the sampling efficiency of the ParaDRAM sampler.
        real(RK)                    :: effective        !<  The computed effective probability of success inferred from fitting the Cyclic Geometric distribution.
        type(PowellMinimum_type)    :: PowellMinimum    !<  An object of class [PowellMinimum_type](@ref optimization_mod::powellminimum_type) containing
                                                        !<  information about the Cyclic Geometric fit to process contribution data.
    end type SuccessProb_type

    !> The `ForkJoin_type` class.
    type :: ForkJoin_type
        type(UniqueProcess_type)    :: UniqueProcess
        type(Contribution_type)     :: Contribution     !<  Contains information similar to `UniqueProcess`, but including all processes, even zero contributors.
        type(SuccessProb_type)      :: SuccessProb
        type(Speedup_type)          :: Speedup
        type(Err_type)              :: Err
    contains
        procedure   , nopass        :: getForkJoinSpeedup
    end type ForkJoin_type

    interface ForkJoin_type
        module procedure :: constructForkJoin
    end interface ForkJoin_type

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

contains

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    !> \brief
    !> Return the statistics of the parallel processors available, depending on the type of parallelism requested.
    !> This is a dynamic member of the [Image_type](@ref image_type) class.
    !>
    !> @param[inout]    Image       :   An object of class [Image_type](@ref image_type).
    !>                                  On output, all properties of `Image` will be reset,
    !>                                  except the attributes `isLeader` and `isRooter`.
    !>
    !> \warning
    !> This routine must not contain any synchronization statements under any circumstances.
    subroutine queryImage(Image)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: queryImage
#endif
        use Constants_mod, only: RK, IK
        use String_mod, only: num2str
        implicit none
        class(Image_type), intent(inout) :: Image

        ! setup general processor / coarray image variables

#if defined CAF_ENABLED
        Image%id             = this_image()
        Image%count          = num_images()
#elif defined MPI_ENABLED
        block
            use mpi
            integer(IK) :: ierrMPI
            logical     :: isInitialized
            call mpi_initialized( isInitialized, ierrMPI )
            if (.not. isInitialized) call mpi_init(ierrMPI)
            call mpi_comm_rank(mpi_comm_world, Image%id, ierrMPI)
            call mpi_comm_size(mpi_comm_world, Image%count, ierrMPI)
            Image%id = Image%id + 1_IK ! make the ranks consistent with Fortran coarray indexing conventions
        end block
#else
        Image%id            = 1_IK
        Image%count         = 1_IK
#endif
        Image%name          = "@process(" // num2str(Image%id) // ")"
        Image%isFirst       = Image%id==1_IK
        Image%isNotFirst    = Image%id/=1_IK
        !Image%isLeader      = .false.   ! ATTN: this is to be set by the user at runtime, depending on the type of parallelism.
        !Image%isRooter      = .false.   ! ATTN: this is to be set by the user at runtime, depending on the type of parallelism.

    end subroutine queryImage

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    !> \brief
    !> Synchronize all existing parallel images and return nothing.
    !> This is a static member of the [Image_type](@ref image_type) class.
    !>
    !> \warning
    !> This routine global Coarray and MPI synchronization barriers and
    !> therefore, must be called by all processes in the current simulation.
    subroutine syncImages()
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: syncImages
#endif
        implicit none
#if defined MPI_ENABLED
        block
            use mpi
            integer(IK) :: ierrMPI
            !logical     :: isFinalized
            !call mpi_finalized( isFinalized, ierrMPI )
            !if (.not. isFinalized) then
                call mpi_barrier(mpi_comm_world,ierrMPI)
            !end if
        end block
#elif defined CAF_ENABLED
        sync all
#endif
    end subroutine syncImages

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    !> \brief
    !> Finalize the current parallel simulation and return nothing.
    !> This is a static member of the [Image_type](@ref image_type) class.
    !>
    !> \warning
    !> MPI communications will be shut down upon calling this routine and further interprocess communications will be impossible.
    subroutine finalizeImages() ! LCOV_EXCL_LINE
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: finalizeImages
#endif
#if defined CAF_ENABLED
        sync all
#elif defined MPI_ENABLED
        use mpi
        implicit none
        integer(IK) :: ierrMPI
        logical     :: isFinalized
        call mpi_finalized( isFinalized, ierrMPI )
        if (.not. isFinalized) then
            call mpi_barrier(mpi_comm_world,ierrMPI)
            call mpi_finalize(ierrMPI)
        end if
#endif
    end subroutine finalizeImages

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    !> \brief
    !> This is the constructor of the [ForkJoin_type](@ref forkjoin_type) class.
    !> Return the predicted speedup of the parallel simulation given the input characteristics and timing information of the simulation.
    !>
    !> @param[in]   processCount    :   The number of processes in the simulation.
    !> @param[in]   lenProcessID    :   The length of `ProcessID` vector.
    !> @param[in]   ProcessID       :   The vector of process IDs.
    !> @param[in]   successProb     :   The success probability (the effective acceptance rate per objective function call).
    !> @param[in]   seqSecTime      :   The timing of the sequential sections of the code.
    !> @param[in]   parSecTime      :   The timing of the parallel sections of the code.
    !> @param[in]   comSecTime      :   The timing of the communication sections of the code.
    !>
    !> \return
    !> `ForkJoin` : An object of class [ForkJoin_type](@ref forkjoin_type) containing the parallelization speedup.
    function constructForkJoin(processCount, lenProcessID, ProcessID, successProb, seqSecTime, parSecTime, comSecTime) result(ForkJoin) ! nonpure
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: constructForkJoin
#endif
        use GeoCyclicFit_mod, only: fitGeoCyclicLogPDF ! LCOV_EXCL_LINE
        use Constants_mod, only: IK, RK, SQRT_EPS_RK, NEGINF_RK
        use String_mod, only: num2str
        use Misc_mod, only: findUnique
        use Sort_mod, only: indexArray

        implicit none

        integer(IK) , intent(in)    :: processCount, lenProcessID, ProcessID(lenProcessID)
        real(RK)    , intent(in)    :: successProb, seqSecTime, parSecTime, comSecTime
        type(ForkJoin_type)         :: ForkJoin

        character(*), parameter     :: PROCEDURE_NAME = MODULE_NAME // "@constructForkJoin()"
        real(RK)    , parameter     :: LOG_HALF = -0.693147180559945_RK

        integer(IK) , allocatable   :: Indx(:)
        integer(IK)                 :: i

        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        ! check for invalid input
        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        ForkJoin%Err%occurred = .false.

        if (successProb<-SQRT_EPS_RK .or. successProb>1._RK+SQRT_EPS_RK) then
            ForkJoin%Err%occurred = .true.
            ForkJoin%Err%msg = PROCEDURE_NAME // ": successProb must be a number between zero and one. The input value is: " // num2str(successProb)
            return
        end if

        if (processCount<1_IK) then
            ForkJoin%Err%occurred = .true.
            ForkJoin%Err%msg = PROCEDURE_NAME // ": processCount cannot be less than one. The input value is: " // num2str(processCount)
            return
        end if

        if (processCount==1_IK .or. successProb <= SQRT_EPS_RK .or. successProb >= 1._RK+SQRT_EPS_RK) then
            ForkJoin%Speedup%Maximum%value = 1._RK
            ForkJoin%Speedup%Maximum%nproc = 1_IK
            ForkJoin%Speedup%current = 1._RK
            ForkJoin%Speedup%count = 1_IK
            ForkJoin%Speedup%Scaling = [1._RK]
            ForkJoin%Contribution%count = processCount
            ForkJoin%Contribution%Identity = [( i, i = 1, processCount )]
            allocate(ForkJoin%Contribution%Frequency(processCount), source = 0_IK)
            ForkJoin%Contribution%Frequency(1) = lenProcessID
            allocate(ForkJoin%Contribution%LogFrequency(processCount), source = NEGINF_RK)
            ForkJoin%Contribution%LogFrequency(1) = log(real(ForkJoin%Contribution%Frequency(1),kind=RK))
            ForkJoin%UniqueProcess%count = 1_IK
            ForkJoin%UniqueProcess%Identity = [1_IK]
            ForkJoin%UniqueProcess%Frequency = [lenProcessID]
            ForkJoin%SuccessProb%current = successProb
            ForkJoin%SuccessProb%effective = successProb
            return
        end if

        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        ! compute the effective successProb (MCMC efficiency) from the Process contributions
        !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        call findUnique ( lenVector = lenProcessID &
                        , Vector = ProcessID &
                        , UniqueValue = ForkJoin%UniqueProcess%Identity &
                        , UniqueCount = ForkJoin%UniqueProcess%Frequency &
                        , lenUnique = ForkJoin%UniqueProcess%count &
                        )
        !maxContributorProcessID = ForkJoin%UniqueProcess%Identity(maxloc(ForkJoin%UniqueProcess%Frequency))

        ! Sort ascending. This rearrangement may not be necessary anymore, since fitGeoCyclicLogPDF() is now order-agnostic.
        ! From the IO report perspective, however, it is important to have it ordered.

        allocate(Indx(ForkJoin%UniqueProcess%count))
        call indexArray( n = ForkJoin%UniqueProcess%count, Array = ForkJoin%UniqueProcess%Identity, Indx = Indx, Err = ForkJoin%Err )
        if (ForkJoin%Err%occurred) then
            ! LCOV_EXCL_START
            ForkJoin%Err%msg = PROCEDURE_NAME // ForkJoin%Err%msg
            return
            ! LCOV_EXCL_STOP
        end if

        ! get all processes contributions

        ForkJoin%UniqueProcess%Identity(:) = ForkJoin%UniqueProcess%Identity(Indx)
        ForkJoin%UniqueProcess%Frequency(:) = ForkJoin%UniqueProcess%Frequency(Indx)

        deallocate(Indx)

        ForkJoin%Contribution%count = processCount
        if (allocated(ForkJoin%Contribution%LogFrequency)) deallocate(ForkJoin%Contribution%LogFrequency);
        allocate(ForkJoin%Contribution%LogFrequency(ForkJoin%Contribution%count), source = 0._RK)
        if (allocated(ForkJoin%Contribution%Frequency)) deallocate(ForkJoin%Contribution%Frequency);
        allocate(ForkJoin%Contribution%Frequency(ForkJoin%Contribution%count), source = 0_IK)
        if (allocated(ForkJoin%Contribution%Identity)) deallocate(ForkJoin%Contribution%Identity);
        allocate(ForkJoin%Contribution%Identity(ForkJoin%Contribution%count))

        ForkJoin%Contribution%Frequency(ForkJoin%UniqueProcess%Identity) = ForkJoin%UniqueProcess%Frequency
        ForkJoin%Contribution%Identity = [(i, i = 1, ForkJoin%Contribution%count)]
        if (ForkJoin%Contribution%Frequency(1)==0_IK) then
            ! LCOV_EXCL_START
            ForkJoin%Err%occurred = .true.
            ForkJoin%Err%msg = PROCEDURE_NAME//": The contribution of the first process to the simulation is zero. This is highly unusual and requires further investigation."
            return
            ! LCOV_EXCL_STOP
        end if
        do i = 1, ForkJoin%Contribution%count
            ! TODO: xxx a better solution instead of this ad-hoc approach would be to fit for the GeoCyclicCDF. Amir.
            if (ForkJoin%Contribution%Frequency(i)/=0_IK) then
                ForkJoin%Contribution%LogFrequency(i) = log(real(ForkJoin%Contribution%Frequency(i),kind=RK))
            else
                ForkJoin%Contribution%LogFrequency(i) = ForkJoin%Contribution%LogFrequency(i-1) + LOG_HALF
            end if
        end do

        ForkJoin%SuccessProb%current = successProb

        ! fit for the effective successProb

        ForkJoin%SuccessProb%PowellMinimum = fitGeoCyclicLogPDF ( maxNumTrial = processCount & ! LCOV_EXCL_LINE
                                                                , numTrial = ForkJoin%Contribution%count & ! LCOV_EXCL_LINE
                                                                , SuccessStep = ForkJoin%Contribution%Identity & ! LCOV_EXCL_LINE
                                                                , LogCount = ForkJoin%Contribution%LogFrequency & ! LCOV_EXCL_LINE
                                                                )
        if (ForkJoin%SuccessProb%PowellMinimum%Err%occurred) then
            ! LCOV_EXCL_START
            ForkJoin%Err = ForkJoin%SuccessProb%PowellMinimum%Err
            ForkJoin%Err%msg = PROCEDURE_NAME // ForkJoin%Err%msg
            return
            ! LCOV_EXCL_STOP
        end if

        ForkJoin%SuccessProb%effective = ForkJoin%SuccessProb%PowellMinimum%xmin(1)

        ! fit for the effective successProb

        call getForkJoinSpeedup ( successProb = ForkJoin%SuccessProb%effective &                ! max(mcmcSamplingEfficiency, ForkJoin%SuccessProb%PowellMinimum%xmin(1)) & ! avoid unstable estimates of the effective efficiency.
                                , seqSecTime = epsilon(seqSecTime) &                            ! time cost of the sequential section of the code, which is negligible here.
                                , parSecTime = parSecTime &                                     ! the serial runtime for the parallel section of the code.
                                , comSecTimePerProc = comSecTime / (processCount - 1) &         ! the communication overhead for each additional image beyond the leader.
                                , minMaxNumProc = 2 * processCount &                            ! speedup will be computed at least up to this process count, if not more.
                                , Speedup = ForkJoin%Speedup%Scaling &                          ! returned Speedup.
                                , lenSpeedup = ForkJoin%Speedup%count &                         ! length of the returned `Speedup%Scaling` vector.
                                , maxSpeedup = ForkJoin%Speedup%Maximum%value &                 ! returned maximum speedup.
                                , maxSpeedupNumProc = ForkJoin%Speedup%Maximum%nproc &          ! returned number of processes for maximum speedup.
                                , Err = ForkJoin%Err &
                                )
        if (ForkJoin%Err%occurred) then
            ! LCOV_EXCL_START
            ForkJoin%Err%msg = PROCEDURE_NAME // ForkJoin%Err%msg
            return
            ! LCOV_EXCL_STOP
        end if
        ForkJoin%Speedup%current = ForkJoin%Speedup%Scaling(processCount)

    end function constructForkJoin

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    !> \brief
    !> Predict the parallel simulation speedup for a range of possible processor counts.
    !>
    !> @param[in]   successProb         :   The success probability (the effective acceptance rate per objective function call).
    !> @param[in]   seqSecTime          :   The timing of the sequential sections of the code.
    !> @param[in]   parSecTime          :   The timing of the parallel sections of the code.
    !> @param[in]   comSecTimePerProc   :   The timing of the communication sections of the code per processor.
    !> @param[in]   minMaxNumProc       :   The minimum number of processes for which the speedup will be computed. It must be at least 1.
    !> @param[out]  Speedup             :   The vector of speedup values for different counts of processes.
    !> @param[out]  lenSpeedup          :   The length of the vector `Speedup`.
    !> @param[out]  maxSpeedupNumProc   :   The number of processes at which the maximum speedup happens.
    !> @param[out]  maxSpeedup          :   The maximum speedup for any number of processors.
    !> @param[out]  Err                 :   An object of [Err_type](@ref err_mod::err_type) class indicating whether and error has occurred upon return.
    subroutine getForkJoinSpeedup(successProb, seqSecTime, parSecTime, comSecTimePerProc, minMaxNumProc, Speedup, lenSpeedup, maxSpeedupNumProc, maxSpeedup, Err)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: getForkJoinSpeedup
#endif
        use Statistics_mod, only: getLogProbGeoCyclic
        use Constants_mod, only: IK, RK
        use String_mod, only: num2str
        use Misc_mod, only: resize

        implicit none

        real(RK)        , intent(in)                :: successProb, seqSecTime, parSecTime, comSecTimePerProc
        integer(IK)     , intent(in)                :: minMaxNumProc ! must be at least 1
        real(RK)        , intent(out), allocatable  :: Speedup(:)
        integer(IK)     , intent(out)               :: lenSpeedup
        real(RK)        , intent(out)               :: maxSpeedup
        integer(IK)     , intent(out)               :: maxSpeedupNumProc
        type(Err_type)  , intent(out), optional     :: Err

        character(*)    , parameter                 :: PROCEDURE_NAME = MODULE_NAME // "@constructForkJoin()"
        integer(IK)     , parameter                 :: ABS_MAX_NUM_PROC = 1000000_IK

        real(RK)                                    :: FirstImageContribution(1)
        real(RK)                                    :: serialRuntime, parSecTimePerProc, comSecTime
        integer(IK)                                 :: numProc, lenSpeedupNew !, maxNumProc
        logical                                     :: maxSpeedupFound, isPresentErr
        integer(IK), parameter                      :: SuccessStep(*) = [1_IK]

        isPresentErr = present(Err)
        if (isPresentErr) Err%occurred = .false.

        lenSpeedup = minMaxNumProc
        allocate(Speedup(lenSpeedup))

        numProc = 1
        Speedup(numProc) = 1._RK
        maxSpeedup = Speedup(numProc)
        maxSpeedupNumProc = numProc
        serialRuntime = seqSecTime + parSecTime ! serial runtime of the code per function call

        maxSpeedupFound = .false.
        loopOverProc: do

            numProc = numProc + 1

            ! resize allocation if needed

            if (numProc>lenSpeedup) then
                if (maxSpeedupFound) exit loopOverProc
                lenSpeedupNew = 2 * lenSpeedup
                call resize(Speedup, from = lenSpeedup, to = lenSpeedupNew)
                lenSpeedup = lenSpeedupNew
            end if

            ! compute the fraction of work done by the first image

            if (successProb==0._RK) then
                FirstImageContribution(1) = 1._RK / numProc
            else

                FirstImageContribution = exp(getLogProbGeoCyclic( successProb = successProb &
                                                                , maxNumTrial = numProc &
                                                                , numTrial = 1_IK &
                                                                , SuccessStep = SuccessStep &
                                                                ) )
            end if

            ! effective runtime of the parallel-section of the code, when executed in parallel on numProc processes

            parSecTimePerProc = parSecTime * FirstImageContribution(1)

            ! communication time grows linearly with the number of processes.

            comSecTime = (numProc-1_IK) * comSecTimePerProc

            Speedup(numProc) = serialRuntime / (seqSecTime + parSecTimePerProc + comSecTime)
            if (maxSpeedup < Speedup(numProc)) then
                maxSpeedup = Speedup(numProc)
                maxSpeedupNumProc = numProc
            else
                maxSpeedupFound = .true.
            end if

            if (numProc<ABS_MAX_NUM_PROC) cycle loopOverProc ! LCOV_EXCL_LINE

            ! LCOV_EXCL_START
            if (isPresentErr) then
                Err%occurred = .true.
                Err%msg =   PROCEDURE_NAME // &
                            ": Failed to find the number of processes with which the maximum speedup occurs. The search continued up to " // &
                            num2str(ABS_MAX_NUM_PROC) // " processes."
            end if
            return
            ! LCOV_EXCL_STOP

        end do loopOverProc

    end subroutine getForkJoinSpeedup

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end module Parallelism_mod ! LCOV_EXCL_LINE

max_stars_repo_name: cdslaborg/paramonte

max_issues_repo_name: cdslaborg/paramonte

max_forks_repo_name: cdslaborg/paramonte


!
!  Copyright 2016 ARTED developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
!This file is "current.f90"
!This file contain one subroutine.
!SUBROUTINE current
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
Subroutine current
  use Global_Variables
  implicit none
  integer :: ib,ik,i,j,ix,iy,iz,ia,ilma
  real(8) :: jx,jy,jz,jxt,jyt,jzt,x,y,z,kr,jav_l(3)
  complex(8) :: uVpsi,uVpsix,uVpsiy,uVpsiz,zs(3)

  jx=0.d0; jy=0.d0; jz=0.d0
  do ik=NK_s,NK_e
    do ib=1,NBoccmax
      jxt=0d0
      do i=1,NL
        jxt=jxt+real(zu(i,ib,ik))**2+aimag(zu(i,ib,ik))**2
      end do
      jyt=jxt ; jzt=jxt
      jx=jx+occ(ib,ik)*kAc(ik,1)*jxt
      jy=jy+occ(ib,ik)*kAc(ik,2)*jyt
      jz=jz+occ(ib,ik)*kAc(ik,3)*jzt

      jxt=0d0;jyt=0d0;jzt=0d0
      select case(Nd)
      case(4)
      do i=1,NL
        zs(1)=nabx(1)*(zu(ifdx(1,i),ib,ik)-zu(ifdx(-1,i),ib,ik))&
            &+nabx(2)*(zu(ifdx(2,i),ib,ik)-zu(ifdx(-2,i),ib,ik))&
            &+nabx(3)*(zu(ifdx(3,i),ib,ik)-zu(ifdx(-3,i),ib,ik))&
            &+nabx(4)*(zu(ifdx(4,i),ib,ik)-zu(ifdx(-4,i),ib,ik))
        zs(2)=naby(1)*(zu(ifdy(1,i),ib,ik)-zu(ifdy(-1,i),ib,ik))&
            &+naby(2)*(zu(ifdy(2,i),ib,ik)-zu(ifdy(-2,i),ib,ik))&
            &+naby(3)*(zu(ifdy(3,i),ib,ik)-zu(ifdy(-3,i),ib,ik))&
            &+naby(4)*(zu(ifdy(4,i),ib,ik)-zu(ifdy(-4,i),ib,ik))
        zs(3)=nabz(1)*(zu(ifdz(1,i),ib,ik)-zu(ifdz(-1,i),ib,ik))&
            &+nabz(2)*(zu(ifdz(2,i),ib,ik)-zu(ifdz(-2,i),ib,ik))&
            &+nabz(3)*(zu(ifdz(3,i),ib,ik)-zu(ifdz(-3,i),ib,ik))&
            &+nabz(4)*(zu(ifdz(4,i),ib,ik)-zu(ifdz(-4,i),ib,ik))
        jxt=jxt+imag(conjg(zu(i,ib,ik))*zs(1))
        jyt=jyt+imag(conjg(zu(i,ib,ik))*zs(2))
        jzt=jzt+imag(conjg(zu(i,ib,ik))*zs(3))
      enddo
      case default
        call err_finalize('Nd /= 4')
      end select
      jx=jx+jxt*occ(ib,ik)
      jy=jy+jyt*occ(ib,ik)
      jz=jz+jzt*occ(ib,ik)

    enddo
  enddo

  jx=jx*Hxyz/aLxyz
  jy=jy*Hxyz/aLxyz
  jz=jz*Hxyz/aLxyz

  jxt=0d0;jyt=0d0;jzt=0d0
  do ik=NK_s,NK_e
    do ia=1,NI
      do j=1,Mps(ia)
        i=Jxyz(j,ia); ix=Jxx(j,ia); iy=Jyy(j,ia); iz=Jzz(j,ia)
        kr=kAc(ik,1)*(Lx(i)*Hx-ix*aLx)+kAc(ik,2)*(Ly(i)*Hy-iy*aLy)+kAc(ik,3)*(Lz(i)*Hz-iz*aLz)
        ekr(j,ia)=exp(zI*kr)
      enddo
    end do
    do ib=1,NBoccmax
      do ilma=1,Nlma
        ia=a_tbl(ilma)
        uVpsi=0.d0; uVpsix=0.d0; uVpsiy=0.d0; uVpsiz=0.d0
        do j=1,Mps(ia)
          i=Jxyz(j,ia); ix=Jxx(j,ia); iy=Jyy(j,ia); iz=Jzz(j,ia)
          x=Lx(i)*Hx-ix*aLx
          y=Ly(i)*Hy-iy*aLy
          z=Lz(i)*Hz-iz*aLz
          uVpsi =uVpsi +uV(j,ilma)*ekr(j,ia)  *zu(i,ib,ik)
          uVpsix=uVpsix+uV(j,ilma)*ekr(j,ia)*x*zu(i,ib,ik)
          uVpsiy=uVpsiy+uV(j,ilma)*ekr(j,ia)*y*zu(i,ib,ik)
          uVpsiz=uVpsiz+uV(j,ilma)*ekr(j,ia)*z*zu(i,ib,ik)
        end do
        uVpsi =uVpsi *Hxyz*iuV(ilma)
        uVpsix=uVpsix*Hxyz
        uVpsiy=uVpsiy*Hxyz
        uVpsiz=uVpsiz*Hxyz
        jxt=jxt+occ(ib,ik)/aLxyz*2*imag(conjg(uVpsix)*uVpsi)
        jyt=jyt+occ(ib,ik)/aLxyz*2*imag(conjg(uVpsiy)*uVpsi)
        jzt=jzt+occ(ib,ik)/aLxyz*2*imag(conjg(uVpsiz)*uVpsi)
      enddo
    enddo
  enddo

  jx=jx+jxt
  jy=jy+jyt
  jz=jz+jzt

  jav_l(1)=jx
  jav_l(2)=jy
  jav_l(3)=jz

  call MPI_ALLREDUCE(jav_l,jav,3,MPI_REAL8,MPI_SUM,NEW_COMM_WORLD,ierr)

  return
end Subroutine current
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
Subroutine current_GS
  use Global_Variables
  implicit none
  integer :: ib,ik,i,j,ix,iy,iz,ia,ilma
  real(8) :: jx(NL),jy(NL),jz(NL),x,y,z,kr,jav_l(3)
  complex(8) :: uVpsi,uVpsix,uVpsiy,uVpsiz,zs(3)

  jx=0.d0; jy=0.d0; jz=0.d0
  do ik=NK_s,NK_e
    do ib=1,NBoccmax
      do i=1,NL
        jx(i)=jx(i)+occ(ib,ik)*kAc(ik,1)*abs(zu_GS(i,ib,ik))**2
        jy(i)=jy(i)+occ(ib,ik)*kAc(ik,2)*abs(zu_GS(i,ib,ik))**2
        jz(i)=jz(i)+occ(ib,ik)*kAc(ik,3)*abs(zu_GS(i,ib,ik))**2
      enddo

      select case(Nd)
      case(4)
      do i=1,NL
        zs(1)=nabx(1)*(zu_GS(ifdx(1,i),ib,ik)-zu_GS(ifdx(-1,i),ib,ik))&
            &+nabx(2)*(zu_GS(ifdx(2,i),ib,ik)-zu_GS(ifdx(-2,i),ib,ik))&
            &+nabx(3)*(zu_GS(ifdx(3,i),ib,ik)-zu_GS(ifdx(-3,i),ib,ik))&
            &+nabx(4)*(zu_GS(ifdx(4,i),ib,ik)-zu_GS(ifdx(-4,i),ib,ik))
        zs(2)=naby(1)*(zu_GS(ifdy(1,i),ib,ik)-zu_GS(ifdy(-1,i),ib,ik))&
            &+naby(2)*(zu_GS(ifdy(2,i),ib,ik)-zu_GS(ifdy(-2,i),ib,ik))&
            &+naby(3)*(zu_GS(ifdy(3,i),ib,ik)-zu_GS(ifdy(-3,i),ib,ik))&
            &+naby(4)*(zu_GS(ifdy(4,i),ib,ik)-zu_GS(ifdy(-4,i),ib,ik))
        zs(3)=nabz(1)*(zu_GS(ifdz(1,i),ib,ik)-zu_GS(ifdz(-1,i),ib,ik))&
            &+nabz(2)*(zu_GS(ifdz(2,i),ib,ik)-zu_GS(ifdz(-2,i),ib,ik))&
            &+nabz(3)*(zu_GS(ifdz(3,i),ib,ik)-zu_GS(ifdz(-3,i),ib,ik))&
            &+nabz(4)*(zu_GS(ifdz(4,i),ib,ik)-zu_GS(ifdz(-4,i),ib,ik))
        jx(i)=jx(i)+occ(ib,ik)*imag(conjg(zu_GS(i,ib,ik))*zs(1))
        jy(i)=jy(i)+occ(ib,ik)*imag(conjg(zu_GS(i,ib,ik))*zs(2))
        jz(i)=jz(i)+occ(ib,ik)*imag(conjg(zu_GS(i,ib,ik))*zs(3))
      enddo
      case default
        call err_finalize('Nd /= 4')
      end select

    enddo
  enddo

  jav_l(1)=sum(jx)*Hxyz/aLxyz
  jav_l(2)=sum(jy)*Hxyz/aLxyz
  jav_l(3)=sum(jz)*Hxyz/aLxyz

  do ik=NK_s,NK_e
    do ia=1,NI
      do j=1,Mps(ia)
        i=Jxyz(j,ia); ix=Jxx(j,ia); iy=Jyy(j,ia); iz=Jzz(j,ia)
        kr=kAc(ik,1)*(Lx(i)*Hx-ix*aLx)+kAc(ik,2)*(Ly(i)*Hy-iy*aLy)+kAc(ik,3)*(Lz(i)*Hz-iz*aLz)
        ekr(j,ia)=exp(zI*kr)
      enddo
    enddo
    do ib=1,NBoccmax
      do ilma=1,Nlma
        ia=a_tbl(ilma)
        uVpsi=0.d0; uVpsix=0.d0; uVpsiy=0.d0; uVpsiz=0.d0
        do j=1,Mps(ia)
          i=Jxyz(j,ia); ix=Jxx(j,ia); iy=Jyy(j,ia); iz=Jzz(j,ia)
          x=Lx(i)*Hx-ix*aLx
          y=Ly(i)*Hy-iy*aLy
          z=Lz(i)*Hz-iz*aLz
          uVpsi =uVpsi +uV(j,ilma)*ekr(j,ia)  *zu_GS(i,ib,ik)
          uVpsix=uVpsix+uV(j,ilma)*ekr(j,ia)*x*zu_GS(i,ib,ik)
          uVpsiy=uVpsiy+uV(j,ilma)*ekr(j,ia)*y*zu_GS(i,ib,ik)
          uVpsiz=uVpsiz+uV(j,ilma)*ekr(j,ia)*z*zu_GS(i,ib,ik)
        end do
        uVpsi =uVpsi *Hxyz*iuV(ilma)
        uVpsix=uVpsix*Hxyz
        uVpsiy=uVpsiy*Hxyz
        uVpsiz=uVpsiz*Hxyz
        jav_l(1)=jav_l(1)+occ(ib,ik)/aLxyz*2*imag(conjg(uVpsix)*uVpsi)
        jav_l(2)=jav_l(2)+occ(ib,ik)/aLxyz*2*imag(conjg(uVpsiy)*uVpsi)
        jav_l(3)=jav_l(3)+occ(ib,ik)/aLxyz*2*imag(conjg(uVpsiz)*uVpsi)
      enddo
    enddo
  enddo
  call MPI_ALLREDUCE(jav_l,jav,3,MPI_REAL8,MPI_SUM,NEW_COMM_WORLD,ierr)

  return
End Subroutine Current_GS
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
#ifdef ARTED_CURRENT_OPTIMIZED
subroutine current_omp_KB_ST(ib,ik,A)
  use Global_Variables
  use opt_variables
  implicit none
  integer,intent(in)    :: ib,ik
  complex(8),intent(in) :: A(0:NL-1)
  real(8) :: nabt(12),zx,zy,zz

  nabt( 1: 4) = nabx(1:4)
  nabt( 5: 8) = naby(1:4)
  nabt( 9:12) = nabz(1:4)

  zx = 0.d0
  zy = 0.d0
  zz = 0.d0
  call current_stencil(nabt,A,zx,zy,zz)

  zcx(ib,ik)=zx
  zcy(ib,ik)=zy
  zcz(ib,ik)=zz
end subroutine
#endif
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
Subroutine current_GS_omp_KB
  use Global_Variables
  use timelog
  implicit none
  real(8) :: jx,jy,jz

  call timelog_begin(LOG_CURRENT)
  call current_GS_omp_KB_impl(zu_GS,jx,jy,jz)
  call current_result_impl(jx,jy,jz)
  call timelog_end(LOG_CURRENT)
end Subroutine
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
Subroutine current_omp_KB
  use Global_Variables
  use timelog
  implicit none
  real(8) :: jx,jy,jz

  call timelog_begin(LOG_CURRENT)
  call current_omp_KB_impl(zu,jx,jy,jz)
  call current_result_impl(jx,jy,jz)
  call timelog_end(LOG_CURRENT)
end Subroutine
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
Subroutine current_init_impl(ik,ib,zutmp,jx,jy,jz)
  use Global_Variables
  implicit none
  integer,intent(in)    :: ik,ib
  complex(8),intent(in) :: zutmp(0:NL-1)
  real(8),intent(out)   :: jx,jy,jz

  integer :: i
  real(8) :: jt

  jt=0.d0
!dir$ vector aligned
  do i=0,NL-1
    jt=jt+real(zutmp(i))**2+imag(zutmp(i))**2
  end do

  jx=occ(ib,ik)*kAc(ik,1)*jt
  jy=occ(ib,ik)*kAc(ik,2)*jt
  jz=occ(ib,ik)*kAc(ik,3)*jt
end Subroutine
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
Subroutine current_stencil_impl(ik,ib,zutmp,nabt,jx,jy,jz)
  use Global_Variables
  use opt_variables
  implicit none
  integer,intent(in)    :: ik,ib
  complex(8),intent(in) :: zutmp(0:NL-1)
  real(8),intent(in)    :: nabt(12)
  real(8),intent(inout) :: jx,jy,jz
  real(8) :: jxt,jyt,jzt

  jxt=0.d0
  jyt=0.d0
  jzt=0.d0
  call current_stencil(nabt,zutmp,jxt,jyt,jzt)

  jx=jx+jxt*occ(ib,ik)
  jy=jy+jyt*occ(ib,ik)
  jz=jz+jzt*occ(ib,ik)
end Subroutine
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
#ifdef ARTED_CURRENT_OPTIMIZED
Subroutine current_RT_stencil_impl(ik,ib,jx,jy,jz)
  use Global_Variables
  use opt_variables
  implicit none
  integer,intent(in)    :: ik,ib
  real(8),intent(inout) :: jx,jy,jz

  jx=jx+zcx(ib,ik)*occ(ib,ik)
  jy=jy+zcy(ib,ik)*occ(ib,ik)
  jz=jz+zcz(ib,ik)*occ(ib,ik)
end Subroutine
#endif
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
Subroutine current_pseudo_impl(ik,ib,zutmp,IaLxyz,jx,jy,jz)
  use Global_Variables
  use omp_lib
  use opt_variables
  use timelog
  implicit none
  integer,intent(in)    :: ik,ib
  complex(8),intent(in) :: zutmp(NL)
  real(8),intent(in)    :: IaLxyz
  real(8),intent(out)   :: jx,jy,jz

  integer    :: ilma,ia,j,i,ix,iy,iz
  real(8)    :: x,y,z
  real(8)    :: jxt,jyt,jzt
  complex(8) :: uVpsi,uVpsix,uVpsiy,uVpsiz

  jxt=0d0
  jyt=0d0
  jzt=0d0

  do ilma=1,Nlma
    ia=a_tbl(ilma)
    uVpsi=0.d0; uVpsix=0.d0; uVpsiy=0.d0; uVpsiz=0.d0
    do j=1,Mps(ia)
      i=Jxyz(j,ia)

      ix=Jxx(j,ia); x=Lx(i)*Hx-ix*aLx
      iy=Jyy(j,ia); y=Ly(i)*Hy-iy*aLy
      iz=Jzz(j,ia); z=Lz(i)*Hz-iz*aLz

      uVpsi =uVpsi +uV(j,ilma)*ekr_omp(j,ia,ik)  *zutmp(i)
      uVpsix=uVpsix+uV(j,ilma)*ekr_omp(j,ia,ik)*x*zutmp(i)
      uVpsiy=uVpsiy+uV(j,ilma)*ekr_omp(j,ia,ik)*y*zutmp(i)
      uVpsiz=uVpsiz+uV(j,ilma)*ekr_omp(j,ia,ik)*z*zutmp(i)
    end do
    uVpsi =uVpsi *Hxyz*iuV(ilma)
    uVpsix=uVpsix*Hxyz
    uVpsiy=uVpsiy*Hxyz
    uVpsiz=uVpsiz*Hxyz
    jxt=jxt+occ(ib,ik)*IaLxyz*2*imag(conjg(uVpsix)*uVpsi)
    jyt=jyt+occ(ib,ik)*IaLxyz*2*imag(conjg(uVpsiy)*uVpsi)
    jzt=jzt+occ(ib,ik)*IaLxyz*2*imag(conjg(uVpsiz)*uVpsi)
  enddo

  jx=jx*Hxyz*IaLxyz+jxt
  jy=jy*Hxyz*IaLxyz+jyt
  jz=jz*Hxyz*IaLxyz+jzt
end Subroutine
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
Subroutine current_GS_omp_KB_impl(zutmp,jxs,jys,jzs)
  use Global_Variables
  implicit none
  complex(8),intent(in) :: zutmp(0:NL-1,NB,NK_s:NK_e)
  real(8),intent(out)   :: jxs,jys,jzs

  integer :: ikb,ib,ik,i,j,ix,iy,iz,ia
  real(8) :: kr,jx,jy,jz,IaLxyz
  real(8) :: nabt(12)

  nabt( 1: 4) = nabx(1:4)
  nabt( 5: 8) = naby(1:4)
  nabt( 9:12) = nabz(1:4)

  IaLxyz = 1.0 / aLxyz

  jxs=0.d0
  jys=0.d0
  jzs=0.d0
!$omp parallel reduction(+:jxs,jys,jzs)

!Constructing nonlocal part
!$omp do private(ik,ia,j,i,ix,iy,iz,kr) collapse(2)
  do ik=NK_s,NK_e
  do ia=1,NI
  do j=1,Mps(ia)
    i=Jxyz(j,ia); ix=Jxx(j,ia); iy=Jyy(j,ia); iz=Jzz(j,ia)
    kr=kAc(ik,1)*(Lx(i)*Hx-ix*aLx)+kAc(ik,2)*(Ly(i)*Hy-iy*aLy)+kAc(ik,3)*(Lz(i)*Hz-iz*aLz)
    ekr_omp(j,ia,ik)=exp(zI*kr)
  end do
  end do
  end do
!$omp end do

!$omp do private(ikb,ik,ib,jx,jy,jz)
  do ikb=1,NKB
    ik=ik_table(ikb)
    ib=ib_table(ikb)

    call current_init_impl(ik,ib,zutmp(:,ib,ik),jx,jy,jz)
    call current_stencil_impl(ik,ib,zutmp(:,ib,ik),nabt,jx,jy,jz)
    call current_pseudo_impl(ik,ib,zutmp(:,ib,ik),IaLxyz,jx,jy,jz)

    jxs=jxs+jx
    jys=jys+jy
    jzs=jzs+jz
  enddo
!$omp end parallel
end Subroutine
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
Subroutine current_omp_KB_impl(zutmp,jxs,jys,jzs)
  use Global_Variables
  implicit none
  complex(8),intent(in) :: zutmp(0:NL-1,NBoccmax,NK_s:NK_e)
  real(8),intent(out)   :: jxs,jys,jzs

  integer :: ikb,ib,ik,i,j,ix,iy,iz,ia
  real(8) :: kr,jx,jy,jz,IaLxyz
  real(8) :: nabt(12)

  nabt( 1: 4) = nabx(1:4)
  nabt( 5: 8) = naby(1:4)
  nabt( 9:12) = nabz(1:4)

  IaLxyz = 1.0 / aLxyz

  jxs=0.d0
  jys=0.d0
  jzs=0.d0
!$omp parallel reduction(+:jxs,jys,jzs)

!Constructing nonlocal part
!$omp do private(ik,ia,j,i,ix,iy,iz,kr) collapse(2)
  do ik=NK_s,NK_e
  do ia=1,NI
  do j=1,Mps(ia)
    i=Jxyz(j,ia); ix=Jxx(j,ia); iy=Jyy(j,ia); iz=Jzz(j,ia)
    kr=kAc(ik,1)*(Lx(i)*Hx-ix*aLx)+kAc(ik,2)*(Ly(i)*Hy-iy*aLy)+kAc(ik,3)*(Lz(i)*Hz-iz*aLz)
    ekr_omp(j,ia,ik)=exp(zI*kr)
  end do
  end do
  end do
!$omp end do

!$omp do private(ikb,ik,ib,jx,jy,jz)
  do ikb=1,NKB
    ik=ik_table(ikb)
    ib=ib_table(ikb)

    call current_init_impl(ik,ib,zutmp(:,ib,ik),jx,jy,jz)
#ifdef ARTED_CURRENT_OPTIMIZED
    call current_RT_stencil_impl(ik,ib,jx,jy,jz)
#else
    call current_stencil_impl(ik,ib,zutmp(:,ib,ik),nabt,jx,jy,jz)
#endif
    call current_pseudo_impl(ik,ib,zutmp(:,ib,ik),IaLxyz,jx,jy,jz)

    jxs=jxs+jx
    jys=jys+jy
    jzs=jzs+jz
  enddo
!$omp end parallel
end Subroutine
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
Subroutine current_result_impl(jx,jy,jz)
  use Global_Variables
  use timelog
  implicit none
  real(8),intent(in) :: jx,jy,jz

  real(8) :: jav_l(3)

  call timelog_begin(LOG_ALLREDUCE)
  jav_l(1)=jx
  jav_l(2)=jy
  jav_l(3)=jz

  call MPI_ALLREDUCE(jav_l,jav,3,MPI_REAL8,MPI_SUM,NEW_COMM_WORLD,ierr)
  call timelog_end(LOG_ALLREDUCE)
end Subroutine


text:
Armando Figueroa was a candidate for ASUCD ASUCD Senate Senate in the Fall 2012 ASUCD Election running on the SMART slate, winning with 669 votes.  As a senator he proved to be both a fierce defender of the communities he came from while working with the exec to compromise on the budget and other issues.  He is now a candidate for ASUCD President with Maxwell Kappes as his VicePresidential candidate.  They are running on a SMART/Independent fusion ticket.

Experience:
  ASUCD Intern Under Tatiana Bush 20102011

  UC Davis Silver Volunteer Service Award 2011

  Educational Opportunity Program (EOP) Peer Advising Counselor 20122014

  Movimiento Estudiantil Chicano de Aztlán Movimiento Estudiantil Chican@ de Aztlán (MEChA) Political Chair 20122013

  Intern at the SRRC Student Recruitment & Retention Center with African Diaspora Cultivating Education (ACE) 20102013

  Casa Cuauhtémoc Cultural Programmer 20112012

 Presidential Candidate Statement

Peace Community,

I’m Armando Figueroa and I previously served as ASUCD senator, but more importantly, a student leader serving all underrepresented progressive student communities, working with various departments such EOP, Student Housing and the SRRC, as well as many student organizations on campus.

I’m Maxwell Kappes, and I previously served as your ASUCD Senator as well, representing the more than six hundred clubs on campus with my experiences from leading Davis Urban Gaming Group, Humans vs. Zombies, and serving as an officer on several other groups.

Together we are running for President and Vice President of ASUCD. We are a coalition of SMART and Independent, looking to serve all students.

Make the Transition Out of the Dorms Easier

When student housing finally welcomes first years and transfers into the community, it can be difficult for students to manage new parts of their lives alongside increasingly challenging schoolwork. We intend to expand Counseling and Psychology Services (CAPS) to include conflict resolution services, specifically to assist with roommate issues and dealing with the challenges of a group living dynamic that apartment dwelling frequently presents. Additionally, we want to expand the Student Farm’s relationship with the Pantry, so students new on their own know what foods to eat and how to cook them.

Advocate for More Study Spaces, 24Hour Computer Labs

As senators, both of us were involved in expanding student involvement in campus renovations and construction. Using that experience, we will serve as fierce advocates for an expansion of roundtheclock studying spaces and computer labs. By the end of our term, when someone mentions “the 24hour room,” you’ll have to ask them to clarify which one.

Bring More Student Services to Everyone

ASUCD is remarkable in that it provides more distinct services than any other UC student government. However, many of these services, such as The Pantry, a food bank for undergraduates, and Aggie Legal Services, which offers free consultation from a lawyer, and underutilized and obscured from the bureaucracy that threatens to reign over them. We want to give all of our units a better opportunity to actually serve students, as well as increase the services provided to all Aggies, including affordable senior portraits and professional audio recording.

Create SMARTer Connections with Student Organizations and Units

Student organizations frame most Aggies’ UC Davis experience, but frequently, it can be difficult to find your niche. The Activities Fair, currently the easiest way to interact with the most amount of clubs, is nearly halfway into fall quarter, long after most rush events and welcome weeks. By working with student housing and the Center for Student Involvement to move this key event to week zero, first year Aggies can get a better understanding of all of the opportunities provided by their peers. We intend to involve ASUCD’s resources in club life beyond the first week of school as well, providing more collaborative opportunities for student organizations.

Active Campus, Active Aggies

Campus life should not end before the buses stop running. After a midterm, or during a study break, being able to relax and interact with your fellow Aggies are among the best ways to maintain a high level of mental health. Our proposed activities range from weekly trivia nights at the oncampus pub to rentable sporting goods at the quad, with frequent open mic nights and karaoke events. We want to see UC Davis as a center of student engagement, not just the place you go to take classes.

Promote Student Advocacy to Administration, not the other way around

Over the last decade, communications between high level administration and the student body has decreased to the point that only the executive office has frequent contact with those in positions that can enact real change. As your executives, we want to open these communications up to all students, as well as serving as the loudest voice in the room. With us, administration won’t be able to try to raise tuition or cut more sports because there will too many students at the table.

Vote for a better, SMARTer future. VOTE FIGUEROA/KAPPES for ASUCD President and Vice President and VOTE S.M.A.R.T. OR Independent for ASUCD Senate.

 Senate Candidate Statement

Peace Community,

My name is Armando Figueroa and I am running for ASUCD Senate this quarter with the Students Matter Activism Retention Teamwork Slate. S.M.A.R.T is a group of socially and politically enlightened students who care about the communities that support them, and work to make decisions for the university that benefit everyone as a whole. I have decided to run because I want to give proper representation to all student groups, and I genuinely care about improving student life on this campus.

As an ASUCD senator, I will address the following and more:

UNIONIZE THE ASUCD COFFEE HOUSE

The first way I will advocate for student organizations is by making facilities more accessible. New campus policy dictates that campus organizations need to pay to utilize the space the university provides, including classrooms and conference rooms. I will pass a bill that will allow students to apply to have ASUCD waive the charge for student organizations who lack the funds to cover it.

Emiliano Zapata once said “the land belongs to those who work it.” I firmly believe that without organization, student workers are in vulnerable positions. So to better support student workers, I will organize outlets for us to dialogue about pressing issues. Furthermore, I will initiate communication with various campus departments to find ways in which we can holistically support all student workers on this campus. I will create an effective medium for student workers to network as well. In order to truly be a model institution, we must look after those who put in the extra hours without pay, who innovate on a daily basis to make our campus attractive and inclusive. Most importantly, we must look after those who enable our university to meet it’s full potential.

CLOSE THE DISCONNECT BETWEEN ADMINISTRATION AND STUDENTS

I will not only work to improve the communication with administration and students, but ensure that it is effective by prioritizing student concerns. I will create more dialogues and programs through various effective avenues, such as working with the student undergraduate advisory board. This will allow for greater student interaction with administration. After the incidents of Nov. 18, I want to ensure students that administration accountability can and will happen, as long as we all engage in the simple collaboration that is needed and wanted from both parties.

VOTE TO PRODUCE THE RESULTS WE WILL BRING ABOUT TOGETHER. VOTE ARMANDO FIGUEROA FOR YOUR #1. VOTE S.M.A.R.T FOR YOUR #2#5

20131009 11:46:40 nbsp The bromance is strong with this one.  Users/MaxwellKappes

Oh god

20140225 16:39:39 nbsp oh god Users/ArmandoFigueroa

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!##############################################################################
!# ****************************************************************************
!# <name> main_program </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module solves an optimal control problem for the stationary and
!# nonstationary Navier-Stokes optimal control problem
!#
!#  $$ min J(y,u) = 1/2||y-z||_{L^2} + \gamma/2||y(T)-z(T)||_{L^2} + \alphga/2||u||^2 $$
!#
!#  $$- \nu Delta(y) + y*\Nabla(y) + \Nabla p = f $$
!#  $$ \Nabla \cdot y = 0$$
!#  $$- \nu Delta(\lambda) - y*\Nabla(\lambda) + \lambda\Nabla y + \Nabla \xi = y-z $$
!#  $$ \Nabla \cdot \lambda = 0$$
!#
!#
!# on a 2D domain for a 2D function $y=(y_1,y_2)$, a pressure $p$,
!# a dual velocity $\lambda$ and a dual pressure $\xi$. $u$ is the control
!# and $z$ a desired target function field.
!#
!# The routine splits up the tasks of reading the domain, creating
!# triangulations, discretisation, solving, postprocessing and creanup into
!# different subroutines. The communication between these subroutines
!# is done using an application-specific structure saving problem data
!# as well as a collection structure for the communication with callback
!# routines.
!#
!# For the nonlinearity, the nonlinear solver is invoked. The
!# defect that is setted up there is preconditioned by a linear Multigrid
!# solver with a simple-VANKA smoother/preconditioner for
!# 2D saddle point problems, Jacobi-Type. As coarse grid solver,
!# UMFPACK is used.
!# </purpose>
!##############################################################################

module main_program

  use fsystem
  use genoutput
  use storage
  use linearsolver
  use boundary
  use bilinearformevaluation
  use linearformevaluation
  use cubature
  use matrixfilters
  use vectorfilters
  use bcassembly
  use triangulation
  use spatialdiscretisation
  use coarsegridcorrection
  use spdiscprojection
  use nonlinearsolver
  use paramlist
  use fparser
  use statistics

  use collection
  use convection

  use externalstorage
  use paramlist

  use spacetimevectors
  use analyticsolution

  use structurespostproc
  use structuresmain

  use structuresoptflow
  use kktsystem
  use kktsystemhierarchy
  use initoptflow
  use newtoniteration

  use postprocessing

  implicit none

!<globals>

  ! Directory containing the data files.
  character(LEN=SYS_STRLEN), save :: DIR_DATA = "./data";

!</globals>

contains

  ! ***************************************************************************

!<subroutine>

  subroutine cc2doptc_getLogFiles (slogfile,serrorfile)

!<description>
  ! Temporarily reads the output DAT file to get the names of the output
  ! files.
!</description>

!<output>
  ! Name of the message log file.
  character(LEN=*), intent(OUT) :: slogfile

  ! Name of the error log file.
  character(LEN=*), intent(OUT) :: serrorfile
!</output>

!</subroutine>

    type(t_parlist) :: rparlist

    ! Init parameter list that accepts parameters for output files
    call parlst_init (rparlist)

    ! Read parameters that configure the output
    call parlst_readfromfile (rparlist, trim(DIR_DATA)//"/output.dat")

    ! Now the real initialisation of the output including log file stuff!
    call parlst_getvalue_string (rparlist,"GENERALOUTPUT",&
        "smsgLog",slogfile,"",bdequote=.true.)

    call parlst_getvalue_string (rparlist,"GENERALOUTPUT",&
        "serrorLog",serrorfile,"",bdequote=.true.)

    ! That temporary parameter list is not needed anymore.
    call parlst_done (rparlist)

    end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine cc2doptc_evalParameters ()

!<description>
  ! Evaluates command line parameters.
  ! In the current implementation, command line parameters are passed as
  ! text file. This routine searches in the main directory for a file
  ! "cmdline.dat". If this file is found, is it opened and evaluated.
  ! Every line may contain a command line parameter in the form of
  ! a DAT file (name=value pairs).
  !
  ! Supported command line parameters:
  !   "datdirectory = [Directory, where DAT files can be found]"
!</description>

!</subroutine>

    ! local variables
    type(t_parlist) :: rparamList
    logical :: bexists

    ! Figure out if the file exists.
    inquire(file="./cmdline.dat", exist=bexists)

    if (bexists) then
      ! Read the file
      call parlst_init (rparamList)
      call parlst_readfromfile (rparamList, "./cmdline.dat")

      ! Evaluate parameters
      call parlst_getvalue_string ( &
          rparamList, "","datdirectory", DIR_DATA, DIR_DATA,bdequote=.true.)

      call parlst_done (rparamList)

    end if

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine main_getDat (rparamList)

!<description>
  ! Reads in all DAT files into the parameter list rparlist
!</description>

!<inputoutput>
  ! The parameter list where the values of the DAT files should be stored.
  ! The structure must have been initialised, the parameters are just added
  ! to the list.
  type(t_parlist), intent(INOUT) :: rparamList
!</inputoutput>

!</subroutine>

    logical :: bexists
    character(LEN=SYS_STRLEN) :: smaster

    ! Check if a command line parameter specifies the master.dat file.
    call sys_getcommandLineArg(1,smaster,sdefault=trim(DIR_DATA)//"/master.dat")

    ! Read the file "master.dat".
    ! If that does not exist, try to manually read files with parameters from a
    ! couple of files.
    inquire(file=smaster, exist=bexists)

    if (bexists) then
      ! Read the master file. That either one contains all parameters or
      ! contains references to subfiles with data.
      call parlst_readfromfile (rparamList, trim(smaster),trim(DIR_DATA))
    else
      call output_line("Master file not found: "//trim(smaster),&
          OU_CLASS_WARNING,ssubroutine="main_getDat")
      call output_line("Reading standard parameters.",&
          OU_CLASS_WARNING,ssubroutine="main_getDat")

      ! Each "readfromfile" command adds the parameter of the specified file
      ! to the parameter list.
      call parlst_readfromfile (rparamList, trim(DIR_DATA)//"/main.dat")
      call parlst_readfromfile (rparamList, trim(DIR_DATA)//"/bdconditions.dat")
      call parlst_readfromfile (rparamList, trim(DIR_DATA)//"/discretisation.dat")
      call parlst_readfromfile (rparamList, trim(DIR_DATA)//"/flows.dat")
      call parlst_readfromfile (rparamList, trim(DIR_DATA)//"/linsol.dat")
      call parlst_readfromfile (rparamList, trim(DIR_DATA)//"/optcontrol.dat")
      call parlst_readfromfile (rparamList, trim(DIR_DATA)//"/output.dat")
      call parlst_readfromfile (rparamList, trim(DIR_DATA)//"/paramtriang.dat")
      call parlst_readfromfile (rparamList, trim(DIR_DATA)//"/postprocessing.dat")
      call parlst_readfromfile (rparamList, trim(DIR_DATA)//"/spacetimesolver.dat")
      call parlst_readfromfile (rparamList, trim(DIR_DATA)//"/timediscr.dat")
    end if

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine main_printSolverStatistics (rnewtonStatistics,&
      dtimeInit,dtimeStartVec,dtimePostproc,dtimeTotal)

!<description>
  ! Prints out statistics about the solution process.
!</description>

!<input>
  ! Statistics of the Newton solver
  type(t_newtonitSolverStat), intent(in) :: rnewtonStatistics

  ! Time for initialisation
  real(DP), intent(in) :: dtimeInit

  ! Time for creation of the start vector
  real(DP), intent(in) :: dtimeStartVec

  ! Time for final postprocessing
  real(DP), intent(in) :: dtimePostproc

  ! Total computation time
  real(DP), intent(in) :: dtimeTotal
!</input>

!</subroutine>

    ! Global time
    call output_line ("Total time                                   = "//&
        sys_sdL(dtimeTotal,10))
    call output_lbrk()

    call output_line ("Time for initialisation                      = "//&
        sys_sdL(dtimeInit,10))
    call output_line ("Time for creation of the start vector        = "//&
        sys_sdL(dtimeStartVec,10))
    call output_line ("Time for postprocessing                      = "//&
        sys_sdL(dtimePostproc+rnewtonStatistics%rtimePostprocessing%delapsedReal,10))
    call output_line ("#Iterations nonlinear space-time solver      = "//&
        sys_siL(rnewtonStatistics%niterations,10))
    call output_line ("Time for nonlinear space-time solver         = "//&
        sys_sdL(rnewtonStatistics%rtotalTime%delapsedReal &
               -rnewtonStatistics%rtimePostprocessing%delapsedReal,10))
    call output_lbrk()

    ! More detailed statistics of the nonlinear solver
    call output_line ("  Time for defect calculations               = "//&
        sys_sdL(rnewtonStatistics%rtimeDefect%delapsedReal,10))
    call output_line ("    Time for forward simulations             = "//&
        sys_sdL(rnewtonStatistics%rtimeForward%delapsedReal,10))
    call output_line ("    Time for backward simulations            = "//&
        sys_sdL(rnewtonStatistics%rtimeBackward%delapsedReal,10))
    call output_line ("    #Iterations nonlinear solver in space    = "//&
        sys_siL(rnewtonStatistics%rspaceslSolverStat%niterations,10))
    call output_line ("    Time nonlinear solver in space           = "//&
        sys_sdL(rnewtonStatistics%rspaceslSolverStat%rtotalTime%delapsedReal,10))
    call output_line ("      Time for nonlinear defects in space    = "//&
        sys_sdL(rnewtonStatistics%rspaceslSolverStat%rtimeDefect%delapsedReal,10))
    call output_line ("      Time for matrix assembly in space      = "//&
        sys_sdL(rnewtonStatistics%rspaceslSolverStat%rtimeMatrixAssembly%delapsedReal,10))
    call output_line ("      #Iterations linear solver in space     = "//&
        sys_siL(rnewtonStatistics%rspaceslSolverStat%rlssSolverStat%niterations,10))
    call output_line ("      Time linear solver in space            = "//&
        sys_sdL(rnewtonStatistics%rspaceslSolverStat%rlssSolverStat%&
                rtotalTime%delapsedReal,10))
    call output_line ("  Time for interpolation to lower levels     = "//&
        sys_sdL(rnewtonStatistics%rtimeProlRest%delapsedReal,10))
    call output_lbrk()

    ! More detailed statistics about the linear space-time solver
    call output_line ("  #Iterations linear space-time solver       = "//&
        sys_siL(rnewtonStatistics%rnewtonlinSolverStat%niterations,10))
    call output_line ("  Time for linear space-time solver          = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtotalTime%delapsedReal,10))

    call output_line ("    Time for smoothing                       = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtimeSmoothing%delapsedReal,10))
    call output_line ("      Time for smoothing (fine)              = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtimeSmoothingFine%delapsedReal,10))
    call output_line ("    Time for prolongation/restriction        = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtimeProlRest%delapsedReal,10))
    call output_line ("    Time for coarse grid solving             = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtimeSolverCoarse%delapsedReal,10))
    call output_lbrk()

    call output_line ("    Time for defect calculations             = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtimeDefect%delapsedReal,10))
    call output_line ("      Time for forward simulations           = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtimeForwardLin%delapsedReal,10))
    call output_line ("      Time for backward simulations          = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtimeBackwardLin%delapsedReal,10))
    call output_line ("      #Iterations linear solver in space     = "//&
        sys_siL(rnewtonStatistics%rnewtonlinSolverStat%rspaceslSolverStat%&
                rlssSolverStat%niterations,10))
    call output_line ("      Time linear solver in space            = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rspaceslSolverStat%&
                rtotalTime%delapsedReal,10))
    call output_lbrk()

    call output_line ("    Time for defect calculations (fine)      = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtimeDefectFine%delapsedReal,10))
    call output_line ("      Time for forward simulations (fine)    = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtimeForwardLinFine%delapsedReal,10))
    call output_line ("      Time for backward simulations (fine)   = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtimeBackwardLinFine%delapsedReal,10))
    call output_line ("      #Iterations linear solver (sp., fine)  = "//&
        sys_siL(rnewtonStatistics%rnewtonlinSolverStat%rspaceslSolverStatFine%&
                rlssSolverStat%niterations,10))
    call output_line ("      Time linear solver in space (fine)     = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rspaceslSolverStatFine%&
                rtotalTime%delapsedReal,10))
    call output_lbrk()

    call output_line ("    Time for defect calculations (coarse)    = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtimeDefectCoarse%delapsedReal,10))
    call output_line ("      Time for forward simulations (coarse)  = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtimeForwardLinCoarse%delapsedReal,10))
    call output_line ("      Time for backward simulations (coarse) = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rtimeBackwardLinCoarse%delapsedReal,10))
    call output_line ("      #Iterations linear solver (sp, coarse) = "//&
        sys_siL(rnewtonStatistics%rnewtonlinSolverStat%rspaceslSolverStatCoarse%&
                rlssSolverStat%niterations,10))
    call output_line ("      Time linear solver in space (coarse)   = "//&
        sys_sdL(rnewtonStatistics%rnewtonlinSolverStat%rspaceslSolverStatCoarse%rtotalTime%delapsedReal,10))

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine cc2doptccalculate (rsettings,rparlist)

!<description>
  ! This is a "separated" Navier-Stokes solver for solving a Navier-Stokes
  ! problem. The different tasks of the problem are separated into
  ! subroutines. The problem uses a problem-specific structure for the
  ! communication: All subroutines add their generated information to the
  ! structure, so that the other subroutines can work with them.
  ! (This is somehow a cleaner implementation than using a collection!).
  ! For the communication to callback routines of black-box subroutines
  ! (matrix-assembly), a collection is used.
  !
  ! The following tasks are performed by the subroutines:
  !
  ! 1.) Read in parametrisation
  ! 2.) Read in triangulation
  ! 3.) Set up RHS
  ! 4.) Set up matrix
  ! 5.) Create solver structure
  ! 6.) Solve the problem
  ! 7.) Write solution to GMV file
  ! 8.) Release all variables, finish
!</description>

!<input>
  ! Basic program settings
  type(t_settings_main), intent(in) :: rsettings

  ! Parameter list containing the DAT file parameters
  type(t_parlist), intent(in) :: rparlist
!</input>

!</subroutine>

    ! Local variables
    !
    ! The following structure configures our solver.
    ! We allocate this dynamically since the structure may be too large
    ! for the stack.
    type(t_settings_optflow), pointer :: p_rsettingsSolver

    ! Structure for the Newton solver
    type(t_spacetimeNewton) :: rsolver

    ! Structure encapsuling a hierarchy of solutions
    type(t_kktsystemHierarchy) :: rkktsystemHierarchy

    ! Structure holding the solution on the maximum level
    type(t_kktSystem) :: rsolution

    ! Some timers
    type(t_timer) :: rtotalTime,rinitTime,rsolverTime,rtimePostProc,rstartvectime
    type(t_newtonitSolverStat) :: rnewtonStatistics

    ! Ok, let us start.
    !
    ! Initialise the external storage management.

    call exstor_init (999,100)
    !CALL exstor_attachDirectory("./ff2storage")

    ! Ok, parameters are read in.
    ! Print the parameters to the terminal.
    if (rsettings%routput%ioutputInit .ge. 2) then
      call parlst_info (rparlist)
      call output_separator (OU_SEP_EQUAL)
    end if

    ! Measure the total time.
    call stat_clearTimer (rtotalTime)
    call stat_startTimer (rtotalTime)

!    ! Initialise the settings of the solver,
!    ! allocate all template matrices etc.
    allocate(p_rsettingsSolver)
!    allocate(p_rnlstsolver)

    call stat_clearTimer (rinitTime)
    call stat_startTimer (rinitTime)

    if (rsettings%routput%ioutputInit .ge. 1) then
      call output_lbrk()
      call output_line ("Initialising solver structures.")
    end if
    call init_initOptFlow (rparlist,rsettings,p_rsettingsSolver,rsettings%routput%ioutputInit)

    ! Create a solution
    call kkth_initHierarchy (rkktsystemHierarchy,&
        p_rsettingsSolver%roperatorAsmHier,&
        p_rsettingsSolver%rspaceTimeHierPrimal,&
        p_rsettingsSolver%rspaceTimeHierDual,&
        p_rsettingsSolver%rspaceTimeHierControl,.false.,&
        p_rsettingsSolver%roptcBDC)

    ! Allocate a solution
    call kkth_initKKTSystem (rsolution,rkktsystemHierarchy,&
        rkktsystemHierarchy%nlevels,p_rsettingsSolver%roperatorAsmHier,&
        p_rsettingsSolver%roptcBDC)

    if (rsettings%routput%ioutputInit .ge. 1) then
      call output_lbrk()
      call output_line ("Initialising initial iterate.")
    end if

    ! Initialise the initial iterate
    call init_initInitialIterate (rparlist,rsettings%ssectionOptControl,&
        p_rsettingsSolver%rphysics,p_rsettingsSolver%rsettingsSpaceDiscr,&
        p_rsettingsSolver%rtriaCoarse,p_rsettingsSolver%rrefinementSpace,&
        p_rsettingsSolver%rfeHierarchyControl,p_rsettingsSolver%rtimeHierarchy,&
        p_rsettingsSolver%rboundary,rsolution%p_rcontrol%p_rvector)
    call sptivec_invalidateVecInPool (rsolution%p_rcontrol%p_rvectorAccess)

    if (rsettings%routput%ioutputInit .ge. 1) then
      call output_lbrk()
      call output_line ("Initialising the descent algorithm loop.")
    end if
    call newtonit_init (rsolver,p_rsettingsSolver,rsettings%ssectionSpaceTimeSolver,rparlist)
    call newtonit_initStructure (rsolver,rkktsystemHierarchy,&
        p_rsettingsSolver%rprjHierSpaceTimePrimal,p_rsettingsSolver%rprjHierSpaceTimeDual,&
        p_rsettingsSolver%rprjHierSpaceTimeControl)

!    ! Create a temp vector
!    call sptivec_initVector (rtemp,&
!        p_rsettingsSolver%rtimeHierarchy%p_rtimeLevels(p_rsettingsSolver%rtimeHierarchy%nlevels),&
!        p_rsettingsSolver%rfeHierPrimalDual% &
!        p_rfeSpaces(p_rsettingsSolver%rfeHierPrimalDual%nlevels)%p_rdiscretisation)
!
    call stat_stopTimer (rinitTime)

    call output_separator (OU_SEP_EQUAL)
    call output_line ("Time for initialisation            = "//&
        sys_sdL(rinitTime%delapsedReal,10))
    call output_separator (OU_SEP_EQUAL)

    ! Solve the system
    call stat_clearTimer (rsolverTime)
    call stat_startTimer (rsolverTime)
    call newtonit_solve (rsolver,rsolution,rnewtonStatistics)
    call stat_stopTimer (rsolverTime)

    call output_separator (OU_SEP_EQUAL)

    ! Pipe the solution through our postprocessing routines
    call output_line ("Postprocessing of the final solution...")
    call output_lbrk()

    call stat_clearTimer (rtimePostProc)
    call stat_startTimer (rtimePostProc)
    call optcpp_postprocessing (p_rsettingsSolver%rpostproc,&
        p_rsettingsSolver%rspaceTimeHierPrimal%p_rfeHierarchy%nlevels,&
        p_rsettingsSolver%rspaceTimeHierPrimal%p_rtimeHierarchy%nlevels,&
        rsolution,rsolver%rkktSubsolvers,p_rsettingsSolver)
!    call optcpp_postprocessSpaceTimeVec (rpostproc,rsolution,rrhsDiscrete,&
!        p_rsettingsSolver%rsettingsOptControl,p_rsettingsSolver)
    call stat_stopTimer (rtimePostProc)

!    ! Release all data
    call newtonit_doneStructure (rsolver)
    call newtonit_done (rsolver)

    ! Release the solution
    call kkth_doneHierarchy (rkktSystemHierarchy)
    call kkt_doneKKTsystem (rsolution)

    ! Release solver structures
    call init_doneOptFlow (p_rsettingsSolver)

    call stat_stopTimer (rtotalTime)

    call output_separator (OU_SEP_EQUAL)
    call main_printSolverStatistics (&
        rnewtonStatistics,rinitTime%delapsedReal,rstartvectime%delapsedReal,&
        rtimePostProc%delapsedReal,rtotalTime%delapsedReal)

    call output_separator (OU_SEP_EQUAL)

    ! Information about external storage usage
    call output_lbrk ()
    call exstor_info (bprintHandles=.true.)

    ! Clean up the external storage management
    call exstor_done ()

  end subroutine

  ! ***************************************************************************

  subroutine cc2doptcmain

    ! Program parameters
    type(t_parlist) :: rparlist
    type(t_settings_main) :: rsettings

    ! The very first thing in every application:
    ! Initialise system-wide settings:
    call sys_init()

    ! Read the program parameters.
    call parlst_init (rparlist)
    call main_getDat (rparlist)

    ! Read the basic parameter settings from the "MAIN" section.
    call smain_initMainParams (rsettings,rparlist,"MAIN")

    ! Initialise log file for output.
    call output_init (rsettings%routput%smsgLog,rsettings%routput%serrorLog)
    OU_LINE_LENGTH = 132
    cdefaultDateTimeLogPolicy = OU_DTP_ADDDATETIME
    cdatetimeLogFormat = 1

    ! Now we can really start!
    !
    ! Initialise the storage management:
    call storage_init(999, 100)

    ! Initialise the parser
    call fparser_init ()

    ! Call the problem to solve.
    call output_lbrk ()
    call output_line ("Calculating cc2doptc-Problem")
    call output_separator (OU_SEP_MINUS)

    call cc2doptccalculate (rsettings,rparlist)

    ! Release the parser
    call fparser_done ()

    ! Print out heap statistics - just to check if everything
    ! is cleaned up.
    ! This should display "Handles in use=0" and "Memory in use=0"!
    call output_lbrk ()
    call storage_info(.true.)

    ! Clean up the storage management, parameter list, finish
    call storage_done()
    call parlst_done (rparlist)

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
Windmere Apartments is a rental housing community developed and managed by Community Housing Opportunities Corporation CHOC for very low, low and moderate income households. Sterling Asset Management Company also provides a resident services program which provides residents with professional staff to assist them in enhancing the quality of life within their residential communities and achieving individual personal goals.

Windmere is located on land dedicated by the City of Davis within Mace Ranch and offers twobedroom apartments and townhomes. Townhouse units feature ninefoot ceilings downstairs, and all units provide a private patio/balcony and exterior storage units. Computer facilities are available for resident use in the community room. All of the apartments have a refrigerator, electric stove, garbage disposal, dishwasher, central heating and air conditioning.  There are two playgrounds in the complex, plus a basketball halfcourt and a foursquare court.  There are no washer and drier hookups in the units, but multiple laundry rooms are conveniently distributed throughout the complex.  No pets are permitted (although birds, fish, and small caged critters seem to get overlooked).

In particular Windmere caters to lowincome families and people with disabilities.  Its an appealing option for lowincome parents whod rather live where their young kids will have other kids next door to play with, instead of being the only children in a complex designed for college students.  Many of the families in Windmere have resided there for years, and there is a strong positive feeling of community among the residents.

To learn more about rental housing in Davis, check out our Housing Guide as well as our Apartments pages.

20051220 13:21:54 nbsp I have three kids under 10; Ive lived in Windmere II since it opened in 1999, and have been very happy with the complex.  Its hard for young/single parents to find affordable housing in Davis, and Windmere has been a great place for me and my family.  I highly recommend it. Users/ErinManley

20110120 13:57:02 nbsp A friend of mine lived in a 2bd 1 bath, and it was BIG with a separate storage unit. She also paid like $650...nice huh? Users/ucdsunshine

20110328 13:49:25 nbsp Affordable and spacious, but WAY too many kids!!! oodles and oodles. Also, some of the residents seem kind of sketchy. This one unit has about 3545 bikes in the patio/stairway. All sizes/brands...looks sketchy to me, especially when ive seen new ones everytime I visit my sister. Users/cora