text:
program vin
  integer :: val2(2,5), i, j, val3(2,5,10)

  do i = 1,5
    do j = 1,2
      val2(j,i) = mod(i * j, i + j)
    end do
  end do

  print *, minloc(val2,mask=val2>2)
  print *, minloc(val2,mask=val2<3)
  print *, minloc(val2(:,1) + val2(:,2) ,mask=val2(:,1)<3)

  print *, maxloc(val2 + val2 ,mask=val2>2)
  print *, maxloc(val2,mask=val2<3)
  print *, maxloc(val2(:,1),mask=val2(:,1)<3)

end program vin

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
*> \brief \b DLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLARRJ + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST,
*                          RTOL, OFFSET, W, WERR, WORK, IWORK,
*                          PIVMIN, SPDIAM, INFO )
*
*       .. Scalar Arguments ..
*       INTEGER            IFIRST, ILAST, INFO, N, OFFSET
*       DOUBLE PRECISION   PIVMIN, RTOL, SPDIAM
*       ..
*       .. Array Arguments ..
*       INTEGER            IWORK( * )
*       DOUBLE PRECISION   D( * ), E2( * ), W( * ),
*      $                   WERR( * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> Given the initial eigenvalue approximations of T, DLARRJ
*> does  bisection to refine the eigenvalues of T,
*> W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
*> guesses for these eigenvalues are input in W, the corresponding estimate
*> of the error in these guesses in WERR. During bisection, intervals
*> [left, right] are maintained by storing their mid-points and
*> semi-widths in the arrays W and WERR respectively.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix.
*> \endverbatim
*>
*> \param[in] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (N)
*>          The N diagonal elements of T.
*> \endverbatim
*>
*> \param[in] E2
*> \verbatim
*>          E2 is DOUBLE PRECISION array, dimension (N-1)
*>          The Squares of the (N-1) subdiagonal elements of T.
*> \endverbatim
*>
*> \param[in] IFIRST
*> \verbatim
*>          IFIRST is INTEGER
*>          The index of the first eigenvalue to be computed.
*> \endverbatim
*>
*> \param[in] ILAST
*> \verbatim
*>          ILAST is INTEGER
*>          The index of the last eigenvalue to be computed.
*> \endverbatim
*>
*> \param[in] RTOL
*> \verbatim
*>          RTOL is DOUBLE PRECISION
*>          Tolerance for the convergence of the bisection intervals.
*>          An interval [LEFT,RIGHT] has converged if
*>          RIGHT-LEFT.LT.RTOL*MAX(|LEFT|,|RIGHT|).
*> \endverbatim
*>
*> \param[in] OFFSET
*> \verbatim
*>          OFFSET is INTEGER
*>          Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET
*>          through ILAST-OFFSET elements of these arrays are to be used.
*> \endverbatim
*>
*> \param[in,out] W
*> \verbatim
*>          W is DOUBLE PRECISION array, dimension (N)
*>          On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
*>          estimates of the eigenvalues of L D L^T indexed IFIRST through
*>          ILAST.
*>          On output, these estimates are refined.
*> \endverbatim
*>
*> \param[in,out] WERR
*> \verbatim
*>          WERR is DOUBLE PRECISION array, dimension (N)
*>          On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
*>          the errors in the estimates of the corresponding elements in W.
*>          On output, these errors are refined.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (2*N)
*>          Workspace.
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (2*N)
*>          Workspace.
*> \endverbatim
*>
*> \param[in] PIVMIN
*> \verbatim
*>          PIVMIN is DOUBLE PRECISION
*>          The minimum pivot in the Sturm sequence for T.
*> \endverbatim
*>
*> \param[in] SPDIAM
*> \verbatim
*>          SPDIAM is DOUBLE PRECISION
*>          The spectral diameter of T.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          Error flag.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date September 2012
*
*> \ingroup auxOTHERauxiliary
*
*> \par Contributors:
*  ==================
*>
*> Beresford Parlett, University of California, Berkeley, USA \n
*> Jim Demmel, University of California, Berkeley, USA \n
*> Inderjit Dhillon, University of Texas, Austin, USA \n
*> Osni Marques, LBNL/NERSC, USA \n
*> Christof Voemel, University of California, Berkeley, USA
*
*  =====================================================================
      SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST,
     $                   RTOL, OFFSET, W, WERR, WORK, IWORK,
     $                   PIVMIN, SPDIAM, INFO )
*
*  -- LAPACK auxiliary routine (version 3.4.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     September 2012
*
*     .. Scalar Arguments ..
      INTEGER            IFIRST, ILAST, INFO, N, OFFSET
      DOUBLE PRECISION   PIVMIN, RTOL, SPDIAM
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   D( * ), E2( * ), W( * ),
     $                   WERR( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TWO, HALF
      PARAMETER        ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
     $                   HALF = 0.5D0 )
      INTEGER   MAXITR
*     ..
*     .. Local Scalars ..
      INTEGER            CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT,
     $                   OLNINT, P, PREV, SAVI1
      DOUBLE PRECISION   DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH
*
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX
*     ..
*     .. Executable Statements ..
*
      INFO = 0
*
      MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
     $           LOG( TWO ) ) + 2
*
*     Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
*     The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
*     Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
*     for an unconverged interval is set to the index of the next unconverged
*     interval, and is -1 or 0 for a converged interval. Thus a linked
*     list of unconverged intervals is set up.
*

      I1 = IFIRST
      I2 = ILAST
*     The number of unconverged intervals
      NINT = 0
*     The last unconverged interval found
      PREV = 0
      DO 75 I = I1, I2
         K = 2*I
         II = I - OFFSET
         LEFT = W( II ) - WERR( II )
         MID = W(II)
         RIGHT = W( II ) + WERR( II )
         WIDTH = RIGHT - MID
         TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )

*        The following test prevents the test of converged intervals
         IF( WIDTH.LT.RTOL*TMP ) THEN
*           This interval has already converged and does not need refinement.
*           (Note that the gaps might change through refining the
*            eigenvalues, however, they can only get bigger.)
*           Remove it from the list.
            IWORK( K-1 ) = -1
*           Make sure that I1 always points to the first unconverged interval
            IF((I.EQ.I1).AND.(I.LT.I2)) I1 = I + 1
            IF((PREV.GE.I1).AND.(I.LE.I2)) IWORK( 2*PREV-1 ) = I + 1
         ELSE
*           unconverged interval found
            PREV = I
*           Make sure that [LEFT,RIGHT] contains the desired eigenvalue
*
*           Do while( CNT(LEFT).GT.I-1 )
*
            FAC = ONE
 20         CONTINUE
            CNT = 0
            S = LEFT
            DPLUS = D( 1 ) - S
            IF( DPLUS.LT.ZERO ) CNT = CNT + 1
            DO 30 J = 2, N
               DPLUS = D( J ) - S - E2( J-1 )/DPLUS
               IF( DPLUS.LT.ZERO ) CNT = CNT + 1
 30         CONTINUE
            IF( CNT.GT.I-1 ) THEN
               LEFT = LEFT - WERR( II )*FAC
               FAC = TWO*FAC
               GO TO 20
            END IF
*
*           Do while( CNT(RIGHT).LT.I )
*
            FAC = ONE
 50         CONTINUE
            CNT = 0
            S = RIGHT
            DPLUS = D( 1 ) - S
            IF( DPLUS.LT.ZERO ) CNT = CNT + 1
            DO 60 J = 2, N
               DPLUS = D( J ) - S - E2( J-1 )/DPLUS
               IF( DPLUS.LT.ZERO ) CNT = CNT + 1
 60         CONTINUE
            IF( CNT.LT.I ) THEN
               RIGHT = RIGHT + WERR( II )*FAC
               FAC = TWO*FAC
               GO TO 50
            END IF
            NINT = NINT + 1
            IWORK( K-1 ) = I + 1
            IWORK( K ) = CNT
         END IF
         WORK( K-1 ) = LEFT
         WORK( K ) = RIGHT
 75   CONTINUE

      SAVI1 = I1
*
*     Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
*     and while (ITER.LT.MAXITR)
*
      ITER = 0
 80   CONTINUE
      PREV = I1 - 1
      I = I1
      OLNINT = NINT

      DO 100 P = 1, OLNINT
         K = 2*I
         II = I - OFFSET
         NEXT = IWORK( K-1 )
         LEFT = WORK( K-1 )
         RIGHT = WORK( K )
         MID = HALF*( LEFT + RIGHT )

*        semiwidth of interval
         WIDTH = RIGHT - MID
         TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )

         IF( ( WIDTH.LT.RTOL*TMP ) .OR.
     $      (ITER.EQ.MAXITR) )THEN
*           reduce number of unconverged intervals
            NINT = NINT - 1
*           Mark interval as converged.
            IWORK( K-1 ) = 0
            IF( I1.EQ.I ) THEN
               I1 = NEXT
            ELSE
*              Prev holds the last unconverged interval previously examined
               IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
            END IF
            I = NEXT
            GO TO 100
         END IF
         PREV = I
*
*        Perform one bisection step
*
         CNT = 0
         S = MID
         DPLUS = D( 1 ) - S
         IF( DPLUS.LT.ZERO ) CNT = CNT + 1
         DO 90 J = 2, N
            DPLUS = D( J ) - S - E2( J-1 )/DPLUS
            IF( DPLUS.LT.ZERO ) CNT = CNT + 1
 90      CONTINUE
         IF( CNT.LE.I-1 ) THEN
            WORK( K-1 ) = MID
         ELSE
            WORK( K ) = MID
         END IF
         I = NEXT

 100  CONTINUE
      ITER = ITER + 1
*     do another loop if there are still unconverged intervals
*     However, in the last iteration, all intervals are accepted
*     since this is the best we can do.
      IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
*
*
*     At this point, all the intervals have converged
      DO 110 I = SAVI1, ILAST
         K = 2*I
         II = I - OFFSET
*        All intervals marked by '0' have been refined.
         IF( IWORK( K-1 ).EQ.0 ) THEN
            W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
            WERR( II ) = WORK( K ) - W( II )
         END IF
 110  CONTINUE
*

      RETURN
*
*     End of DLARRJ
*
      END

ext: f

max_stars_repo_name: zerothi/lapack

max_issues_repo_name: zerothi/lapack

max_forks_repo_name: zerothi/lapack


text:
      double precision function DCSHMM (X)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c>> 1994-10-20 DCSHMM Krogh  Changes to use M77CON
c>> 1993-05-07 DCSHMM WVSnyder JPL Initial code
c
c     Compute cosh(x) - 1 - x**2 using a rational approximation when
c     abs(x) is less than 2.7, else use the Fortran intrinsic function.
c
c--D replaces "?": ?CSHMM
      double precision ZP3, ZP2, ZP1, ZQ4, ZQ3, ZQ2, ZQ1, X, XS
      parameter ( ZP3 = 5.59297116264720D-07 )
      parameter ( ZP2 = 1.77943488030894D-04 )
      parameter ( ZP1 = 1.69800461894792D-02 )
      parameter ( ZQ4 = 1.33412535492375D-09 * 24.0d0 )
      parameter ( ZQ3 = -5.80858944138663D-07 * 24.0d0 )
      parameter ( ZQ2 = 1.27814964403863D-04 * 24.0d0 )
      parameter ( ZQ1 = -1.63532871439181D-02 * 24.0d0 )
c     DATA ZP3/5.59297116264720D-07/,
c    *     ZP2/1.77943488030894D-04/,
c    *     ZP1/1.69800461894792D-02/,
c    *     ZQ4/1.33412535492375D-09/,
c    *     ZQ3/-5.80858944138663D-07/,
c    *     ZQ2/1.27814964403863D-04/,
c    *     ZQ1/-1.63532871439181D-02/
      xs = x * x
      if (xs .lt. 7.29d0) then
         dcshmm = ((((zp3*xs+zp2)*xs+zp1)*xs+1.0d0)*xs*xs)/
     *            ((((zq4*xs+zq3)*xs+zq2)*xs+zq1)*xs+24.0d0)
      else
         dcshmm = cosh(x) - 1.0d0 - 0.5d0*xs
      end if
      return
      end


text:
c
c ---------------------------------------------------------
c
      subroutine cleanup(nvar,naux)
c
c :::::::::::::::::::::: CLEANUP ::::::::::::::::::::::::::::::::;
!> this is just a check to make sure all storage was accounted for.
!! routine is called after all the data has been checkpointed.
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::;

      use amr_module
      implicit double precision (a-h,o-z)
c
c      ## clean up storage to double check that everything taken care of
c      ## done after the checkpoint so pointers sitll work on restart
       do  120 level = 1, lfine
         call putsp(1,level,nvar,naux)
         mptr =  lstart(level)
 110        nx      = node(ndihi,mptr) - node(ndilo,mptr) + 1
            ny      = node(ndjhi,mptr) - node(ndjlo,mptr) + 1
            mitot  = nx + 2*nghost
            mjtot  = ny + 2*nghost
            nwords  = mitot*mjtot*nvar
            call reclam(node(store1, mptr), nwords)
            if (level .lt. mxnest)
     .         call reclam(node(store2, mptr), nwords)
            if (naux .gt. 0)
     .         call reclam(node(storeaux, mptr), mitot*mjtot*naux)
        mptr = node(levelptr, mptr)
        if (mptr .ne. 0) go to 110
120    continue

      return
      end

max_stars_repo_name: navravi/amrclaw

max_issues_repo_name: navravi/amrclaw

max_forks_repo_name: BrisaDavis/amrclaw


text:
!=======================================================================
! This is part of the 2DECOMP&FFT library
!
! 2DECOMP&FFT is a software framework for general-purpose 2D (pencil)
! decomposition. It also implements a highly scalable distributed
! three-dimensional Fast Fourier Transform (FFT).
!
! Copyright (C) 2009-2012 Ning Li, the Numerical Algorithms Group (NAG)
!
!=======================================================================

