text:
c-----------------------------------------------------------------------------
c
c Copyright (c) 2000 - 2014, Lawrence Livermore National Security, LLC
c Produced at the Lawrence Livermore National Laboratory
c LLNL-CODE-442911
c All rights reserved.
c
c This file is  part of VisIt. For  details, see https://visit.llnl.gov/.  The
c full copyright notice is contained in the file COPYRIGHT located at the root
c of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
c
c Redistribution  and  use  in  source  and  binary  forms,  with  or  without
c modification, are permitted provided that the following conditions are met:
c
c  - Redistributions of  source code must  retain the above  copyright notice,
c    this list of conditions and the disclaimer below.
c  - Redistributions in binary form must reproduce the above copyright notice,
c    this  list of  conditions  and  the  disclaimer (as noted below)  in  the
c    documentation and/or other materials provided with the distribution.
c  - Neither the name of  the LLNS/LLNL nor the names of  its contributors may
c    be used to endorse or promote products derived from this software without
c    specific prior written permission.
c
c THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT  HOLDERS AND CONTRIBUTORS "AS IS"
c AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO, THE
c IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR  PURPOSE
c ARE  DISCLAIMED. IN  NO EVENT  SHALL LAWRENCE  LIVERMORE NATIONAL  SECURITY,
c LLC, THE  U.S.  DEPARTMENT OF  ENERGY  OR  CONTRIBUTORS BE  LIABLE  FOR  ANY
c DIRECT,  INDIRECT,   INCIDENTAL,   SPECIAL,   EXEMPLARY,  OR   CONSEQUENTIAL
c DAMAGES (INCLUDING, BUT NOT  LIMITED TO, PROCUREMENT OF  SUBSTITUTE GOODS OR
c SERVICES; LOSS OF  USE, DATA, OR PROFITS; OR  BUSINESS INTERRUPTION) HOWEVER
c CAUSED  AND  ON  ANY  THEORY  OF  LIABILITY,  WHETHER  IN  CONTRACT,  STRICT
c LIABILITY, OR TORT  (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY  WAY
c OUT OF THE  USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
c DAMAGE.
c
c-----------------------------------------------------------------------------

c-----------------------------------------------------------------
c Program: main
c
c Programmer: Brad Whitlock
c Date:       Thu Mar 11 12:13:06 PST 2010
c
c Modifications:
c
c-----------------------------------------------------------------
      program main
      implicit none
      include "visitfortransimV2interface.inc"
ccc   local variables
      integer err

      call simulationarguments()
      err = visitsetupenv()
      err = visitinitializesim("fmaterial", 9,
     . "Demonstrates material data access function", 42,
     . "/no/useful/path", 15,
     . VISIT_F77NULLSTRING, VISIT_F77NULLSTRINGLEN,
     . VISIT_F77NULLSTRING, VISIT_F77NULLSTRINGLEN,
     . VISIT_F77NULLSTRING, VISIT_F77NULLSTRINGLEN)
      call mainloop()
      stop
      end

c-----------------------------------------------------------------
c mainloop
c-----------------------------------------------------------------
      subroutine mainloop()
      implicit none
      include "visitfortransimV2interface.inc"
ccc   local variables
      integer visitstate, result, blocking
ccc   SIMSTATE common block
      integer runflag, simcycle
      double precision simtime
      common /SIMSTATE/ runflag, simcycle, simtime
      save /SIMSTATE/

c     main loop
      runflag = 1
      simcycle = 0
      simtime = 0.
      do 10
          if(runflag.eq.1) then
              blocking = 0
          else
              blocking = 1
          endif

          visitstate = visitdetectinput(blocking, -1)

          if (visitstate.lt.0) then
              goto 1234
          elseif (visitstate.eq.0) then
              call simulate_one_timestep()
          elseif (visitstate.eq.1) then
              runflag = 0
              result = visitattemptconnection()
              if (result.eq.1) then
                  write (6,*) 'VisIt connected!'
              else
                  write (6,*) 'VisIt did not connect!'
              endif
          elseif (visitstate.eq.2) then
              runflag = 0
              if (visitprocessenginecommand().eq.0) then
                  result = visitdisconnect()
                  runflag = 1
              endif
          endif
10    continue
1234  end

      subroutine simulate_one_timestep()
ccc   SIMSTATE common block
      integer runFlag, simcycle
      double precision simtime
      common /SIMSTATE/ runflag, simcycle, simtime
c Simulate one time step
      simcycle = simcycle + 1
      simtime = simtime + 0.0134
      write (6,*) 'Simulating time step: cycle=',simcycle, ' time=', simtime
      call sleep(1)
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c These functions must be defined to satisfy the visitfortransimV2interface lib.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c---------------------------------------------------------------------------
c visitcommandcallback
c---------------------------------------------------------------------------
      subroutine visitcommandcallback (cmd, lcmd, args, largs)
      implicit none
      character*8 cmd, args
      integer     lcmd, largs
      include "visitfortransimV2interface.inc"
ccc   SIMSTATE common block
      integer runflag, simcycle
      double precision simtime
      common /SIMSTATE/ runflag, simcycle, simtime

c     Handle the commands that we define in visitgetmetadata.
      if(visitstrcmp(cmd, lcmd, "halt", 4).eq.0) then
          runflag = 0
      elseif(visitstrcmp(cmd, lcmd, "step", 4).eq.0) then
          call simulate_one_timestep()
      elseif(visitstrcmp(cmd, lcmd, "run", 3).eq.0) then
          runflag = 1
      endif
      end

c---------------------------------------------------------------------------
c visitbroadcastintfunction
c---------------------------------------------------------------------------
      integer function visitbroadcastintfunction(value, sender)
      implicit none
      integer value, sender
c     REPLACE WITH MPI COMMUNICATION IF SIMULATION IS PARALLEL
      visitbroadcastintfunction = 0
      end

c---------------------------------------------------------------------------
c visitbroadcaststringfunction
c---------------------------------------------------------------------------
      integer function visitbroadcaststringfunction(str, lstr, sender)
      implicit none
      character*8 str
      integer     lstr, sender
c     REPLACE WITH MPI COMMUNICATION IF SIMULATION IS PARALLEL
      visitbroadcaststringfunction = 0
      end

c---------------------------------------------------------------------------
c visitslaveprocesscallback
c---------------------------------------------------------------------------
      subroutine visitslaveprocesscallback ()
      implicit none
c     REPLACE WITH MPI COMMUNICATION IF SIMULATION IS PARALLEL
      end

c---------------------------------------------------------------------------
c visitactivatetimestep
c---------------------------------------------------------------------------
      integer function visitactivatetimestep()
      implicit none
      include "visitfortransimV2interface.inc"
      visitactivatetimestep = VISIT_OKAY
      end

c---------------------------------------------------------------------------
c visitgetmetadata
c---------------------------------------------------------------------------
      integer function visitgetmetadata()
      implicit none
      include "visitfortransimV2interface.inc"
ccc   SIMSTATE common block
      integer runflag, simcycle
      double precision simtime
      common /SIMSTATE/ runflag, simcycle, simtime
      integer md, mmd, mat, cmd, err

      if(visitmdsimalloc(md).eq.VISIT_OKAY) then
          err = visitmdsimsetcycletime(md, simcycle, simtime)
          if(runflag.eq.1) then
              err = visitmdsimsetmode(md, VISIT_SIMMODE_RUNNING)
          else
              err = visitmdsimsetmode(md, VISIT_SIMMODE_STOPPED)
          endif

c     Add a 2D rectilinear mesh
          err = visitmdmeshalloc(mmd);
          if(err.eq.VISIT_OKAY) then
              err = visitmdmeshsetname(mmd, "mesh2d", 6)
              err = visitmdmeshsetmeshtype(mmd,
     .            VISIT_MESHTYPE_RECTILINEAR)
              err = visitmdmeshsettopologicaldim(mmd, 2)
              err = visitmdmeshsetspatialdim(mmd, 2)
              err = visitmdmeshsetnumdomains(mmd, 1)
              err = visitmdmeshsetdomaintitle(mmd, "Domains", 7)
              err = visitmdmeshsetdomainpiecename(mmd, "domain", 6)
              err = visitmdmeshsetxunits(mmd, "cm", 2)
              err = visitmdmeshsetyunits(mmd, "cm", 2)
              err = visitmdmeshsetxlabel(mmd, "Width", 5)
              err = visitmdmeshsetylabel(mmd, "Height", 6)
              err = visitmdmeshsetcellorigin(mmd, 1)
              err = visitmdmeshsetnodeorigin(mmd, 1)
              err = visitmdsimaddmesh(md, mmd)
          endif

c     Add a material
          err = visitmdmatalloc(mat)
          if(err.eq.VISIT_OKAY) then
              err = visitmdmatsetname(mat, "mat", 3)
              err = visitmdmatsetmeshname(mat, "mesh2d", 6)
              err = visitmdmataddmaterialname(mat, "Water", 5)
              err = visitmdmataddmaterialname(mat,
     .           "Membrane", 8)
              err = visitmdmataddmaterialname(mat,
     .           "Air", 3)
              err = visitmdsimaddmaterial(md, mat)
          endif

c     Add simulation commands
          err = visitmdcmdalloc(cmd)
          if(err.eq.VISIT_OKAY) then
              err = visitmdcmdsetname(cmd, "halt", 4)
              err = visitmdsimaddgenericcommand(md, cmd)
          endif
          err = visitmdcmdalloc(cmd)
          if(err.eq.VISIT_OKAY) then
              err = visitmdcmdsetname(cmd, "step", 4)
              err = visitmdsimaddgenericcommand(md, cmd)
          endif
          err = visitmdcmdalloc(cmd)
          if(err.eq.VISIT_OKAY) then
              err = visitmdcmdsetname(cmd, "run", 3)
              err = visitmdsimaddgenericcommand(md, cmd)
          endif
      endif
      visitgetmetadata = md
      end

c---------------------------------------------------------------------------
c visitgetmesh
c---------------------------------------------------------------------------
      integer function visitgetmesh(domain, name, lname)
      implicit none
      character*8 name
      integer     domain, lname
      include "visitfortransimV2interface.inc"
ccc   RECTMESH
      integer NX, NY
      parameter (NX = 5)
      parameter (NY = 4)

      real rmx(NX), rmy(NY)
      data rmx/0., 1., 2., 3., 4./
      data rmy/0., 1., 2., 3./
ccc   local variables
      integer h, x, y, err

      h = VISIT_INVALID_HANDLE
      if(visitstrcmp(name, lname, "mesh2d", 6).eq.0) then
c Create a rectilinear mesh here
          if(visitrectmeshalloc(h).eq.VISIT_OKAY) then
              err = visitvardataalloc(x)
              err = visitvardataalloc(y)
              err = visitvardatasetf(x,VISIT_OWNER_COPY,1,NX,rmx)
              err = visitvardatasetf(y,VISIT_OWNER_COPY,1,NY,rmy)
              err = visitrectmeshsetcoordsxy(h, x, y)
          endif
      endif
      visitgetmesh = h
      end

c---------------------------------------------------------------------------
c visitgetvariable
c---------------------------------------------------------------------------
      integer function visitgetvariable(domain, name, lname)
      implicit none
      character*8 name
      integer     domain, lname
      include "visitfortransimV2interface.inc"
      visitgetvariable = VISIT_INVALID_HANDLE
      end

c---------------------------------------------------------------------------
c visitgetcurve
c---------------------------------------------------------------------------
      integer function visitgetcurve(name, lname)
      implicit none
      character*8 name
      integer     lname
      include "visitfortransimV2interface.inc"
      visitgetcurve = VISIT_INVALID_HANDLE
      end

c---------------------------------------------------------------------------
c visitgetdomainlist
c---------------------------------------------------------------------------
      integer function visitgetdomainlist(name, lname)
      implicit none
      character*8 name
      integer     lname
      include "visitfortransimV2interface.inc"
      visitgetdomainlist = VISIT_INVALID_HANDLE
      end

c---------------------------------------------------------------------------
c visitgetdomainbounds
c---------------------------------------------------------------------------
      integer function visitgetdomainbounds(name, lname)
      implicit none
      character*8 name
      integer     lname
      include "visitfortransimV2interface.inc"
      visitgetdomainbounds = VISIT_INVALID_HANDLE
      end

c---------------------------------------------------------------------------
c visitgetdomainnesting
c---------------------------------------------------------------------------
      integer function visitgetdomainnesting(name, lname)
      implicit none
      character*8 name
      integer     lname
      include "visitfortransimV2interface.inc"
      visitgetdomainnesting = VISIT_INVALID_HANDLE
      end

c---------------------------------------------------------------------------
c visitgetmaterial
c---------------------------------------------------------------------------
      integer function visitgetmaterial(domain, name, lname)
      implicit none
      character*8 name
      integer     domain, lname
      include "visitfortransimV2interface.inc"
ccc   RECTMESH
      integer NX, NY
      parameter (NX = 5)
      parameter (NY = 4)

c The matlist table indicates the material numbers that are found in
c each cell. Every 3 numbers indicates the material numbers in a cell.
c A material number of 0 means that the material entry is not used.
      integer matlist(3, NX-1, NY-1)
      data matlist/3,0,0,2,3,0,1,2,0,1,0,0,
     .             3,0,0,2,3,0,1,2,0,1,0,0,
     .             3,0,0,2,3,0,1,2,3,1,2,0/
c The mat_vf table indicates the material volume fractions that are
c found in a cell.
      real mat_vf(3, NX-1, NY-1)
      data mat_vf/1.,0.,0.,0.75,0.25,0.,0.8125,0.1875, 0.,1.,0.,0.,
     . 1.,0.,0.,0.625,0.375,0.,0.5625,0.4375,0.,1.,0.,0.,
     . 1.,0.,0.,0.3,0.7,0.,0.2,0.4,0.4,0.55,0.45,0./

c Size the material object so it has the right dimensions (equal to
c the number of cells in the mesh. Unstructured grids would have
c /ncells,1,1/ size for the material. That means that the cellid array
c would have one linear index in its first element.
      integer err, I, J, m, cellid, nmats, h
      integer matno(3), cellmat(3)
      real cellmatvf(3)

      h = VISIT_INVALID_HANDLE
      err = visitmatdataalloc(h)
      err = visitmatdataappendcells(h, (NX-1) * (NY-1))
      err = visitmatdataaddmat(h, "Water", 5, matno(1))
      err = visitmatdataaddmat(h, "Membrane", 8, matno(2))
      err = visitmatdataaddmat(h, "Air", 3, matno(3))

      cellid = 0
      do 2020 J=1,NY-1
          do 2010 I=1,NX-1
              nmats = 0
              do 2000 m=1,3
                  if(matlist(m,I,J).gt.0) then
                      nmats = nmats + 1
                      cellmat(nmats) = matno(matlist(m,I,J))
                      cellmatvf(nmats) = mat_vf(m,I,J)
                  endif
2000          continue
              if(nmats.gt.1) then
                  err = visitmatdataaddmixedcell(h,
     . cellid, cellmat, cellmatvf, nmats)
              else
                  err = visitmatdataaddcleancell(h,
     . cellid, cellmat)
              endif
              cellid = cellid + 1
2010      continue
2020  continue

      visitgetmaterial = h
      end


text:
C Copyright(C) 2011-2017 National Technology & Engineering Solutions
C of 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
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 -*- Mode: fortran -*-
C=======================================================================
      SUBROUTINE NEWXYZ (XN, YN, XN3, YN3, ZN3, ATRIB, A)
C=======================================================================

C   --*** NEWXYZ *** (GENSHELL) Calculate 3D coordinates
C   --
C   --NEWXYZ calculates the coordinate array for the 3D database.
C   --
C   --Parameters:
C   --   XN, YN - IN - the 2D coordinates, destroyed
C   --   XN3, YN3, ZN3 - OUT - the 3D coordinates
C   --
C   --Common Variables:
C   --   Uses NDIM, NUMNP of /DBNUMS/
C   --   Uses NDIM3, NUMNP3 of /DBNUM3/
C   --   Uses ITRANT, NNREPL, DIM3, NRTRAN, D3TRAN, ZGRAD,
C   --      CENTER, NUMCOL, NUMROW of /PARAMS/
C   --   Uses XOFFS, YOFFS, ZOFFS of /XYZOFF/
C   --   Uses ROT3D, ROTMAT of /XYZROT/

      INCLUDE 'gs_dbnums.blk'
      INCLUDE 'gs_dbnum3.blk'
      INCLUDE 'gs_params.blk'
      INCLUDE 'gs_xyzoff.blk'
      INCLUDE 'gs_xyzrot.blk'
      INCLUDE 'gs_xyzmir.blk'
      INCLUDE 'gs_xyzero.blk'
      INCLUDE 'gs_xyzscl.blk'
      INCLUDE 'gs_splxyz.blk'

      REAL XN(NUMNP), YN(NUMNP),
     &   XN3(NUMNP3), YN3(NUMNP3), ZN3(NUMNP3)
      REAL ATRIB(*)
      REAL A(*)

      IF (ITRANT .EQ. 1) THEN
         CALL TRNXYZ (XN, YN, XN3, YN3, ZN3, ATRIB)
C$$$      ELSE IF (ITRANT .EQ. 2) THEN
C$$$         CALL ARCXYZ (XN, YN, XN3, YN3, ZN3, ATRIB)
      ELSE IF (ITRANT .EQ. 4) THEN
         CALL WRPXYZ (XN, YN, XN3, YN3, ZN3, ATRIB)
C$$$      ELSE IF (ITRANT .EQ. 8) THEN
C$$$         CALL TWIXYZ (XN, YN, XN3, YN3, ZN3, ATRIB)
c$$$      ELSE IF (ITRANT .EQ. 16) THEN
c$$$         CALL PROXYZ (XN, YN, XN3, YN3, ZN3, ATRIB)
c$$$      ELSE IF (ITRANT .EQ. 32) THEN
c$$$         CALL EXPARC (XN, YN, XN3, YN3, ZN3, ATRIB)
      ELSE IF (ITRANT .EQ. 64) THEN
         CALL SPLXYZ (XN, YN, XN3, YN3, ZN3, ATRIB,
     $      NSPL(1), A(KRSPLA), A(KZSPLA), A(KSPL2A), A(KDISTA),
     $        A(KSCRA), SLLFT, SLRGT, RDTHET, SWEEP, NOSCAL)
c$$$      ELSE IF (ITRANT .EQ. 128) THEN
c$$$         CALL SPTXYZ (XN, YN, XN3, YN3, ZN3, ATRIB, NSPLT,
c$$$     $        A(KZSPL), A(KXSPL), A(KXSPL2), A(KYSPL), A(KYSPL2),
c$$$     $        A(KSCR))
      END IF

C   --Revolve 3D mesh, if needed

      IF (ROT3D) THEN
         DO 30 JNP = 1, NUMNP3
            X = XN3(JNP) - ROTCEN(1)
            Y = YN3(JNP) - ROTCEN(2)
            Z = ZN3(JNP) - ROTCEN(3)
            XN3(JNP) = X*ROTMAT(1,1) + Y*ROTMAT(2,1) + Z*ROTMAT(3,1)
     &         + ROTCEN(1)
            YN3(JNP) = X*ROTMAT(1,2) + Y*ROTMAT(2,2) + Z*ROTMAT(3,2)
     &         + ROTCEN(2)
            ZN3(JNP) = X*ROTMAT(1,3) + Y*ROTMAT(2,3) + Z*ROTMAT(3,3)
     &         + ROTCEN(3)
   30    CONTINUE
      END IF

C   --Add offset, if any, to coordinates

      IF (XOFFS .NE. 0.0) THEN
         DO 40 JNP = 1, NUMNP3
            XN3(JNP) = XN3(JNP) + XOFFS
   40    CONTINUE
      END IF
      IF (YOFFS .NE. 0.0) THEN
         DO 50 JNP = 1, NUMNP3
            YN3(JNP) = YN3(JNP) + YOFFS
   50    CONTINUE
      END IF
      IF (ZOFFS .NE. 0.0) THEN
         DO 60 JNP = 1, NUMNP3
            ZN3(JNP) = ZN3(JNP) + ZOFFS
   60    CONTINUE
      END IF

C   --Mirror coordinates if any specified

      IF (XMIRR .LT. 0.) THEN
         DO 70 JNP = 1, NUMNP3
            XN3(JNP) = -1.0 * XN3(JNP)
   70    CONTINUE
      END IF
      IF (YMIRR .LT. 0.) THEN
         DO 80 JNP = 1, NUMNP3
            YN3(JNP) = -1.0 * YN3(JNP)
   80    CONTINUE
      END IF
      IF (ZMIRR .LT. 0.) THEN
         DO 90 JNP = 1, NUMNP3
            ZN3(JNP) = -1.0 * ZN3(JNP)
   90    CONTINUE
      END IF

C     --Randomize coordinates if any specified

      IDUM = 1
      IF (XRAND .NE. 0.) THEN
         DO 61 JNP = 1, NUMNP
            XN3(JNP) = (2.0*RAN1(IDUM)-1.0) * XRAND + XN3(JNP)
 61      CONTINUE
      END IF
      IF (YRAND .NE. 0.) THEN
         DO 71 JNP = 1, NUMNP
            YN3(JNP) = (2.0*RAN1(IDUM)-1.0) * YRAND + YN3(JNP)
 71      CONTINUE
      END IF
      IF (ZRAND .NE. 0.0) THEN
         DO 81 JNP = 1, NUMNP
            ZN3(JNP) = (2.0*RAN1(IDUM)-1.0) * ZRAND + ZN3(JNP)
 81      CONTINUE
      END IF

C --- Zero coordinates if *ZERO is not equal to zero

      IF (XZERO .NE. 0.) THEN
         DO 100 JNP = 1, NUMNP3
            XN3(JNP) = (0.5+SIGN( 0.5, ABS(XN3(JNP))-XZERO)) * XN3(JNP)
  100    CONTINUE
      END IF
      IF (YZERO .NE. 0.) THEN
         DO 110 JNP = 1, NUMNP3
            YN3(JNP) = (0.5+SIGN( 0.5, ABS(YN3(JNP))-YZERO)) * YN3(JNP)
  110    CONTINUE
      END IF
      IF (ZZERO .NE. 0.) THEN
         DO 120 JNP = 1, NUMNP3
            ZN3(JNP) = (0.5+SIGN( 0.5, ABS(ZN3(JNP))-ZZERO)) * ZN3(JNP)
  120    CONTINUE
      END IF

C --- Scale the coordinates if any Scaled

      IF (XSCAL .NE. 1.) THEN
         DO 130 JNP = 1, NUMNP
            XN3(JNP) = XSCAL * XN3(JNP)
  130    CONTINUE
      END IF
      IF (YSCAL .NE. 1.) THEN
         DO 140 JNP = 1, NUMNP
            YN3(JNP) = YSCAL * YN3(JNP)
  140    CONTINUE
      END IF
      IF (ZSCAL .NE. 1. .AND. NDIM .EQ. 3) THEN
         DO 150 JNP = 1, NUMNP
            ZN3(JNP) = ZSCAL * ZN3(JNP)
  150    CONTINUE
      END IF

      RETURN
      END

ext: f

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


! { dg-do compile }
! PR fortran/32669
!
! Contributed by Janus Weil <jaydub66@gmail.com>
!
program tfe
implicit none

real,dimension(-1:1) ::  w
real,dimension(1:4) ::  x
real,dimension(0:3) ::  y
real,dimension(-1:2) ::  z

call sub(x(:))
call sub(y(:))
call sub(z(:))
call sub(w(:)) ! { dg-warning "too few elements" }

contains
  subroutine sub(a)
    implicit none
    real,dimension(1:4) :: a
  end subroutine sub
end program tfe


text:
! Test of fix for PR18878
!
! Based on example in PR by Steve Kargl
!
module a
  integer, parameter :: b = kind(1.d0)
  real(b)            :: z
end module a
program d
  use a, only : e => b, f => b, u => z, v => z
  real(e) x
  real(f) y
  x = 1.e0_e
  y = 1.e0_f
  u = 99.0
  if (kind(x).ne.kind(y)) call abort ()
  if (v.ne.u) call abort ()
end program d

! { dg-final { cleanup-modules "a" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
program driver
  !
  !-----------------------------------------------------------
  ! 2014/06/13: Original (Ryuji Yoshida)
  !
  !-----------------------------------------------------------

  use mpi !external
  use mod_vars
  use mod_prof
  use mod_comm

  implicit none

  integer :: i, j, k, l, t
  integer :: value
  integer :: ierr

  namelist / PARAM_INTERCOMM / &
     NX,             &
     NY,             &
     NZ,             &
     NT,             &
     DOMAIN_NUM,     &
     IAM_PARENT,     &
     IAM_CHILD,      &
     LAUNCH_PRC
  !-----------------------------------------------------------

  call get_command_argument(1,fnamelist)

  ! default settings
  NX         = 8
  NY         = 8
  NZ         = 40
  NT         = 5
  DOMAIN_NUM = 0
  LAUNCH_PRC = 0
  IAM_PARENT = .false.
  IAM_CHILD  = .false.

  value      = 1

  ! read namelist
  open  ( FID_NML, file=trim(fnamelist), status='old', delim='apostrophe' )
  read  ( FID_NML, nml=PARAM_INTERCOMM )
  close ( FID_NML )

  ! initialize MPI
  call MPI_INIT( ierr )
  call MPI_COMM_RANK( MPI_COMM_WORLD, myproc, ierr )
  call MPI_COMM_SIZE( MPI_COMM_WORLD, nprocs, ierr )
  call prof_initialize
  call prof_start ( 'init' )

  ! open log file
  LFID = DOMAIN_NUM*1000 + myproc
  write(domnum,'(I2.2)') DOMAIN_NUM
  write(prcnum,'(I3.3)') myproc
  flogfile = "LOG.d"//domnum//".pe"//prcnum
  open  ( LFID, file=trim(flogfile), status='replace' )

  write ( LFID, * ) "---------------------------------------------------"
  write ( LFID, * ) "     START MINI TEST OF INTER COMMUNICATION"
  write ( LFID, * ) "---------------------------------------------------"
  write ( LFID, * ) ""
  write ( LFID, * ) "NAMELIST: DOMAIN_NUM = ", DOMAIN_NUM
  write ( LFID, * ) "NAMELIST: NX = ", NX
  write ( LFID, * ) "NAMELIST: NY = ", NY
  write ( LFID, * ) "NAMELIST: NZ = ", NZ
  write ( LFID, * ) "NAMELIST: NT = ", NT
  write ( LFID, * ) "NAMELIST: IAM_PARENT = ", IAM_PARENT
  write ( LFID, * ) "NAMELIST: IAM_CHILD = ", IAM_CHILD
  if ( IAM_PARENT ) then
     write ( LFID, * ) "NAMELIST: LAUNCH_PRC = ", LAUNCH_PRC
  endif
  write ( LFID, * ) ""
  write ( LFID, * ) "---------------------------------------------------"
  write ( LFID, * ) "MPI INITIALIZE: TOTAL PROCESSES = ", nprocs
  write ( LFID, * ) "MPI INITIALIZE: MY PROCESS NUMBER = ", myproc

  ! memory allocate
  allocate ( var  ( NX, NY, NZ ) )
  allocate ( dummy( 1,  1,  1  ) )

  ! launch child process
  write ( LFID, * ) ""
  write ( LFID, * ) "---------------------------------------------------"
  call prof_start ( 'launch' )
  if ( IAM_PARENT ) then
     call comm_spawn ( DOMAIN_NUM, LAUNCH_PRC )
  endif
  if ( IAM_CHILD ) then
     call comm_get_parent
  endif
  call prof_end ( 'launch' )

  ! initialize variables
  if ( IAM_PARENT ) then
     do k = 1, nz
     do j = 1, ny
     do i = 1, nx
        var(i,j,k) = float( myproc+1 )*1000 + float( DOMAIN_NUM )*100 + float(value)
     enddo
     enddo
     enddo
  else
     var(:,:,:) = 0.0D0
  endif

  ! setup inter-domain communication
  write ( LFID, * ) ""
  write ( LFID, * ) "---------------------------------------------------"
  if ( IAM_PARENT ) then
     call prof_start ( 'make yp' )
     call comm_makeyp ( nprocs, LAUNCH_PRC, parent )
     call prof_end ( 'make yp' )
  endif
  if ( IAM_CHILD ) then
     call prof_start ( 'make yp' )
     call comm_makeyp ( PARENT_SIZE, nprocs, child )
     call prof_end ( 'make yp' )
  endif

  sleep_loop = DOMAIN_NUM

  if ( IAM_PARENT ) then
     call prof_start ( 'barrier' )
     call comm_barrier ( parent )
     call prof_end ( 'barrier' )
  endif
  if ( IAM_CHILD ) then
     call prof_start ( 'barrier' )
     call comm_barrier ( child )
     call prof_end ( 'barrier' )
  endif

  call prof_end ( 'init' )
  call prof_start ( 'main' )

  !----------------------------------------------------------->>>
  do t = 1, NT
  write ( LFID, * ) ""
  write ( LFID, '(1X,A,I2,A,I2)' ) "time step: ", t, "/", NT
  write ( LFID, * ) "---------------------------------------------------"
  write ( LFID, '(1X,A,F9.3,A,F9.3)' ) &
        "before comm; var check: max = ", maxval(var), "  min = ", minval(var)

  ! issue send or recv
  if ( IAM_PARENT ) then
     call prof_start ( 'send' )
     call comm_send ( var, parent )
     call prof_end ( 'send' )
  endif
  if ( IAM_CHILD ) then
     call prof_start ( 'recv' )
     call comm_recv ( child )
     call prof_end ( 'recv' )
  endif

  ! dummy main routine
  do l = 1, sleep_loop
     write ( LFID, '(3X,A,I2,A,I2)' ) "DO DUMMY JOB: ", l, "/", sleep_loop
     call sleep(1)
  enddo

  ! update a value
  if ( IAM_PARENT ) then
     value = value + 1
     do k = 1, nz
     do j = 1, ny
     do i = 1, nx
        var(i,j,k) = float( myproc+1 )*1000 + float( DOMAIN_NUM )*100 + float(value)
     enddo
     enddo
     enddo
  endif

  ! issue wait
  if ( IAM_PARENT ) then
     call prof_start ( 'wait' )
     call comm_wait ( dummy, parent )
     call prof_end ( 'wait' )
  endif
  if ( IAM_CHILD ) then
     call prof_start ( 'wait' )
     call comm_wait ( var, child )
     call prof_end ( 'wait' )
  endif

  ! check variables
  write ( LFID, '(1X,A,F9.3,A,F9.3)' ) &
        "after comm; var check: max = ", maxval(var), "  min = ", minval(var)

  enddo
  !-----------------------------------------------------------<<<

  call prof_end ( 'main' )
  call prof_start ( 'finaliz' )

  ! finalize the process
  write ( LFID, * ) ""
  write ( LFID, * ) "---------------------------------------------------"
  write ( LFID, * ) "finalizing the process"

  if ( IAM_PARENT ) then
     call prof_start ( 'barrier' )
     call comm_barrier ( parent )
     call prof_end ( 'barrier' )
     call prof_start ( 'free' )
     call comm_free ( parent )
     call prof_end ( 'free' )
  endif
  if ( IAM_CHILD ) then
     call prof_start ( 'barrier' )
     call comm_barrier ( child )
     call prof_end ( 'barrier' )
     call prof_start ( 'free' )
     call comm_free ( child )
     call prof_end ( 'free' )
  endif
  call prof_end ( 'finaliz' )

  ! report profile
  write ( LFID, * ) ""
  write ( LFID, * ) "---------------------------------------------------"
  call prof_finalize

  call MPI_FINALIZE(ierr)
  write ( LFID, * ) ""
  write ( LFID, * ) "---------------------------------------------------"
  write ( LFID, * ) "finished"
  close ( LFID )
  stop

  !-----------------------------------------------------------
end program driver

max_stars_repo_name: slayoo/scale

max_issues_repo_name: slayoo/scale

max_forks_repo_name: slayoo/scale


text:
!
! Updated on 6/3/2013
! by Layne Watson
!

MODULE LINEAR_SHEPARD_MOD
INTERFACE
   SUBROUTINE DGELSS(M, N, NRHS, A, LDA, B, LDB, SIGMA, RCOND, RANK, &
        WORK, LWORK, INFO)
     USE REAL_PRECISION
     INTEGER, INTENT(IN) :: M, N, NRHS, LDA, LDB, LWORK
     INTEGER, INTENT(OUT) :: RANK, INFO
     REAL(KIND = R8), INTENT(IN) :: RCOND
     REAL(KIND = R8), DIMENSION(LDA, *), INTENT(INOUT) :: A
     REAL(KIND = R8), DIMENSION(LDB, *), INTENT(INOUT) :: B
     REAL(KIND = R8), DIMENSION(N), INTENT(OUT) :: SIGMA
     REAL(KIND = R8), DIMENSION(LWORK), INTENT(OUT) :: WORK
   END SUBROUTINE DGELSS
END INTERFACE
PRIVATE
PUBLIC LSHEP, LSHEPVAL, RIPPLE
CONTAINS

SUBROUTINE LSHEP ( M, N, X, F, A, RW, IER, RLSHEP )
! This subroutine computes a set of parameters defining a function that
! interpolates N data values F(i) at scattered nodes X(i) in M dimensions. The
! interpolant may be evaluated at an arbitrary point by the function LSHEPVAL.
! The interpolation scheme is a modified linear Shepard method, and
! can use either the Shepard weighted least squares fit or robust
! M-estimation for each local approximation.
!
! Input parameters:
!   M is the dimension of the data.
!   N is the number of nodes.
!   X(M, N) contains the coordinates of the nodes.
!   F(N) contains the function values of the nodes.
! Output parameters:
!   A(M, N) contains the coefficients of the local fits.
!   RW(N) is an array containing the radius of influence for each point.
!   IER
!   = 0, if no errors were encountered.
!   = 1, if N is too small relative to M.
!   = 2, if any least squares problem is rank deficient, in which
!        case an SVD based minimum norm solution is returned; the
!        local fit should still be reasonable.
!   = 3, if the IRLS subroutine returns an error, meaning the IRLS
!        iteration has failed, and the returned local fit may
!        be bad.
! Optional input:
!   RLSHEP specifies by its presence that robust M-estimation is to be used to
!          compute the local approximation.
USE REAL_PRECISION
IMPLICIT NONE

INTEGER, INTENT(IN) :: M, N
REAL(KIND=R8), DIMENSION(M,N), INTENT(IN) :: X
REAL(KIND=R8), DIMENSION(N), INTENT(IN) :: F
REAL(KIND=R8), DIMENSION(M,N), INTENT(OUT) :: A
REAL(KIND=R8), DIMENSION(N), INTENT(OUT) :: RW
INTEGER, INTENT(OUT) :: IER
LOGICAL, OPTIONAL :: RLSHEP

! Local variables.
INTEGER :: FUNC  ! IRLS influence function switch, 1 = Huber, 2 = Tukey's
                 ! bisquare.
INTEGER :: IER1  ! error indicator for IRLS subroutine.
INTEGER :: INFO  ! error indicator for DGELSS.
INTEGER :: ITER  ! number of iterations the IRLS algorithm will perform.
INTEGER :: I, J, K  ! loop control variables.
INTEGER :: LWORK ! length of the work array WORK for DGELSS.
INTEGER :: NP    ! number of nodes used for the local fit.
INTEGER :: RANK  ! rank of the array (used for DGELSS).
INTEGER, DIMENSION(N-1) :: IDIST ! index array for DIST.
REAL(KIND=R8) :: DIAM  ! diameter of the whole date set.
REAL(KIND=R8) :: MAR   ! the median absolute residual.
REAL(KIND=R8) :: RCOND ! reciprocal condition number used in DGELSS.
REAL(KIND=R8) :: RP    ! least squares fit radius R_p.
REAL(KIND=R8), DIMENSION(MIN(N-1, (3*M+1)/2)) :: BETA_IN
               ! variable used in DGELSS.
REAL(KIND=R8), DIMENSION(N-1) :: DIST ! array that stores all the distances.
REAL(KIND=R8), DIMENSION(MIN(N-1, (3*M+1)/2)) :: OMEGA
               ! weights used for linear least squares fit.
REAL(KIND=R8), DIMENSION(MIN(N-1, (3*M+1)/2)) :: RES
               ! the residuals from the linear fit.
REAL(KIND=R8), DIMENSION(M) :: SIGMA   ! array used for DGELSS.
REAL(KIND=R8), DIMENSION(3*M+2*MIN(N-1, (3*M+1)/2)) :: WORK
               ! work array of length LWORK for DGELSS.
REAL(KIND=R8), DIMENSION(MIN(N-1, (3*M+1)/2), M) :: ALPHA, ALPHA_IN
REAL(KIND=R8), DIMENSION(MIN(N-1, (3*M+1)/2), 1) :: BETA
               ! variables used in DGELSS.

IER = 0
! Check the validity of M and N.
IF ( ( M < 1 ) .OR. ( N <= (M + 1) ) ) THEN
   IER = 1
   RETURN
END IF
! Set the number of nodes used in local linear fit.
NP = MIN( N - 1,  (3 * M + 1) / 2 )
! Set the RCOND parameter used in DGELSS.
RCOND = NP * EPSILON( 1.0_R8 )
! Set the LWORK parameter used in DGELSS.
LWORK = 3 * M + 2 * NP
! Calculate RW and A.
DIAM = 0.0_R8
DO K = 1, N
   J = 0
   DO I = 1, N
      IF (I /= K) THEN
         J = J + 1
         DIST(J) = DOT_PRODUCT( X(:, I) - X(:, K), X(:, I) - X(:, K) )
         IDIST(J) = I
      END IF
   END DO
   DIAM = MAX( DIAM, MAXVAL( DIST(1:N-1) ) )
   CALL SORT( DIST, IDIST, NP, N - 1 )
   DIST(1:NP) = SQRT( DIST(1:NP) )
   BETA_IN(1:NP) = F(IDIST(1:NP)) - F(K)
   DO I = 1, NP
      ALPHA_IN(I, 1:M) = X(1:M, IDIST(I)) - X(1:M, K)
   END DO
   RW(K) = DIST(NP)
   IF ( PRESENT(RLSHEP) ) THEN
      OMEGA(1:NP) = 1.0_R8
   ELSE
      RP = 1.1_R8 * RW(K)
      OMEGA(1:NP) =  (RP - DIST(1:NP)) / (RP * DIST(1:NP))
   END IF
   BETA(1:NP,1) = OMEGA(1:NP) * BETA_IN(1:NP)
   DO I = 1, NP
      ALPHA(I, 1:M) = OMEGA(I) * ALPHA_IN(I, 1:M)
   END DO
   CALL DGELSS( NP, M, 1, ALPHA(1:NP, 1:M), NP, BETA(1:NP,1), NP, SIGMA, &
        RCOND, RANK, WORK, LWORK, INFO )
   A(1:M, K) = BETA(1:M,1)
   IF ( RANK < M )  IER = 2
   ROBUST: IF ( PRESENT(RLSHEP) ) THEN
      RES(1:NP) = MATMUL( ALPHA_IN(1:NP, 1:M), BETA(1:M,1) ) - BETA_IN(1:NP)
      MAR = MAD( RES(1:NP), NP )
      IF ( MAR > EPSILON( 1.0_R8 ) ) THEN
         FUNC = 1
         ITER = 5
         CALL IRLS(M, N, MAR, ITER, FUNC, OMEGA, BETA, BETA_IN, ALPHA,&
              ALPHA_IN, RES, NP, IER1 )
         IF ( IER1 /= 0 ) THEN
            IER = 3
            CYCLE
         ELSE
            A(1:M, K) = BETA(1:M,1)
         END IF
         MAR = MAD( RES(1:NP), NP )
         IF ( MAR > EPSILON( 1.0_R8 ) ) THEN
            FUNC = 2
            ITER = 5
            CALL IRLS(M, N, MAR, ITER, FUNC, OMEGA, BETA, BETA_IN, &
                 ALPHA, ALPHA_IN, RES, NP, IER1 )
            IF ( IER1 /= 0 ) THEN
               IER = 3
               CYCLE
            ELSE
               A(1:M, K) = BETA(1:M,1)
            END IF
         END IF
      END IF
! Adjust RW(K) according to weights.
      DO I = 1, NP
         IF ( OMEGA(I) < 0.8_R8 ) THEN
            IF ( I == 1 ) THEN
               RW(K) = DIST(I) / 2.0_R8
            ELSE
               RW(K) = ( DIST(I - 1) + DIST(I) ) / 2.0_R8
            END IF
            EXIT
         END IF
      END DO
   END IF ROBUST
END DO
RW(1:N) = MIN( SQRT( DIAM ) / 2.0_R8, RW(1:N) )
RETURN
END SUBROUTINE LSHEP

SUBROUTINE IRLS( M, N,  S, ITER, FUNC, OMEGA, BETA, BETA_IN, ALPHA, &
           ALPHA_IN, RES, NP, IER )
! IRLS returns the coefficients of the linear least squares fit at the current
! node using robust M-estimation. For Tukey's bisquare influence function,
! there is an additional convergence test.
!
! Input parameters:
!  M    is the dimension of the problem.
!  N    is the number of points.
!  S    is the scale constant.
!  ITER is the number of iterations to be used in the algorithm.
!  FUNC identifies which IRLS influence function to use here, 1 = Huber,
!       2 = Tukey's bisquare.
!  BETA_IN holds distances.
!  ALPHA, ALPHA_IN are used for DGELSS calls.
!  RES is the residuals from the linear fit.
! Output parameters:
!  OMEGA stores the final weights.
!  BETA stores the coefficients of the linear fit using robust M-estimation.
!  NP is the size of OMEGA(:) and BETA(:,1).
!  IER
!   = 0, if no errors were encountered.
!   = 1, if any least squares problem is rank deficient.
!   = 2, if the bisquare IRLS did not converge.
!   = 3, if the input parameters are invalid.
USE REAL_PRECISION
IMPLICIT NONE
INTEGER, INTENT(IN) :: M, N
REAL(KIND=R8), INTENT(INOUT) :: S
INTEGER, INTENT(IN) :: ITER, FUNC
REAL(KIND=R8), INTENT(OUT), DIMENSION(:) :: OMEGA
REAL(KIND=R8), INTENT(OUT), DIMENSION(:,:) :: BETA
REAL(KIND=R8), INTENT(IN), DIMENSION(:) :: BETA_IN
REAL(KIND=R8), INTENT(OUT), DIMENSION(:,:) :: ALPHA
REAL(KIND=R8), INTENT(IN), DIMENSION(:,:) :: ALPHA_IN
REAL(KIND=R8), INTENT(INOUT), DIMENSION(:) :: RES
INTEGER, INTENT(IN) :: NP
INTEGER, INTENT(OUT) :: IER

! Local variables.
INTEGER :: I,J     ! loop control variables.
INTEGER:: INFO   ! error code from DGELSS.
INTEGER :: LWORK ! size of work array for DGELSS.
INTEGER :: RANK  ! rank of array (used for DGELSS).
REAL(KIND=R8) :: RCOND ! reciprocal condition number used in DGELSS.
REAL(KIND=R8), DIMENSION(M) :: SIGMA ! array used for DGELSS.
REAL(KIND=R8) :: ERR_NEW      ! objective function value after bisquare IRLS.
REAL(KIND=R8) :: ERR_OLD      ! objective function value before bisquare IRLS.
REAL(KIND=R8), DIMENSION(3*M+2*MIN(N-1, (5*M+1)/2)) :: WORK
REAL(KIND=R8), PARAMETER, DIMENSION(2) :: TUNING = (/1.0_R8, 3.0_R8/)
               ! tuning parameters for influence functions.
IF ((FUNC /= 1) .AND. (FUNC /= 2)) THEN
    IER=3
    RETURN
END IF
S = TUNING(FUNC) * S
! Set the RCOND parameter used in DGELSS.
RCOND = NP * EPSILON( 1.0_R8 )
! Set the LWORK parameter used in DGELSS.
LWORK = 3 * M + 2 * NP
IF ( FUNC ==  2 ) THEN
   ERR_OLD = SUM( 1.0_R8 - (MAX( 0.0_R8, 1.0_R8 - (RES(1:NP) /  S)**2 ))**3 )
END IF
DO I = 1, ITER
! Update the weight from the previous linear fit residual.
   SELECT CASE ( FUNC )
! Huber function
   CASE (1)
      DO J = 1, NP
         IF ( RES(J) == 0.0_R8 ) THEN
            OMEGA(J) = 1.0_R8
         ELSE
            OMEGA(J) = SQRT( MIN( 1.0_R8, S / ABS( RES(J) ) ) )
         END IF
      END DO
! Bisquare function
   CASE(2)
      OMEGA(1:NP) = MAX( 0.0_R8, 1.0_R8 - (RES(1:NP) / S)**2 )
   END SELECT
! Compute the linear fit with a new weight.
   BETA(1:NP,1) = OMEGA(1:NP) * BETA_IN(1:NP)
   DO J = 1, NP
      ALPHA(J, 1:M) = OMEGA(J) * ALPHA_IN(J, 1:M)
   END DO
   CALL DGELSS( NP, M, 1, ALPHA(1:NP,1:M), NP, BETA(1:NP,1), NP, SIGMA, &
        RCOND, RANK, WORK, LWORK, INFO )
   IF ( RANK < M ) THEN
      IER = 1
      RETURN
   END IF
   RES(1:NP) = MATMUL( ALPHA_IN(1:NP,1:M), BETA(1:M,1) ) - BETA_IN(1:NP)
END DO
IF ( FUNC ==  2 ) THEN
   ERR_NEW = SUM( 1.0_R8 - (MAX( 0.0_R8, 1.0_R8 - (RES(1:NP) / S)**2 ))**3 )
   IF ( ERR_NEW > ERR_OLD ) THEN
      IER = 2
      RETURN
   END IF
END IF
IER = 0
RETURN
END SUBROUTINE IRLS

FUNCTION MAD( RES, NP )
! The function MAD computes the median absolute deviation (MAD) scaled
! estimate for robust M-estimation from the residual array RES.
USE REAL_PRECISION
IMPLICIT NONE
REAL(KIND=R8), DIMENSION(:), INTENT(IN) :: RES
INTEGER, INTENT(IN)::NP
REAL(KIND=R8) :: MAD

! Local variables.
INTEGER :: NMID
REAL(KIND=R8) :: FACTOR, MED

FACTOR = 1.0_R8 / 0.6745_R8
NMID = NP / 2
IF ( MOD(NP, 2) == 0 ) THEN
   MED = (QUANTILE( RES(:), NMID ) + QUANTILE( RES(:), NMID + 1 )) / 2.0_R8
   MAD = FACTOR * (QUANTILE( ABS( RES(:) - MED ), NMID ) + &
        QUANTILE( ABS( RES(:) - MED ), NMID + 1 )) / 2.0_R8
ELSE
   MED = QUANTILE( RES(:), NMID + 1 )
   MAD = FACTOR * QUANTILE( ABS( RES(:) - MED ), NMID + 1 )
END IF
RETURN
END FUNCTION MAD

RECURSIVE FUNCTION QUANTILE( A, K ) RESULT( VALUE )
! Recursive function QUANTILE returns the K-th smallest element in array A.
USE REAL_PRECISION
IMPLICIT NONE
REAL(KIND=R8), DIMENSION (:), INTENT(IN) :: A
INTEGER, INTENT(IN) :: K
REAL(KIND=R8) :: VALUE

! Local variables.
INTEGER :: J
REAL(KIND=R8) :: AK

AK = A(K)
J = COUNT( A(:) < AK )
IF ( J >= K ) THEN
   VALUE = QUANTILE( PACK( A(:), A(:) < AK ), K )
ELSE
   J = COUNT( A(:) > AK ) + K - SIZE( A(:) )
   IF ( J > 0 ) THEN
      VALUE = QUANTILE( PACK( A(:), A(:) > AK ), J )
   ELSE
      VALUE = AK
   END IF
END IF
RETURN
END FUNCTION QUANTILE

SUBROUTINE SORT( DIST, IDIST, NUM, LENGTH )
! The subroutine SORT sorts the real array DIST of length LENGTH in ascending
! order for the smallest NUM elements. IDIST stores the original label for
! each element in array DIST.
! Local variables.
USE REAL_PRECISION
IMPLICIT NONE
REAL(KIND=R8), INTENT(INOUT), DIMENSION(:)::DIST
INTEGER, INTENT(INOUT), DIMENSION(:)::IDIST
INTEGER :: I, ITEMP, J, LENGTH, NUM
REAL(KIND=R8) :: TEMP

DO I = 2, LENGTH
   DO J = 1, MIN( I-1, NUM )
      IF ( DIST(I) < DIST(J) ) THEN
         TEMP = DIST(I)
         ITEMP = IDIST(I)
         DIST(J+1:I) = DIST(J:I-1)
         IDIST(J+1:I) = IDIST(J:I-1)
         DIST(J) = TEMP
         IDIST(J) = ITEMP
         EXIT
      END IF
   END DO
END DO
END SUBROUTINE SORT

FUNCTION LSHEPVAL( XP, M, N, X, F, A, RW, IER )
! LSHEPVAL returns the linear Shepard approximation at the point XP, using the
! local linear approximations computed by LSHEP.
!
! Input parameters:
!  XP is the point at which the linear Shepard interpolant function
!     approximation function is to be evaluated.
!  M is the dimension of the data.
!  N is the number of interpolation points.
!  X contains the interpolation nodes, by column.
!  F contains the function values of the interpolation nodes.
!  A is an M by N matrix containing the coefficients for linear nodal functions
!    returned by LSHEP.
!  RW contains the radius of influence about each interpolation node returned
!    by LSHEP.
! Output parameter:
!  IER
!   = 0, normal returns.
!   = 1, if the point XP is outside the radius of influence RW(i) for all
!        nodes, in which case LSHEPVAL is computed using the original Shepard
!        algorithm with the M+1 closest points.
!   = 2, if the hyperplane normals of the local approximations with positive
!        weights are significantly different. For a nonlinear underlying
!        function f(x), e.g., quadratic f(x), very different normals are
!        typical. For a piecewise linear underlying function f(x), IER = 2
!        signals a potentially large error in LSHEPVAL, since local
!        approximations from different facets of f(x) have been used.
USE REAL_PRECISION
IMPLICIT NONE

INTEGER, INTENT(IN)  :: M, N
REAL(KIND=R8), DIMENSION(M), INTENT(IN) :: XP
REAL(KIND=R8), DIMENSION(M, N), INTENT(IN) :: X
REAL(KIND=R8), DIMENSION(N), INTENT(IN) :: F
REAL(KIND=R8), DIMENSION(M, N), INTENT(IN) :: A
REAL(KIND=R8), DIMENSION(N), INTENT(IN) :: RW
INTEGER, INTENT(OUT):: IER
REAL(KIND=R8) :: LSHEPVAL

! Local variables.
INTEGER :: I   ! number of nodal functions used in Shepard approximation.
INTEGER :: J   ! temporary integer variable.
INTEGER :: K   ! loop control variable.
INTEGER, DIMENSION(M+1) :: ID ! data indices in the original Shepard scheme.
REAL(KIND=R8) :: DIST ! distance between XP and the interpolation nodes.
REAL(KIND=R8) :: TEMP ! temporary real variable.
REAL(KIND=R8) :: W    ! weight value, depends on D.
REAL(KIND=R8), DIMENSION(M+1) :: D ! weights in the original Shepard scheme.
REAL(KIND=R8), DIMENSION(M+1) :: SLOPE ! the hyperplane normal of the Shepard
               ! approximation function.
REAL(KIND=R8), DIMENSION(M+1, N) :: SSLOPE ! the hyperplane normals of the
               ! local approximations with positive weights.
REAL(KIND=R8), DIMENSION(N) :: SW   ! weight values.
REAL(KIND=R8), DIMENSION(N) :: SWC  ! nodal (linear) function values.

IER = 0
I = 0
J = 1
D(1:M+1) = 0.0_R8
TEMP = 0.0_R8
DO K = 1, N
   DIST = SQRT( DOT_PRODUCT( XP(:) - X(:, K), XP(:) - X(:, K) ) )
   IF ( DIST < RW(K) ) THEN
      IF ( RW(K) - DIST == RW(K) ) THEN
         LSHEPVAL = F(K)
         RETURN
      END IF
      I = I + 1
      W = ((RW(K) - DIST) / (RW(K) * DIST))**2
      SW(I) = W
      SWC(I) = DOT_PRODUCT( A(:, K), (XP(:) - X(:, K)) ) + F(K)
      SSLOPE(1:M, I) = A(1:M, K)
      SSLOPE(M+1, I) = -1.0_R8
   ELSE
      W = 1.0_R8 / DIST
      IF ( W > TEMP ) THEN
         D(J) = W
         ID(J) = K
         J = MINLOC( D(1:M+1), DIM = 1 )
         TEMP = D(J)
      END IF
   END IF
END DO
! I = 0 iff the point XP is not within the radius RW(K) of X(:,K) for all K;
! IER = 1.
IF ( I == 0 ) THEN
   D = D**2
   LSHEPVAL = DOT_PRODUCT( D / SUM( D ), F(ID) )
   IER = 1
   RETURN
ELSE
   SW(1:I) = SW(1:I) / SUM( SW(1:I) )
   LSHEPVAL = DOT_PRODUCT( SW(1:I), SWC(1:I) )
! Return IER = 2 iff the angle between two local approximation hyperplanes is
! too large.
   SLOPE(1:M+1) = MATMUL( SSLOPE(1:M+1, 1:I), SW(1:I) )
   SLOPE(:) = SLOPE(:) / SQRT( DOT_PRODUCT( SLOPE(:), SLOPE(:) ) )
   DO K = 1, I
      SSLOPE(:, K) = SSLOPE(:, K) / SQRT( DOT_PRODUCT( SSLOPE(:, K), &
           SSLOPE(:, K) ) )
      IF ( DOT_PRODUCT( SLOPE(1:M+1), SSLOPE(1:M+1, K) ) < 0.9_R8 ) THEN
         IER = 2
         EXIT
      END IF
   END DO
END IF
RETURN
END FUNCTION LSHEPVAL

SUBROUTINE RIPPLE( M, N, X, F, A, RW, IER )
! The subroutine RIPPLE, residual initiated polynomial-time piecewise linear
! estimation, computes a set of parameters defining a function that
! interpolates N data values F(i) at scattered nodes X(i) in M dimensions.
! The interpolant may be evaluated at an arbitrary point by the function
! LSHEPVAL.
! The interpolation scheme is a modified linear Shepard method.
! The time complexity of the algorithm RIPPLE is between the robust
! M-estimation and LMS estimation.
!
! Input parameters:
!   M is the dimension of the data.
!   N is the number of nodes.
!   X(M, N) contains the coordinates of the nodes.
!   F(N) contains the function values of the nodes.
! Output parameters:
!   A(M, N) contains the coefficients of the local fits.
!   RW(N) is an array containing the radius of influence for each point.
!   IER
!   = 0, if no errors were encountered.
!   = 1, if N is too small relative to M.
!   = 2, if the IRLS subroutine returns an error.
USE REAL_PRECISION
IMPLICIT NONE

INTEGER, INTENT(IN) :: M, N
REAL(KIND=R8), DIMENSION(M,N), INTENT(IN) :: X
REAL(KIND=R8), DIMENSION(N), INTENT(IN) :: F
REAL(KIND=R8), DIMENSION(M,N), INTENT(OUT) :: A
REAL(KIND=R8), DIMENSION(N), INTENT(OUT) :: RW
INTEGER, INTENT(OUT) :: IER

! Local variables
INTEGER :: FUNC ! IRLS influence function switch, 1 = Huber, 2 = Tukey's
                ! bisquare.
INTEGER :: IER1 ! error indicator for IRLS subroutine.
INTEGER :: INFO ! error indicator for DGELSS.
INTEGER :: ITER ! number of iterations for IRLS algorithm.
INTEGER :: I, J, K, L  ! loop control variables.
INTEGER :: LWORK ! length of the work array WORK for DGELSS.
INTEGER :: NP    ! number of nodes used for the local fit.
INTEGER :: RANK  ! rank of array (used for DGELSS).
INTEGER, DIMENSION(N-1) :: IDIST    ! index array for DIST.
INTEGER, DIMENSION(MIN(N-1, (3*M+1)/2), M+3) :: D
                    ! distance matrix of indices.
INTEGER, DIMENSION(N, MIN(N-1, (3*M+5)/2)) :: S
                    ! set of indices close to X(:, K).
REAL(KIND=R8) :: DIAM  ! diameter of the whole date set.
REAL(KIND=R8) :: LSE   ! the least squares error.
REAL(KIND=R8) :: MAR   ! the median absolute residual.
REAL(KIND=R8) :: RCOND ! reciprocal condition number used in DGELSS.
REAL(KIND=R8), DIMENSION(MIN(N-1, (5*M+1)/2)) :: BETA_IN
               ! variables used for DGELSS.
REAL(KIND=R8), DIMENSION(N-1) :: DIST  ! array that stores all distances.
REAL(KIND=R8), DIMENSION(MIN(N-1, (5*M+1)/2)) :: OMEGA
               ! weights used for least squares fit.
REAL(KIND=R8), DIMENSION(MIN(N-1, (5*M+1)/2)) :: RES
               ! the residuals from the linear fit.
REAL(KIND=R8), DIMENSION(M) :: SIGMA ! array used for DGELSS.
REAL(KIND=R8), DIMENSION(3*M+2*MIN(N-1, (5*M+1)/2)) :: WORK
               ! work array for DGELSS.
REAL(KIND=R8), DIMENSION(MIN(N-1, (5*M+1)/2), M):: ALPHA, ALPHA_IN
REAL(KIND=R8), DIMENSION(MIN(N-1, (5*M+1)/2), 1) :: BETA
               ! variables used for DGELSS.

IER = 0
! Check the validity of M and N.
IF ( ( M < 1 ) .OR. ( N <= (M + 3) ) ) THEN
   IER = 1
   RETURN
END IF
! Compute the set of indices S.
DIAM = 0.0_R8
NP = MIN( N - 1, (3 * M + 5) / 2 )
DO K = 1, N
   J = 0
   DO I = 1, N
      IF (I /= K) THEN
         J = J + 1
         DIST(J) = DOT_PRODUCT( X(:, I) - X(:, K), X(:, I) - X(:, K) )
         IDIST(J) = I
      END IF
   END DO
   DIAM = MAX( DIAM, MAXVAL( DIST(1:N-1) ) )
   CALL SORT( DIST, IDIST, NP, N - 1 )
   S(K, 1:NP) = IDIST(1:NP)
END DO
! Compute the local linear approximation.
DO K = 1, N
   NP = MIN( N - 1,  (3 * M + 1) / 2 )
   D(1:NP, 1) = S(K, 1:NP)
   DO J = 2, M + 3
      DO I = 1, NP
         DO L = 1, J
            IF ( ALL( S(D(I, J-1), L) /= D(I, 1:J-1) ) .AND. &
                 S(D(I, J-1), L) /= K )  D(I, J) = S(D(I, J-1), L)
         END DO
      END DO
   END DO
   LSE = HUGE(1.0_R8)
! Set the number of nodes used in searching for the set R.
   DO J = 1, NP
      BETA_IN(1:M+3) = F(D(J, 1:M+3)) - F(K)
      DO I = 1, M + 3
         ALPHA_IN(I, 1:M) = X(1:M, D(J, I)) - X(1:M, K)
      END DO
! Eliminate 2 of the M+3 points.
      RCOND = (M + 1) * EPSILON( 1.0_R8 )
      LWORK = 5 * M + 2
      DO I = 3, M + 3
         DO L = 2, I - 1
            ALPHA(1:L-1, 1:M) = ALPHA_IN(1:L-1, 1:M)
            ALPHA(L:I-2, 1:M) = ALPHA_IN(L+1:I-1, 1:M)
            ALPHA(I-1:M+1, 1:M) = ALPHA_IN(I+1:M+3, 1:M)
            BETA(1:L-1,1) = BETA_IN(1:L-1)
            BETA(L:I-2,1) = BETA_IN(L+1:I-1)
            BETA(I-1:M+1,1) = BETA_IN(I+1:M+3)
            CALL DGELSS( M+1,M,1,ALPHA(1:M+1,1:M),M+1,BETA(1:M+1,1),M+1, &
                 SIGMA, RCOND, RANK, WORK, LWORK, INFO )
            IF ( ABS( BETA(M+1,1) ) < LSE )  THEN
               LSE = ABS( BETA(M+1,1) )
               A(1:M, K) = BETA(1:M,1)
               IDIST(1:L-1) = D(J, 1:L-1)
               IDIST(L:I-2) = D(J, L+1:I-1)
               IDIST(I-1:M+1) = D(J, I+1:M+3)
            END IF
         END DO
      END DO
   END DO
! Compute the scale estimate and distance vector from the optimal set R.
   BETA_IN(1:M+1) = F(IDIST(1:M+1)) - F(K)
   DO I = 1, M + 1
      ALPHA_IN(I, 1:M) = X(1:M, IDIST(I)) - X(1:M, K)
      DIST(I) = DOT_PRODUCT( ALPHA_IN(I, 1:M), ALPHA_IN(I, 1:M) )
   END DO
   RES(1:M+1) = MATMUL( ALPHA_IN(1:M+1, 1:M), A(1:M, K) ) - BETA_IN(1:M+1)
   MAR = MAX( MAD( RES(1:M+1),NP ), 1.0E-6_R8 )
! Determine the set T for the final local approx.
   BETA_IN(1:NP) = F(S(K, 1:NP)) - F(K)
   DO I = 1, NP
      ALPHA_IN(I, 1:M) = X(1:M, S(K, I)) - X(1:M, K)
   END DO
   RW(K) = DOT_PRODUCT( ALPHA_IN(NP, :), ALPHA_IN(NP, :) )
   CALL SORT( DIST, IDIST, M + 1, M + 1 )
   DO I = 1, M + 1
      IF ( DIST(I) > RW(K) .OR. ( DIST(I) == RW(K) .AND. ALL( IDIST(I) &
           /= S(K, 1: MIN( N - 1, (3 * M + 1) / 2 )) ) ) ) THEN
         NP = NP + 1
         BETA_IN(NP) = F(IDIST(I)) - F(K)
         ALPHA_IN(NP, 1:M) = X(1:M, IDIST(I)) - X(1:M, K)
      END IF
   END DO
   RW(K) = SQRT( RW(K) )
   RES(1:NP) = MATMUL( ALPHA_IN(1:NP, 1:M), A(1:M, K) ) - BETA_IN(1:NP)
   FUNC = 1
   ITER = 5
   CALL IRLS(M, N, MAR, ITER, FUNC, OMEGA, BETA, BETA_IN, ALPHA, ALPHA_IN, &
        RES, NP, IER1 )
   IF ( IER1  /= 0 ) THEN
      IER = 2
      CYCLE
   ELSE
      A(1:M, K) = BETA(1:M,1)
   END IF
   FUNC = 2
   ITER = 5
   CALL IRLS(M, N, MAR, ITER, FUNC, OMEGA, BETA, BETA_IN, ALPHA, ALPHA_IN, &
        RES, NP, IER1 )
   IF ( IER1 /= 0 ) THEN
      IER = 2
      CYCLE
   ELSE
      A(1:M, K) = BETA(1:M,1)
   END IF
! Adjust RW(K) according to the weights.
   DO I = 1, NP
      IF ( OMEGA(I) < 0.8_R8 ) THEN
         IF ( I == 1 ) THEN
            RW(K) = SQRT( DOT_PRODUCT( ALPHA_IN(I, :), ALPHA_IN(I, :) ) ) &
                 / 2.0_R8
         ELSE
            RW(K) = (SQRT( DOT_PRODUCT( ALPHA_IN(I-1, :), ALPHA_IN(I-1, :) ) )&
                 + SQRT( DOT_PRODUCT( ALPHA_IN(I, :), ALPHA_IN(I, :) ) )) &
                 / 2.0_R8
         END IF
         EXIT
      END IF
   END DO
END DO
RW(1:N) = MIN( SQRT( DIAM ) / 2.0_R8, RW(1:N) )
RETURN
END SUBROUTINE RIPPLE

END MODULE LINEAR_SHEPARD_MOD

max_stars_repo_name: vtopt/vtmop

max_issues_repo_name: vtopt/vtmop

max_forks_repo_name: vtopt/vtmop


text:
Community Employment Services provides goaloriented support to adults with disabilities, including persons with autism, down syndrome, and physical disabilities. CES employees work with clients to help them form and attain goals in their personal and professional lives.  They have offices in Davis, Citrus Heights, Crescent City and Eureka.

By employing through CES, you can get tax credits.

20081205 19:23:48 nbsp I work here. This is not an employment agency. Please dont call asking for jobs ) Users/EricaMacGregor

20081205 19:24:37 nbsp Ill add to the page soon Users/EricaMacGregor

20090531 19:01:54 nbsp I work here too, it is a great place to work!  I hope that the States Budget Woes dont cause hurt us too much.   Users/HazelWatson

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

!>    @brief info, verbosity level
module mod_linfos

  !> @brief Array of verbosity level switches.
  !> @details
  !> Array of verbosity level switches. \n\n
  !> linfos(1)=\n
  !> - `0`, only necessary outputs\n
  !> - `1`, time steps information\n
  !> - `2`, with initialisation and file reading outputs\n
  !>
  !> linfos(2)=\n
  !> - `0`, 'inversion/simulation' with only necessary outputs\n
  !> - `1`, 'inversion/simulation' with some more outputs\n
  !> - `2`, 'inversion/simulation' with all outputs (optimisation outputs)\n
  !> - `3`, 'inversion/simulation' with all debug outputs\n
  !>
  !> linfos(3)=\n
  !> - `-1`, 'non linear solver' disables some error outputs\n
  !> - `0`, 'non linear solver' with only necessary outputs\n
  !> - `1`, 'non linear solver' with some stage outputs\n
  !> - `2`, 'non linear solver' with all outputs\n
  !>
  !> linfos(4)=\n
  !> - `0`, 'system solver' without any outputs\n
  !> - `1`, 'system solver' with some outputs\n
  !> - `2`, 'system solver' with some debug outputs\n
  !> - `3`, 'system solver' with all debug outputs\n
  integer, dimension (4) :: linfos

end module mod_linfos


! =================================================================================================
! This file is part of the CodeVault project. The project is licensed under Apache Version 2.0.
! CodeVault is part of the EU-project PRACE-4IP (WP7.3.C).
!
! Author(s):
!   Mariusz Uchronski <mariusz.uchronski@pwr.edu.pl>
!
! This example demonstrates the use of Fortran 90 with OpenMP for matrix-matrix multiplication.
! The example is set-up to perform single precision matrix-matrix multiplication.
!
! See [http://openmp.org] for the full OpenMP documentation.
!
! =================================================================================================

program gemm_openmp

  implicit none
  integer,parameter :: seed = 86456
  integer, parameter :: n = 400
  integer, parameter :: m = 600
  integer, parameter :: k = 800
  real, allocatable :: a(:)
  real, allocatable :: b(:)
  real, allocatable :: c(:)

  allocate(a(n*m))
  allocate(b(m*k))
  allocate(c(n*k))

  call fill_random(a, n, m)
  call fill_random(b, m, k)

  call run_gemm_openmp(a, b, c, n, m, k)

  deallocate(a)
  deallocate(b)
  deallocate(c)

end program gemm_openmp

subroutine fill_random(a, n, m)

  implicit none

  integer, intent(in) :: n, m
  real, dimension(n*m), intent(inout) :: a
  integer :: i, j
  do i=1, n
    do j=1, m
      a((i-1)*m+j) = rand()
    end do
  end do

end subroutine fill_random

subroutine print_array(a, n, m)

  implicit none

  integer, intent(in) :: n, m
  real, dimension(n*m), intent(in) :: a
  integer :: i, j
  do i=1, n
    do j=1, m
      print '(F8.4,$)', a((i-1)*m+j)
    end do
    print *, ''
  end do

end subroutine print_array

subroutine run_gemm_openmp(a, b, c, n, m, k)

  integer, intent(in) :: n, m , k
  real, dimension(n*m), intent(in) :: a
  real, dimension(m*k), intent(in) :: b
  real, dimension(n*k), intent(inout) :: c
  integer :: i, j, l
  real :: summ

!$OMP PARALLEL DO DEFAULT (NONE) &
!$OMP SHARED(n,m,k,a,b,c) PRIVATE(i,j,l,summ)
  do i=1, n
    do j=1, k
      summ = 0
      do l=1, m
        summ = summ + a((i-1)*m+l) * b((l-1)*k+j)
      end do
      c((i-1)*k+j) = summ
    end do
  end do
!$OMP END PARALLEL DO

end subroutine


text:
!==== This file contains the box coefficients for the Q(i2) pieces of the following amplitude

!===== (i1)^- (i2)^+ (i3)^- + (i4)^+ + gamma(i5) + gluon(i6)

! where i3 and i4 are the leptons from a W decay.

!=== C. Williams Nov 2013

!================================================================================================

!======= Q2 h5=1 h6=1 LEADING COLOR AMPLITUDES (order matches KC)
!=====  s(34)

      double complex function Bcoeff_q2lcpp_s34(i1,i2,i3,i4,i5,i6,za,zb)
      implicit none
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      integer i1,i2,i3,i4,i5,i6
      double complex zab2,t,ans
!===== statement functions
      zab2(i1,i2,i3,i4)=za(i1,i2)*zb(i2,i4)+za(i1,i3)*zb(i3,i4)
      t(i1,i2,i3)=s(i1,i2)+s(i2,i3)+s(i1,i3)

      Bcoeff_q2lcpp_s34=(-za(i1,i2)*zb(i4,i5)*
     -    (za(i1,i3)/(-s(i3,i4) + t(i1,i2,i6)) +
     -      (za(i1,i5)*za(i4,i3)*zb(i5,i4))/
     -       (2.*(-s(i3,i4) + t(i1,i2,i6))**2)))/
     -  (za(i1,i6)*za(i2,i5)*za(i2,i6))

      return
      end

      double complex function Bcoeff_q2lcpp_s16(i1,i2,i3,i4,i5,i6,za,zb)
      implicit none
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      integer i1,i2,i3,i4,i5,i6
      double complex zab2,t,ans
!===== statement functions
      zab2(i1,i2,i3,i4)=za(i1,i2)*zb(i2,i4)+za(i1,i3)*zb(i3,i4)
      t(i1,i2,i3)=s(i1,i2)+s(i2,i3)+s(i1,i3)

      Bcoeff_q2lcpp_s16=   (-za(i3,i2)*zab2(i2,i3,i4,i5)*zb(i2,i6)*
     -    ((2*za(i1,i3))/(-s(i1,i6) + t(i1,i2,i6)) -
     -      (za(i1,i6)*za(i2,i3)*zb(i6,i2))/
     -       (2.*(-s(i1,i6) + t(i1,i2,i6))**2)))/
     -  (za(i2,i5)*za(i2,i6)*za(i3,i4)*zab2(i5,i3,i4,i5))
      return
      end

!============ SUBLEADING COLORS + +

      double complex function Bcoeff_q2slcpp_s34
     &     (i1,i2,i3,i4,i5,i6,za,zb)
      implicit none
            include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      integer i1,i2,i3,i4,i5,i6
      double complex zab2,t,ans
!===== statement functions
      zab2(i1,i2,i3,i4)=za(i1,i2)*zb(i2,i4)+za(i1,i3)*zb(i3,i4)
      t(i1,i2,i3)=s(i1,i2)+s(i2,i3)+s(i1,i3)

      Bcoeff_q2slcpp_s34=((za(i1,i3)/zab2(i5,i3,i4,i5) -
     -       (za(i3,i4)*za(i5,i1)*zb(i4,i5))/
     -        (2.*zab2(i5,i3,i4,i5)**2))*zb(i5,i4))/
     -   (za(i2,i6)*za(i5,i6)) +
     -  (za(i1,i5)**2*(za(i6,i3)/zab2(i5,i3,i4,i5) -
     -       (za(i3,i4)*za(i5,i6)*zb(i4,i5))/
     -        (2.*zab2(i5,i3,i4,i5)**2))*zb(i5,i4))/
     -   (za(i1,i6)*za(i2,i5)*za(i5,i6)**2) +
     -  (((za(i1,i5)*za(i3,i6))/
     -        (za(i5,i6)*zab2(i6,i3,i4,i6)) +
     -       (za(i3,i4)*za(i6,i1)*zb(i4,i6))/
     -        zab2(i6,i3,i4,i6)**2)*zb(i6,i4))/
     -   (za(i2,i5)*za(i5,i6))

      return
      end

      double complex function Bcoeff_q2slcpp_s16
     &     (i1,i2,i3,i4,i5,i6,za,zb)
      implicit none
            include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      integer i1,i2,i3,i4,i5,i6
      double complex zab2,t,ans
!===== statement functions
      zab2(i1,i2,i3,i4)=za(i1,i2)*zb(i2,i4)+za(i1,i3)*zb(i3,i4)
      t(i1,i2,i3)=s(i1,i2)+s(i2,i3)+s(i1,i3)

       Bcoeff_q2slcpp_s16=  (za(i3,i2)*(-(s(i2,i6)*za(i1,i6)*za(i2,i3))/
     -       (2.*zab2(i2,i1,i6,i2)**2) +
     -      (za(i1,i3)*za(i2,i6) + za(i1,i2)*za(i3,i6))/
     -       zab2(i2,i1,i6,i2))*zab2(i2,i3,i4,i5)*zb(i2,i6))/
     -  (za(i2,i5)*za(i2,i6)**2*za(i3,i4)*zab2(i5,i3,i4,i5))

      return
      end

      double complex function Bcoeff_q2slcpp_s125
     &     (i1,i2,i3,i4,i5,i6,za,zb)
      implicit none
            include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      integer i1,i2,i3,i4,i5,i6
      double complex zab2,t,ans
!===== statement functions
      zab2(i1,i2,i3,i4)=za(i1,i2)*zb(i2,i4)+za(i1,i3)*zb(i3,i4)
      t(i1,i2,i3)=s(i1,i2)+s(i2,i3)+s(i1,i3)

      Bcoeff_q2slcpp_s125= -((za(i1,i3)*za(i3,i5))/
     -     (za(i2,i5)*za(i3,i4)*za(i5,i6)**2)) +
     -  (za(i1,i2)*za(i3,i5)**2*zb(i5,i2))/
     -   (za(i2,i5)*za(i3,i4)*za(i5,i6)**2*zab2(i5,i1,i2,i5))
     -   - (((za(i1,i5)*za(i3,i6))/
     -        (za(i5,i6)*zab2(i6,i3,i4,i6)) +
     -       (za(i3,i4)*za(i6,i1)*zb(i4,i6))/
     -        zab2(i6,i3,i4,i6)**2)*zb(i6,i4))/
     -   (za(i2,i5)*za(i5,i6))

      return
      end

      double complex function Bcoeff_q2slcpp_s12
     &     (i1,i2,i3,i4,i5,i6,za,zb)
      implicit none
            include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      integer i1,i2,i3,i4,i5,i6
      double complex zab2,t,ans
!===== statement functions
      zab2(i1,i2,i3,i4)=za(i1,i2)*zb(i2,i4)+za(i1,i3)*zb(i3,i4)
      t(i1,i2,i3)=s(i1,i2)+s(i2,i3)+s(i1,i3)

        Bcoeff_q2slcpp_s12=-((za(i1,i2)*za(i3,i5)**2*zb(i5,i2))/
     -     (za(i2,i5)*za(i3,i4)*za(i5,i6)**2*
     -       zab2(i5,i1,i2,i5))) +
     -  (za(i2,i1)*za(i3,i6)*zb(i6,i2)*
     -     (((za(i2,i3) - (za(i2,i6)*za(i3,i5))/za(i5,i6))*
     -          zab2(i6,i3,i4,i5))/
     -        (za(i2,i6)*zab2(i6,i1,i2,i6)) +
     -       (t(i3,i4,i5)*za(i3,i6)*zb(i6,i5))/
     -        zab2(i6,i1,i2,i6)**2))/
     -   (za(i2,i6)*za(i3,i4)*za(i5,i6)*zab2(i5,i3,i4,i5))

        return
        end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
module writer_mod
! module provides output routines in raw binary and .nrrd formats
!
!
    use constants, only : wp
    use utils,     only : str

implicit none

    interface nrrd_write
        module procedure write_3d_r8_nrrd
    end interface nrrd_write

    interface raw_write
        module procedure write_3d_r8_raw
    end interface raw_write

    private
    public :: raw_write, nrrd_write, normalise_fluence, write

    contains

        function normalise_fluence(grid, array, nphotons) result(out)
        ! normalise fluence in the Lucy 1999 way

            use gridMod

            implicit none

            type(cart_grid), intent(IN) :: grid
            real(kind=wp),   intent(IN) :: array(:, :, :)
            integer,         intent(IN) :: nphotons

            real(kind=wp), allocatable :: out(:, :, :)

            real(kind=wp) :: xmax, ymax, zmax
            integer       :: nxg, nyg, nzg

            nxg = grid%nxg
            nyg = grid%nyg
            nzg = grid%nzg
            xmax = grid%xmax
            ymax = grid%ymax
            zmax = grid%zmax

            allocate(out(size(array, 1), size(array, 2), size(array, 3)))

            out  = array * ((2._wp*xmax*2._wp*ymax)/(nphotons * (2._wp * xmax / nxg) * (2._wp * ymax / nyg) * (2._wp * zmax / nzg)))

        end function normalise_fluence

        subroutine write(array, filename, dict)
        ! routine automatically selects which way to write ouresults based upon file extension
            use fhash,        only : fhash_tbl_t

            implicit none

            real(kind=wp),          intent(IN) :: array(:,:,:)
            character(*),           intent(IN) :: filename
            type(fhash_tbl_t), optional, intent(INOUT) :: dict

            integer :: pos

            pos = index(filename, ".nrrd")
            if(pos > 0)then
                if(present(dict))then
                    call nrrd_write(array, filename, dict)
                else
                    call nrrd_write(array, filename)
                end if
                return
            end if

            pos = index(filename, ".raw")
            if(pos > 0)then
                call raw_write(array, filename)
                return
            end if

            pos = index(filename, ".dat")
            if(pos > 0)then
                call raw_write(array, filename)
                return
            end if

            error stop "File type not supported!"

        end subroutine write

        subroutine write_3d_r8_raw(array, filename)

            implicit none

            real(kind=wp), intent(IN) :: array(:, :, :)
            character(*),  intent(IN) :: filename

            integer :: u
            character(len=:), allocatable :: file

            if(check_file(filename))then
                file = get_new_file_name(filename)
            else
                file = filename
            end if
            open(newunit=u,file=file,access='stream',status='REPLACE',form='unformatted')
            write(u) array
            close(u)

        end subroutine write_3d_r8_raw

        function get_new_file_name(file) result(res)

            implicit none

            character(len=*), intent(IN) :: file
            character(len=:), allocatable :: res
            integer :: pos, i

            i = 1
            do
                pos = scan(trim(file), ".", back=.true.)
                res = file(1:pos-1)//" ("//str(i)//")"//file(pos:)
                if(.not. check_file(res))exit
                i = i + 1
            end do

        end function get_new_file_name

        logical function check_file(file) result(res)

            implicit none

            character(len=*), intent(IN) :: file

            inquire(file=trim(file), exist=res)

        end function check_file

        subroutine write_hdr(u, sizes, type)

            implicit none

            character(*), intent(IN) :: type
            integer,      intent(IN) :: sizes(:), u

            character(len=100) :: string
            integer :: i

            string = ""
            do i = 1, size(sizes)
                if(i == 1)then
                    string = str(sizes(i))
                else
                    string = trim(string) // " " // str(sizes(i))
                end if
            end do

            write(u,"(A)")"NRRD0004"
            write(u,"(A)")"type: "//type
            write(u,"(A)")"dimension: "//str(size(sizes))
            write(u,"(A)")"sizes: "//trim(string)
            write(u,"(A)")"space dimension: "//str(size(sizes))
            write(u,"(A)")"encoding: raw"
            write(u,"(A)")"endian: little"

        end subroutine write_hdr

        subroutine write_3d_r8_nrrd(array, filename, dict)

            use fhash,           only : fhash_tbl_t, key=>fhash_tbl_t, fhash_key_t
            use iso_fortran_env, only : int32, int64, real32, real64

            implicit none

            character(*),                intent(IN)    :: filename
            real(kind=wp),               intent(IN)    :: array(:, :, :)
            type(fhash_tbl_t), optional, intent(INOUT) :: dict

            class(fhash_key_t), pointer   :: keyd
            character(len=:), allocatable :: key_out, file, type_str
            integer :: u, i, stat

            integer(kind=int32) :: val_i32
            integer(kind=int64) :: val_i64
            real(kind=real32)   :: val_r32
            real(kind=real64)   :: val_r64
            character(len=:), allocatable :: val_char
            logical :: val_bool

            if(check_file(filename))then
                file = get_new_file_name(filename)
            else
                file = filename
            end if

            open(newunit=u,file=file,form="formatted")
            !to do fix precision
            call write_hdr(u, [size(array, 1), size(array, 2), size(array, 3)], "double")

            if(present(dict))then
                do
                    call dict%next(keyd, key_out, type_str, stat)
                    if(stat < 0)exit
                    select case(type_str)
                        case("integer32")
                            call dict%get(keyd, val_i32)
                            write(u, "(A)")trim(key_out)//": "//str(val_i32)
                        case("integer64")
                            call dict%get(keyd, val_i64)
                            write(u, "(A)")trim(key_out)//": "//str(val_i64)
                        case("real32")
                            call dict%get(keyd, val_r32)
                            write(u, "(A)")trim(key_out)//": "//str(val_r32)
                        case("real64")
                            call dict%get(keyd, val_r64)
                            write(u, "(A)")trim(key_out)//": "//str(val_r64)
                        case("character*")
                            call dict%get(keyd, val_char)
                            write(u, "(A)")trim(key_out)//": "//val_char
                        case("logical")
                            call dict%get(keyd, val_bool)
                            write(u, "(A)")trim(key_out)//": "//str(val_bool)
                        case default
                            continue
                    end select
                end do
            end if
            write(u,"(A)")new_line("C")
            close(u)
            open(newunit=u,file=file,access="stream",form="unformatted",position="append")
            write(u)array
            close(u)

        end subroutine write_3d_r8_nrrd
end module writer_mod

max_stars_repo_name: lewisfish/signedMCRT

max_issues_repo_name: lewisfish/signedMCRT

max_forks_repo_name: lewisfish/signedMCRT