text:
      SUBROUTINE iau_FK54Z ( R2000, D2000, BEPOCH,
     :                       R1950, D1950, DR1950, DD1950 )
*+
*  - - - - - - - - - -
*   i a u _ F K 5 4 Z
*  - - - - - - - - - -
*
*  Convert a J2000.0 FK5 star position to B1950.0 FK4, assuming zero
*  proper motion in FK5 and parallax.
*
*  This routine is part of the International Astronomical Union's
*  SOFA (Standards of Fundamental Astronomy) software collection.
*
*  Status:  support routine.
*
*  Given:
*     R2000,D2000      d     J2000.0 FK5 RA,Dec (rad)
*     BEPOCH           d     Besselian epoch (e.g. 1950D0)
*
*  Returned:
*     R1950,D1950      d     B1950.0 FK4 RA,Dec (rad) at epoch BEPOCH
*     DR1950,DD1950    d     B1950.0 FK4 proper motions (rad/trop.yr)
*
*  Notes:
*
*  1) In contrast to the iau_FK524 routine, here the FK5 proper motions,
*     the parallax and the radial velocity are presumed zero.
*
*  2) This routine converts a star position from the IAU 1976 FK5
*     (Fricke) system to the former FK4 (Bessel-Newcomb) system, for
*     cases such as distant radio sources where it is presumed there is
*     zero parallax and no proper motion.  Because of the E-terms of
*     aberration, such objects have (in general) non-zero proper motion
*     in FK4, and the present routine returns those fictitious proper
*     motions.
*
*  3) Conversion from B1950.0 FK4 to J2000.0 FK5 only is provided for.
*     Conversions involving other equinoxes would require additional
*     treatment for precession.
*
*  4) The position returned by this routine is in the B1950.0 FK4
*     reference system but at Besselian epoch BEPOCH.  For comparison
*     with catalogs the BEPOCH argument will frequently be 1950D0.  (In
*     this context the distinction between Besselian and Julian epoch is
*     insignificant.)
*
*  5) The RA component of the returned (fictitious) proper motion is
*     dRA/dt rather than cos(Dec)*dRA/dt.
*
*  Called:
*     iau_ANP      normalize angle into range 0 to 2pi
*     iau_C2S      p-vector to spherical
*     iau_FK524    FK4 to FK5
*     iau_S2C      spherical to p-vector
*
*  This revision:   2018 January 11
*
*  SOFA release 2020-07-21
*
*  Copyright (C) 2020 IAU SOFA Board.  See notes at end.
*
*-----------------------------------------------------------------------

      IMPLICIT NONE

      DOUBLE PRECISION R2000, D2000, BEPOCH,
     :                 R1950, D1950, DR1950, DD1950

      DOUBLE PRECISION R, D, PR, PD, PX, RV, P(3), W, V(3)
      INTEGER I

      DOUBLE PRECISION iau_ANP

* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

*  FK5 equinox J2000.0 to FK4 equinox B1950.0.
      CALL iau_FK524 ( R2000, D2000, 0D0, 0D0, 0D0, 0D0,
     :                 R, D, PR, PD, PX, RV )

*  Spherical to Cartesian.
      CALL iau_S2C ( R, D, P )

*  Fictitious proper motion (radians per year).
      V(1) = - PR*P(2) - PD*COS(R)*SIN(D)
      V(2) =   PR*P(1) - PD*SIN(R)*SIN(D)
      V(3) =             PD*COS(D)

*  Apply the motion.
      W = BEPOCH - 1950D0
      DO 1 I=1,3
         P(I) = P(I) + W*V(I)
 1    CONTINUE

*  Cartesian to spherical.
      CALL iau_C2S ( P, W, D1950 )
      R1950 = iau_ANP ( W )

*  Fictitious proper motion.
      DR1950 = PR
      DD1950 = PD

*  Finished.

*+----------------------------------------------------------------------
*
*  Copyright (C) 2020
*  Standards Of Fundamental Astronomy Board
*  of the International Astronomical Union.
*
*  =====================
*  SOFA Software License
*  =====================
*
*  NOTICE TO USER:
*
*  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING SIX TERMS AND
*  CONDITIONS WHICH APPLY TO ITS USE.
*
*  1. The Software is owned by the IAU SOFA Board ("SOFA").
*
*  2. Permission is granted to anyone to use the SOFA software for any
*     purpose, including commercial applications, free of charge and
*     without payment of royalties, subject to the conditions and
*     restrictions listed below.
*
*  3. You (the user) may copy and distribute SOFA source code to others,
*     and use and adapt its code and algorithms in your own software,
*     on a world-wide, royalty-free basis.  That portion of your
*     distribution that does not consist of intact and unchanged copies
*     of SOFA source code files is a "derived work" that must comply
*     with the following requirements:
*
*     a) Your work shall be marked or carry a statement that it
*        (i) uses routines and computations derived by you from
*        software provided by SOFA under license to you; and
*        (ii) does not itself constitute software provided by and/or
*        endorsed by SOFA.
*
*     b) The source code of your derived work must contain descriptions
*        of how the derived work is based upon, contains and/or differs
*        from the original SOFA software.
*
*     c) The names of all routines in your derived work shall not
*        include the prefix "iau" or "sofa" or trivial modifications
*        thereof such as changes of case.
*
*     d) The origin of the SOFA components of your derived work must
*        not be misrepresented;  you must not claim that you wrote the
*        original software, nor file a patent application for SOFA
*        software or algorithms embedded in the SOFA software.
*
*     e) These requirements must be reproduced intact in any source
*        distribution and shall apply to anyone to whom you have
*        granted a further right to modify the source code of your
*        derived work.
*
*     Note that, as originally distributed, the SOFA software is
*     intended to be a definitive implementation of the IAU standards,
*     and consequently third-party modifications are discouraged.  All
*     variations, no matter how minor, must be explicitly marked as
*     such, as explained above.
*
*  4. You shall not cause the SOFA software to be brought into
*     disrepute, either by misuse, or use for inappropriate tasks, or
*     by inappropriate modification.
*
*  5. The SOFA software is provided "as is" and SOFA makes no warranty
*     as to its use or performance.   SOFA does not and cannot warrant
*     the performance or results which the user may obtain by using the
*     SOFA software.  SOFA makes no warranties, express or implied, as
*     to non-infringement of third party rights, merchantability, or
*     fitness for any particular purpose.  In no event will SOFA be
*     liable to the user for any consequential, incidental, or special
*     damages, including any lost profits or lost savings, even if a
*     SOFA representative has been advised of such damages, or for any
*     claim by any third party.
*
*  6. The provision of any version of the SOFA software under the terms
*     and conditions specified herein does not imply that future
*     versions will also be made available under the same terms and
*     conditions.
*
*  In any published work or commercial product which uses the SOFA
*  software directly, acknowledgement (see www.iausofa.org) is
*  appreciated.
*
*  Correspondence concerning SOFA software should be addressed as
*  follows:
*
*      By email:  sofa@ukho.gov.uk
*      By post:   IAU SOFA Center
*                 HM Nautical Almanac Office
*                 UK Hydrographic Office
*                 Admiralty Way, Taunton
*                 Somerset, TA1 2DN
*                 United Kingdom
*
*-----------------------------------------------------------------------

      END