! This file contains the routines that transpose data from Y to Z pencil

#if defined(_OPENACC)
  subroutine transpose_yt_to_zp(src, dst, yt, zp, normfft, opt_decomp)

    implicit none

    real(mytype), dimension(:,:,:), intent(INOUT) :: src
    real(mytype), dimension(:,:,:), intent(INOUT) :: dst
    real(mytype), dimension(:,:,:), intent(INOUT) :: yt
    real(mytype), dimension(0:,0:,0:), intent(INOUT) :: zp
    real(mytype), intent(IN) :: normfft

    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
    attributes( managed ) :: src, dst
    attributes( device  ) :: yt
    attributes( managed ) :: zp
    TYPE(DECOMP_INFO) :: decomp

    integer :: s1,s2,s3,d1,d2,d3
    integer :: ierror, istat, m, i1, i2, pos
    integer :: iter, dest, sorc, pow2
    integer :: i,j,k,kstart,kend,jstart,jend,decomp_y2distm, idx
    if (present(opt_decomp)) then
       decomp = opt_decomp
    else
       decomp = decomp_main
    end if

    s1 = SIZE(src,1)
    s2 = SIZE(src,2)
    s3 = SIZE(src,3)
    d1 = SIZE(dst,1)
    d2 = SIZE(dst,2)
    d3 = SIZE(dst,3)

#if defined(_EPA2A)

#if !defined(_EPHC)
    !$cuf kernel do(3) <<<*,*>>>
    do k=1,s3
    do i=1,s2
    do j=1,s1
      src(j,i,k) = yt(i,j,k)
    enddo
    enddo
    enddo
#endif

    if(IAND(dims(2),dims(2)-1)==0) then
      pow2 = 1
    else
      pow2 = 0
    endif

    ! rearrange source array as send buffer
    do iter=1,dims(2)-1
       if( pow2 ) then
         dest = IEOR(row_rank,iter)
       else
         dest = mod(row_rank + iter, dims(2))
       endif
       m = dest
       pos = decomp%y2disp(m) + 1
#if !defined(_EPHC)
       istat = cudaMemcpy2DAsync( work1_r_d(pos), s1*(decomp%y2dist(m)), src(1,decomp%y2idx(m),1), s1*s2, s1*(decomp%y2dist(m)), s3, stream=a2a_d2h )
#else
       jstart = decomp%y2idx(m)
       jend   = jstart + decomp%y2dist(m) - 1
       decomp_y2distm = decomp%y2dist(m)
       !$cuf kernel do(3) <<<*,*,stream=a2a_d2h>>>
       do k=1,s3
       do j=jstart, jend
       do i=1,s1
         idx = i + s1*(j-jstart) + s1*decomp_y2distm*(k-1)
         work1_r_d(pos + idx - 1) = yt( j, i, k )
       enddo
       enddo
       enddo
#endif
       istat = cudaEventRecord( a2a_event(iter), a2a_d2h )
    end do

    ! self
    m = row_rank
#if !defined(_EPHC)
    istat = cudaMemcpy2DAsync( dst(1,1,decomp%z2idx(m)), s1*(decomp%y2dist(m)), src(1,decomp%y2idx(m),1), s1*s2, s1*(decomp%y2dist(m)), s3, stream=a2a_comp )
#else
    kstart = decomp%z2idx(m)
    jstart = decomp%y2idx(m)
    !$cuf kernel do(3) <<<*,(8,8,8),stream=a2a_comp>>>
    do k=1,s3
    do j=1,d2
    do i=1,s1
      zp(i,j,k+kstart-1) = yt(j+jstart-1,i,k)*normfft
    enddo
    enddo
    enddo
#endif

    do iter=1,dims(2)-1
      if( pow2 ) then
        sorc = IEOR(row_rank,iter)
      else
        sorc = mod(row_rank - iter + dims(2), dims(2))
      endif
      m = sorc
      call MPI_IRECV( work2_r_d(decomp%z2disp(m)+1), decomp%z2cnts(m), real_type, m, 0, DECOMP_2D_COMM_ROW, a2a_requests(iter), ierror)
    end do

    do iter=1,dims(2)-1
       if( pow2 ) then
          dest = IEOR(row_rank,iter)
          sorc = dest
       else
          dest = mod(row_rank + iter, dims(2))
          sorc = mod(row_rank - iter + dims(2), dims(2))
       endif
       m = dest
       istat = cudaEventSynchronize( a2a_event(iter) )
       call MPI_SEND( work1_r_d(decomp%y2disp(m)+1), decomp%y2cnts(m), real_type, m, 0, DECOMP_2D_COMM_ROW, ierror)
       call MPI_WAIT(a2a_requests(iter), MPI_STATUS_IGNORE, ierror)
       m = sorc
       pos = decomp%z2disp(m) + 1
#if !defined(_EPHC)
       istat = cudaMemcpyAsync( dst(1,1,decomp%z2idx(m)), work2_r_d(pos), decomp%z2cnts(m), a2a_h2d )
#else
       kstart = decomp%z2idx(m)
       kend = kstart + decomp%z2dist(m) - 1

       !$cuf kernel do(3) <<<*,*,stream=a2a_h2d>>>
       do k=kstart, kend
       do j=1,d2
       do i=1,d1
         idx = i + d1*(j-1) + d1*d2*(k-1)
         zp( i, j, k ) = work2_r_d( idx )*normfft
       enddo
       enddo
       enddo
#endif
    end do

#if !defined(_EPHC)
    !$cuf kernel do(3) <<<*,*>>>
    do k=1,d3
    do j=1,d2
    do i=1,d1
      zp(i,j,k) = dst(i,j,k)*normfft
    enddo
    enddo
    enddo
#endif

    istat = cudaEventRecord( a2a_event(0), 0 )
    istat = cudaEventSynchronize( a2a_event(0) )
#else
    ! rearrange source array as send buffer
    do m=0,dims(2)-1
       if (m==0) then
          i1 = 1
          i2 = decomp%y2dist(0)
       else
          i1 = i2+1
          i2 = i1+decomp%y2dist(m)-1
       end if
       pos = decomp%y2disp(m) + 1
       istat = cudaMemcpy2D( work1_r(pos), s1*(i2-i1+1), src(1,i1,1), s1*s2, s1*(i2-i1+1), s3, cudaMemcpyDeviceToHost )
    end do

    call MPI_ALLTOALLV(work1_r, decomp%y2cnts, decomp%y2disp, &
         real_type, work2_r, decomp%z2cnts, decomp%z2disp, &
         real_type, DECOMP_2D_COMM_ROW, ierror)

    istat = cudaMemcpy( dst, work2_r, d1*d2*d3, cudaMemcpyHostToDevice )
#endif

    return
  end subroutine transpose_yt_to_zp

  subroutine transpose_y_to_z_real_d(src, dst, opt_decomp)

    implicit none

    real(mytype), dimension(:,:,:), intent(IN) :: src
    real(mytype), dimension(:,:,:), intent(OUT) :: dst
    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp
    attributes( managed ) :: src, dst
    TYPE(DECOMP_INFO) :: decomp

    integer :: s1,s2,s3,d1,d2,d3
    integer :: ierror, istat, m, i1, i2, pos
    integer :: iter, dest, sorc, pow2
    if (present(opt_decomp)) then
       decomp = opt_decomp
    else
       decomp = decomp_main
    end if

    s1 = SIZE(src,1)
    s2 = SIZE(src,2)
    s3 = SIZE(src,3)
    d1 = SIZE(dst,1)
    d2 = SIZE(dst,2)
    d3 = SIZE(dst,3)

#if defined(_EPA2A)
    if(IAND(dims(2),dims(2)-1)==0) then
      pow2 = 1
    else
      pow2 = 0
    endif

    ! rearrange source array as send buffer
    do iter=1,dims(2)-1
       if( pow2 ) then
         dest = IEOR(row_rank,iter)
       else
         dest = mod(row_rank + iter, dims(2))
       endif
       m = dest
       pos = decomp%y2disp(m) + 1
       istat = cudaMemcpy2DAsync( work1_r_d(pos), s1*(decomp%y2dist(m)), src(1,decomp%y2idx(m),1), s1*s2, s1*(decomp%y2dist(m)), s3, stream=a2a_d2h )
       istat = cudaEventRecord( a2a_event(iter), a2a_d2h )
    end do

    ! self
    m = row_rank
    istat = cudaMemcpy2DAsync( dst(1,1,decomp%z2idx(m)), s1*(decomp%y2dist(m)), src(1,decomp%y2idx(m),1), s1*s2, s1*(decomp%y2dist(m)), s3, stream=a2a_comp )

    do iter=1,dims(2)-1
      if( pow2 ) then
        sorc = IEOR(row_rank,iter)
      else
        sorc = mod(row_rank - iter + dims(2), dims(2))
      endif
      m = sorc
      call MPI_IRECV( work2_r_d(decomp%z2disp(m)+1), decomp%z2cnts(m), real_type, m, 0, DECOMP_2D_COMM_ROW, a2a_requests(iter), ierror)
    end do

    do iter=1,dims(2)-1
       if( pow2 ) then
          dest = IEOR(row_rank,iter)
          sorc = dest
       else
          dest = mod(row_rank + iter, dims(2))
          sorc = mod(row_rank - iter + dims(2), dims(2))
       endif
       m = dest
       istat = cudaEventSynchronize( a2a_event(iter) )
       call MPI_SEND( work1_r_d(decomp%y2disp(m)+1), decomp%y2cnts(m), real_type, m, 0, DECOMP_2D_COMM_ROW, ierror)
       call MPI_WAIT(a2a_requests(iter), MPI_STATUS_IGNORE, ierror)
       m = sorc
       pos = decomp%z2disp(m) + 1
       istat = cudaMemcpyAsync( dst(1,1,decomp%z2idx(m)), work2_r_d(pos), decomp%z2cnts(m), a2a_h2d )
    end do
    istat = cudaEventRecord( a2a_event(0), 0 )
    istat = cudaEventSynchronize( a2a_event(0) )
#else
    ! rearrange source array as send buffer
    do m=0,dims(2)-1
       if (m==0) then
          i1 = 1
          i2 = decomp%y2dist(0)
       else
          i1 = i2+1
          i2 = i1+decomp%y2dist(m)-1
       end if
       pos = decomp%y2disp(m) + 1
       istat = cudaMemcpy2D( work1_r(pos), s1*(i2-i1+1), src(1,i1,1), s1*s2, s1*(i2-i1+1), s3, cudaMemcpyDeviceToHost )
    end do

    call MPI_ALLTOALLV(work1_r, decomp%y2cnts, decomp%y2disp, &
         real_type, work2_r, decomp%z2cnts, decomp%z2disp, &
         real_type, DECOMP_2D_COMM_ROW, ierror)

    istat = cudaMemcpy( dst, work2_r, d1*d2*d3, cudaMemcpyHostToDevice )
#endif

    return
  end subroutine transpose_y_to_z_real_d
#endif

  subroutine transpose_y_to_z_real(src, dst, opt_decomp)

    implicit none

    real(mytype), dimension(:,:,:), intent(IN) :: src
    real(mytype), dimension(:,:,:), intent(OUT) :: dst
    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp

    TYPE(DECOMP_INFO) :: decomp

#if defined(_SHM)
    real(mytype) :: work1(*), work2(*)
    POINTER  (work1_p, work1), (work2_p, work2)  ! Cray pointers
#endif

    integer :: s1,s2,s3,d1,d2,d3
    integer :: ierror

    if (present(opt_decomp)) then
       decomp = opt_decomp
    else
       decomp = decomp_main
    end if

    s1 = SIZE(src,1)
    s2 = SIZE(src,2)
    s3 = SIZE(src,3)
    d1 = SIZE(dst,1)
    d2 = SIZE(dst,2)
    d3 = SIZE(dst,3)

    ! rearrange source array as send buffer
#if defined(_SHM)
    work1_p = decomp%ROW_INFO%SND_P
    call mem_split_yz_real(src, s1, s2, s3, work1, dims(2), &
         decomp%y2dist, decomp)
#else
    call mem_split_yz_real(src, s1, s2, s3, work1_r, dims(2), &
         decomp%y2dist, decomp)
#endif

    ! define receive buffer
#if defined(_SHM)
    work2_p = decomp%ROW_INFO%RCV_P
    call MPI_BARRIER(decomp%ROW_INFO%CORE_COMM, ierror)
#endif

#if defined(_SHM)
    if (decomp%ROW_INFO%CORE_ME==1) THEN
       call MPI_ALLTOALLV(work1, decomp%y2cnts_s, decomp%y2disp_s, &
            real_type, work2, decomp%z2cnts_s, decomp%z2disp_s, &
            real_type, decomp%ROW_INFO%SMP_COMM, ierror)
    end if
#else
#if defined(_EVEN)
    if (decomp%even) then
       call MPI_ALLTOALL(work1_r, decomp%y2count, &
            real_type, dst, decomp%z2count, &
            real_type, DECOMP_2D_COMM_ROW, ierror)
    else
       call MPI_ALLTOALL(work1_r, decomp%y2count, &
            real_type, work2_r, decomp%z2count, &
            real_type, DECOMP_2D_COMM_ROW, ierror)
    end if
#else
    call MPI_ALLTOALLV(work1_r, decomp%y2cnts, decomp%y2disp, &
         real_type, dst, decomp%z2cnts, decomp%z2disp, &
         real_type, DECOMP_2D_COMM_ROW, ierror)
#endif
#endif

    ! rearrange receive buffer
#if defined(_SHM)
    call MPI_BARRIER(decomp%ROW_INFO%CORE_COMM, ierror)
    call mem_merge_yz_real(work2, d1, d2, d3, dst, dims(2), &
         decomp%z2dist, decomp)
#else
#if defined(_EVEN)
    if (.not. decomp%even) then
       call mem_merge_yz_real(work2_r, d1, d2, d3, dst, dims(2), &
         decomp%z2dist, decomp)
    end if
#else
    ! note the receive buffer is already in natural (i,j,k) order
    ! so no merge operation needed
#endif
#endif

    return
  end subroutine transpose_y_to_z_real

#if defined(_OCC)
  subroutine transpose_y_to_z_real_start(handle, src, dst, sbuf, rbuf, &
       opt_decomp)

    implicit none

    integer :: handle
    real(mytype), dimension(:,:,:) :: src, dst, sbuf, rbuf
    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp

    TYPE(DECOMP_INFO) :: decomp

    integer :: s1,s2,s3
    integer :: ierror

    if (present(opt_decomp)) then
       decomp = opt_decomp
    else
       decomp = decomp_main
    end if

    s1 = SIZE(src,1)
    s2 = SIZE(src,2)
    s3 = SIZE(src,3)

    ! rearrange source array as send buffer
    call mem_split_yz_real(src, s1, s2, s3, sbuf, dims(2), &
         decomp%y2dist, decomp)

#if defined(_EVEN)
    call NBC_IALLTOALL(sbuf, decomp%y2count, real_type, &
         rbuf, decomp%z2count, real_type, &
         DECOMP_2D_COMM_ROW, handle, ierror)
#else
    call NBC_IALLTOALLV(sbuf, decomp%y2cnts, decomp%y2disp, real_type, &
         rbuf, decomp%z2cnts, decomp%z2disp, real_type, &
         DECOMP_2D_COMM_ROW, handle, ierror)
#endif

    return
  end subroutine transpose_y_to_z_real_start

  subroutine transpose_y_to_z_real_wait(handle, src, dst, sbuf, rbuf, &
       opt_decomp)

    implicit none

    integer :: handle
    real(mytype), dimension(:,:,:) :: src, dst, sbuf, rbuf
    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp

    TYPE(DECOMP_INFO) :: decomp

    integer :: ierror

    if (present(opt_decomp)) then
       decomp = opt_decomp
    else
       decomp = decomp_main
    end if

    call NBC_WAIT(handle, ierror)

    dst = rbuf

    return
  end subroutine transpose_y_to_z_real_wait
#endif

  subroutine transpose_y_to_z_complex(src, dst, opt_decomp)

    implicit none

    complex(mytype), dimension(:,:,:), intent(IN) :: src
    complex(mytype), dimension(:,:,:), intent(OUT) :: dst
    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp

    TYPE(DECOMP_INFO) :: decomp

#if defined(_SHM)
    complex(mytype) :: work1(*), work2(*)
    POINTER  (work1_p, work1), (work2_p, work2)  ! Cray pointers
#endif

    integer :: s1,s2,s3,d1,d2,d3
    integer :: ierror

    if (present(opt_decomp)) then
       decomp = opt_decomp
    else
       decomp = decomp_main
    end if

    s1 = SIZE(src,1)
    s2 = SIZE(src,2)
    s3 = SIZE(src,3)
    d1 = SIZE(dst,1)
    d2 = SIZE(dst,2)
    d3 = SIZE(dst,3)

    ! rearrange source array as send buffer
#if defined(_SHM)
    work1_p = decomp%ROW_INFO%SND_P_c
    call mem_split_yz_complex(src, s1, s2, s3, work1, dims(2), &
         decomp%y2dist, decomp)
#else
    call mem_split_yz_complex(src, s1, s2, s3, work1_c, dims(2), &
         decomp%y2dist, decomp)
#endif

    ! define receive buffer
#if defined(_SHM)
    work2_p = decomp%ROW_INFO%RCV_P_c
    call MPI_BARRIER(decomp%ROW_INFO%CORE_COMM, ierror)
#endif

#if defined(_SHM)
    if (decomp%ROW_INFO%CORE_ME==1) THEN
       call MPI_ALLTOALLV(work1, decomp%y2cnts_s, decomp%y2disp_s, &
            complex_type, work2, decomp%z2cnts_s, decomp%z2disp_s, &
            complex_type, decomp%ROW_INFO%SMP_COMM, ierror)
    end if
#else
#if defined(_EVEN)
    if (decomp%even) then
       call MPI_ALLTOALL(work1_c, decomp%y2count, &
            complex_type, dst, decomp%z2count, &
            complex_type, DECOMP_2D_COMM_ROW, ierror)
    else
       call MPI_ALLTOALL(work1_c, decomp%y2count, &
            complex_type, work2_c, decomp%z2count, &
            complex_type, DECOMP_2D_COMM_ROW, ierror)
    end if
#else
    call MPI_ALLTOALLV(work1_c, decomp%y2cnts, decomp%y2disp, &
         complex_type, dst, decomp%z2cnts, decomp%z2disp, &
         complex_type, DECOMP_2D_COMM_ROW, ierror)
#endif
#endif

    ! rearrange receive buffer
#if defined(_SHM)
    call MPI_BARRIER(decomp%ROW_INFO%CORE_COMM, ierror)
    call mem_merge_yz_complex(work2, d1, d2, d3, dst, dims(2), &
         decomp%z2dist, decomp)
#else
#if defined(_EVEN)
    if (.not. decomp%even) then
       call mem_merge_yz_complex(work2_c, d1, d2, d3, dst, dims(2), &
         decomp%z2dist, decomp)
    end if
#else
    ! note the receive buffer is already in natural (i,j,k) order
    ! so no merge operation needed
#endif
#endif

    return
  end subroutine transpose_y_to_z_complex

#if defined(_OCC)
  subroutine transpose_y_to_z_complex_start(handle, src, dst, sbuf, &
       rbuf, opt_decomp)

    implicit none

    integer :: handle
    complex(mytype), dimension(:,:,:) :: src, dst, sbuf, rbuf
    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp

    TYPE(DECOMP_INFO) :: decomp

    integer :: s1,s2,s3
    integer :: ierror

    if (present(opt_decomp)) then
       decomp = opt_decomp
    else
       decomp = decomp_main
    end if

    s1 = SIZE(src,1)
    s2 = SIZE(src,2)
    s3 = SIZE(src,3)

    ! rearrange source array as send buffer
    call mem_split_yz_complex(src, s1, s2, s3, sbuf, dims(2), &
         decomp%y2dist, decomp)

#if defined(_EVEN)
    call NBC_IALLTOALL(sbuf, decomp%y2count, &
         complex_type, rbuf, decomp%z2count, &
         complex_type, DECOMP_2D_COMM_ROW, handle, ierror)
#else
    call NBC_IALLTOALLV(sbuf, decomp%y2cnts, decomp%y2disp, &
         complex_type, rbuf,decomp%z2cnts, decomp%z2disp, &
         complex_type, DECOMP_2D_COMM_ROW, handle, ierror)
#endif

    return
  end subroutine transpose_y_to_z_complex_start

  subroutine transpose_y_to_z_complex_wait(handle, src, dst, sbuf, &
       rbuf, opt_decomp)

    implicit none

    integer :: handle
    complex(mytype), dimension(:,:,:) :: src, dst, sbuf, rbuf
    TYPE(DECOMP_INFO), intent(IN), optional :: opt_decomp

    TYPE(DECOMP_INFO) :: decomp

    integer :: ierror

    if (present(opt_decomp)) then
       decomp = opt_decomp
    else
       decomp = decomp_main
    end if

    call NBC_WAIT(handle, ierror)

    dst = rbuf

    return
  end subroutine transpose_y_to_z_complex_wait
#endif

  ! pack/unpack ALLTOALL(V) buffers

  subroutine mem_split_yz_real(in,n1,n2,n3,out,iproc,dist,decomp)

    implicit none

    integer, intent(IN) :: n1,n2,n3
    real(mytype), dimension(n1,n2,n3), intent(IN) :: in
    real(mytype), dimension(*), intent(OUT) :: out
    integer, intent(IN) :: iproc
    integer, dimension(0:iproc-1), intent(IN) :: dist
    TYPE(DECOMP_INFO), intent(IN) :: decomp

    integer :: i,j,k, m,i1,i2,pos

    do m=0,iproc-1
       if (m==0) then
          i1 = 1
          i2 = dist(0)
       else
          i1 = i2+1
          i2 = i1+dist(m)-1
       end if

#if defined(_SHM)
       pos = decomp%y2disp_o(m) + 1
#else
#if defined(_EVEN)
       pos = m * decomp%y2count + 1
#else
       pos = decomp%y2disp(m) + 1
#endif
#endif

       do k=1,n3
          do j=i1,i2
             do i=1,n1
                out(pos) = in(i,j,k)
                pos = pos + 1
             end do
          end do
       end do
    end do

    return
  end subroutine mem_split_yz_real

  subroutine mem_split_yz_complex(in,n1,n2,n3,out,iproc,dist,decomp)

    implicit none

    integer, intent(IN) :: n1,n2,n3
    complex(mytype), dimension(n1,n2,n3), intent(IN) :: in
    complex(mytype), dimension(*), intent(OUT) :: out
    integer, intent(IN) :: iproc
    integer, dimension(0:iproc-1), intent(IN) :: dist
    TYPE(DECOMP_INFO), intent(IN) :: decomp

    integer :: i,j,k, m,i1,i2,pos

    do m=0,iproc-1
       if (m==0) then
          i1 = 1
          i2 = dist(0)
       else
          i1 = i2+1
          i2 = i1+dist(m)-1
       end if

#if defined(_SHM)
       pos = decomp%y2disp_o(m) + 1
#else
#if defined(_EVEN)
       pos = m * decomp%y2count + 1
#else
       pos = decomp%y2disp(m) + 1
#endif
#endif

       do k=1,n3
          do j=i1,i2
             do i=1,n1
                out(pos) = in(i,j,k)
                pos = pos + 1
             end do
          end do
       end do
    end do

    return
  end subroutine mem_split_yz_complex

  subroutine mem_merge_yz_real(in,n1,n2,n3,out,iproc,dist,decomp)

    implicit none

    integer, intent(IN) :: n1,n2,n3
    real(mytype), dimension(*), intent(IN) :: in
    real(mytype), dimension(n1,n2,n3), intent(OUT) :: out
    integer, intent(IN) :: iproc
    integer, dimension(0:iproc-1), intent(IN) :: dist
    TYPE(DECOMP_INFO), intent(IN) :: decomp

    integer :: i,j,k, m,i1,i2, pos

    do m=0,iproc-1
       if (m==0) then
          i1 = 1
          i2 = dist(0)
       else
          i1 = i2+1
          i2 = i1+dist(m)-1
       end if

#if defined(_SHM)
       pos = decomp%z2disp_o(m) + 1
#else
#if defined(_EVEN)
       pos = m * decomp%z2count + 1
#else
       pos = decomp%z2disp(m) + 1
#endif
#endif

       do k=i1,i2
          do j=1,n2
             do i=1,n1
                out(i,j,k) = in(pos)
                pos = pos + 1
             end do
          end do
       end do
    end do

    return
  end subroutine mem_merge_yz_real

  subroutine mem_merge_yz_complex(in,n1,n2,n3,out,iproc,dist,decomp)

    implicit none

    integer, intent(IN) :: n1,n2,n3
    complex(mytype), dimension(*), intent(IN) :: in
    complex(mytype), dimension(n1,n2,n3), intent(OUT) :: out
    integer, intent(IN) :: iproc
    integer, dimension(0:iproc-1), intent(IN) :: dist
    TYPE(DECOMP_INFO), intent(IN) :: decomp

    integer :: i,j,k, m,i1,i2, pos

    do m=0,iproc-1
       if (m==0) then
          i1 = 1
          i2 = dist(0)
       else
          i1 = i2+1
          i2 = i1+dist(m)-1
       end if

#if defined(_SHM)
       pos = decomp%z2disp_o(m) + 1
#else
#if defined(_EVEN)
       pos = m * decomp%z2count + 1
#else
       pos = decomp%z2disp(m) + 1
#endif
#endif

       do k=i1,i2
          do j=1,n2
             do i=1,n1
                out(i,j,k) = in(pos)
                pos = pos + 1
             end do
          end do
       end do
    end do

    return
  end subroutine mem_merge_yz_complex

max_stars_repo_name: fspiga/FluTAS

max_issues_repo_name: fspiga/FluTAS

max_forks_repo_name: fspiga/FluTAS


text:
c
c     This file contains the Stokes FMM wrappers
c
c**********************************************************************
c
c     We take the following conventions for the Stokes kernels
c
c     For a source y and target x, let r_i = x_i-y_i
c     and let r = sqrt(r_1^2 + r_2^2 + r_3^2)
c
c     The Stokeslet, G_{ij}, and its associated pressure tensor, P_j,
c     (without the 1/4pi scaling) are
c
c     G_{ij}(x,y) = (r_i r_j)/(2r^3) + delta_{ij}/(2r)
c     P_j(x,y) = r_j/r^3
c
c     The (Type I) stresslet, T_{ijk}, and its associated pressure
c     tensor, PI_{jk}, (without the 1/4pi scaling) are
c
c     T_{ijk}(x,y) = -3 r_i r_j r_k/ r^5
c     PI_{jk} = -2 delta_{jk} + 6 r_j r_k/r^5
c

      subroutine stfmm3d(nd, eps,
     $                 nsource, source,
     $                 ifstoklet, stoklet, ifstrslet, strslet, strsvec,
     $                 ifppreg, pot, pre, grad, ntarg, targ,
     $                 ifppregtarg, pottarg, pretarg, gradtarg,ier)