ext: for

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
C    @(#)swap_buf.f	20.3 2/13/96
	subroutine swap_buf (m, n)

	common /shellsort/ buffer(1000)
        integer buffer

        junk = buffer(m)
        buffer(m) = buffer(n)
        buffer(n) = junk

        return
        end

ext: f

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
      double complex function BSYA0ggppmp(e1,p2,p3,e4,za,zb,zab,zba)
      implicit none
C-----Authors: John Campbell and Keith Ellis, March 2012
C---- arXiv:1101.5947 [hep-ph], Eq. (A6), fully BSY compliant
C---- (These are twiddle functions, c.f.arXiv:1101.5947[hep-ph],Eq.(91))
      include 'constants.f'
      include 'masses.f'
      include 'zprods_decl.f'
      include 'zabprods_decl.f'
      include 'sprods_com.f'
      include 'qdef.f'
      integer e1,p2,p3,e4
      BSYA0ggppmp=mt*zab(p3,q1,p2)
     & *(za(e1,e4)*zab(p3,q1,p2)-zb(p2,p3)*za(p3,e1)*za(p3,e4))
     & /(s(p2,p3)*zab(p2,q1,p2))

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


      subroutine soiltest_all_init

      use hru_module, only : hru, wfsh, ihru, sol_plt_ini
      use soil_module
      use plant_module
      use maximum_data_module
      use soil_data_module
      use hydrograph_module, only : sp_ob

      implicit none

      integer :: isolt            !           |
      integer :: isol_pl          !           |

      do ihru = 1, sp_ob%hru
        isol_pl = hru(ihru)%dbs%soil_plant_init
        isolt = sol_plt_ini(isol_pl)%nut
        if (isolt > 0) then
          call soiltest_init (ihru, isolt)
        end if
      end do

      return
      end subroutine soiltest_all_init


text:
module mod_monolis_solver_GMRES
  use mod_monolis_prm
  use mod_monolis_com
  use mod_monolis_mat
  use mod_monolis_precond
  use mod_monolis_matvec
  use mod_monolis_linalg
  use mod_monolis_linalg_util
  use mod_monolis_converge

  implicit none

contains

  subroutine monolis_solver_GMRES(monoPRM, monoCOM, monoMAT)
    implicit none
    type(monolis_prm) :: monoPRM
    type(monolis_com) :: monoCOM
    type(monolis_mat) :: monoMAT
    integer(kint), parameter :: R  = 1
    integer(kint), parameter :: ZP = R + 1
    integer(kint), parameter :: ZQ = R + 2
    integer(kint), parameter :: S  = R + 3
    integer(kint), parameter :: W  = S + 1
    integer(kint), parameter :: Y  = W
    integer(kint), parameter :: AV = Y  + 1
    integer(kint), parameter :: V  = AV + 1
    integer(kint) :: N, NP, NDOF, NNDOF, CS, SN, NRK
    integer(kint) :: i, iter, ik, irow, jj, k, kk, nrest
    real(kdouble), pointer :: B(:), X(:), SS(:)
    real(kdouble), pointer :: WW(:,:), H(:,:)
    real(kdouble) :: BNRM2, DNRM2, RNORM, RESID, tol
    real(kdouble) :: coef, val, VCS, VSN, DTEMP, AA, BB, R0, scale, RR
    logical :: is_converge

    N     = monoMAT%N
    NP    = monoMAT%NP
    NDOF  = monoMAT%NDOF
    NNDOF = N*NDOF
    X => monoMAT%X
    B => monoMAT%B
    tol = monoPRM%tol

    if(monoPRM%is_init_x) X = 0.0d0

    NREST = 10
    if (NREST >= NDOF*NP-1) NREST = NDOF*NP-2

    NRK = NREST + 7
    allocate (H (NRK,NRK), source = 0.0d0)
    allocate (WW(NDOF*NP,NRK), source = 0.0d0)
    allocate (SS(NRK), source = 0.0d0)

    CS  = NREST + 1
    SN  = CS    + 1

    call monolis_residual(monoCOM, monoMAT, X, B, WW(:,R), monoPRM%tspmv, monoPRM%tcomm_spmv)
    call monolis_set_converge(monoPRM, monoCOM, monoMAT, B, BNRM2, is_converge, monoPRM%tdotp, monoPRM%tcomm_dotp)
    if(is_converge) return

    ITER = 0
    OUTER:do

      call monolis_inner_product_R(monoCOM, N, NDOF, WW(:,R), WW(:,R), DNRM2, monoPRM%tdotp, monoPRM%tcomm_dotp)
      if (DNRM2 < tol) exit OUTER ! converged

      RNORM= dsqrt(DNRM2)
      coef= 1.0d0/RNORM
      do ik= 1, NNDOF
        WW(ik,V)= WW(ik,R) * coef
      enddo

      WW(1 ,S) = RNORM
      do k = 2, NNDOF
        WW(k,S) = 0.0d0
      enddo

      do I = 1, NREST
        ITER= ITER + 1

        call monolis_precond_apply(monoPRM, monoCOM, monoMAT, WW(:,V+I-1), WW(:,ZQ))

        call monolis_matvec(monoCOM, monoMAT, WW(:,ZQ), WW(:,W), monoPRM%tspmv, monoPRM%tcomm_spmv)

        do K = 1, I
          call monolis_inner_product_R(monoCOM, N, NDOF, WW(:,W), WW(:,V+K-1), val, monoPRM%tdotp, monoPRM%tcomm_dotp)

          do ik= 1, NNDOF
            WW(ik,W)= WW(ik,W) - val * WW(ik,V+K-1)
          enddo
          H(K,I)= val
        enddo

        call monolis_inner_product_R(monoCOM, N, NDOF, WW(:,W), WW(:,W), val, monoPRM%tdotp, monoPRM%tcomm_dotp)

        if (val < tol) exit ! converged

        H(I+1,I)= dsqrt(val)
        coef= 1.0d0 / H(I+1,I)
        do ik= 1, NNDOF
          WW(ik,V+I+1-1)= WW(ik,W) * coef
        enddo

        do k = 1, I-1
          VCS= H(k,CS)
          VSN= H(k,SN)
          DTEMP   = VCS*H(k  ,I) + VSN*H(k+1,I)
          H(k+1,I)= VCS*H(k+1,I) - VSN*H(k  ,I)
          H(k  ,I)= DTEMP
        enddo

        AA = H(I  ,I)
        BB = H(I+1,I)
        R0= BB
        if (dabs(AA).gt.dabs(BB)) R0= AA
        scale= dabs(AA) + dabs(BB)

        if (scale.ne.0.d0) then
          RR= scale * dsqrt((AA/scale)**2+(BB/scale)**2)
          RR= dsign(1.d0,R0)*RR
          H(I,CS)= AA/RR
          H(I,SN)= BB/RR
        else
          H(I,CS)= 1.d0
          H(I,SN)= 0.d0
          RR     = 0.d0
        endif

        VCS= H(I,CS)
        VSN= H(I,SN)
        DTEMP    = VCS*H(I  ,I) + VSN*H(I+1,I)
        H (I+1,I)= VCS*H(I+1,I) - VSN*H(I  ,I)
        H (I  ,I)= DTEMP

        DTEMP    = VCS*WW(I  ,S) + VSN*WW(I+1,S)
        WW(I+1,S)= VCS*WW(I+1,S) - VSN*WW(I  ,S)
        WW(I  ,S)= DTEMP

        RESID = dabs ( WW(I+1,S))/dsqrt(BNRM2)

        if ( RESID.le.TOL ) then
          !C-- [H]{y}= {s_tld}
          do ik= 1, I
            SS(ik)= WW(ik,S)
          enddo
          IROW= I
          WW(IROW,Y)= SS(IROW) / H(IROW,IROW)

          do kk= IROW-1, 1, -1
            do jj= IROW, kk+1, -1
              SS(kk)= SS(kk) - H(kk,jj)*WW(jj,Y)
            enddo
            WW(kk,Y)= SS(kk) / H(kk,kk)
          enddo

          !C-- {x}= {x} + {y}{V}
          do kk= 1, NNDOF
            WW(kk, AV)= 0.d0
          enddo

          jj= IROW
          do jj= 1, IROW
            do kk= 1, NNDOF
              WW(kk,AV)= WW(kk,AV) + WW(jj,Y)*WW(kk,V+jj-1)
            enddo
          enddo

          call monolis_precond_apply(monoPRM, monoCOM, monoMAT, WW(:,AV), WW(:,ZQ))

          do kk= 1, NNDOF
            X(kk)= X(kk) + WW(kk,ZQ)
          enddo

          exit OUTER
        endif

        if(iter > monoPRM%maxiter) exit OUTER
      end do

      !C-- [H]{y}= {s_tld}
      do ik= 1, NREST
        SS(ik)= WW(ik,S)
      enddo

      IROW = NREST
      if(H(IROW,IROW) /= 0.0d0) WW(IROW,Y)= SS(IROW) / H(IROW,IROW)

      do kk= IROW-1, 1, -1
        do jj= IROW, kk+1, -1
          SS(kk)= SS(kk) - H(kk,jj)*WW(jj,Y)
        enddo
        if(H(kk,kk) /= 0.0d0)  WW(kk,Y)= SS(kk) / H(kk,kk)
      enddo

      !C-- {x}= {x} + {y}{V}
      do kk= 1, NNDOF
        WW(kk, AV)= 0.d0
      enddo

      jj= IROW
      do jj= 1, IROW
        do kk= 1, NNDOF
          WW(kk,AV)= WW(kk,AV) + WW(jj,Y)*WW(kk,V+jj-1)
        enddo
      enddo

      call monolis_precond_apply(monoPRM, monoCOM, monoMAT, WW(:,AV), WW(:,ZQ))

      do kk= 1, NNDOF
        X(kk)= X(kk) + WW(kk,ZQ)
      enddo

      call monolis_residual(monoCOM, monoMAT, X, B, WW(:,R), monoPRM%tspmv, monoPRM%tcomm_spmv)

      call monolis_inner_product_R(monoCOM, N, NDOF, WW(:,R), WW(:,R), DNRM2, monoPRM%tdotp, monoPRM%tcomm_dotp)

      WW(I+1,S)= dsqrt(DNRM2/BNRM2)
      RESID    = WW( I+1,S )

      call monolis_check_converge(monoPRM, monoCOM, monoMAT, WW(:,R), BNRM2, iter, is_converge, monoPRM%tdotp, monoPRM%tcomm_dotp)
      if(is_converge) exit
      if(iter > monoPRM%maxiter) exit OUTER
    end do OUTER

    call monolis_update_R(monoCOM, NDOF, X, monoPRM%tcomm_spmv)

  end subroutine monolis_solver_GMRES

end module mod_monolis_solver_GMRES

max_stars_repo_name: nqomorita/monolis

max_issues_repo_name: nqomorita/monolis

max_forks_repo_name: nqomorita/monolis


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> Solve the matrix equation AX = B
MODULE LinearSolversModule
  USE DataTypesModule, ONLY : NTREAL
  USE DMatrixModule, ONLY : DestructMatrix
  USE LoadBalancerModule, ONLY : PermuteMatrix, UndoPermuteMatrix
  USE LoggingModule, ONLY : EnterSubLog, ExitSubLog, WriteElement, &
       & WriteHeader, WriteListElement
  USE PMatrixMemoryPoolModule, ONLY : MatrixMemoryPool_p, &
       & DestructMatrixMemoryPool
  USE PSMatrixAlgebraModule, ONLY : IncrementMatrix, MatrixNorm, &
       & MatrixMultiply, MatrixTrace, ScaleMatrix
  USE PSMatrixModule, ONLY : Matrix_ps, ConstructEmptyMatrix, &
       & TransposeMatrix, DestructMatrix, ConjugateMatrix, CopyMatrix, &
       & FillMatrixIdentity, PrintMatrixInformation
  USE SolverParametersModule, ONLY : SolverParameters_t, PrintParameters, &
       & DestructSolverParameters
  IMPLICIT NONE
  PRIVATE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  PUBLIC :: CGSolver
CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Solve the matrix equation AX = B using the conjugate gradient method.
  SUBROUTINE CGSolver(AMat, XMat, BMat, solver_parameters_in)
    !> The matrix A, must be hermitian, positive definite.
    TYPE(Matrix_ps), INTENT(IN)  :: AMat
    !> The solved for matrix X.
    TYPE(Matrix_ps), INTENT(INOUT) :: XMat
    !> The right hand side.
    TYPE(Matrix_ps), INTENT(IN)  :: BMat
    !> Parameters for the solver
    TYPE(SolverParameters_t), INTENT(IN), OPTIONAL :: solver_parameters_in
    !! Handling Optional Parameters
    TYPE(SolverParameters_t) :: solver_parameters
    !! Local Variables
    TYPE(Matrix_ps) :: Identity
    TYPE(Matrix_ps) :: ABalanced
    TYPE(Matrix_ps) :: BBalanced
    TYPE(Matrix_ps) :: RMat, PMat, QMat
    TYPE(Matrix_ps) :: RMatT, PMatT
    TYPE(Matrix_ps) :: TempMat
    !! Temporary Variables
    INTEGER :: outer_counter
    REAL(NTREAL) :: norm_value
    TYPE(MatrixMemoryPool_p) :: pool
    REAL(NTREAL) :: top, bottom, new_top, step_size

    !! Optional Parameters
    IF (PRESENT(solver_parameters_in)) THEN
       solver_parameters = solver_parameters_in
    ELSE
       solver_parameters = SolverParameters_t()
    END IF

    !! Print out parameters
    IF (solver_parameters%be_verbose) THEN
       CALL WriteHeader("Linear Solver")
       CALL EnterSubLog
       CALL WriteElement(key="Method", VALUE="CG")
       CALL PrintParameters(solver_parameters)
    END IF

    !! Setup all the matrices
    CALL ConstructEmptyMatrix(Identity, AMat)
    CALL FillMatrixIdentity(Identity)
    CALL ConstructEmptyMatrix(ABalanced, AMat)
    CALL ConstructEmptyMatrix(BBalanced, AMat)
    CALL ConstructEmptyMatrix(RMat, AMat)
    CALL ConstructEmptyMatrix(PMat, AMat)
    CALL ConstructEmptyMatrix(QMat, AMat)
    CALL ConstructEmptyMatrix(TempMat, AMat)

    !! Load Balancing Step
    IF (solver_parameters%do_load_balancing) THEN
       CALL PermuteMatrix(Identity, Identity, &
            & solver_parameters%BalancePermutation, memorypool_in=pool)
       CALL PermuteMatrix(AMat, ABalanced, &
            & solver_parameters%BalancePermutation, memorypool_in=pool)
       CALL PermuteMatrix(BMat, BBalanced, &
            & solver_parameters%BalancePermutation, memorypool_in=pool)
    ELSE
       CALL CopyMatrix(AMat,ABalanced)
       CALL CopyMatrix(BMat,BBalanced)
    END IF

    !! Initial Matrix Values
    CALL CopyMatrix(Identity, XMat)
    !! Compute residual
    CALL MatrixMultiply(ABalanced, Xmat, TempMat, &
         & threshold_in=solver_parameters%threshold, memory_pool_in=pool)
    CALL CopyMatrix(BBalanced,RMat)
    CALL IncrementMatrix(TempMat, RMat, -1.0_NTREAL)
    CALL CopyMatrix(RMat,PMat)

    !! Iterate
    IF (solver_parameters%be_verbose) THEN
       CALL WriteHeader("Iterations")
       CALL EnterSubLog
    END IF
    norm_value = solver_parameters%converge_diff + 1.0_NTREAL
    DO outer_counter = 1,solver_parameters%max_iterations
       IF (solver_parameters%be_verbose .AND. outer_counter .GT. 1) THEN
          CALL WriteListElement(key="Round", VALUE=outer_counter-1)
          CALL EnterSubLog
          CALL WriteListElement(key="Convergence", VALUE=norm_value)
          CALL ExitSubLog
       END IF
       IF (norm_value .LE. solver_parameters%converge_diff) THEN
          EXIT
       END IF

       !! Compute the Step Size
       CALL MatrixMultiply(ABalanced, PMat, QMat, &
            & threshold_in=solver_parameters%threshold, memory_pool_in=pool)

       CALL TransposeMatrix(RMat,RMatT)
       IF (RMatT%is_complex) THEN
          CALL ConjugateMatrix(RMatT)
       END IF
       CALL MatrixMultiply(RMatT, RMat, TempMat, &
            & threshold_in=solver_parameters%threshold, memory_pool_in=pool)
       CALL MatrixTrace(TempMat, top)
       CALL TransposeMatrix(PMat,PMatT)
       IF (PMatT%is_complex) THEN
          CALL ConjugateMatrix(PMatT)
       END IF
       CALL MatrixMultiply(PMatT, QMat, TempMat, &
            & threshold_in=solver_parameters%threshold, memory_pool_in=pool)
       CALL MatrixTrace(TempMat, bottom)
       step_size = top/bottom

       !! Update
       CALL IncrementMatrix(PMat, XMat, alpha_in=step_size)
       norm_value = ABS(step_size*MatrixNorm(PMat))
       CALL IncrementMatrix(QMat, RMat, alpha_in=-1.0_NTREAL*step_size)

       !! Update PMat
       CALL TransposeMatrix(RMat,RMatT)
       IF (RMatT%is_complex) THEN
          CALL ConjugateMatrix(RMatT)
       END IF
       CALL MatrixMultiply(RMatT, RMat, TempMat, &
            & threshold_in=solver_parameters%threshold, memory_pool_in=pool)
       CALL MatrixTrace(TempMat, new_top)
       step_size = new_top / top
       CALL ScaleMatrix(PMat, step_size)
       CALL IncrementMatrix(RMat, PMat)

    END DO
    IF (solver_parameters%be_verbose) THEN
       CALL ExitSubLog
       CALL WriteElement(key="Total_Iterations", VALUE=outer_counter-1)
       CALL PrintMatrixInformation(XMat)
    END IF

    !! Undo Load Balancing Step
    IF (solver_parameters%do_load_balancing) THEN
       CALL UndoPermuteMatrix(XMat,XMat, &
            & solver_parameters%BalancePermutation, memorypool_in=pool)
    END IF

    !! Cleanup
    IF (solver_parameters%be_verbose) THEN
       CALL ExitSubLog
    END IF
    CALL DestructMatrix(TempMat)
    CALL DestructMatrix(RMat)
    CALL DestructMatrix(PMat)
    CALL DestructMatrix(QMat)
    CALL DestructMatrix(Identity)
    CALL DestructMatrix(ABalanced)
    CALL DestructMatrix(BBalanced)
    CALL DestructMatrixMemoryPool(pool)
    CALL DestructSolverParameters(solver_parameters)
  END SUBROUTINE CGSolver
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END MODULE LinearSolversModule


text:
      MODULE GwfResModule

        public

        INTEGER, SAVE,POINTER   ::NRES,IRESCB,NRESOP,IRESPT,NPTS
        INTEGER, SAVE, DIMENSION(:,:), POINTER ::IRES
        INTEGER, SAVE, DIMENSION(:,:), POINTER ::IRESL
        double precision,    SAVE, DIMENSION(:,:), POINTER ::BRES
        double precision,    SAVE, DIMENSION(:,:), POINTER ::CRES
        double precision,    SAVE, DIMENSION(:,:), POINTER ::BBRES
        double precision,    SAVE, DIMENSION(:),   POINTER ::HRES
        double precision,    SAVE, DIMENSION(:,:), POINTER ::HRESSE
      TYPE GWFRESTYPE
        INTEGER,POINTER   ::NRES,IRESCB,NRESOP,IRESPT,NPTS
        INTEGER,  DIMENSION(:,:), POINTER ::IRES
        INTEGER,  DIMENSION(:,:), POINTER ::IRESL
        double precision,     DIMENSION(:,:), POINTER ::BRES
        double precision,     DIMENSION(:,:), POINTER ::CRES
        double precision,     DIMENSION(:,:), POINTER ::BBRES
        double precision,     DIMENSION(:),   POINTER ::HRES
        double precision,     DIMENSION(:,:), POINTER ::HRESSE
      END TYPE
      TYPE(GWFRESTYPE), SAVE :: GWFRESDAT(10)

      contains

      SUBROUTINE GWF2RES7DA(IGRID)
C  Deallocate RES MEMORY
C
        CALL SGWF2RES7PNT(IGRID)
        DEALLOCATE(NRES)
        DEALLOCATE(IRESCB)
        DEALLOCATE(NRESOP)
        DEALLOCATE(IRESPT)
        DEALLOCATE(NPTS)
        DEALLOCATE(IRES)
        DEALLOCATE(IRESL)
        DEALLOCATE(BRES)
        DEALLOCATE(CRES)
        DEALLOCATE(BBRES)
        DEALLOCATE(HRES)
        DEALLOCATE(HRESSE)
C
      RETURN
      END SUBROUTINE GWF2RES7DA

      SUBROUTINE SGWF2RES7PNT(IGRID)
C  Change RES data to a different grid.
C
        NRES=>GWFRESDAT(IGRID)%NRES
        IRESCB=>GWFRESDAT(IGRID)%IRESCB
        NRESOP=>GWFRESDAT(IGRID)%NRESOP
        IRESPT=>GWFRESDAT(IGRID)%IRESPT
        NPTS=>GWFRESDAT(IGRID)%NPTS
        IRES=>GWFRESDAT(IGRID)%IRES
        IRESL=>GWFRESDAT(IGRID)%IRESL
        BRES=>GWFRESDAT(IGRID)%BRES
        CRES=>GWFRESDAT(IGRID)%CRES
        BBRES=>GWFRESDAT(IGRID)%BBRES
        HRES=>GWFRESDAT(IGRID)%HRES
        HRESSE=>GWFRESDAT(IGRID)%HRESSE
C
      RETURN
      END SUBROUTINE SGWF2RES7PNT

      SUBROUTINE SGWF2RES7PSV(IGRID)
C  Save RES data for a grid.
C
        GWFRESDAT(IGRID)%NRES=>NRES
        GWFRESDAT(IGRID)%IRESCB=>IRESCB
        GWFRESDAT(IGRID)%NRESOP=>NRESOP
        GWFRESDAT(IGRID)%IRESPT=>IRESPT
        GWFRESDAT(IGRID)%NPTS=>NPTS
        GWFRESDAT(IGRID)%IRES=>IRES
        GWFRESDAT(IGRID)%IRESL=>IRESL
        GWFRESDAT(IGRID)%BRES=>BRES
        GWFRESDAT(IGRID)%CRES=>CRES
        GWFRESDAT(IGRID)%BBRES=>BBRES
        GWFRESDAT(IGRID)%HRES=>HRES
        GWFRESDAT(IGRID)%HRESSE=>HRESSE
C
      RETURN
      END SUBROUTINE SGWF2RES7PSV

      END MODULE GwfResModule


text: #redirect Users/yzhelfman

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! { dg-do run }
! Tests the patch to implement the array version of the TRANSFER
! intrinsic (PR17298).
! Contributed by Paul Thomas  <pault@gcc.gnu.org>

! Bigendian test posted by Perseus in comp.lang.fortran on 4 July 2005.
! Original had parameter but this fails, at present, if is_gimple_var with -Ox, x>0

   LOGICAL :: bigend
   integer :: icheck = 1

   character(8) :: ch(2) = (/"lmnoPQRS","LMNOpqrs"/)

   bigend = IACHAR(TRANSFER(icheck,"a")) == 0

! tests numeric transfers other than original testscase.

   call test1 ()

! tests numeric/character transfers.

   call test2 ()

! Test dummies, automatic objects and assumed character length.

   call test3 (ch, ch, ch, 8)

contains

   subroutine test1 ()
     real(4) :: a(4, 4)
     integer(2) :: it(4, 2, 4), jt(32)

! Check multi-dimensional sources and that transfer works as an actual
! argument of reshape.

     a = reshape ((/(rand (), i = 1, 16)/), (/4,4/))
     jt = transfer (a, it)
     it = reshape (jt, (/4, 2, 4/))
     if (any (reshape (transfer (it, a), (/4,4/)) .ne. a)) STOP 1

   end subroutine test1

   subroutine test2 ()
     integer(4) :: y(4), z(2)
     character(4) :: ch(4)

! Allow for endian-ness
     if (bigend) then
       y = (/(i + 3 + ishft (i + 2, 8) + ishft (i + 1, 16) &
                + ishft (i, 24), i = 65, 80 , 4)/)
     else
       y = (/(i + ishft (i + 1, 8) + ishft (i + 2, 16) &
                + ishft (i + 3, 24), i = 65, 80 , 4)/)
     end if

! Check source array sections in both directions.

     ch = "wxyz"
     ch(1:2) = transfer (y(2:4:2), ch)
     if (any (ch(1:2) .ne. (/"EFGH","MNOP"/))) STOP 2
     ch = "wxyz"
     ch(1:2) = transfer (y(4:2:-2), ch)
     if (any (ch(1:2) .ne. (/"MNOP","EFGH"/))) STOP 3

! Check that a complete array transfers with size absent.

     ch = transfer (y, ch)
     if (any (ch .ne. (/"ABCD","EFGH","IJKL","MNOP"/))) STOP 4

! Check that a character array section is OK

     z = transfer (ch(2:3), y)
     if (any (z .ne. y(2:3))) STOP 5

! Check dest array sections in both directions.

     ch = "wxyz"
     ch(3:4) = transfer (y, ch, 2)
     if (any (ch(3:4) .ne. (/"ABCD","EFGH"/))) STOP 6
     ch = "wxyz"
     ch(3:2:-1) = transfer (y, ch, 2)
     if (any (ch(2:3) .ne. (/"EFGH","ABCD"/))) STOP 7

! Make sure that character to numeric is OK.

     ch = "wxyz"
     ch(1:2) = transfer (y, ch, 2)
     if (any (ch(1:2) .ne. (/"ABCD","EFGH"/))) STOP 8

     z = transfer (ch, y, 2)
     if (any (y(1:2) .ne. z)) STOP 9

   end subroutine test2

   subroutine test3 (ch1, ch2, ch3, clen)
     integer clen
     character(8) :: ch1(:)
     character(*) :: ch2(2)
     character(clen) :: ch3(2)
     character(8) :: cntrl(2) = (/"lmnoPQRS","LMNOpqrs"/)
     integer(8) :: ic(2)
     ic = transfer (cntrl, ic)

! Check assumed shape.

     if (any (ic .ne. transfer (ch1, ic))) STOP 10

! Check assumed character length.

     if (any (ic .ne. transfer (ch2, ic))) STOP 11

! Check automatic character length.

     if (any (ic .ne. transfer (ch3, ic))) STOP 12

  end subroutine test3

end


text:
      module comrxnb
!     comrxnb
!***********************************************************************
!  Copyright, 1993, 2004,  The  Regents of the University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
!**********************************************************************
!D1
!D1 PURPOSE
!D1
!D1 Include file containing passed parameters and pointers related to
!D1 solute module.
!D1
!**********************************************************************
!D2
!D2 REVISION HISTORY
!D2
!D2 Revision                    ECD
!D2 Date         Programmer     Number  Comments
!D2
!D2 09-13-93     B. Robinson    00022   Initial Implementation
!D2
!D2 $Log:   /pvcs.config/fehm90/src/comrxnb.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:42:40   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 08:59:12   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:07:38   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:22:40   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 11:58:12   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:39:56 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
!D2
!D2    Rev 1.1   03/18/94 16:23:22   gaz
!D2 Added solve_new and cleaned up memory management.
!D2
!D2    Rev 1.0   01/20/94 10:22:18   pvcs
!D2 original version in process of being certified
!D2
!**********************************************************************
!D3
!D3 INTERFACES
!D3
!D3 None
!D3
!**********************************************************************
!D4
!D4 GLOBAL OBJECTS
!D4
!D4 Global Constants
!D4
!D4   None
!D4
!D4 Global Types
!D4
!D4   None
!D4
!D4 Global Variables
!D4
!D4                                 COMMON
!D4   Identifier           Type     Block    Description
!D4
!D4   ***** COMMON Block comrxnb pointers and associated variables *****
!D4   ipstored_derivative  POINTER  comrxnb  Pointer for array stored_derivative
!D4   ipstored_residual    POINTER  comrxnb  Pointer for array stored_residual
!D4   stored_derivatives   REAL*8   comrxnb  Array of derivatives for each
!D4                                            solute at each node
!D4   stored_residual      REAL*8   comrxnb  Array of residuals for each
!D4                                            solute at each node
!D4
!**********************************************************************
!D5
!D5 LOCAL IDENTIFIERS
!D5
!D5 None
!D5
!**********************************************************************
!D6
!D6 FUNCTIONAL DESCRIPTION
!D6
!D6 None
!D6
!**********************************************************************
!D7
!D7 ASSUMPTIONS AND LIMITATIONS
!D7
!D7 None
!D7
!**********************************************************************
!D8
!D8 SPECIAL COMMENTS
!D8
!D8 None
!D8
!**********************************************************************
!D9
!D9 REQUIREMENTS TRACEABILITY
!D9
!D9 None
!D9
!**********************************************************************
!DA
!DA REFERENCES
!DA
!DA None
!DA
!**********************************************************************
!PS
!PS PSEUDOCODE
!PS
!PS None
!PS
!**********************************************************************

      real*8, allocatable ::  stored_derivative(:)
      real*8, allocatable ::  stored_residual(:)

      end module comrxnb

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


!=================================================
! SUBROUTINES
!	BIBSP
!
!=================================================

module mod_bibsp

	implicit none

contains
!==========
subroutine BIBSP(nd2, nbdeg, ndd2, nbd2, lc_bi)

	use prm_matrix
	use mod_bspline

	integer, intent(in)  :: nd2, nbdeg, ndd2, nbd2
	real(8), intent(out) :: lc_bi(Ndd,16) ! bi bspline

	integer :: norder, m1, m0, lx, ly, mx, my, j
	real(8) :: bx(4), bs(Nd,4)
	real(8) :: px

	norder = nbdeg + 1

	do m1 = 1, nd2

		px = (dble(m1) - 0.5d0) / dble(nd2)
		CALL BSPLINE(px,bx,nbdeg)

		do j = 1, norder
			bs(m1,j) = bx(j)
		end do

	end do

	do m0 = 1, nbd2

		lx = m0 - int((m0-1) / norder) * norder
		ly = 1  + int((m0-1) / norder)

		do m1 = 1, ndd2
			mx = m1 - int((m1-1) / nd2) * nd2
			my = 1  + int((m1-1) / nd2)

			lc_bi(m1,m0) = bs(mx,lx) * bs(my,ly)
		end do

	end do

return
end subroutine BIBSP

!==========

end module mod_bibsp


text:
! ###################################################################
! Copyright (c) 2013-2019, Marc De Graef Research Group/Carnegie Mellon University
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
!     - Redistributions of source code must retain the above copyright notice, this list
!        of conditions and the following disclaimer.
!     - Redistributions in binary form must reproduce the above copyright notice, this
!        list of conditions and the following disclaimer in the documentation and/or
!        other materials provided with the distribution.
!     - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
!        of its contributors may be used to endorse or promote products derived from
!        this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################

!--------------------------------------------------------------------------
! EMsoft:EMmkxtal.f90
!--------------------------------------------------------------------------
!
! PROGRAM: mkxtal
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief create a crystal structure file (very simple program)
!
!> @date   11/23/13 MDG 1.0 original
!> @date   08/10/14 MDG 1.1 added GenerateSymmetry call
!> @date   03/30/15 MDG 2.0 modified output to HDF format
!> @date   10/07/15 MDG 2.1 minor cleanup in preparation for release 3.0
!> @date   09/06/16 MDG 3.0 added Wyckoff position entry as an option
!> @date   04/20/18 MDG 4.0 start incorporation of 2D and 3D quasi-crystal structures
!--------------------------------------------------------------------------
program EMmkxtal

use local
use typedefs
use HDFsupport
use io
use crystal
use symmetry
use files

IMPLICIT NONE

type(unitcell), pointer         :: cell
character(fnlen)                :: progname, progdesc, fname, source
integer(kind=irg)               :: numarg, i
integer(kind=irg)               :: iargc        !< external function for command line
character(fnlen)                :: arg          !< to be read from the command line
logical                         :: useWyckoff

 progname = 'EMmkxtal.f90'
 progdesc = 'Create an HDF crystal structure file and place it in the XtalFolder'

 call EMsoft(progname, progdesc)

 useWyckoff = .FALSE.
 numarg = iargc()
 if (numarg.gt.0) then ! there is at least one argument
  do i=1,numarg
    call getarg(i,arg)
!    mess = 'Found the following argument: '//trim(arg); call Message("(/A/)")
! does the argument start with a '-' character?
    if (arg(1:1).eq.'-') then
        if (trim(arg).eq.'-h') then
         call Message(' Program should be called as follows: ', frm = "(/A)")
         call Message('        '//trim(progname)//' [-h] [-w] ', frm = "(A)")
         call Message(' To produce this message, type '//trim(progname)//' -h', frm = "(A)")
         call Message(' To use Wyckoff positions to enter atom coordinates, use -w option', frm = "(A)")
        end if
        if (trim(arg).eq.'-w') then
! with this option the GetAsymPosWyckoff routine will ask the user for Wyckoff Positions instead of regular cordinate strings
         useWyckoff = .TRUE.
        end if
    end if
  end do
 end if

 allocate(cell)

 cell%SYM_SGset=0
 call GetLatParm(cell)
 call GetSpaceGroup(cell)
 call GenerateSymmetry(cell,.TRUE.)
 if (useWyckoff) then
  call GetAsymPosWyckoff(cell)
 else
  call GetAsymPos(cell)
 end if
 call ReadValue('Enter output file name (*.xtal) ', fname)

 call ReadValue('Enter the source for this data [max. 512 characters] ', source)
 cell%source = trim(source)
 cell%fname = fname

 write (*,*) 'Source = '//trim(cell%source)
 call SaveDataHDF(cell)

end program EMmkxtal

max_stars_repo_name: ZachClayburn/EMsoft

max_issues_repo_name: ZachClayburn/EMsoft

max_forks_repo_name: ZachClayburn/EMsoft


text:
tinyalu_pkg.sv
tinyalu_bfm.sv
top.sv
+incdir+./tb_classes

max_stars_repo_name: OpenEDF/uvmprimer

max_issues_repo_name: OpenEDF/uvmprimer

max_forks_repo_name: albasili/uvmprimer


text:
      logical*4 function tflogi(word,exist)
      use ffs
      implicit none
      integer*4 i
      character*(*) ,intent(in):: word
      character*32 fname1
      logical*4 ,intent(out):: exist
      tflogi=.false.
      exist=.true.
      if(word .eq. ' ')then
        exist=.false.
        tflogi=.false.
      endif
      do i=1,nflag
        if(word .eq. fname(i))then
          tflogi=flags(i)
          return
        else
          fname1='NO'//fname(i)
          if(word .eq. fname1)then
            tflogi=.not. flags(i)
            return
          elseif(word .eq. sino(i))then
            tflogi=.not. flags(i)
            return
          else
            fname1='NO'//sino(i)
            if(word .eq. fname1)then
              tflogi=flags(i)
              return
            endif
          endif
        endif
      enddo
      exist=.false.
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


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

module Batse_mod

    use Constants_mod, only: IK, RK, PI, LN10

    implicit none

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

    integer(IK) , parameter :: NLGRB = 1366_IK, NSGRB = 565_IK, NVAR = 4_IK
    real(RK)    , parameter :: MAX_LOG10PH53_4_LOGPBOLZERO  = 6.318167895318538_RK
    real(RK)    , parameter :: MIN_LOG10PH53_4_LOGPBOLZERO  = 4.92_RK
    real(RK)    , parameter :: MAX_LOGPH53_4_LOGPBOLZERO    = MAX_LOG10PH53_4_LOGPBOLZERO * LN10
    real(RK)    , parameter :: MIN_LOGPH53_4_LOGPBOLZERO    = MIN_LOG10PH53_4_LOGPBOLZERO * LN10
    real(RK)    , parameter :: DIF_LOGPH53_4_LOGPBOLZERO    = MAX_LOGPH53_4_LOGPBOLZERO - MIN_LOGPH53_4_LOGPBOLZERO

    ! Parameters for the effective peak flux of SGRBs (P64ms conversion to P1024):

!   ! The scale of the change in BATSE efficiency for different GRB durations.
!   real(RK)    , parameter :: THRESH_ERFC_BASE             = +0.146314238936889_RK
!
!   ! The scale of the change in BATSE efficiency for different GRB durations.
!   real(RK)    , parameter :: THRESH_ERFC_AMP              = +0.570285374263156_RK
!
!   ! Mean duration in the Error function used to model the connection between the peak fluxes in 64 and 1024 ms.
!   real(RK)    , parameter :: THRESH_ERFC_AVG              = -0.480811530417719_RK
!
!   ! scale of the duration in the Error function used to model the connection between the peak fluxes in 64 and 1024 ms.
!   real(RK)    , parameter :: THRESH_ERFC_STD              = +1.292443569094922_RK

    ! The scale of the change in BATSE efficiency for different GRB durations.
    real(RK)    , parameter :: THRESH_ERFC_BASE             = +0.146314238936889_RK * LN10

    ! The scale of the change in BATSE efficiency for different GRB durations.
    real(RK)    , parameter :: THRESH_ERFC_AMP              = +0.282313526464596_RK * LN10  ! Serfc

    ! Mean duration in the Error function used to model the connection between the peak fluxes in 64 and 1024 ms.
    real(RK)    , parameter :: THRESH_ERFC_AVG              = -0.483553339256463_RK * LN10  ! meandur

    ! scale of the duration in the Error function used to model the connection between the peak fluxes in 64 and 1024 ms.
    real(RK)    , parameter :: THRESH_ERFC_STD              = 1.0514698984694800_RK * LN10  ! scaledur

    ! inverse scale of the duration in the Error function used to model the connection between the peak fluxes in 64 and 1024 ms.
    real(RK)    , parameter :: THRESH_ERFC_STD_INV          = 1._RK / THRESH_ERFC_STD       ! inverse scaledur

    ! The height of the ERFC function.
    real(RK)    , parameter :: THRESH_ERFC_HEIGHT           = 2_IK * THRESH_ERFC_AMP

    ! correction that must be added to logPbol64ms to convert it to effective peak flux.
    ! Effective LogPbol limit above which trigger efficiency is 100%, for any Log(Epk) and Log(dur). It is equivalent to maximum Log(Pbol) at very long durations.
    real(RK)    , parameter :: THRESH_LOGPBOL64_CORRECTION  = DIF_LOGPH53_4_LOGPBOLZERO - MIN_LOGPH53_4_LOGPBOLZERO + THRESH_ERFC_HEIGHT ! equivalent to lpb_correction

    ! GRB attributes
    type :: Event_type
        real(RK) :: logPbol, logEpk, logSbol, logT90, logPF53
    end type Event_type

    type :: GRB_type
        integer(IK) :: count
        type(Event_type), allocatable :: Event(:)
        !type(Event_type) :: Event(NLGRB)
    end type GRB_type

#if defined CAF_ENABLED
    type(GRB_type) :: GRB[*]
#else
    type(GRB_type) :: GRB
#endif

    integer(IK), allocatable :: Trigger(:)
    !integer(IK) :: Trigger(NLGRB)
    !integer(IK) :: TriggerSGRB(NSGRB)

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

contains

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

    ! returns all log of data in natural (Neper) base.
    subroutine readDataGRB(inFilePath,outFilePath,isLgrb)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: readDataGRB
#endif

        use Constants_mod, only: IK, RK
        implicit none
        character(*), intent(in)            :: inFilePath, outFilePath
        integer(IK)                         :: inFileUnit, outFileUnit, igrb
        logical     , intent(in)            :: isLgrb

        if (isLgrb) then
            GRB%count = NLGRB
        else
            GRB%count = NSGRB
        end if

        allocate(GRB%Event(GRB%count), Trigger(GRB%count))

#if defined CAF_ENABLED
        if (this_image()==1) then
#endif
            open(newunit=inFileUnit,file=inFilePath,status="old")
            open(newunit=outFileUnit,file=outFilePath,status="replace")
            write(outFileUnit,"(9a30)"  ) "trigger"             &
                                        , "logPbol_1eV_20MeV"   &
                                        , "logSbol_1eV_20MeV"   &
                                        , "logEpk"              &
                                        , "logEPR1024"          &
                                        , "logEFR"              &
                                        , "logFPR1024"          &
                                        , "logT90"              &
                                        , "logEffPF53"

            ! skip the header row in the input file
            read(inFileUnit,*)

            ! read BATSE GRB data
            do igrb = 1, GRB%count

                read(inFileUnit,*   ) Trigger(igrb)             &
                                    , GRB%Event(igrb)%logPF53   &
                                    , GRB%Event(igrb)%logEpk    &
                                    , GRB%Event(igrb)%logSbol   &
                                    , GRB%Event(igrb)%logT90

                ! convert all values to logarithm in base Neper

                GRB%Event(igrb)%logPF53 = LN10 * GRB%Event(igrb)%logPF53
                GRB%Event(igrb)%logEpk  = LN10 * GRB%Event(igrb)%logEpk
                GRB%Event(igrb)%logSbol = LN10 * GRB%Event(igrb)%logSbol
                GRB%Event(igrb)%logT90  = LN10 * GRB%Event(igrb)%logT90

                ! convert photon count data to energy in units of ergs

                GRB%Event(igrb)%logPbol = getLogPbol( GRB%Event(igrb)%logEpk, GRB%Event(igrb)%logPF53 )
                if (isLgrb) then
                    GRB%Event(igrb)%logSbol = getLogPbol( GRB%Event(igrb)%logEpk, GRB%Event(igrb)%logSbol )
                else
                    GRB%Event(igrb)%logPF53 = GRB%Event(igrb)%logPF53 - THRESH_ERFC_AMP * erfc( (GRB%Event(igrb)%logT90-THRESH_ERFC_AVG) * THRESH_ERFC_STD_INV )
                end if

                ! write the converted data to output file
                write(outFileUnit,"(i30,8E30.6)") Trigger(igrb)                                     &
                                                , GRB%Event(igrb)%logPbol                           &
                                                , GRB%Event(igrb)%logSbol                           &
                                                , GRB%Event(igrb)%logEpk                            &
                                                , GRB%Event(igrb)%logEpk-GRB%Event(igrb)%logPbol    &
                                                , GRB%Event(igrb)%logEpk-GRB%Event(igrb)%logSbol    &
                                                , GRB%Event(igrb)%logSbol-GRB%Event(igrb)%logPbol   &
                                                , GRB%Event(igrb)%logT90                            &
                                                , GRB%Event(igrb)%logPF53

            end do
        close(outFileUnit)
        close(inFileUnit)

#if defined CAF_ENABLED
        sync images(*)

    else

        sync images(*)
        do igrb = 1, GRB%count
            GRB%Event(igrb)%logPbol = GRB[1]%Event(igrb)%logPbol
            GRB%Event(igrb)%logSbol = GRB[1]%Event(igrb)%logSbol
            GRB%Event(igrb)%logPF53 = GRB[1]%Event(igrb)%logPF53
            GRB%Event(igrb)%logEpk  = GRB[1]%Event(igrb)%logEpk
            GRB%Event(igrb)%logT90  = GRB[1]%Event(igrb)%logT90
        end do

    end if
#endif

    end subroutine readDataGRB

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

    pure function getLog10PF53(log10epk,log10pbol) result(log10PF53)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: getLog10PF53
#endif
        ! Given Log10(Epk [KeV]) of an LGRB and its bolometric (0.0001-20000 KeV) peak flux [in units of Ergs/s], Log(Pbol),
        ! this function calculates the corresponding Log10(peak photon flux) in the BATSE detection energy range [50,300] KeV.
        ! Amir Shahmoradi, Wednesday June 27, 2012, 7:15 PM, IFS, The University of Texas at Austin.
        ! Amir Shahmoradi, Wednesday June 27, 2012, 9:29 PM, IFS, The University of Texas at Austin.
        use Constants_mod, only: RK
        implicit none
        real(RK), intent(in)    :: log10epk,log10pbol
        real(RK)                :: log10PF53
        if (log10epk<-2.915056638230699_RK) then
            log10PF53 = log10pbol    + 4.9200000000000000_RK
        elseif (log10epk<1.5_RK) then
            log10PF53 = log10pbol   + 5.7361200000000000_RK + log10epk * &
                                    ( 0.3093600000000000_RK + log10epk * &
                                    ( 0.0045600000000000_RK + log10epk * &
                                    ( 0.0015900000000000_RK + log10epk * &
                                    ( 0.0001533360000000_RK - log10epk   &
                                    * 0.0003574800000000_RK ))))
        elseif (log10epk<2.5_RK) then
            log10PF53 = log10pbol   + 1.9112800000000000_RK + log10epk * &
                                    ( 39.710390000000000_RK - log10epk * &
                                    ( 96.606280000000000_RK - log10epk * &
                                    ( 109.24696000000000_RK - log10epk * &
                                    ( 67.271800000000000_RK - log10epk * &
                                    ( 23.402390000000000_RK - log10epk * &
                                    ( 4.3454400000000000_RK - log10epk   &
                                    * 0.3360600000000000_RK ))))))
        elseif (log10epk<4._RK) then
            log10PF53 = log10pbol   + 2.8020600000000000_RK + log10epk * &
                                    ( 4.5690700000000000_RK - log10epk * &
                                    ( 1.9277200000000000_RK - log10epk * &
                                    ( 0.2938100000000000_RK - log10epk   &
                                    * 0.0148900000000000_RK )))
        elseif (log10epk<5.4093868613659435_RK) then
            log10PF53 = log10pbol   - 10.465330000000000_RK + log10epk * &
                                    ( 26.706370000000000_RK - log10epk * &
                                    ( 14.476310000000000_RK - log10epk * &
                                    ( 3.5404100000000000_RK - log10epk * &
                                    ( 0.4095700000000000_RK - log10epk   &
                                    * 0.0183100000000000_RK ))))
        else ! if (log10epk>=5.4093868613659435_RK) then
            log10PF53 = log10pbol   + 4.9200000000000000_RK
        end if
    end function getLog10PF53

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

    pure function getLogPF53(logEpk,logPbol) result(logPF53)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: getLogPF53
#endif
        ! Given Log(Epk [KeV]) of an LGRB and its bolometric (0.0001-20000 KeV) peak flux [in units of Ergs/s], Log(Pbol),
        ! this function calculates the corresponding Log(peak photon flux) in the BATSE detection energy range [50,300] KeV.
        ! Amir Shahmoradi, Wednesday June 27, 2012, 7:15 PM, IFS, The University of Texas at Austin.
        ! Amir Shahmoradi, Wednesday June 27, 2012, 9:29 PM, IFS, The University of Texas at Austin.
        use Constants_mod, only: RK
        implicit none
        real(RK), intent(in)    :: logEpk,logPbol
        real(RK)                :: logPF53
        if (logEpk<-6.712165960423344_RK) then
            logPF53 = logPbol    + 11.32871865753070600_RK
        elseif (logEpk<3.453877639491069_RK) then
            logPF53 = logPbol   + 13.20790440362500600_RK + logEpk * &
                                ( 0.309360000000000000_RK + logEpk * &
                                ( 0.001980382837478830_RK + logEpk * &
                                ( 0.000299892598248466_RK + logEpk * &
                                ( 1.25602147173493e-05_RK - logEpk   &
                                * 1.27171265917873e-05_RK ))))
        elseif (logEpk<5.756462732485115_RK) then
            logPF53 = logPbol   + 4.400884836537660000_RK + logEpk * &
                                ( 39.71039000000000000_RK - logEpk * &
                                ( 41.95557432120050000_RK - logEpk * &
                                ( 20.60525451895990000_RK - logEpk * &
                                ( 5.510436247342930000_RK - logEpk * &
                                ( 0.832525333390336000_RK - logEpk * &
                                ( 0.067135977132248900_RK - logEpk   &
                                * 0.002254876138523550_RK ))))))
        elseif (logEpk<9.210340371976184_RK) then
            logPF53 = logPbol   + 6.451981585674900000_RK + logEpk * &
                                ( 4.569070000000000000_RK - logEpk * &
                                ( 0.837198158654537000_RK - logEpk * &
                                ( 0.055416002698982300_RK - logEpk   &
                                * 0.001219684856402480_RK )))
        elseif (logEpk<12.455573549219071_RK) then
            logPF53 = logPbol   - 24.09731285126340000_RK + logEpk * &
                                ( 26.70637000000000000_RK - logEpk * &
                                ( 6.286981551320860000_RK - logEpk * &
                                ( 0.667762738216888000_RK - logEpk * &
                                ( 0.033549115287895400_RK - logEpk   &
                                * 0.000651366755890191_RK ))))
        else
            logPF53 = logPbol   + 11.32871865753070600_RK
        end if
        !write(*,"(*(g0.13,:,', '))") "logEpk, logPbol, logPF53<0.0: ", logEpk, logPbol, logPF53
    end function getLogPF53

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

    pure function getLogPbol(logEpk,logPF53) result(logPbol)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: getLogPbol
#endif
        ! Given Log10(Epk [KeV]) of an LGRB and its bolometric (0.0001-20000 KeV) peak flux [in units of Ergs/s], Log(Pbol),
        ! this function calculates the corresponding Log10(peak photon flux) in the BATSE detection energy range [50,300] KeV.
        ! Amir Shahmoradi, Wednesday June 27, 2012, 7:15 PM, IFS, The University of Texas at Austin.
        ! Amir Shahmoradi, Wednesday June 27, 2012, 9:29 PM, IFS, The University of Texas at Austin.
        use Constants_mod, only: RK
        implicit none
        real(RK), intent(in)    :: logEpk, logPF53
        real(RK)                :: logPbol
        logPbol = logPF53 - getLogPF53(logEpk,0._RK)
    end function getLogPbol

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

    pure function getLogEffectivePeakPhotonFlux(logPeakPhotonFlux64ms,logT90) result(logEffectivePeakPhotonFlux)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: getEffectivePeakPhotonFlux
#endif
        ! Converts an input natural-log peak photon flux in 64ms timescale to an effective triggering peak photon flux.
        ! To do so, the observed T90 duration of the event is also necessary as input.
        ! Reference: Eqn A4 of Shahmoradi and Nemiroff 2015, MNRAS, Short versus long gamma-ray bursts.
        use, intrinsic :: iso_fortran_env, only: real32
        use Constants_mod, only: RK
        implicit none
        real(RK), intent(in)    :: logPeakPhotonFlux64ms, logT90
        real(RK)                :: logEffectivePeakPhotonFlux
        logEffectivePeakPhotonFlux  = logPeakPhotonFlux64ms - getLogEffectivePeakPhotonFluxCorrection(logT90)
    end function getLogEffectivePeakPhotonFlux

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

    pure function getLogEffectivePeakPhotonFluxCorrection(logT90) result(logEffectivePeakPhotonFluxCorrection)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: getEffectivePeakPhotonFlux
#endif
        ! Converts an input natural-log peak photon flux in 64ms timescale to an effective triggering peak photon flux.
        ! To do so, the observed T90 duration of the event is also necessary as input.
        ! Reference: Eqn A4 of Shahmoradi and Nemiroff 2015, MNRAS, Short versus long gamma-ray bursts.
        use, intrinsic :: iso_fortran_env, only: real32
        use Constants_mod, only: RK
        implicit none
        real(RK), intent(in)    :: logT90
        real(RK)                :: logEffectivePeakPhotonFluxCorrection
        logEffectivePeakPhotonFluxCorrection    = THRESH_ERFC_AMP * erfc(real((logT90-THRESH_ERFC_AVG)/THRESH_ERFC_STD,kind=real32))
                                              ! + THRESH_ERFC_BASE ! adding this term will make the effective peak flux equivalent to PF1024ms
    end function getLogEffectivePeakPhotonFluxCorrection

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

end module Batse_mod


text:
module mm_factory_module
    use :: util_api, only : &
            string, &
            dictionary, &
            dictionary_converter

    use :: matrix_multiplication_module, only : matrix_multiplication

    implicit none
    private

    public :: mm_factory

    type, abstract :: mm_factory
    contains
        procedure :: get => get
        procedure :: create => create
        procedure(create_interface), deferred :: create_from_key
        procedure(get_available_interface), deferred :: get_available_mm_drivers
    end type mm_factory

    abstract interface
        subroutine create_interface(this, mm, key)
            import :: mm_factory
            import :: matrix_multiplication
            import :: string

            class(mm_factory), intent(in) :: this
            class(matrix_multiplication), allocatable, intent(inout) :: mm
            type(string), intent(in) :: key
        end subroutine create_interface

        function get_available_interface(this) result(drivers)
            import :: mm_factory
            import :: string

            class(mm_factory), intent(in) :: this
            type(string), dimension(:), allocatable :: drivers
        end function get_available_interface
    end interface

    character(len=*), parameter :: mm_driver_key = "mm_driver"
    character(len=*), parameter :: default_mm_driver = "intrinsic"
contains
    function get(this, driver, options, priorities) result(mm)
        class(mm_factory), intent(in) :: this
        character(len=*), intent(in), optional :: driver
        type(dictionary), intent(in), optional :: options
        type(string), dimension(:), intent(in), optional :: priorities
        class(matrix_multiplication), allocatable :: mm

        call this%create(mm, driver, options, priorities)
    end function get

    subroutine create(this, mm, driver, options, priorities)
        class(mm_factory), intent(in) :: this
        class(matrix_multiplication), allocatable, intent(inout) :: mm
        character(len=*), intent(in), optional :: driver
        type(dictionary), intent(in), optional :: options
        type(string), dimension(:), intent(in), optional :: priorities

        type(string) :: key
        type(dictionary_converter) :: conv

        if ( present(driver) ) then
            key = driver
        else
            key = conv%to_string(mm_driver_key, options, priorities, default_mm_driver)
        end if
        call this%create_from_key(mm, key)
    end subroutine create
end module mm_factory_module

max_stars_repo_name: cheshyre/ntcl-algorithms

max_issues_repo_name: cheshyre/ntcl-algorithms

max_forks_repo_name: cheshyre/ntcl-algorithms


!     Output the open bnd nodes in hgrid.gr3 as build pts for viz
!     Input:
!           (1) hgrid.gr3
!     Output: openbnd.bp (the depth is segment #)

!     ifort -Bstatic -assume byterecl -O3 -o check_openbnd check_openbnd.f90
!     pgf90 -mcmodel=medium  -Bstatic -o check_openbnd check_openbnd.f90

      implicit real*8(a-h,o-z)
      parameter(nbyte=4)
      allocatable :: xnd(:),ynd(:),dp(:),eta(:,:),nond(:),iond(:,:)

      open(14,file='hgrid.gr3',status='old')
      open(12,file='openbnd.bp',status='replace')
      read(14,*); read(14,*)ne,np
      allocate(xnd(np),ynd(np),dp(np),stat=istat)
      if(istat/=0) stop 'Failed to alloc (1)'
      do i=1,np
        read(14,*)j,xnd(i),ynd(i),dp(i)
      enddo !i
      do i=1,ne
        read(14,*)j,k !,elnode(1:3,i)
      enddo !i
      read(14,*)nope
      read(14,*)neta
      nlines_b4=2+np+ne+2 !for rewind and re-read the open bnd part
      allocate(nond(nope))
      do i=1,nope
        read(14,*)nond(i)
        do j=1,nond(i)
          read(14,*) !iond(j)
        enddo !j
      enddo !i
      mnond=maxval(nond(:))
      nond_all=sum(nond(:))
      allocate(iond(nope,mnond))
      rewind(14)
      do l=1,nlines_b4
        read(14,*)
      enddo !l

      write(12,*)'Open bnd nodes'
      write(12,*)nond_all
      icount=0
      do i=1,nope
        read(14,*) !nond(i)
        do j=1,nond(i)
          read(14,*) iond(i,j)
          icount=icount+1
          write(12,*)icount,xnd(iond(i,j)),ynd(iond(i,j)),i
        enddo !j
      enddo !i
      close(14)

      if(icount/=nond_all) stop 'Mismatch'

      stop
      end


text:
C Copyright(C) 2009-2017 National Technology & Engineering Solutions of
C Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C=======================================================================
      SUBROUTINE HIDMID (NLNKF, LINKF1, XN, YN, ZN, LINSET,
     &   IPSET, TVHMAX, ICROSS, XP, YP, ZP, NPART)
C=======================================================================

C   --*** HIDMID *** (MESH) Hide nodes beneath a face
C   --   Written by Amy Gilkey - revised 02/24/88
C   --
C   --HIDMID marks a node as hidden if it is beneath a face.
C   --
C   --Parameters:
C   --   NLNKF - IN - the number of nodes per face
C   --   LINKF1 - IN - the connectivity for the face
C   --   XN, YN, ZN - IN - the nodal coordinates
C   --   LINSET - IN - the sorted line set
C   --   IPSET - IN/OUT - the indices of the partial line set
C   --   NPART - IN/OUT - the number of lines in the partial line set
C   --
C   --Common Variables:

      include 'debug.blk'
      include 'dbnums.blk'
      include 'd3nums.blk'

      INTEGER LINKF1(NLNKF)
      REAL XN(*), YN(*), ZN(*)
      INTEGER LINSET(LLNSET,*)
      INTEGER IPSET(*)
      REAL TVHMAX(*)
      INTEGER ICROSS(*)
      REAL XP(*), YP(*), ZP(*)

      REAL X(4), Y(4), Z(4)
      REAL XNORM(4), YNORM(4), ZNORM(4)
      LOGICAL FIRSIN, FIRSTZ, FIRST(4)

C   --Calculate enclosing X-Y-Z box for face

      DO 100 ILINK = 1, NLNKF
         X(ILINK) = XN(LINKF1(ILINK))
         Y(ILINK) = YN(LINKF1(ILINK))
         Z(ILINK) = ZN(LINKF1(ILINK))
  100 CONTINUE

      XMIN = MIN (X(1), X(2), X(3), X(4))
      XMAX = MAX (X(1), X(2), X(3), X(4))
      YMIN = MIN (Y(1), Y(2), Y(3), Y(4))
      YMAX = MAX (Y(1), Y(2), Y(3), Y(4))
      ZMIN = MIN (Z(1), Z(2), Z(3), Z(4))
      ZMAX = MAX (Z(1), Z(2), Z(3), Z(4))

      FIRSIN = .TRUE.
      FIRSTZ = .TRUE.

C   --Check for nodes behind the face, mark as hidden

      IP = 1
  110 CONTINUE
      IF (IP .LE. NPART) THEN

C      --Check if outside of box enclosing face

         IF (XMAX .LT. XP(IP)) GOTO 120
         IF (XMIN .GT. XP(IP)) GOTO 120
         IF (YMAX .LT. YP(IP)) GOTO 120
         IF (YMIN .GT. YP(IP)) GOTO 120
         IF (ZMAX .LT. ZP(IP)) GOTO 120

         IV = LINSET(1,IPSET(IP))
         IH = LINSET(2,IPSET(IP))
         IF (((IV .EQ. LINKF1(1)) .OR. (IV .EQ. LINKF1(2)) .OR.
     &        (IV .EQ. LINKF1(3)) .OR. (IV .EQ. LINKF1(4))) .AND.
     &       ((IH .EQ. LINKF1(1)) .OR. (IH .EQ. LINKF1(2)) .OR.
     &        (IH .EQ. LINKF1(3)) .OR. (IH .EQ. LINKF1(4)))) GOTO 120

         IF (FIRSIN) THEN
            FIRSIN = .FALSE.

C         --Pre-compute factors for "within" test

            X2X1 = X(2) - X(1)
            Y2Y1 = Y(2) - Y(1)
            XY21 = X(1)*Y(2) - X(2)*Y(1)
            X3X2 = X(3) - X(2)
            Y3Y2 = Y(3) - Y(2)
            XY32 = X(2)*Y(3) - X(3)*Y(2)
            X4X3 = X(4) - X(3)
            Y4Y3 = Y(4) - Y(3)
            XY43 = X(3)*Y(4) - X(4)*Y(3)
            X1X4 = X(1) - X(4)
            Y1Y4 = Y(1) - Y(4)
            XY14 = X(4)*Y(1) - X(1)*Y(4)
         END IF

C      --Check if all 4 triangles formed by 2 nodes of the face
C      --and the node have a positive area

         X0 = XP(IP)
         Y0 = YP(IP)

         A12 = X2X1*Y0 - X0*Y2Y1 + XY21
         IF (A12 .LT. 0.0) GOTO 120
         A23 = X3X2*Y0 - X0*Y3Y2 + XY32
         IF (A23 .LT. 0.0) GOTO 120
         A34 = X4X3*Y0 - X0*Y4Y3 + XY43
         IF (A34 .LT. 0.0) GOTO 120
         A41 = X1X4*Y0 - X0*Y1Y4 + XY14
         IF (A41 .LT. 0.0) GOTO 120

C      --Check if the node's Z is behind the face's Z at that point

         IF (FIRSTZ) THEN
            FIRSTZ = .FALSE.
            FIRST(1) = .TRUE.
            FIRST(2) = .TRUE.
            FIRST(3) = .TRUE.
            FIRST(4) = .TRUE.

C         --Calculate epsilon based on length of "standard" side

            A = A12 + A23 + A34 + A41
            EPS = A * SQRT (A) * .001

C         --Calculate center of face

            XCEN = 0.25 * (X(1) + X(2) + X(3) + X(4))
            YCEN = 0.25 * (Y(1) + Y(2) + Y(3) + Y(4))
            ZCEN = 0.25 * (Z(1) + Z(2) + Z(3) + Z(4))
         END IF

C      --Check if the node's Z is behind the face's minimum Z

         IF ((ZMIN - ZP(IP)) .LT. EPS) THEN

C         --Calculate normal for center of face with two nodes
C         --that have the smallest area with the given node

            AMIN = MIN (A12, A23, A34, A41)

            IF (AMIN .EQ. A12) THEN
               ISIDE = 1
            ELSE IF (AMIN .EQ. A23) THEN
               ISIDE = 2
            ELSE IF (AMIN .EQ. A34) THEN
               ISIDE = 3
            ELSE IF (AMIN .EQ. A41) THEN
               ISIDE = 4
            END IF
            IF (FIRST(ISIDE)) THEN
               FIRST(ISIDE) = .FALSE.
               IF (ISIDE .LT. 4) THEN
                  N2 = ISIDE + 1
               ELSE
                  N2 = 1
               END IF
               AX = X(ISIDE) - XCEN
               AY = Y(ISIDE) - YCEN
               AZ = Z(ISIDE) - ZCEN
               BX = X(N2) - XCEN
               BY = Y(N2) - YCEN
               BZ = Z(N2) - ZCEN
               XNORM(ISIDE) = 0.5 * (AY*BZ - BY*AZ)
               YNORM(ISIDE) = 0.5 * (AZ*BX - BZ*AX)
               ZNORM(ISIDE) = 0.5 * (AX*BY - BX*AY)
            END IF

C         --Check if the node's Z is behind the face's Z at that point

C#????d           z0c = xnorm(iside)*(x0-xcen) + ynorm(iside)*(y0-ycen) !#VAX
C#????d    &         + znorm(iside)*(zp(ip)-zcen) !#VAX
C#????d           write (99, 9001, iostat=idum) inp, LINKF1, !#VAX
C#????d    &         a12, a23, a34, a41, eps, z0c, (z0c.ge.-eps) !#VAX
C#????d9001        format (5i5, 5f8.4, 1x, f8.4, (2x, l1)) !#VAX

            IF ((XNORM(ISIDE)*(X0-XCEN) + YNORM(ISIDE)*(Y0-YCEN)
     &         + ZNORM(ISIDE)*(ZP(IP)-ZCEN)) .GE. -EPS) GOTO 120
         END IF

C      --Hide the entire line (as is), and move the line to hidden section

C#????d        write (99, *) 'invisible line', (linset(i,ipset(ip)), i=1,3) !#VAX
         I = IPSET(IP)
         IPSET(IP) = IPSET(NPART)
         IPSET(NPART) = I
         T = TVHMAX(IP)
         TVHMAX(IP) = TVHMAX(NPART)
         TVHMAX(NPART) = T
         I = ICROSS(IP)
         ICROSS(IP) = ICROSS(NPART)
         ICROSS(NPART) = I
         XP(IP) = XP(NPART)
         YP(IP) = YP(NPART)
         ZP(IP) = ZP(NPART)
         NPART = NPART - 1
         IP = IP - 1

  120    CONTINUE
         IP = IP + 1
         GOTO 110
      END IF

      RETURN
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


program hello_world

  print *, 'Hello world!'

end program