c
c     Stokes FMM in R^{3}: evaluate all pairwise particle
c     interactions (ignoring self-interactions) and
c     interactions with targs.
c
c     This routine computes sums of the form
c
c       u(x) = sum_m G_{ij}(x,y^{(m)}) sigma^{(m)}_j
c                + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k
c
c     where sigma^{(m)} is the Stokeslet charge, mu^{(m)} is the
c     stresslet charge, and nu^{(m)} is the stresslet orientation
c     (note that each of these is a 3 vector per source point y^{(m)}).
c     For x a source point, the self-interaction in the sum is omitted.
c
c     Optionally, the associated pressure p(x) and gradient grad u(x)
c     are returned
c
c       p(x) = sum_m P_j(x,y^m) sigma^{(m)}_j
c          + sum_m T_{ijk}(x,y^{(m)}) PI_{jk} mu^{(m)}_j nu^{(m)}_k
c
c       grad u(x) = grad[sum_m G_{ij}(x,y^m) sigma^{(m)}_j
c                + sum_m T_{ijk}(x,y^{(m)}) mu^{(m)}_j nu^{(m)}_k]
c
c-----------------------------------------------------------------------
c     INPUT PARAMETERS:
c
c   nd:    in: integer
c              number of densities
c
c   eps:   in: double precision
c              requested precision
c
c   nsource in: integer
c               number of sources
c
c   source  in: double precision (3,nsource)
c               source(k,j) is the kth component of the jth
c               source locations
c
c   ifstoklet  in: integer
c               Stokeslet charge computation flag
c               ifstoklet = 1   =>  include Stokeslet contribution
c                                   otherwise do not
c
c   stoklet in: double precision (nd,3,nsource)
c               Stokeslet charge strengths (sigma vectors above)
c
c   ifstrslet in: integer
c               stresslet computation flag
c               ifstrslet = 1   =>  include standard stresslet
c                                   (type I)
c
c            NOT YET IMPLEMENTED
c
c               ifstrslet = 2   =>  include symmetric stresslet
c                                   (type II)
c               ifstrslet = 3   =>  include rotlet
c               ifstrslet = 4   =>  include Stokes doublet
c                      otherwise do not include
c
c   strslet  in: double precision (nd,3,nsource)
c               stresslet strengths (mu vectors above)
c
c   strsvec  in: double precision (nd,3,nsource)
c               stresslet orientations (nu vectors above)
c
c     ifppreg    in: integer
c               flag for evaluating potential, gradient, and pressure
c               at the sources
c               ifppreg = 1, only potential
c               ifppreg = 2, potential and pressure
c         GRADIENT NOT IMPLEMENTED
c               ifppreg = 3, potential, pressure, and gradient
c
c   ntarg   in: integer
c              number of targs
c
c   targ    in: double precision (3,ntarg)
c             targ(k,j) is the kth component of the jth
c             targ location
c
c   ifppregtarg in: integer
c                flag for evaluating potential, gradient, and pressure
c                at the targets
c                ifppregtarg = 1, only potential
c                ifppregtarg = 2, potential and pressure
c                ifppregtarg = 3, potential, pressure, and gradient
c
c-----------------------------------------------------------------------
c
c   OUTPUT parameters:
c
c   pot   out: double precision(nd,3,nsource)
c           velocity at the source locations
c
c   pre   out: double precision(nd,nsource)
c           pressure at the source locations
c
c         GRADIENT NOT IMPLEMENTED
c   grad   out: double precision(nd,3,3,nsource)
c              gradient of velocity at the source locations
c              grad(l,i,j,k) is the ith component of the
c              gradient of the jth component of the velocity
c              for the lth density at the kth source location
c
c   pottarg   out: double precision(nd,3,ntarg)
c               velocity at the targets
c
c   pretarg   out: double precision(nd,ntarg)
c               pressure at the targets
c
c   gradtarg   out: double precision(nd,3,3,ntarg)
c               gradient of velocity at the targets
c               gradtarg(l,i,j,k) is the ith component of the
c               gradient of the jth component of the velocity
c               for the lth density at the kth target
c     ier     out: integer
c               error flag
c
c     TODO: implement other stresslet options and gradient
c------------------------------------------------------------------
      implicit none
      integer nd, ifstoklet, ifstrslet, ntarg
      double precision eps
      integer nsource, ifppreg, ifppregtarg
      double precision source(3, *), targ(3, *)
      double precision stoklet(nd, 3, *), strslet(nd, 3, *)
      double precision strsvec(nd, 3, *)
      double precision pot(nd, 3, *), pre(nd,*), grad(nd, 3, 3, *)
      double precision pottarg(nd, 3, *), pretarg(nd,*),
     1     gradtarg(nd, 3, 3, *)

c     local
      double precision, allocatable :: charge(:,:,:), dipvec(:,:,:,:)
      double precision, allocatable :: potl(:,:,:), gradl(:,:,:,:),
     1     hessl(:,:,:,:), pottargl(:,:,:), gradtargl(:,:,:,:),
     2     hesstargl(:,:,:,:)
      double precision :: pt(3), gl(3), hl(6), vel(3), velgrad(3,3)
      double precision :: press, pl, pv, dmu(3), dnu(3), sigma(3)

      integer ndl, ifchargel, ifdipolel, ifpghl, ifpghtargl
      integer ndper

      integer i, j, ii, ifppreg1, l, npt, ier,iper

      ndper = 0

      if (ifstrslet .eq. 1 .or. ifstoklet .eq. 1) then
         ndper = 4
      endif

      ifdipolel = 0
      ifchargel = 0

      if (ifstoklet .eq. 1) ifchargel = 1
      if (ifstrslet .eq. 1) ifdipolel = 1

      ndper = 4
      ndl = ndper*nd

      ifpghl = 3
      ifpghtargl = 3

c     allocate necessary arrays

      allocate(charge(ndper,nd,nsource),dipvec(ndper,nd,3,nsource),
     1     potl(ndper,nd,nsource),pottargl(ndper,nd,ntarg),
     2     gradl(ndper,nd,3,nsource),gradtargl(ndper,nd,3,ntarg),
     3     hessl(ndper,nd,6,nsource),hesstargl(ndper,nd,6,ntarg),
     4     stat=ier)
      if(ier .ne. 0) then
         print *, "In stfmm3d: cannot allocate Laplace call storage"
         print *, "ndper =",ndper
         print *, "nd =",nd
         print *, "nsource =",nsource
         print *, "ntarg =",ntarg
         stop
      endif

c     set-up appropriate vector charge and dipole arrays

c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,l,sigma,dmu,dnu)
c$OMP$ PRIVATE(pl,pv)
      do i = 1,nsource

         do j = 1,nd
            do l = 1,ndper
               charge(l,j,i) = 0
               dipvec(l,j,1,i) = 0
               dipvec(l,j,2,i) = 0
               dipvec(l,j,3,i) = 0
            enddo
         enddo

         do j = 1,nd
            if (ifstoklet .eq. 1) then
               sigma(1) = stoklet(j,1,i)
               sigma(2) = stoklet(j,2,i)
               sigma(3) = stoklet(j,3,i)
            endif
            if (ifstrslet .ge. 1) then
               dmu(1) = strslet(j,1,i)
               dmu(2) = strslet(j,2,i)
               dmu(3) = strslet(j,3,i)
               dnu(1) = strsvec(j,1,i)
               dnu(2) = strsvec(j,2,i)
               dnu(3) = strsvec(j,3,i)
            endif

            do l = 1,3

               if (ifstoklet .eq. 1) then
                  charge(l,j,i) = charge(l,j,i) + sigma(l)/2
               endif
               if (ifstrslet .eq. 1) then
                  dipvec(l,j,1,i) = dipvec(l,j,1,i) - (dmu(l)*dnu(1) +
     1                 dmu(1)*dnu(l))/2
                  dipvec(l,j,2,i) = dipvec(l,j,2,i) - (dmu(l)*dnu(2) +
     1                 dmu(2)*dnu(l))/2
                  dipvec(l,j,3,i) = dipvec(l,j,3,i) - (dmu(l)*dnu(3) +
     1                 dmu(3)*dnu(l))/2
               endif
            enddo

            l = 4

            if (ifstoklet .eq. 1) then
               pl = sigma(1)*source(1,i) + sigma(2)*source(2,i) +
     1              sigma(3)*source(3,i)
               charge(l,j,i) = charge(l,j,i) + pl/2
            endif
            if (ifstrslet .eq. 1) then
               pl = dmu(1)*source(1,i) + dmu(2)*source(2,i) +
     1              dmu(3)*source(3,i)
               pv = dnu(1)*source(1,i) + dnu(2)*source(2,i) +
     1              dnu(3)*source(3,i)

               dipvec(l,j,1,i) = dipvec(l,j,1,i) -
     1              (dmu(1)*pv + dnu(1)*pl)/2
               dipvec(l,j,2,i) = dipvec(l,j,2,i) -
     1              (dmu(2)*pv + dnu(2)*pl)/2
               dipvec(l,j,3,i) = dipvec(l,j,3,i) -
     1              (dmu(3)*pv + dnu(3)*pl)/2
            endif

         enddo

      enddo
c$OMP END PARALLEL DO

c     call Laplace FMM

      iper = 0
      ier = 0
      call lfmm3d(ndl,eps,nsource,source,ifchargel,charge,
     1     ifdipolel,dipvec,iper,ifpghl,potl,gradl,hessl,ntarg,
     2     targ,ifpghtargl,pottargl,gradtargl,hesstargl,ier)

      npt = ntarg + nsource

c     unpack stacked Laplace FMM calls

c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,l,ifppreg1,ii)
c$OMP$ PRIVATE(pt,pl,gl,hl,vel,velgrad,press)
      do i = 1,npt

         do j = 1,nd

            vel(1) = 0
            vel(2) = 0
            vel(3) = 0
            velgrad(1,1) = 0
            velgrad(2,1) = 0
            velgrad(3,1) = 0
            velgrad(1,2) = 0
            velgrad(2,2) = 0
            velgrad(3,2) = 0
            velgrad(1,3) = 0
            velgrad(2,3) = 0
            velgrad(3,3) = 0
            press = 0

            do l = 1,ndper

               if (i .gt. ntarg) then
                  ifppreg1 = ifppreg
                  ii = i-ntarg
                  pt(1) = source(1,ii)
                  pt(2) = source(2,ii)
                  pt(3) = source(3,ii)
                  pl = potl(l,j,ii)
                  gl(1) = gradl(l,j,1,ii)
                  gl(2) = gradl(l,j,2,ii)
                  gl(3) = gradl(l,j,3,ii)
                  hl(1) = hessl(l,j,1,ii)
                  hl(2) = hessl(l,j,2,ii)
                  hl(3) = hessl(l,j,3,ii)
                  hl(4) = hessl(l,j,4,ii)
                  hl(5) = hessl(l,j,5,ii)
                  hl(6) = hessl(l,j,6,ii)
               else
                  ifppreg1 = ifppregtarg
                  ii = i
                  pt(1) = targ(1,ii)
                  pt(2) = targ(2,ii)
                  pt(3) = targ(3,ii)
                  pl = pottargl(l,j,ii)
                  gl(1) = gradtargl(l,j,1,ii)
                  gl(2) = gradtargl(l,j,2,ii)
                  gl(3) = gradtargl(l,j,3,ii)
                  hl(1) = hesstargl(l,j,1,ii)
                  hl(2) = hesstargl(l,j,2,ii)
                  hl(3) = hesstargl(l,j,3,ii)
                  hl(4) = hesstargl(l,j,4,ii)
                  hl(5) = hesstargl(l,j,5,ii)
                  hl(6) = hesstargl(l,j,6,ii)
               endif

               if (l .ge. 1 .and. l .le. 3) then

                  vel(l) = vel(l) + pl
                  vel(1) = vel(1) - pt(l)*gl(1)
                  vel(2) = vel(2) - pt(l)*gl(2)
                  vel(3) = vel(3) - pt(l)*gl(3)
                  press = press - gl(l)*2
                  if (ifppreg1 .eq. 3) then
                     velgrad(1,l) =  velgrad(1,l) + gl(1)
                     velgrad(2,l) =  velgrad(2,l) + gl(2)
                     velgrad(3,l) =  velgrad(3,l) + gl(3)

                     velgrad(l,1) = velgrad(l,1) - gl(1)
                     velgrad(l,2) = velgrad(l,2) - gl(2)
                     velgrad(l,3) = velgrad(l,3) - gl(3)

c     confirm hessian ordering convention...
                     velgrad(1,1) = velgrad(1,1) - pt(l)*hl(1)
                     velgrad(2,1) = velgrad(2,1) - pt(l)*hl(4)
                     velgrad(3,1) = velgrad(3,1) - pt(l)*hl(5)
                     velgrad(1,2) = velgrad(1,2) - pt(l)*hl(4)
                     velgrad(2,2) = velgrad(2,2) - pt(l)*hl(2)
                     velgrad(3,2) = velgrad(3,2) - pt(l)*hl(6)
                     velgrad(1,3) = velgrad(1,3) - pt(l)*hl(5)
                     velgrad(2,3) = velgrad(2,3) - pt(l)*hl(6)
                     velgrad(3,3) = velgrad(3,3) - pt(l)*hl(3)
                  endif

               else if (l .eq. 4) then

                  vel(1) = vel(1) + gl(1)
                  vel(2) = vel(2) + gl(2)
                  vel(3) = vel(3) + gl(3)
                  if (ifppreg1 .eq. 3) then
c     confirm hessian ordering convention...
                     velgrad(1,1) = velgrad(1,1) + hl(1)
                     velgrad(2,1) = velgrad(2,1) + hl(4)
                     velgrad(3,1) = velgrad(3,1) + hl(5)
                     velgrad(1,2) = velgrad(1,2) + hl(4)
                     velgrad(2,2) = velgrad(2,2) + hl(2)
                     velgrad(3,2) = velgrad(3,2) + hl(6)
                     velgrad(1,3) = velgrad(1,3) + hl(5)
                     velgrad(2,3) = velgrad(2,3) + hl(6)
                     velgrad(3,3) = velgrad(3,3) + hl(3)
                  endif

               endif
            enddo

            if (i .gt. ntarg) then
               if (ifppreg1 .ge. 1) then
                  pot(j,1,ii) = vel(1)
                  pot(j,2,ii) = vel(2)
                  pot(j,3,ii) = vel(3)
               endif
               if (ifppreg1 .ge. 2) then
                  pre(j,ii) = press
               endif
               if (ifppreg1 .ge. 3) then
                  grad(j,1,1,ii) = velgrad(1,1)
                  grad(j,2,1,ii) = velgrad(2,1)
                  grad(j,3,1,ii) = velgrad(3,1)
                  grad(j,1,2,ii) = velgrad(1,2)
                  grad(j,2,2,ii) = velgrad(2,2)
                  grad(j,3,2,ii) = velgrad(3,2)
                  grad(j,1,3,ii) = velgrad(1,3)
                  grad(j,2,3,ii) = velgrad(2,3)
                  grad(j,3,3,ii) = velgrad(3,3)
               endif
            else
               if (ifppreg1 .ge. 1) then
                  pottarg(j,1,ii) = vel(1)
                  pottarg(j,2,ii) = vel(2)
                  pottarg(j,3,ii) = vel(3)
               endif
               if (ifppreg1 .ge. 2) then
                  pretarg(j,ii) = press
               endif
               if (ifppreg1 .ge. 3) then
                  gradtarg(j,1,1,ii) = velgrad(1,1)
                  gradtarg(j,2,1,ii) = velgrad(2,1)
                  gradtarg(j,3,1,ii) = velgrad(3,1)
                  gradtarg(j,1,2,ii) = velgrad(1,2)
                  gradtarg(j,2,2,ii) = velgrad(2,2)
                  gradtarg(j,3,2,ii) = velgrad(3,2)
                  gradtarg(j,1,3,ii) = velgrad(1,3)
                  gradtarg(j,2,3,ii) = velgrad(2,3)
                  gradtarg(j,3,3,ii) = velgrad(3,3)
               endif
            endif
         enddo
      enddo
c$OMP END PARALLEL DO

      return
      end


text:
      SUBROUTINE SOLVE_LAPLACE
     &     (BCD,     NSD,
     &      NNW_SD,  NGP_SD,
     &      NW_IPAR, NW_RPAR,
     &      CRD,     PHI,     PHN,
     &      S0,      S1,      D0,      D1)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'tme_params.inc'
      INCLUDE 'chk_params.inc'
      INCLUDE 'bct_params.inc'
C
      INTEGER(KIND=IK) BCD
      INTEGER(KIND=IK) NSD
      INTEGER(KIND=IK) NNW_SD (*)
      INTEGER(KIND=IK) NGP_SD (*)
      INTEGER(KIND=IK) NW_IPAR (N_IP, *)
      REAL   (KIND=RK) NW_RPAR (N_RP, *)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) S0 (*)
      REAL   (KIND=RK) S1 (*)
      REAL   (KIND=RK) D0 (*)
      REAL   (KIND=RK) D1 (*)
C
      IF ((GRID_TIME_DEPENDENT .OR. FIRST_TIME)
     &    .AND. BCD .NE. BCD_PHI_T) THEN
C        Depending on the type of boundary condition some grids must be
C        recomputed (e.g. R&F bc's)
         CALL REGRID
     &        (NSD, NNW_SD, NGP_SD, NW_IPAR, NW_RPAR, CRD)
C        Compute the source and dipole coefficients only if the grid is
C        time-dependent or if this is the first time
         CALL COEFFICIENTS
     &        (NSD, NNW_SD, NGP_SD, NW_IPAR, CRD, S0, S1, D0, D1)
      END IF
C     Compute the boundary conditions for those networks that use an
C     analytic solution.
      CALL BOUNDARY_CONDITIONS
     &     (BCD, NSD, NNW_SD, NW_IPAR, NW_RPAR, CRD, PHI, PHN)
C     Check if the solution satisfies the BIE
      IF (CHK_COEFS_PRE) CALL CHECK_COEFS
     &     (NSD, NGP_SD, PHI, PHN, S0, S1, D0, D1)
C     Solve the Boundary Integral equations
      CALL SOLVE_BIE
     &     (NSD,
     &      NNW_SD, NGP_SD,  NW_IPAR,
     &      S0,     S1,
     &      D0,     D1,
     &      CRD,    PHI,     PHN)
C     Check if the solution satisfies the BIE
      IF (CHK_COEFS_POST) CALL CHECK_COEFS
     &     (NSD, NGP_SD, PHI, PHN, S0, S1, D0, D1)
C     The first boundary value problem has been solved
      FIRST_TIME = .FALSE.
C
      RETURN
      END

max_stars_repo_name: meggermo/waterwaves

max_issues_repo_name: meggermo/waterwaves

max_forks_repo_name: meggermo/waterwaves


text:
      subroutine fzero(a,ns,ne)
c     implicit double precision (a-h,o-z)                                    dp
      dimension a(*)
      do 10 i=ns,ne
      a(i)=0.
   10 continue
      return
      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


text:
module Mod_RpMeshAllocator
   use typre
   use Mod_Memor
   use Mod_Mesh
   use Mod_Element
   implicit none

contains

   subroutine AllocR1P(Memor,Mesh,r1p_array)
      type(MemoryMan) :: Memor
      class(FemMesh) :: Mesh
      type(r1p), allocatable :: r1p_array(:)

      class(FiniteElement), pointer :: e => NULL()
      integer(ip) :: ielem,nelem

      call Mesh%GetNelem(nelem)
      call Memor%alloc(nelem,r1p_array,'DPostprocess','MD_SetupPostprocess')

      call Mesh%ElementAlloc(e,Memor,'ForceClosedRule','plcd_memall')
      elements : do ielem = 1,nelem
         !Load Element
         call Mesh%ElementLoad(ielem,e)
         call Memor%palloc(e%pgaus,r1p_array(ielem)%a,'DPostprocess','plcd_memall')
      enddo elements
      call Mesh%ElementdeAlloc(e,Memor,'ForceClosedRule','plcd_memall')

   end subroutine

   subroutine DeAllocR1P(Memor,Mesh,r1p_array)
      type(MemoryMan) :: Memor
      class(FemMesh) :: Mesh
      type(r1p), allocatable :: r1p_array(:)

      class(FiniteElement), pointer :: e => NULL()
      integer(ip) :: ielem,nelem

      call Mesh%GetNelem(nelem)
      call Mesh%ElementAlloc(e,Memor,'ForceClosedRule','plcd_memall')
      elements : do ielem = 1,nelem
         !Load Element
         call Mesh%ElementLoad(ielem,e)
         call Memor%pdealloc(e%pgaus,r1p_array(ielem)%a,'DPostprocess','plcd_memall')
      enddo elements
      call Mesh%ElementdeAlloc(e,Memor,'ForceClosedRule','plcd_memall')

      call Memor%dealloc(nelem,r1p_array,'DPostprocess','MD_SetupPostprocess')

   end subroutine

   subroutine AllocR2P(Memor,Mesh,ndime1,r1p_array)
      type(MemoryMan) :: Memor
      class(FemMesh) :: Mesh
      type(r2p), allocatable :: r1p_array(:)
      integer(ip) :: ndime1

      class(FiniteElement), pointer :: e => NULL()
      integer(ip) :: ielem,nelem

      call Mesh%GetNelem(nelem)
      call Memor%alloc(nelem,r1p_array,'DPostprocess','MD_SetupPostprocess')

      call Mesh%ElementAlloc(e,Memor,'ForceClosedRule','plcd_memall')
      elements : do ielem = 1,nelem
         !Load Element
         call Mesh%ElementLoad(ielem,e)
         call Memor%palloc(ndime1,e%pgaus,r1p_array(ielem)%a,'DPostprocess','plcd_memall')
      enddo elements
      call Mesh%ElementdeAlloc(e,Memor,'ForceClosedRule','plcd_memall')

   end subroutine

   subroutine DeAllocR2P(Memor,Mesh,ndime1,r1p_array)
      type(MemoryMan) :: Memor
      class(FemMesh) :: Mesh
      type(r2p), allocatable :: r1p_array(:)
      integer(ip) :: ndime1

      class(FiniteElement), pointer :: e => NULL()
      integer(ip) :: ielem,nelem

      call Mesh%GetNelem(nelem)
      call Mesh%ElementAlloc(e,Memor,'ForceClosedRule','plcd_memall')
      elements : do ielem = 1,nelem
         !Load Element
         call Mesh%ElementLoad(ielem,e)
         call Memor%pdealloc(ndime1,e%pgaus,r1p_array(ielem)%a,'DPostprocess','plcd_memall')
      enddo elements
      call Mesh%ElementdeAlloc(e,Memor,'ForceClosedRule','plcd_memall')

      call Memor%dealloc(nelem,r1p_array,'DPostprocess','MD_SetupPostprocess')

   end subroutine

   subroutine AllocR3P(Memor,Mesh,ndime1,ndime2,r1p_array)
      type(MemoryMan) :: Memor
      class(FemMesh) :: Mesh
      type(r3p), allocatable :: r1p_array(:)
      integer(ip) :: ndime1,ndime2

      class(FiniteElement), pointer :: e => NULL()
      integer(ip) :: ielem,nelem

      call Mesh%GetNelem(nelem)
      call Memor%alloc(nelem,r1p_array,'DPostprocess','MD_SetupPostprocess')

      call Mesh%ElementAlloc(e,Memor,'ForceClosedRule','plcd_memall')
      elements : do ielem = 1,nelem
         !Load Element
         call Mesh%ElementLoad(ielem,e)
         call Memor%palloc(ndime1,ndime2,e%pgaus,r1p_array(ielem)%a,'DPostprocess','plcd_memall')
      enddo elements
      call Mesh%ElementdeAlloc(e,Memor,'ForceClosedRule','plcd_memall')

   end subroutine

   subroutine DeAllocR3P(Memor,Mesh,ndime1,ndime2,r1p_array)
      type(MemoryMan) :: Memor
      class(FemMesh) :: Mesh
      type(r3p), allocatable :: r1p_array(:)
      integer(ip) :: ndime1,ndime2

      class(FiniteElement), pointer :: e => NULL()
      integer(ip) :: ielem,nelem

      elements : do ielem = 1,size(r1p_array,1)
         call Memor%pdealloc(size(r1p_array(ielem)%a,1),size(r1p_array(ielem)%a,2),size(r1p_array(ielem)%a,3),r1p_array(ielem)%a,'DPostprocess','plcd_memall')
      enddo elements

      call Memor%dealloc(size(r1p_array,1),r1p_array,'DPostprocess','MD_SetupPostprocess')

   end subroutine

end module Mod_RpMeshAllocator

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
! THIS VERSION: 25/06/2002 AT 14:00:00 PM.
! Updated 25/06/2002: additional warning information added

!-*-*-*-*-*  L A N C E L O T  -B-  DUMMY AD02_FORWARD  M O D U L E S *-*-*-*

!  Nick Gould, for GALAHAD productions
!  Copyright reserved
!  June 28th 1996

MODULE HSL_AD02_FORWARD_DOUBLE

      IMPLICIT NONE
      PRIVATE
      INTEGER, PARAMETER :: WP = KIND(1D0)
      PUBLIC :: AD02_INITIALIZE

!  Dummy HSL_AD02_FORWARD_DOUBLE module

      TYPE, PUBLIC :: AD02_REAL
        PRIVATE
        INTEGER :: P
      END TYPE AD02_REAL

      TYPE, PUBLIC :: AD02_DATA
        PRIVATE
        INTEGER :: LEVEL
      END TYPE AD02_DATA

CONTAINS

      SUBROUTINE AD02_INITIALIZE(DEGREE,A,VALUE,DATA,FULL_THRESHOLD)
        INTEGER, INTENT (IN) :: DEGREE
!       TYPE (AD02_REAL), INTENT (OUT) :: A
        TYPE (AD02_REAL) :: A
        INTEGER, OPTIONAL, INTENT (IN) :: FULL_THRESHOLD
        TYPE (AD02_DATA), POINTER :: DATA
        REAL (WP), INTENT (IN) :: VALUE

!  Dummy subroutine available with LANCELOT

        WRITE ( 6, 2000 )
        STOP

!  Non-executable statements

 2000    FORMAT( /, ' We regret that the solution options that you have ', /, &
                    ' chosen are not all freely available with LANCELOT.', //,&
                    ' If you have HSL (formerly the Harwell Subroutine',      &
                    ' Library), this ', /,                                    &
                    ' option may be enabled by replacing the dummy ', /,      &
                    ' module HSL_AD02_FORWARD_DOUBLE with its H...', /,       &
                    ' namesake and dependencies. See', /,                     &
                    '   $GALAHAD/src/makedefs/packages for details.', //,     &
                    ' *** EXECUTION TERMINATING *** ', / )

      END SUBROUTINE AD02_INITIALIZE

END MODULE HSL_AD02_FORWARD_DOUBLE

!  THIS VERSION: 28/06/1996 AT 09:00:00 AM

!-*-*-*-*-*  L A N C E L O T  -B-  DUMMY AD02_BACKWARD  M O D U L E S *-*-*-*

!  Nick Gould, for GALAHAD productions
!  Copyright reserved
!  June 28th 1996

MODULE HSL_AD02_BACKWARD_DOUBLE

      IMPLICIT NONE
      PRIVATE
      INTEGER, PARAMETER :: WP = KIND(1D0)
      PUBLIC :: AD02_INITIALIZE

!  Dummy HSL_AD02_BACKWARD_DOUBLE module

      TYPE, PUBLIC :: AD02_REAL
        PRIVATE
        INTEGER :: P
      END TYPE AD02_REAL

      TYPE, PUBLIC :: AD02_DATA
        PRIVATE
        INTEGER :: LEVEL
      END TYPE AD02_DATA

CONTAINS

      SUBROUTINE AD02_INITIALIZE(DEGREE,A,VALUE,DATA,FULL_THRESHOLD)
        INTEGER, INTENT (IN) :: DEGREE
!       TYPE (AD02_REAL), INTENT (OUT) :: A
        TYPE (AD02_REAL) :: A
        INTEGER, OPTIONAL, INTENT (IN) :: FULL_THRESHOLD
        TYPE (AD02_DATA), POINTER :: DATA
        REAL (WP), INTENT (IN) :: VALUE

!  Dummy subroutine available with LANCELOT

        WRITE ( 6, 2000 )
        STOP

!  Non-executable statements

 2000    FORMAT( /, ' We regret that the solution options that you have ', /, &
                    ' chosen are not all freely available with LANCELOT.', //,&
                    ' If you have HSL (formerly the Harwell Subroutine',      &
                    ' Library), this ', /,                                    &
                    ' option may be enabled by replacing the dummy ', /,      &
                    ' module HSL_AD02_BACKWARD_DOUBLE with its HSL ',/,       &
                    ' namesake and dependencies. See', /,                     &
                    '   $GALAHAD/src/makedefs/packages for details.', //,     &
                    ' *** EXECUTION TERMINATING *** ', / )

      END SUBROUTINE AD02_INITIALIZE

END MODULE HSL_AD02_BACKWARD_DOUBLE

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


text:
Auburn is a suburb of Sacramento.

40 minutes NorthEast of Davis via I 80 the super highway.

Here are two websites about Auburn, CA:
   http://www.auburn.ca.gov/
   http://www.auburnca.com/

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20090707 03:15:02 nbsp Welcome to the Wiki Howdy and Welcome to the Wiki!    My names Evan, pleased to meet you!  Other editors are removing your edits because they... well, quite frankly, they look and smell like spam.  Spamming isnt nice... its basically like going down to your local community park and spray painting your link on the slide.  Somebody has to clean it up, and its just hurting your neighbors.  Not a nice thing to do at all. Users/JabberWokky Evan JabberWokky Edwards

20090707 00:22:06 nbsp Also, please dont create sockpuppet multiple accounts.  Youre a person, and you are actually quite welcome here... but act like a neighbor, not like somebody trying to squeeze every advantage out of the volunteer work of the people you live next to.  Try being yourself, not an invented Users/Davissushi.blogspot.com advertising gimmick.  Users/JabberWokky


text:
#include "scalar.fpp"

!=====================================================================!
! Module that contains vector type definition and procedures to add
! and get entries from the vector. This class also defines the
! interface to add elements into the vector.
!
! Author: Komahan Boopathy (komahan@gatech.edu)
!=====================================================================!

module vector_interface

  implicit none

  private
  public :: vector

  type, abstract :: vector

     ! type options
     private

     ! attributes
     integer :: size ! the number of columns in a vector

   contains

     ! defined procedures
     procedure :: get_size
     procedure :: set_size

     ! deferred procedures
     procedure(add_entry_interface), deferred :: add_entry
     procedure(get_entry_interface), deferred :: get_entry

  end type vector

  ! Define interfaces to all abstract procedures

  abstract interface

     !================================================================!
     ! Interface to adding entryies into the vector
     !================================================================!

     subroutine add_entry_interface(this, idx, val)

       import vector

       class(vector) , intent(inout) :: this
       integer       , intent(in)    :: idx
       type(scalar)  , intent(in)    :: val

     end subroutine add_entry_interface

     type(scalar) function get_entry_interface(this, idx)

       import vector

       class(vector) , intent(in) :: this
       type(integer) , intent(in) :: idx

     end function get_entry_interface

  end interface

contains

  !===================================================================!
  ! Returns the size of the vector
  !===================================================================!

  pure type(integer) function get_size(this)

    class(vector), intent(in) :: this

    get_size = this % size

  end function get_size

  !===================================================================!
  ! Sets the size of the vector
  !===================================================================!

  pure subroutine set_size(this, size)

    class(vector), intent(inout) :: this
    type(integer), intent(in)    :: size

    this % size = size

  end subroutine set_size

end module vector_interface

max_stars_repo_name: komahanb/time-marching-adjoint

max_issues_repo_name: komahanb/diffusion

max_forks_repo_name: komahanb/diffusion


text:
MODULE jmpMod
  USE iso_Fortran_env, ONLY: wp => real64
  USE mpi
  USE MPIMod
  USE Param
  USE NL
  USE sim
  IMPLICIT NONE

  ! IO:
  CHARACTER (LEN=*), PARAMETER :: ioDir = "./"
  CHARACTER (LEN=*), PARAMETER :: outDir = ioDir // "results/"
  CHARACTER (LEN=*), PARAMETER :: shockDir = ioDir // "shocks/"
  CHARACTER, PARAMETER:: tab = CHAR(9)

  ! Grids
  REAL(wp), DIMENSION(ySz) :: y, EdCre
  REAL(wp), DIMENSION(ySz, ySz) :: exoTran
  REAL(wp), DIMENSION(bSSz) :: bS
  REAL(wp), DIMENSION(bLSz) :: bL

  ! ArrayIxs: values, prices, policies
  REAL(wp), DIMENSION(ySz, bSSz, bLSz) :: V0, V1
  REAL(wp), DIMENSION(ySz, bSSz, bLSz) :: Va0, Va1

  REAL(wp), DIMENSION(ySz) :: Vd0, Vd1
  REAL(wp), DIMENSION(ySz) :: hy, uhy, Vaut

  REAL(wp), DIMENSION(ySz, bSSz, bLSz) :: qS0, qS1, qL0, qL1
  REAL(wp), DIMENSION(ySz, bSSz, bLSz) :: qSa0, qSa1, qLa0, qLa1
  REAL(wp), DIMENSION(ySz, bSSz, bLSz) :: xiS, xiL, xiSkeep, xiLkeep

  REAL(wp), DIMENSION(ySz, bSSz, bLSz) :: defPr
  REAL(wp), DIMENSION(:, :, :, :), ALLOCATABLE :: chPr
  INTEGER, DIMENSION(:, :, :, :), ALLOCATABLE :: chS, chL

  REAL(wp), DIMENSION(:, :, :, :), ALLOCATABLE :: chPrA
  INTEGER, DIMENSION(:, :, :, :), ALLOCATABLE :: chSa, chLa

  REAL(wp), DIMENSION(ySz, chPrSz) :: gammaPr
  INTEGER, DIMENSION(ySz, chPrSz) :: gammaS, gammaL

  TYPE stateSp
    INTEGER :: yIx, bSIx, bLIx
  END TYPE stateSp

  TYPE(stateSp), DIMENSION(:), ALLOCATABLE :: toSS

CONTAINS
  !
  !
  !

  PURE FUNCTION uu(cc) RESULT(uuval)
    IMPLICIT NONE
    REAL(wp), INTENT(IN) :: cc
    REAL(wp) :: uuval

    IF (cc <= 0.0_wp) THEN
      uuval = veryNeg
    ELSEIF (crra == 1.0) THEN
      uuval = LOG(cc)
    ELSEIF (crra == 2.0) THEN
      uuval = 1.0_wp - 1.0_wp / cc
    ELSE
      uuval = (cc**(1.0_wp - crra) - 1.0_wp) / (1.0_wp - crra)
    END IF
  END FUNCTION uu

  SUBROUTINE loadShocks()
    IMPLICIT NONE

    OPEN(unit=100, file=shockDir // "y.tab", status="old", action="read", position="rewind")
    READ (100, *) y
    CLOSE(100)

    OPEN(unit=100, file=shockDir // "PI.tab", status="old", action="read", position="rewind")
    READ (100, *) exoTran
    CLOSE(100)
    exoTran = TRANSPOSE(exoTran)
    exoTran(:, ySz) = MAX( 0.0_wp , 1.0_wp - SUM(exoTran(:, 1:(ySz-1)), 2) )
  END SUBROUTINE loadShocks

  SUBROUTINE generateGrids()
    IMPLICIT NONE
    INTEGER :: ix, yIx, bSIx, bLIx

    ALLOCATE(toSS(spaceSz))
    ix = 1
    DO bLIx = 1,bLSz
    DO bSIx = 1,bSSz
    DO yIx = 1,ySz
      toSS(ix)%yIx = yIx
      toSS(ix)%bSIx = bSIx
      toSS(ix)%bLIx = bLIx
      ix = ix + 1
    END DO
    END DO
    END DO

    CALL linspace(bS(1:bSsplit), bSmin, bSmid, bSsplit)
    CALL linspace(bS(bSsplit:bSSz), bSmid, bSmax, bSSz - bSsplit + 1)

    CALL linspace(bL(1:bLsplit), bLmin, bLmid, bLsplit)
    CALL linspace(bL(bLsplit:bLSz), bLmid, bLmax, bLSz - bLsplit + 1)

    DO yIx = 1,ySz
      hy(yIx) = y(yIx) - MAX( 0.0_wp, lbd0 * y(yIx) + lbd1 * y(yIx) * y(yIx) )
      ! hy(yIx) = MIN( y(yIx), 1.0_wp - lbd )
      uhy(yIx) = uu( hy(yIx) )
    END DO
  END SUBROUTINE generateGrids

  SUBROUTINE findVaut()
    IMPLICIT NONE
    INTEGER :: yIx
    REAL(wp), DIMENSION(ySz) :: Vaut0
    REAL(wp) :: errV

    Vaut = 0.0_wp

    errV = 1.0_wp
    DO WHILE (errV > epsV)
      Vaut0 = Vaut
      Vaut = (1.0_wp - beta) * uhy + beta * MATMUL( exoTran, Vaut0 )
      errV = MAXVAL(ABS(Vaut - Vaut0))
    END DO
  END SUBROUTINE

  SUBROUTINE saveEquilibrium()
    IMPLICIT NONE
    OPEN(unit=100, file=outDir // "parameters.txt")
    WRITE (100, "(I4)") ySz
    WRITE (100, "(I4)") bSSz
    WRITE (100, "(I4)") bLSz
    WRITE (100, "(I4)") chPrSz
    WRITE (100, "(F9.6)") crra
    WRITE (100, "(F9.6)") beta
    WRITE (100, "(F9.6)") rf
    WRITE (100, "(F9.6)") eta
    WRITE (100, "(F9.6)") etaA
    WRITE (100, "(F9.6)") deltaS
    WRITE (100, "(F9.6)") deltaL
    WRITE (100, "(F9.6)") alpha
    WRITE (100, "(F9.6)") muS
    WRITE (100, "(F9.6)") lbd0
    WRITE (100, "(F9.6)") lbd1
    WRITE (100, "(F15.10)") rho
    WRITE (100, "(F15.10)") rhoDef
    WRITE (100, "(F15.10)") rhoNash
    WRITE (100, "(F9.6)") adj
    WRITE (100, "(F9.6)") bSshare
    CLOSE(100)

    OPEN(unit=100, file=outDir // "bS.txt")
    WRITE (100, "(F9.6)") bS
    CLOSE(100)

    OPEN(unit=100, file=outDir // "bL.txt")
    WRITE (100, "(F9.6)") bL
    CLOSE(100)

    OPEN(unit=100, file=outDir // "V.txt")
    WRITE (100, "(F16.6)") V1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "Va.txt")
    WRITE (100, "(F16.6)") Va1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "qL.txt")
    WRITE (100, "(F9.6)") qL1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "qS.txt")
    WRITE (100, "(F9.6)") qS1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "qLa.txt")
    WRITE (100, "(F9.6)") qLa1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "qSa.txt")
    WRITE (100, "(F9.6)") qSa1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "xiL.txt")
    WRITE (100, "(F9.6)") xiL
    CLOSE(100)

    OPEN(unit=100, file=outDir // "xiS.txt")
    WRITE (100, "(F9.6)") xiS
    CLOSE(100)

    OPEN(unit=100, file=outDir // "defPr.txt")
    WRITE (100, "(F9.6)") defPr
    CLOSE(100)

    OPEN(unit=100, file=outDir // "Vd.txt")
    WRITE (100, "(F16.6)") Vd1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "Vaut.txt")
    WRITE (100, "(F16.6)") Vaut
    CLOSE(100)

    OPEN(unit=100, file=outDir // "gammaPr.txt")
    WRITE (100, "(F9.6)") gammaPr
    CLOSE(100)

    OPEN(unit=100, file=outDir // "gammaS.txt")
    WRITE (100, "(I4)") gammaS
    CLOSE(100)

    OPEN(unit=100, file=outDir // "gammaL.txt")
    WRITE (100, "(I4)") gammaL
    CLOSE(100)

    OPEN(unit=100, file=outDir // "chPr.txt")
    WRITE (100, "(F9.6)") chPr
    CLOSE(100)

    OPEN(unit=100, file=outDir // "chS.txt")
    WRITE (100, "(I4)") chS
    CLOSE(100)

    OPEN(unit=100, file=outDir // "chL.txt")
    WRITE (100, "(I4)") chL
    CLOSE(100)

    OPEN(unit=100, file=outDir // "chPrA.txt")
    WRITE (100, "(F9.6)") chPrA
    CLOSE(100)

    OPEN(unit=100, file=outDir // "chSa.txt")
    WRITE (100, "(I4)") chSa
    CLOSE(100)

    OPEN(unit=100, file=outDir // "chLa.txt")
    WRITE (100, "(I4)") chLa
    CLOSE(100)

    OPEN(unit=100, file=outDir // "EdCre.txt")
    WRITE (100, "(F9.6)") EdCre
    CLOSE(100)
  END SUBROUTINE saveEquilibrium

  ! Load arrays from .txt files
  SUBROUTINE loadEquilibrium()
    IMPLICIT NONE

    IF (workerId == 0) WRITE (*, *) "Loading policies and values from files."
    IF (workerId == 0) WRITE (*, *) "Make sure the grids and parameters are the same."

    OPEN(unit=100, file=outDir // "V.txt", status="old", position="rewind")
    READ (100, "(F16.6)") V1
    V0 = V1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "Va.txt", status="old", position="rewind")
    READ (100, "(F16.6)") Va1
    Va0 = Va1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "qL.txt", status="old", position="rewind")
    READ (100, "(F9.6)") qL1
    qL0 = qL1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "qS.txt", status="old", position="rewind")
    READ (100, "(F9.6)") qS1
    qS0 = qS1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "qLa.txt", status="old", position="rewind")
    READ (100, "(F9.6)") qLa1
    qLa0 = qLa1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "qSa.txt", status="old", position="rewind")
    READ (100, "(F9.6)") qSa1
    qSa0 = qSa1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "xiL.txt", status="old", position="rewind")
    READ (100, "(F9.6)") xiL
    CLOSE(100)

    OPEN(unit=100, file=outDir // "xiS.txt", status="old", position="rewind")
    READ (100, "(F9.6)") xiS
    CLOSE(100)

    OPEN(unit=100, file=outDir // "defPr.txt", status="old", position="rewind")
    READ (100, "(F9.6)") defPr
    CLOSE(100)

    OPEN(unit=100, file=outDir // "Vd.txt", status="old", position="rewind")
    READ (100, "(F16.6)") Vd1
    Vd0 = Vd1
    CLOSE(100)

    OPEN(unit=100, file=outDir // "gammaPr.txt", status="old", position="rewind")
    READ (100, "(F9.6)") gammaPr
    CLOSE(100)

    OPEN(unit=100, file=outDir // "gammaS.txt", status="old", position="rewind")
    READ (100, "(I4)") gammaS
    CLOSE(100)

    OPEN(unit=100, file=outDir // "gammaL.txt", status="old", position="rewind")
    READ (100, "(I4)") gammaL
    CLOSE(100)

    OPEN(unit=100, file=outDir // "chPr.txt", status="old", position="rewind")
    READ (100, "(F9.6)") chPr
    CLOSE(100)

    OPEN(unit=100, file=outDir // "chS.txt", status="old", position="rewind")
    READ (100, "(I4)") chS
    CLOSE(100)

    OPEN(unit=100, file=outDir // "chL.txt", status="old", position="rewind")
    READ (100, "(I4)") chL
    CLOSE(100)

    OPEN(unit=100, file=outDir // "chPrA.txt", status="old", position="rewind")
    READ (100, "(F9.6)") chPrA
    CLOSE(100)

    OPEN(unit=100, file=outDir // "chSa.txt", status="old", position="rewind")
    READ (100, "(I4)") chSa
    CLOSE(100)

    OPEN(unit=100, file=outDir // "chLa.txt", status="old", position="rewind")
    READ (100, "(I4)") chLa
    CLOSE(100)

    OPEN(unit=100, file=outDir // "EdCre.txt", status="old", position="rewind")
    READ (100, "(F9.6)") EdCre
    CLOSE(100)
  END SUBROUTINE loadEquilibrium

  ! Simulate
  SUBROUTINE simulate()
    IMPLICIT NONE
    INTEGER, ALLOCATABLE, DIMENSION(:) :: yPath, bSpath, bLpath, dPath
    REAL(wp), ALLOCATABLE, DIMENSION(:) :: cPath, gdpPath, qSpath, qLpath
    INTEGER, ALLOCATABLE, DIMENSION(:) :: yEv, yDefEv, bSev, bLev, gammaSev, gammaLev
    INTEGER :: tt, evCnt, ix, drawCh, yDefaulted, returnS, returnL
    REAL(wp) :: rno, ddPr
    LOGICAL :: inDef, inArr, justReturn

    ALLOCATE(yPath(simSz), bSpath(simSz), bLpath(simSz), dPath(simSz))
    ALLOCATE(cPath(simSz), gdpPath(simSz), qSpath(simSz), qLpath(simSz))
    ALLOCATE(yEv(simSz), yDefEv(simSz),bSev(simSz), bLev(simSz), gammaSev(simSz), gammaLev(simSz))

    CALL fixSeed()

    yPath = 1
    bSpath = 1
    bLpath = 1
    dPath = 0
    cPath = 0.0_wp
    gdpPath = 0.0_wp
    qSpath = 0.0_wp
    qLpath = 0.0_wp

    yDefaulted = 1
    evCnt = 1

    DO tt = 2,simSz-1

      CALL simMarkov(ySz, exoTran, yPath(tt), yPath(tt-1))
      justReturn = .FALSE.
      returnS = 1
      returnL = 1

      IF ( dPath(tt-1) == 2 ) THEN ! Prev in arr
        CALL simUniform(rno)
        IF (rno < etaA) THEN ! To market
          dPath(tt) = 0
          ! Check Default
          ! Use chPr later
        ELSE ! Stay
          dPath(tt) = 2
          ! Use chPrA
          CALL simDiscrete(chPrSz, chPrA(yPath(tt), bSpath(tt-1), bLpath(tt-1), :), drawCh)
          bSpath(tt) = chSa(yPath(tt), bSpath(tt-1), bLpath(tt-1), drawCh)
          bLpath(tt) = chLa(yPath(tt), bSpath(tt-1), bLpath(tt-1), drawCh)
        END IF
      ELSE IF ( dPath(tt-1) == 1 ) THEN ! Prev in default
        CALL simUniform(rno)
        IF (rno < eta) THEN ! To arr
          dPath(tt) = 2
          ! Use gamma-s
          CALL simDiscrete(chPrSz, gammaPr(yPath(tt), :), drawCh)
          ! bSpath(tt) = gammaS(yPath(tt), drawCh)
          ! bLpath(tt) = gammaL(yPath(tt), drawCh)
          returnS = gammaS(yPath(tt), drawCh)
          returnL = gammaL(yPath(tt), drawCh)

          yEv(evCnt) = yPath(tt)
          yDefEv(evCnt) = yDefaulted
          bSev(evCnt) = bSpath(tt-1)
          bLev(evCnt) = bLpath(tt-1)
          gammaSev(evCnt) = returnS
          gammaLev(evCnt) = returnL
          evCnt = evCnt + 1

          justReturn = .TRUE.

          CALL simDiscrete(chPrSz, chPrA(yPath(tt), returnS, returnL, :), drawCh)
          bSpath(tt) = chSa(yPath(tt), returnS, returnL, drawCh)
          bLpath(tt) = chLa(yPath(tt), returnS, returnL, drawCh)
        ELSE ! Stay
          dPath(tt) = 1
          bSpath(tt) = bSpath(tt-1)
          bLpath(tt) = bLpath(tt-1)
        END IF
      ELSE
        ! CALL simUniform(rno)
        dPath(tt) = 0
        ! Check Default
        ! Use chPr later
      END IF

      IF (dPath(tt) == 0) THEN
          CALL simUniform(rno)
          ddPr = defPr(yPath(tt), bSpath(tt-1), bLpath(tt-1))
          IF (rho < ddPr) THEN
            ! Choose default
            dPath(tt) = 1
            yDefaulted = yPath(tt)
            bSpath(tt) = bSpath(tt-1)
            bLpath(tt) = bLpath(tt-1)
          ELSE
            CALL simDiscrete(chPrSz, chPr(yPath(tt), bSpath(tt-1), bLpath(tt-1), :) / (1.0_wp - ddPr), drawCh)
            bSpath(tt) = chS(yPath(tt), bSpath(tt-1), bLpath(tt-1), drawCh)
            bLpath(tt) = chL(yPath(tt), bSpath(tt-1), bLpath(tt-1), drawCh)
          END IF
      END IF

      IF (dPath(tt) == 1) THEN ! Default
        gdpPath(tt) = hy(yPath(tt))
        cPath(tt) = gdpPath(tt)
        qSpath(tt) = xiS(yPath(tt), bSpath(tt), bLpath(tt))
        qLpath(tt) = xiL(yPath(tt), bSpath(tt), bLpath(tt))
      ELSE IF (dPath(tt) == 2) THEN ! Arrears
        gdpPath(tt) = y(yPath(tt))
        qSpath(tt) = qSa1(yPath(tt), bSpath(tt), bLpath(tt))
        qLpath(tt) = qLa1(yPath(tt), bSpath(tt), bLpath(tt))
        IF (justReturn) THEN
          cPath(tt) = gdpPath(tt) - kS * bS(returnS) - kL * bL(returnL) &
            + qSpath(tt) * ( bS(bSpath(tt)) - (1.0_wp - deltaS) * bS(returnS) ) &
            + qLpath(tt) * ( bL(bLpath(tt)) - (1.0_wp - deltaL) * bL(returnL) )
        ELSE
          cPath(tt) = gdpPath(tt) - kS * bS(bSpath(tt-1)) - kL * bL(bLpath(tt-1)) &
            + qSpath(tt) * ( bS(bSpath(tt)) - (1.0_wp - deltaS) * bS(bSpath(tt-1)) ) &
            + qLpath(tt) * ( bL(bLpath(tt)) - (1.0_wp - deltaL) * bL(bLpath(tt-1)) )
          !IF ( bS(bSpath(tt)) + bL(bLpath(tt)) > 0.0_wp ) THEN
          !    cPath(tt) = cPath(tt) - adj * ( &
          !      bS(bSpath(tt)) / (bS(bSpath(tt)) + bL(bLpath(tt))) - bSshare )**2
          !END IF
        END IF
      ELSE
        gdpPath(tt) = y(yPath(tt))
        qSpath(tt) = qS1(yPath(tt), bSpath(tt), bLpath(tt))
        qLpath(tt) = qL1(yPath(tt), bSpath(tt), bLpath(tt))
        cPath(tt) = gdpPath(tt) - kS * bS(bSpath(tt-1)) - kL * bL(bLpath(tt-1)) &
          + qSpath(tt) * ( bS(bSpath(tt)) - (1.0_wp - deltaS) * bS(bSpath(tt-1)) ) &
          + qLpath(tt) * ( bL(bLpath(tt)) - (1.0_wp - deltaL) * bL(bLpath(tt-1)) )
        IF ( bS(bSpath(tt)) + bL(bLpath(tt)) > 0.0_wp ) THEN
            cPath(tt) = cPath(tt) - adj * ( &
              bS(bSpath(tt)) / (bS(bSpath(tt)) + bL(bLpath(tt))) - bSshare )**2
        END IF
      END IF
    END DO

    OPEN(unit=100, file=outDir // "events.tab")
    DO ix = 1,evCnt-1
      WRITE (100, "(7(I10,A))") ix, TAB, yEv(ix), TAB, bSev(ix), TAB, &
        bLev(ix), TAB, gammaSev(ix), TAB, gammaLev(ix), TAB, yDefEv(ix), TAB
    END DO
    CLOSE(100)

    OPEN(unit=100, file=outDir // "simulation.tab")
    WRITE (100, "(A9,A,3(A3,A,A10,A),A1,A,4(A10,A))") "T", TAB, "yIx", TAB, "y", TAB, "bSix", TAB, &
      "bS", TAB, "bLix", TAB, "bL", TAB, "d", TAB, "qS", TAB, "qL", TAB, "gdp", TAB, "c", TAB
    DO ix = 5000,simSz-500
      WRITE (100, "(I9,A,3(I3,A,F10.6,A),I1,A,4(F10.6,A))") ix, TAB, yPath(ix), TAB, y(yPath(ix)), TAB, &
        bSpath(ix), TAB, bS(bSpath(ix)), TAB, bLpath(ix), TAB, bL(bLpath(ix)), TAB, &
        dPath(ix), TAB, qSpath(ix), TAB, qLpath(ix), TAB, gdpPath(ix), TAB, cPath(ix), TAB
    END DO
    CLOSE(100)

!    OPEN(unit=100, file=outDir // "qSpath.txt")
!    WRITE (100, "(F9.6)") qSpath(10:simSz-1)
!    CLOSE(100)
!
!    OPEN(unit=100, file=outDir // "qLpath.txt")
!    WRITE (100, "(F9.6)") qLpath(10:simSz-1)
!    CLOSE(100)
!
!    OPEN(unit=100, file=outDir // "cPath.txt")
!    WRITE (100, "(F9.6)") cPath(10:simSz-1)
!    CLOSE(100)
!
!    OPEN(unit=100, file=outDir // "gdpPath.txt")
!    WRITE (100, "(F9.6)") gdpPath(10:simSz-1)
!    CLOSE(100)
!
!    OPEN(unit=100, file=outDir // "yPath.txt")
!    WRITE (100, "(I4)") yPath(10:simSz-1)
!    CLOSE(100)
!
!    OPEN(unit=100, file=outDir // "bSpath.txt")
!    WRITE (100, "(I4)") bSpath(10:simSz-1)
!    CLOSE(100)
!
!    OPEN(unit=100, file=outDir // "bLpath.txt")
!    WRITE (100, "(I4)") bLpath(10:simSz-1)
!    CLOSE(100)
!
!    OPEN(unit=100, file=outDir // "dPath.txt")
!    WRITE (100, "(I4)") dPath(10:simSz-1)
!    CLOSE(100)
  END SUBROUTINE simulate

END MODULE jmpMod

max_stars_repo_name: gabrielgggg/MaturityExtension

max_issues_repo_name: gabrielgggg/MaturityExtension

max_forks_repo_name: gabrielgggg/MaturityExtension


text:
!==============================================================================!
  subroutine Cgns_Mod_Initialize_Counters
!------------------------------------------------------------------------------!
!   Initializes counters used while CGNS file is being read.                   !
!------------------------------------------------------------------------------!
  implicit none
!==============================================================================!

  cnt_nodes    = 0  ! number of nodes (except boundary and interface)
  cnt_cells    = 0  ! number of cells (except boundary and interface)
  cnt_blocks   = 0  ! number of block
  cnt_bnd_cond = 0  ! number of boundary

  cnt_hex = 0  ! number of hexahedral  cells
  cnt_pyr = 0  ! number of pyramid     cells
  cnt_wed = 0  ! number of wedge/prism cells
  cnt_tet = 0  ! number of tetrahedral cells
  cnt_mix = 0  ! number of mixed       cells

  cnt_bnd_tri = 0  ! number of triangles defining boundary conditions
  cnt_bnd_qua = 0  ! number of quads     defining boundary conditions

  cnt_int       = 0  ! number of unique interfaces
  cnt_int_tri   = 0  ! number of triangles defining block interface
  cnt_int_qua   = 0  ! number of quads     defining block interface
  cnt_int_cells = 0  ! number of cells on interfaces

  pos_of_last_3d_cell = 0  ! position of last 3d cell

  end subroutine

max_stars_repo_name: EdinSmartLab/T-Flows

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


program jsonio_test
  use jsonio, IK => json_IK, RK => json_RK
  use iso_fortran_env, only: int64
  implicit none

  logical, parameter :: debug = .false.

  character(len=128), parameter :: jsonfile = 'jsonio_test.json'
  integer(int64), parameter :: large_int64 = (2_8)**32
  integer, parameter :: endian_flag = 1
  integer, parameter :: ndim = 2
  integer, parameter :: nx = 16
  integer, parameter :: ny = 32
  integer, parameter :: ns = 2
  integer, parameter :: shape_emf(ndim+1) = (/6_4, ny, nx/)
  integer, parameter :: shape_mom(ndim+1) = (/10_4, ny, nx/)
  real(8), parameter :: mass(ns) = (/1.0_8, 100.0_8/)
  real(8), parameter :: charge(ns) = (/-1.0_8, +1.0_8/)

  integer, parameter :: num_data_max = 32
  character(len=16)  :: json_name(num_data_max)
  character(len=16)  :: json_datatype(num_data_max)
  integer(int64)     :: json_offset(num_data_max)
  integer(int64)     :: json_size(num_data_max)
  integer(4)         :: json_ndim(num_data_max)
  integer(4)         :: json_shape(ndim+1,num_data_max)

  call write_json(jsonfile)
  call read_json(jsonfile)

contains

  ! store metadata
  subroutine store_metadata(n, name, datatype, offset, dsize, ndim, dshape)
    integer, intent(in)          :: n
    character(len=*), intent(in) :: name
    character(len=*), intent(in) :: datatype
    integer(int64), intent(in)   :: offset
    integer(int64), intent(in)   :: dsize
    integer(4), intent(in)       :: ndim
    integer(4), intent(in)       :: dshape(ndim)

    json_name(n)         = trim(name)
    json_datatype(n)     = trim(datatype)
    json_offset(n)       = offset
    json_size(n)         = dsize
    json_ndim(n)         = ndim
    json_shape(1:ndim,n) = dshape

  end subroutine store_metadata

  ! check metadata
  subroutine check_metadata(n, name, datatype, offset, dsize, ndim, dshape)
    integer, intent(in)          :: n
    character(len=*), intent(in) :: name
    character(len=*), intent(in) :: datatype
    integer(int64), intent(in)   :: offset
    integer(int64), intent(in)   :: dsize
    integer(4), intent(in)       :: ndim
    integer(4), intent(in)       :: dshape(ndim)

    integer :: i, errcnt = 0

    if( .not. trim(json_datatype(n)) == trim(datatype) ) then
       errcnt = errcnt + 1
    end if

    if( .not. json_offset(n) == offset ) then
       errcnt = errcnt + 1
    end if

    if( .not. json_size(n) == dsize ) then
       errcnt = errcnt + 1
    end if

    if( .not. json_ndim(n) == ndim ) then
       errcnt = errcnt + 1
    end if

    if( .not. all(json_shape(1:ndim,n) == dshape) ) then
       errcnt = errcnt + 1
    end if

    !
    ! report error
    !
    if( errcnt > 0 ) then
       write(0,'("*** Some errors detected for ", A)') trim(name)
       write(0,'(A10, " : ", A10, " <===> ", A10)') &
            & 'datatype', trim(json_datatype(n)), trim(datatype)
       write(0,'(A10, " : ", i10, " <===> ", i10)') &
            & 'offset', json_offset(n), offset
       write(0,'(A10, " : ", i10, " <===> ", i10)') &
            & 'size', json_size(n), dsize
       write(0,'(A10, " : ", i10, " <===> ", i10)') &
            & 'ndim', json_ndim(n), ndim
       ! shape
       write(0,'(A10, " : ")', advance='no') 'shape'
       write(0,'(A)', advance='no') '['
       do i = 1, ndim
          write(0,'(i6)', advance='no') json_shape(i,n)
       end do
       write(0,'(A)', advance='no') ' ] <===> ['
       do i = 1, ndim
          write(0,'(i6)', advance='no') dshape(i)
       end do
       write(0,'(A)', advance='no') ' ]'
       write(0,*)
    else if ( debug ) then
       write(0,'("No error detected for ", A)') trim(name)
    end if

  end subroutine check_metadata

  ! write a json file
  subroutine write_json(filename)
    implicit none
    character(len=*), intent(in) :: filename

    type(json_core) :: json
    type(json_value), pointer :: root, p

    integer(int64)     :: offset, dsize, large_int
    integer            :: numdata, endian, ndim, dshape(3)
    character(len=128) :: desc

    endian  = endian_flag
    numdata = 1

    call json%initialize()
    call json%create_object(root, 'root')

    !
    ! meta
    !
    call json%create_object(p, 'meta')
    call jsonio_check_error(json, 'write_json')

    call json%add(root, p)
    call jsonio_check_error(json, 'write_json')

    call json%add(p, 'endian', endian)
    call jsonio_check_error(json, 'write_json')

    call json%add(p, 'rawfile', 'test.raw')
    call jsonio_check_error(json, 'write_json')

    call json%add(p, 'info', 'some information')
    call jsonio_check_error(json, 'write_json')

    ! check for 64bit integer
    large_int = large_int64
    call json%add(p, 'large_int', large_int)

    nullify(p)

    !
    ! attributes
    !
    call json%create_object(p, 'attribute')
    call jsonio_check_error(json, 'write_json')

    call json%add(root, p)
    call jsonio_check_error(json, 'write_json')

    ! nx
    offset = 0
    dsize  = 4
    ndim   = 1
    dshape = (/1, 0, 0/)
    desc   = '# grid in x'
    call jsonio_put_attribute(json, p, nx, 'nx', offset, desc)
    call store_metadata(numdata, 'nx', 'i4', offset, dsize, ndim, dshape)
    numdata = numdata + 1

    ! ny
    offset = offset + dsize
    dsize  = 4
    ndim   = 1
    dshape = (/1, 0, 0/)
    desc   = '# grid in y'
    call jsonio_put_attribute(json, p, ny, 'ny', offset, desc)
    call store_metadata(numdata, 'ny', 'i4', offset, dsize, ndim, dshape)
    numdata = numdata + 1

    ! ns
    offset = offset + dsize
    dsize  = 4
    ndim   = 1
    dshape = (/1, 0, 0/)
    desc   = '# species'
    call jsonio_put_attribute(json, p, ns, 'ns', offset, desc)
    call store_metadata(numdata, 'ns', 'i4', offset, dsize, ndim, dshape)
    numdata = numdata + 1

    ! mass
    offset = offset + dsize
    dsize  = 8 * 2
    ndim   = 1
    dshape = (/2, 0, 0/)
    desc   = 'mass'
    call jsonio_put_attribute(json, p, mass, 'mass', offset, desc)
    call store_metadata(numdata, 'mass', 'f8', offset, dsize, ndim, dshape)
    numdata = numdata + 1

    ! charge
    offset = offset + dsize
    dsize  = 8 * 2
    ndim   = 1
    dshape = (/2, 0, 0/)
    desc   = 'charge'
    call jsonio_put_attribute(json, p, charge, 'charge', offset, desc)
    call store_metadata(numdata, 'charge', 'f8', offset, dsize, ndim, dshape)
    numdata = numdata + 1

    nullify(p)

    !
    ! datasets
    !
    call json%create_object(p, 'dataset')
    call jsonio_check_error(json, 'write_json')

    call json%add(root, p)
    call jsonio_check_error(json, 'write_json')

    ! emf
    offset = offset + dsize
    dsize  = product(shape_emf) * 8
    ndim   = 3
    dshape = shape_emf
    desc   = 'electromagnetic fields'
    call jsonio_put_metadata(json, p, 'emf', 'f8', offset, &
         & dsize, ndim, dshape, desc)
    call store_metadata(numdata, 'emf', 'f8', offset, dsize, &
         & ndim, dshape)
    numdata = numdata + 1

    ! mom
    offset = offset + dsize
    dsize  = product(shape_mom) * 8
    ndim   = 3
    dshape = shape_mom
    desc   = 'moments'
    call jsonio_put_metadata(json, p, 'mom', 'f8', offset, &
         & dsize, ndim, dshape, desc)
    call store_metadata(numdata, 'mom', 'f8', offset, dsize, &
         & ndim, dshape)
    numdata = numdata + 1

    nullify(p)

    ! output
    call json%print(root, filename)

    nullify(root)
    call json%destroy()

  end subroutine write_json

  ! read and check a json file
  subroutine read_json(filename)
   implicit none
    character(len=*), intent(in) :: filename

    type(json_file) :: file
    type(json_core) :: json
    type(json_value), pointer :: root, p

    integer(int64)     :: offset, dsize, large_int
    integer            :: numdata, endian, ndim, dshape(3)
    integer            :: nx, ny, ns
    real(8)            :: mass(2), charge(2)

    numdata = 1

    call json%initialize()
    call jsonio_check_error(json, 'read_json')

    call file%initialize()
    call jsonio_check_error(file, 'read_json')

    call file%load(filename)
    call jsonio_check_error(file, 'read_json')

    ! get root
    call file%get(root)
    call jsonio_check_error(file, 'read_json')

    !
    ! meta
    !
    call json%get(root, 'meta', p)
    call json%get(p, 'endian', endian)
    call json%get(p, 'large_int', large_int)

    if( large_int /= large_int64) then
       write(0,*) 'Error in reading 64bit integer: ', large_int
    end if

    !
    ! attributes
    !
    call json%get(root, 'attribute', p)
    call jsonio_check_error(json, 'read_json')

    call jsonio_get_attribute(json, p, 'nx', offset, nx)
    call jsonio_get_metadata(json, p, 'nx', offset, dsize, ndim, dshape)
    call check_metadata(numdata, 'nx', 'i4', offset, dsize, ndim, dshape)
    numdata = numdata + 1

    call jsonio_get_attribute(json, p, 'ny', offset, ny)
    call jsonio_get_metadata(json, p, 'ny', offset, dsize, ndim, dshape)
    call check_metadata(numdata, 'ny', 'i4', offset, dsize, ndim, dshape)
    numdata = numdata + 1

    call jsonio_get_attribute(json, p, 'ns', offset, ns)
    call jsonio_get_metadata(json, p, 'ns', offset, dsize, ndim, dshape)
    call check_metadata(numdata, 'ns', 'i4', offset, dsize, ndim, dshape)
    numdata = numdata + 1

    call jsonio_get_attribute(json, p, 'mass', offset, mass)
    call jsonio_get_metadata(json, p, 'mass', offset, dsize, ndim, dshape)
    call check_metadata(numdata, 'mass', 'f8', offset, dsize, ndim, dshape)
    numdata = numdata + 1

    call jsonio_get_attribute(json, p, 'charge', offset, charge)
    call jsonio_get_metadata(json, p, 'charge', offset, dsize, ndim, dshape)
    call check_metadata(numdata, 'charge', 'f8', offset, dsize, ndim, dshape)
    numdata = numdata + 1

    nullify(p)

    !
    ! datasets
    !
    call json%get(root, 'dataset', p)
    call jsonio_check_error(json, 'read_json')

    call jsonio_get_metadata(json, p, 'emf', offset, dsize, ndim, dshape)
    call check_metadata(numdata, 'emf', 'f8', offset, dsize, ndim, dshape)
    numdata = numdata + 1

    call jsonio_get_metadata(json, p, 'mom', offset, dsize, ndim, dshape)
    call check_metadata(numdata, 'mom', 'f8', offset, dsize, ndim, dshape)
    numdata = numdata + 1

    nullify(root)
    nullify(p)

    call json%destroy()
    call file%destroy()

  end subroutine read_json

end program jsonio_test


text:
module a

integer::some = 1

contains
  subroutine b0(val1)
    integer,intent(in)::val1
    print *,some
    print *, val1
  end subroutine b0

  subroutine b1(val)
    print *, val
    call b2(val)

    contains

      subroutine b2(val2)
        print *,val2
      end subroutine b2

  end subroutine b1

end module a

program vin

use a

call internal(13.413)

contains

 subroutine internal(val)
  call b0(10)
  call b1(val)
 end subroutine internal

end program vin