text:
!===============================================================================
!===============================================================================
MODULE opacity_mod
  USE io_mod
  implicit none
  private
  type, public:: opacity_t
  contains
    procedure:: init
  end type
  type (opacity_t), public:: opacity
CONTAINS

!===============================================================================
!> Read in the data from a specific table file
!===============================================================================
SUBROUTINE init (self)
  class (opacity_t):: self
END SUBROUTINE init

END MODULE opacity_mod

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch


text:
      recursive
     *subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin)
      integer n,ldr
      double precision alpha
      double precision r(ldr,n),w(n),b(n),cos(n),sin(n)
c     **********
c
c     subroutine rwupdt
c
c     given an n by n upper triangular matrix r, this subroutine
c     computes the qr decomposition of the matrix formed when a row
c     is added to r. if the row is specified by the vector w, then
c     rwupdt determines an orthogonal matrix q such that when the
c     n+1 by n matrix composed of r augmented by w is premultiplied
c     by (q transpose), the resulting matrix is upper trapezoidal.
c     the matrix (q transpose) is the product of n transformations
c
c           g(n)*g(n-1)* ... *g(1)
c
c     where g(i) is a givens rotation in the (i,n+1) plane which
c     eliminates elements in the (n+1)-st plane. rwupdt also
c     computes the product (q transpose)*c where c is the
c     (n+1)-vector (b,alpha). q itself is not accumulated, rather
c     the information to recover the g rotations is supplied.
c
c     the subroutine statement is
c
c       subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin)
c
c     where
c
c       n is a positive integer input variable set to the order of r.
c
c       r is an n by n array. on input the upper triangular part of
c         r must contain the matrix to be updated. on output r
c         contains the updated triangular matrix.
c
c       ldr is a positive integer input variable not less than n
c         which specifies the leading dimension of the array r.
c
c       w is an input array of length n which must contain the row
c         vector to be added to r.
c
c       b is an array of length n. on input b must contain the
c         first n elements of the vector c. on output b contains
c         the first n elements of the vector (q transpose)*c.
c
c       alpha is a variable. on input alpha must contain the
c         (n+1)-st element of the vector c. on output alpha contains
c         the (n+1)-st element of the vector (q transpose)*c.
c
c       cos is an output array of length n which contains the
c         cosines of the transforming givens rotations.
c
c       sin is an output array of length n which contains the
c         sines of the transforming givens rotations.
c
c     subprograms called
c
c       fortran-supplied ... dabs,dsqrt
c
c     argonne national laboratory. minpack project. march 1980.
c     burton s. garbow, dudley v. goetschel, kenneth e. hillstrom,
c     jorge j. more
c
c     **********
      integer i,j,jm1
      double precision cotan,one,p5,p25,rowj,tan,temp,zero
      data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/
c
      do 60 j = 1, n
         rowj = w(j)
         jm1 = j - 1
c
c        apply the previous transformations to
c        r(i,j), i=1,2,...,j-1, and to w(j).
c
         if (jm1 .lt. 1) go to 20
         do 10 i = 1, jm1
            temp = cos(i)*r(i,j) + sin(i)*rowj
            rowj = -sin(i)*r(i,j) + cos(i)*rowj
            r(i,j) = temp
   10       continue
   20    continue
c
c        determine a givens rotation which eliminates w(j).
c
         cos(j) = one
         sin(j) = zero
         if (rowj .eq. zero) go to 50
         if (dabs(r(j,j)) .ge. dabs(rowj)) go to 30
            cotan = r(j,j)/rowj
            sin(j) = p5/dsqrt(p25+p25*cotan**2)
            cos(j) = sin(j)*cotan
            go to 40
   30    continue
            tan = rowj/r(j,j)
            cos(j) = p5/dsqrt(p25+p25*tan**2)
            sin(j) = cos(j)*tan
   40    continue
c
c        apply the current transformation to r(j,j), b(j), and alpha.
c
         r(j,j) = cos(j)*r(j,j) + sin(j)*rowj
         temp = cos(j)*b(j) + sin(j)*alpha
         alpha = -sin(j)*b(j) + cos(j)*alpha
         b(j) = temp
   50    continue
   60    continue
      return
c
c     last card of subroutine rwupdt.
c
      end

ext: f

max_stars_repo_name: Ennosigaeon/scipy

max_issues_repo_name: Ennosigaeon/scipy

max_forks_repo_name: Ennosigaeon/scipy


text:
c     Calculate LU decomposition of a matrix
c     Author: Pearu Peterson, March 2002
c
c     prefixes: d,z,s,c   (double,complex double,float,complex float)
c     suffixes: _c,_r     (column major order,row major order)

      subroutine dlu_c(p,l,u,a,m,n,k,piv,info,permute_l,m1)
      integer m,n,piv(k),i,j,k,permute_l,m1,info
      double precision l(m,k),u(k,n),a(m,n)
      double precision p(m1,m1)

cf2py intent(in,copy) :: a
cf2py intent(out) :: info
cf2py integer intent(hide,cache),depend(k),dimension(k) :: piv
cf2py integer intent(hide),depend(a) :: m = shape(a,0)
cf2py integer intent(hide),depend(a) :: n = shape(a,1)
cf2py integer intent(hide),depend(m,n) :: k = (m<n?m:n)
cf2py intent(out) :: p,l,u
cf2py integer optional,intent(in):: permute_l = 0
cf2py integer intent(hide),depend(permute_l,m) :: m1 = (permute_l?1:m)
cf2py depend(m1) :: p

cf2py callprotoargument double*,double*,double*,double*,int*,int*,int*,int*,int*,int*,int*

      external dgetrf,dlaswp
      call dgetrf(m,n,a,m,piv,info)
      if (info.lt.0) then
         return
      endif
      do 20 j=1,k
         l(j,j) = 1d0
         do 10, i=j+1,m
            l(i,j) = a(i,j)
 10      continue
 20   continue
      do 40 j=1,k
         do 30, i=1,j
            u(i,j) = a(i,j)
 30      continue
 40   continue
      do 60 j=k+1,n
         do 50, i=1,k
            u(i,j) = a(i,j)
 50      continue
 60   continue
      if (permute_l.ne.0) then
         call dlaswp(k,l,m,1,k,piv,-1)
      else
         do 65 i=1,m
            p(i,i)=1d0
 65       continue
         call dlaswp(m,p,m,1,k,piv,-1)
      endif
      end

      subroutine zlu_c(p,l,u,a,m,n,k,piv,info,permute_l,m1)
      integer m,n,piv(k),i,j,k,permute_l,m1,info
      complex*16 l(m,k),u(k,n),a(m,n)
      double precision p(m1,m1)

cf2py intent(in,copy) :: a
cf2py intent(out) :: info
cf2py integer intent(hide,cache),depend(k),dimension(k) :: piv
cf2py integer intent(hide),depend(a) :: m = shape(a,0)
cf2py integer intent(hide),depend(a) :: n = shape(a,1)
cf2py integer intent(hide),depend(m,n) :: k = (m<n?m:n)
cf2py intent(out) :: p,l,u
cf2py integer optional,intent(in):: permute_l = 0
cf2py integer intent(hide),depend(permute_l,m) :: m1 = (permute_l?1:m)
cf2py depend(m1) :: p

cf2py callprotoargument double*,complex_double*,complex_double*,complex_double*,int*,int*,int*,int*,int*,int*,int*

      external zgetrf,zlaswp,dlaswp
      call zgetrf(m,n,a,m,piv,info)
      if (info.lt.0) then
         return
      endif
      do 20 j=1,k
         l(j,j) = 1d0
         do 10, i=j+1,m
            l(i,j) = a(i,j)
 10      continue
 20   continue
      do 40 j=1,k
         do 30, i=1,j
            u(i,j) = a(i,j)
 30      continue
 40   continue
      do 60 j=k+1,n
         do 50, i=1,k
            u(i,j) = a(i,j)
 50      continue
 60   continue
      if (permute_l.ne.0) then
         call zlaswp(k,l,m,1,k,piv,-1)
      else
         do 65 i=1,m
            p(i,i)=1d0
 65       continue
         call dlaswp(m,p,m,1,k,piv,-1)
      endif
      end

      subroutine slu_c(p,l,u,a,m,n,k,piv,info,permute_l,m1)
      integer m,n,piv(k),i,j,k,permute_l,m1,info
      real l(m,k),u(k,n),a(m,n)
      real p(m1,m1)

cf2py intent(in,copy) :: a
cf2py intent(out) :: info
cf2py integer intent(hide,cache),depend(k),dimension(k) :: piv
cf2py integer intent(hide),depend(a) :: m = shape(a,0)
cf2py integer intent(hide),depend(a) :: n = shape(a,1)
cf2py integer intent(hide),depend(m,n) :: k = (m<n?m:n)
cf2py intent(out) :: p,l,u
cf2py integer optional,intent(in):: permute_l = 0
cf2py integer intent(hide),depend(permute_l,m) :: m1 = (permute_l?1:m)
cf2py depend(m1) :: p

cf2py callprotoargument float*,float*,float*,float*,int*,int*,int*,int*,int*,int*,int*

      external sgetrf,slaswp
      call sgetrf(m,n,a,m,piv,info)
      if (info.lt.0) then
         return
      endif
      do 20 j=1,k
         l(j,j) = 1d0
         do 10, i=j+1,m
            l(i,j) = a(i,j)
 10      continue
 20   continue
      do 40 j=1,k
         do 30, i=1,j
            u(i,j) = a(i,j)
 30      continue
 40   continue
      do 60 j=k+1,n
         do 50, i=1,k
            u(i,j) = a(i,j)
 50      continue
 60   continue
      if (permute_l.ne.0) then
         call slaswp(k,l,m,1,k,piv,-1)
      else
         do 65 i=1,m
            p(i,i)=1e0
 65       continue
         call slaswp(m,p,m,1,k,piv,-1)
      endif
      end

      subroutine clu_c(p,l,u,a,m,n,k,piv,info,permute_l,m1)
      integer m,n,piv(k),i,j,k,permute_l,m1,info
      complex l(m,k),u(k,n),a(m,n)
      real p(m1,m1)

cf2py intent(in,copy) :: a
cf2py intent(out) :: info
cf2py integer intent(hide,cache),depend(k),dimension(k) :: piv
cf2py integer intent(hide),depend(a) :: m = shape(a,0)
cf2py integer intent(hide),depend(a) :: n = shape(a,1)
cf2py integer intent(hide),depend(m,n) :: k = (m<n?m:n)
cf2py intent(out) :: p,l,u
cf2py integer optional,intent(in):: permute_l = 0
cf2py integer intent(hide),depend(permute_l,m) :: m1 = (permute_l?1:m)
cf2py depend(m1) :: p

cf2py callprotoargument float*,complex_float*,complex_float*,complex_float*,int*,int*,int*,int*,int*,int*,int*

      external cgetrf,claswp,slaswp
      call cgetrf(m,n,a,m,piv,info)
      if (info.lt.0) then
         return
      endif
      do 20 j=1,k
         l(j,j) = 1d0
         do 10, i=j+1,m
            l(i,j) = a(i,j)
 10      continue
 20   continue
      do 40 j=1,k
         do 30, i=1,j
            u(i,j) = a(i,j)
 30      continue
 40   continue
      do 60 j=k+1,n
         do 50, i=1,k
            u(i,j) = a(i,j)
 50      continue
 60   continue
      if (permute_l.ne.0) then
         call claswp(k,l,m,1,k,piv,-1)
      else
         do 65 i=1,m
            p(i,i)=1e0
 65       continue
         call slaswp(m,p,m,1,k,piv,-1)
      endif
      end


text:
	SUBROUTINE RU_MMXW  ( report, lenr, wnknot, above, ipoint, data,
     +			      nlev, iret )
C************************************************************************
C* RU_MMXW								*
C*									*
C* This subroutine decodes data for the max wind level(s) from the TTAA *
C* reports.  The output data are ordered PRES DRCT SPED.		*
C*									*
C* RU_MMXW  ( REPORT, LENR, WNKNOT, ABOVE, IPOINT, DATA, NLEV, IRET )	*
C*									*
C* Input parameters:							*
C*	REPORT		CHAR*		Station report			*
C*	LENR		INTEGER		Length of report		*
C*	WNKNOT		LOGICAL		Flag for speed in knots		*
C*	ABOVE		LOGICAL		Above 100 mb flag               *
C*									*
C* Input and output parameters:						*
C*	IPOINT		INTEGER		Pointer within report		*
C*									*
C* Output parameters:							*
C*	DATA (3,NLEV)	REAL		Max wind station data		*
C*	NLEV		INTEGER		Number of levels		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* D. Kidwell/NCEP	 2/01						*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'rucmn.cmn'
C*
	CHARACTER*(*)	report
	LOGICAL		wnknot, above
	REAL		data (3,*)
C*
	CHARACTER	pp*2, field*10
	LOGICAL		done
C*
	LOGICAL		ENDRPT
C*
	INCLUDE		'ERMISS.FNC'
C
C*	Function to check for end of message.
C
	ENDRPT ( field, lenf, ier ) =
     +				( ( lenf .gt. MAXFLN )  .or.
     +				  ( ier  .ne. 0 )  .or.
     +				  ( field (1:5) .eq. '51515' ) )
C------------------------------------------------------------------------
	iret = 0
	nlev = 0
	done = .false.
C
C*	Read the max wind data level by level.
C
	DO WHILE ( .not. done )
C
C*	    Get the first field which contains the pressure.
C
	    CALL RU_GFLD  ( report, lenr, ipoint, field, lenf, ier )
C
C*	    Check for the end of message.
C
	    IF  ( ENDRPT ( field, lenf, ier ) )  THEN
		ipoint = lenr + 1
		RETURN
	      ELSE
		pp = field ( 1:2 )
	    END IF
C
C*	    Check if this is a max wind level and for a pressure value.
C
	    IF ( pp .eq. '77' .or. pp .eq. '66' ) THEN
      	        IF ( field ( 3:lenf ) .ne. '999' ) THEN
C
C*		    This was the expected level.  Parse it.
C
		    CALL ST_CRNM ( field ( 3:lenf ), pres, iret )
     		    IF ( .not. ERMISS ( pres ) ) THEN
		        CALL RU_GFLD  ( report, lenr, ipoint, field,
     +				        lenf, ier )
		        CALL RU_WIND  ( field, wnknot, dir, spd, ier )
C
C*		        Check for missing speed and direction.
C
		        IF ( ERMISS ( dir ) .and. ERMISS ( spd ) ) THEN
C
C*			    Do nothing.
C
		          ELSE
		            nlev = nlev + 1
			    IF ( above ) pres = pres / 10.
			    data ( 1, nlev ) = pres
			    data ( 2, nlev ) = dir
			    data ( 3, nlev ) = spd
			    IF ( nlev .eq. LLMXLV ) done = .true.
		        END IF
		    END IF
		  ELSE IF ( pp .eq. '77' ) THEN
		    RETURN
		END IF
	    END IF
	END DO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
! ###################################################################
! Copyright (c) 2013-2020, 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:Indexingmod.f90
!--------------------------------------------------------------------------
!
! MODULE: Indexingmod
!
!> @author Saransh Singh, Carnegie Mellon University
!
!> @brief everything to do with indexing and such
!
!> @date 11/17/15 SS  1.0 original
!> @date 11/02/19 MDG 1.1 moved EBSDDISubroutine from EBSDDImod to this module (due to OpenCL usage)
!---------------------------------------------------------------------------
module Indexingmod

use local
use NameListTypedefs

IMPLICIT NONE

contains

!--------------------------------------------------------------------------
!
! SUBROUTINE:EBSDDISubroutine
!
!> @author Saransh Singh, Carnegie Mellon University
!
!> @brief Master subroutine to control dictionary generation, inner products computations, sorting values
!> and indexing of points, all in parallel using OpenCL/openMP
!
!> @param dinl ped indexing namelist pointer
!> @param mcnl Monte Carlo namelist structure
!> @param mpnl Master PAttern namelist structure
!> @param EBSDMCdata Monte Carlo data arrays
!> @param EBSDMPdata master pattern data arrays
!> @param EBSDdetector detector arrays
!> @param progname name of the program
!> @param nmldeffile namelist filename
!
!> @date 03/30/15  SS 1.0 original
!> @date 05/05/15 MDG 1.1 removed getenv() call; replaced by global path strings
!> @date 02/07/16 MDG 1.2 added image quality computation
!> @date 02/08/16 MDG 1.3 added confidence index output to HDF5 file
!> @date 02/10/16 MDG 1.4 added average dot product map to HDF5 file
!> @date 02/23/16 MDG 1.5 converted program to CLFortran instead of fortrancl.
!> @date 06/07/16 MDG 1.6 added tmpfile for variable temporary data file name
!> @date 11/14/16 MDG 1.7 added code to read h5 file instead of compute on-the-fly (static mode)
!> @date 07/24/17 MDG 1.8 temporary code to change the mask layout for fast DI tests
!> @date 08/30/17 MDG 1.9 added option to read custom mask from file
!> @date 11/13/17 MDG 2.0 moved OpenCL code from InnerProdGPU routine to main code
!> @date 01/09/18 MDG 2.1 first attempt at OpenMP for pattern pre-processing
!> @date 02/13/18 MDG 2.2 added support for multiple input formats for experimental patterns
!> @date 02/22/18 MDG 2.3 added support for Region-of-Interest (ROI) selection
!> @date 04/04/18 MDG 3.0 revised name list use as well as MC and MP data structures
!> @date 11/10/18 MDG 3.1 added EDAX/TSL .ang output format
!> @date 05/19/19 MDG 3.2 removes several averaging routines that were never used
!> @date 06/19/19 MDG 3.3 add support for overlap master patterns generated by EMEBSDoverlap.f90
!--------------------------------------------------------------------------
subroutine EBSDDISubroutine(dinl, mcnl, mpnl, EBSDMCdata, EBSDMPdata, EBSDdetector, progname, nmldeffile)
!DEC$ ATTRIBUTES DLLEXPORT :: EBSDDISubroutine

  use local
  use typedefs
  use NameListTypedefs
  use NameListHandlers
  use files
  use dictmod
  use patternmod
  use Lambert
  use others
  use crystal
  use initializersHDF
  use gvectors
  use filters
  use error
  use io
  use diffraction
  use symmetry
  use quaternions
  use constants
  use rotations
  use so3
  use math
  use EBSDmod
  use EBSDDImod
  use clfortran
  use CLsupport
  use omp_lib
  use HDF5
  use h5im
  use h5lt
  use HDFsupport
  use EMh5ebsd
  use commonmod
  use EBSDiomod
  use NameListHDFwriters
  use ECPmod, only: GetPointGroup
  use ISO_C_BINDING
  use notifications
  use timing

  IMPLICIT NONE

  type(EBSDIndexingNameListType),INTENT(INOUT)        :: dinl
!f2py intent(in,out) ::  dinl
  type(MCCLNameListType),INTENT(INOUT)                :: mcnl
!f2py intent(in,out) ::  mcnl
  type(EBSDMasterNameListType),INTENT(INOUT)          :: mpnl
!f2py intent(in,out) ::  mpnl
  type(EBSDMCdataType),INTENT(INOUT)                  :: EBSDMCdata
!f2py intent(in,out) ::  EBSDMCdata
  type(EBSDMPdataType),INTENT(INOUT)                  :: EBSDMPdata
!f2py intent(in,out) ::  EBSDMPdata
  type(EBSDDetectorType),INTENT(INOUT)                :: EBSDdetector
!f2py intent(in,out) ::  EBSDdetector
  character(fnlen),INTENT(IN)                         :: progname
  character(fnlen),INTENT(IN)                         :: nmldeffile

  type(unitcell)                                      :: cell
  type(DynType)                                       :: Dyn
  type(gnode)                                         :: rlp
  logical                                             :: verbose

  integer(c_intptr_t),allocatable, target             :: platform(:)
  integer(c_intptr_t),allocatable, target             :: device(:)
  integer(c_intptr_t),target                          :: context
  integer(c_intptr_t),target                          :: command_queue
  integer(c_intptr_t),target                          :: cl_expt,cl_dict
  character(len = 50000), target                      :: source
  integer(kind=irg), parameter                        :: source_length = 50000
  integer(kind=irg), target                           :: source_l
  character(len=source_length, KIND=c_char),TARGET    :: csource
  type(c_ptr), target                                 :: psource
  integer(c_int32_t)                                  :: ierr2, pcnt
  integer(c_intptr_t),target                          :: prog
  integer(c_intptr_t),target                          :: kernel
  integer(c_size_t)                                   :: cnum
  character(9),target                                 :: kernelname
  character(10, KIND=c_char),target                   :: ckernelname

  integer(kind=irg)                                   :: num,ierr,irec,istat, jpar(7), SGnum, nlines
  integer(kind=irg),parameter                         :: iunit = 40
  integer(kind=irg),parameter                         :: iunitexpt = 41
  integer(kind=irg),parameter                         :: iunitdict = 42
  character(fnlen)                                    :: info ! info about the GPU
  real(kind=dbl),parameter                            :: nAmpere = 6.241D+18   ! Coulomb per second

  integer(kind=irg)                                   :: Ne,Nd,L,totnumexpt,numdictsingle,numexptsingle,imght,imgwd,nnk, &
                                                         recordsize, fratio, cratio, fratioE, cratioE, iii, itmpexpt, hdferr,&
                                                         recordsize_correct, patsz, tickstart, tickstart2, tock, npy, sz(3), jjj
  integer(kind=8)                                     :: size_in_bytes_dict,size_in_bytes_expt
  real(kind=sgl),pointer                              :: dict(:), T0dict(:)
  real(kind=sgl),allocatable,TARGET                   :: dict1(:), dict2(:)
  !integer(kind=1),allocatable                         :: imageexpt(:),imagedict(:)
  real(kind=sgl),allocatable                          :: imageexpt(:),imagedict(:), mask(:,:),masklin(:), exptIQ(:), &
                                                         exptCI(:), exptFit(:), exppatarray(:), tmpexppatarray(:)
  real(kind=sgl),allocatable                          :: imageexptflt(:),binned(:,:),imagedictflt(:),imagedictfltflip(:), &
                                                         tmpimageexpt(:), OSMmap(:,:)
  real(kind=sgl),allocatable, target                  :: results(:),expt(:),dicttranspose(:),resultarray(:),&
                                                         eulerarray(:,:),eulerarray2(:,:),resultmain(:,:),resulttmp(:,:)
  integer(kind=irg),allocatable                       :: acc_array(:,:), ppend(:), ppendE(:)
  integer*4,allocatable                               :: iexptCI(:,:), iexptIQ(:,:)
  real(kind=sgl),allocatable                          :: meandict(:),meanexpt(:),wf(:),mLPNH(:,:,:),mLPSH(:,:,:),accum_e_MC(:,:,:)
  real(kind=sgl),allocatable                          :: mLPNH_simple(:,:), mLPSH_simple(:,:), eangle(:)
  real(kind=sgl),allocatable                          :: EBSDpattern(:,:), FZarray(:,:), dpmap(:), lstore(:,:), pstore(:,:)
  real(kind=sgl),allocatable                          :: EBSDpatternintd(:,:), lp(:), cp(:), EBSDpat(:,:)
  integer(kind=irg),allocatable                       :: EBSDpatterninteger(:,:), EBSDpatternad(:,:), EBSDpint(:,:)
  character(kind=c_char),allocatable                  :: EBSDdictpat(:,:,:)
  real(kind=sgl),allocatable                          :: EBSDdictpatflt(:,:)
  real(kind=dbl),allocatable                          :: rdata(:,:), fdata(:,:), rrdata(:,:), ffdata(:,:), ksqarray(:,:)
  complex(kind=dbl),allocatable                       :: hpmask(:,:)
  complex(C_DOUBLE_COMPLEX),allocatable               :: inp(:,:), outp(:,:)
  real(kind=dbl)                                      :: w, Jres
  integer(kind=irg)                                   :: dims(2)
  character(11)                                       :: dstr
  character(15)                                       :: tstrb
  character(15)                                       :: tstre
  character(3)                                        :: vendor
  character(fnlen, KIND=c_char),allocatable,TARGET    :: stringarray(:)
  character(fnlen)                                    :: groupname, dataset, fname, clname, ename, sourcefile, &
                                                         datagroupname, dictfile, attname
  integer(hsize_t)                                    :: expwidth, expheight
  integer(hsize_t),allocatable                        :: iPhase(:), iValid(:)
  integer(c_size_t),target                            :: slength
  integer(c_int)                                      :: numd, nump
  type(C_PTR)                                         :: planf, HPplanf, HPplanb
  integer(HSIZE_T)                                    :: dims2(2), offset2(2), dims3(3), offset3(3)

  integer(kind=irg)                                   :: i,j,ii,jj,kk,ll,mm,pp,qq
  integer(kind=irg)                                   :: FZcnt, pgnum, io_int(4), ncubochoric, pc
  type(FZpointd),pointer                              :: FZlist, FZtmp
  integer(kind=irg),allocatable                       :: indexlist(:),indexarray(:),indexmain(:,:),indextmp(:,:)
  real(kind=sgl)                                      :: dmin,voltage,scl,ratio, mi, ma, ratioE, io_real(2), tstart, tmp, &
                                                         totnum_el, vlen, tstop, ttime
  real(kind=dbl)                                      :: prefactor
  character(fnlen)                                    :: xtalname
  integer(kind=irg)                                   :: binx,biny,TID,nthreads,Emin,Emax, iiistart, iiiend, jjend
  real(kind=sgl)                                      :: sx,dx,dxm,dy,dym,rhos,x,projweight, dp, mvres, nel, emult
  real(kind=sgl)                                      :: dc(3),quat(4),ixy(2),bindx
  integer(kind=irg)                                   :: nix,niy,nixp,niyp
  real(kind=sgl)                                      :: euler(3)
  integer(kind=irg)                                   :: indx
  integer(kind=irg)                                   :: correctsize
  logical                                             :: f_exists, init, ROIselected

  integer(kind=irg)                                   :: ipar(10)

  character(fnlen),ALLOCATABLE                        :: MessageLines(:)
  integer(kind=irg)                                   :: NumLines
  character(fnlen)                                    :: TitleMessage, exectime
  character(100)                                      :: c
  character(1000)                                     :: charline
  character(3)                                        :: stratt

  type(HDFobjectStackType)                            :: HDF_head

  call timestamp(datestring=dstr, timestring=tstrb)

  if (trim(dinl%indexingmode).eq.'static') then
      !
      ! Initialize FORTRAN interface.
      CALL h5open_EMsoft(hdferr)

      ! get the full filename
      dictfile = trim(EMsoft_getEMdatapathname())//trim(dinl%dictfile)
      dictfile = EMsoft_toNativePath(dictfile)

      call Message('-->  '//'Opening HDF5 dictionary file '//trim(dinl%dictfile))
      nullify(HDF_head%next)

      hdferr =  HDF_openFile(dictfile, HDF_head)
      if (hdferr.ne.0) call HDF_handleError(hdferr,'HDF_openFile ')

      ! we need the point group number (derived from the space group number)
      ! if EBSDMPdata%newSGnumber is set to 2, then pgnum must be set to 1 for
      ! overlap master patterns  [ added by MDG, 06/19/19 ]
      if (EBSDMPdata%AveragedMP.eqv..TRUE.) then
          pgnum = EBSDMPdata%newPGnumber
          io_int = pgnum
          call WriteValue(' Setting point group number to ',io_int,1)
      else
          groupname = SC_CrystalData
          hdferr = HDF_openGroup(groupname, HDF_head)
          if (hdferr.ne.0) call HDF_handleError(hdferr,'HDF_openGroup:CrystalData')

          dataset = SC_SpaceGroupNumber
          call HDF_readDatasetInteger(dataset, HDF_head, hdferr, SGnum)
          if (hdferr.ne.0) call HDF_handleError(hdferr,'HDF_readDatasetInteger:SpaceGroupNumber')
          call HDF_pop(HDF_head)
      ! get the point group number
          if (SGnum.ge.221) then
            pgnum = 32
          else
            i=0
            do while (SGPG(i+1).le.SGnum)
              i = i+1
            end do
            pgnum = i
          end if
          io_int = pgnum
          call WriteValue(' Setting point group number to ',io_int,1)
      end if

      ! then read some more data from the EMData group
      groupname = SC_EMData
      hdferr = HDF_openGroup(groupname, HDF_head)
      if (hdferr.ne.0) call HDF_handleError(hdferr,'HDF_openGroup:EMData')

      datagroupname = 'EBSD'
      hdferr = HDF_openGroup(datagroupname, HDF_head)
      if (hdferr.ne.0) call HDF_handleError(hdferr,'HDF_openGroup:EBSD')

  ! test the HDF_FileVersion to make sure that the dictionary file is recent enough
      attname = 'HDF_FileVersion'
      hdferr = HDF_getStringAttributeFromGroup(attname, stratt, 3_SIZE_T, HDF_head)

      if (stratt.eq.'4.0') then
          call Message('The dictionary file was created with an older version of the EMEBSD program.')
          call Message('This file can not be used by the present program; must be version 4.1 or higher.')
          call Message('')
          call FatalError('EBSDDISubroutine','Incompatible dictionary file; please rerun the EMEBSD program.')
      end if

      ! we already have the xtalname string from the Monte Carlo name list
      xtalname = trim(mcnl%xtalname)

      ! number of Eulerangles numangles
      dataset = SC_numangles
      call HDF_readDatasetInteger(dataset, HDF_head, hdferr, FZcnt)
      if (hdferr.ne.0) call HDF_handleError(hdferr,'HDF_readDatasetInteger:numangles')

      ! euler angle list Eulerangles
      dataset = SC_Eulerangles
      call HDF_readDatasetFloatArray2D(dataset, dims2, HDF_head, hdferr, eulerarray2)
      eulerarray2 = eulerarray2 * 180.0/sngl(cPi)
      if (hdferr.ne.0) call HDF_handleError(hdferr,'HDF_readDatasetFloatArray2D:Eulerangles')

      ! we leave this file open since we still need to read all the patterns...
      !=====================================================
      call Message('-->  completed initial reading of dictionary file ')
  end if

  if (sum(dinl%ROI).ne.0) then
    ROIselected = .TRUE.
    iiistart = dinl%ROI(2)
    iiiend = dinl%ROI(2)+dinl%ROI(4)-1
    jjend = dinl%ROI(3)
  else
    ROIselected = .FALSE.
    iiistart = 1
    iiiend = dinl%ipf_ht
    jjend = dinl%ipf_wd
  end if

  verbose = .FALSE.
  init = .TRUE.
  Ne = dinl%numexptsingle
  Nd = dinl%numdictsingle
  L = dinl%numsx*dinl%numsy/dinl%binning**2
  if (ROIselected.eqv..TRUE.) then
      totnumexpt = dinl%ROI(3)*dinl%ROI(4)
  else
      totnumexpt = dinl%ipf_wd*dinl%ipf_ht
  end if
  imght = dinl%numsx/dinl%binning
  imgwd = dinl%numsy/dinl%binning
  dims = (/imght, imgwd/)
  nnk = dinl%nnk
  ncubochoric = dinl%ncubochoric
  recordsize = L*4
  itmpexpt = 43
  w = dinl%hipassw
  source_l = source_length

  ! these will eventually need to be read from an experimental data file but we'll set defaults values here.
  dinl%WD = 10.0

  ! nullify the dict and T0dict pointers
  nullify(dict,T0dict)

  ! make sure that correctsize is a multiple of 16; if not, make it so
  if (mod(L,16) .ne. 0) then
      correctsize = 16*ceiling(float(L)/16.0)
  else
      correctsize = L
  end if

  ! determine the experimental and dictionary sizes in bytes
  size_in_bytes_dict = Nd*correctsize*sizeof(correctsize)
  size_in_bytes_expt = Ne*correctsize*sizeof(correctsize)
  recordsize_correct = correctsize*4
  patsz              = correctsize

  if (trim(dinl%indexingmode).eq.'dynamic') then
      !=====================================================
      ! EXTRACT POINT GROUP NUMBER FROM CRYSTAL STRUCTURE FILE
      !=====================================================
      write (*,*) 'reading from xtalfile '//trim(mcnl%xtalname)

      if (EBSDMPdata%AveragedMP.eqv..TRUE.) then
          pgnum = EBSDMPdata%newPGnumber
      else
          pgnum = GetPointGroup(mcnl%xtalname)
      end if

      !=====================================================
      ! make sure the minimum energy is set smaller than the maximum
      !=====================================================
      if (dinl%energymin.gt.dinl%energymax) then
          call Message('Minimum energy is larger than maximum energy; please correct input file')
          stop
      end if

      !=====================================================
      ! get the indices of the minimum and maximum energy
      !=====================================================
      Emin = nint((dinl%energymin - mcnl%Ehistmin)/mcnl%Ebinsize) +1
      if (Emin.lt.1)  Emin=1
      if (Emin.gt.EBSDMCdata%numEbins)  Emin=EBSDMCdata%numEbins

      Emax = nint((dinl%energymax - mcnl%Ehistmin)/mcnl%Ebinsize) + 1
      if (Emax .lt. 1) Emax = 1
      if (Emax .gt. EBSDMCdata%numEbins) Emax = EBSDMCdata%numEbins

      sz = shape(EBSDMPdata%mLPNH)
      dinl%nE = sz(3)

      ! intensity prefactor
      nel = float(mcnl%totnum_el) * float(EBSDMCdata%multiplier)
      emult = nAmpere * 1e-9 / nel  ! multiplicative factor to convert MC data to an equivalent incident beam of 1 nanoCoulomb
      ! intensity prefactor  (redefined by MDG, 3/23/18)
      prefactor = emult * dinl%beamcurrent * dinl%dwelltime * 1.0D-6
  end if

  !====================================
  ! init a bunch of parameters
  !====================================
  ! binned pattern array
  binx = dinl%numsx/dinl%binning
  biny = dinl%numsy/dinl%binning
  bindx = 1.0/float(dinl%binning)**2

  ! allocate the square-Lambert arrays
  npy = mpnl%npx
  if (trim(dinl%indexingmode).eq.'dynamic') then
    allocate(mLPNH(-mpnl%npx:mpnl%npx,-npy:npy,EBSDMCdata%numEbins))
    allocate(mLPSH(-mpnl%npx:mpnl%npx,-npy:npy,EBSDMCdata%numEbins))
    allocate(accum_e_MC(EBSDMCdata%numEbins,dinl%numsx,dinl%numsy),stat=istat)
    accum_e_MC = EBSDdetector%accum_e_detector
    mLPNH = EBSDMPdata%mLPNH
    mLPSH = EBSDMPdata%mLPSH
  end if

  !=====================================================
  ! SAMPLING OF RODRIGUES FUNDAMENTAL ZONE
  !=====================================================
  ! if eulerfile is not defined, then we use the standard RFZ sampling;
  ! if it is defined, then we read the Eulerangle triplets from the file
  ! and generate the FZlist here... this can be useful to index patterns that
  ! have only a small misorientation range with respect to a known orientation,
  ! so that it is not necessary to scan all of orientation space.
  if (trim(dinl%indexingmode).eq.'dynamic') then
      nullify(FZlist)
      FZcnt = 0
      if (trim(dinl%eulerfile).eq.'undefined') then
        call Message('Orientation space sampling mode set to RFZ')
        io_int(1) = pgnum
        io_int(2) = ncubochoric
        call WriteValue('Point group number and number of cubochoric sampling points : ',io_int,2,"(I4,',',I5)")

        call sampleRFZ(ncubochoric, pgnum, 0, FZcnt, FZlist)
      else
      ! read the euler angle file and create the linked list
        call getEulersfromFile(dinl%eulerfile, FZcnt, FZlist)
        call Message('Orientation space sampling mode set to MIS')
        io_int(1) = pgnum
        io_int(2) = FZcnt
        call WriteValue('Point group number and number of sampling points : ',io_int,2,"(I4,',',I5)")
      end if

      ! allocate and fill FZarray for OpenMP parallelization
      allocate(FZarray(4,FZcnt),stat=istat)
      FZarray = 0.0

      FZtmp => FZlist
      do ii = 1,FZcnt
          FZarray(1:4,ii) = FZtmp%rod(1:4)
          FZtmp => FZtmp%next
      end do
      io_int(1) = FZcnt
      call WriteValue(' Number of unique orientations sampled =        : ', io_int, 1, "(I8)")
  end if

  !================================
  ! INITIALIZATION OF OpenCL DEVICE
  !================================
  call Message('--> Initializing OpenCL device')

  call CLinit_PDCCQ(platform, nump, dinl%platid, device, numd, dinl%devid, info, context, command_queue)

  ! read the cl source file
  sourcefile = 'DictIndx.cl'
  call CLread_source_file(sourcefile, csource, slength)

  ! allocate device memory for experimental and dictionary patterns
  cl_expt = clCreateBuffer(context, CL_MEM_READ_WRITE, size_in_bytes_expt, C_NULL_PTR, ierr)
  call CLerror_check('EBSDDISubroutine:clCreateBuffer', ierr)

  cl_dict = clCreateBuffer(context, CL_MEM_READ_WRITE, size_in_bytes_dict, C_NULL_PTR, ierr)
  call CLerror_check('EBSDDISubroutine:clCreateBuffer', ierr)

  !================================
  ! the following lines were originally in the InnerProdGPU routine, but there is no need
  ! to execute them each time that routine is called so we move them here...
  !================================
  ! create the program
  pcnt = 1
  psource = C_LOC(csource)
  !prog = clCreateProgramWithSource(context, pcnt, C_LOC(psource), C_LOC(source_l), ierr)
  prog = clCreateProgramWithSource(context, pcnt, C_LOC(psource), C_LOC(slength), ierr)
  call CLerror_check('InnerProdGPU:clCreateProgramWithSource', ierr)

  ! build the program
  ierr = clBuildProgram(prog, numd, C_LOC(device), C_NULL_PTR, C_NULL_FUNPTR, C_NULL_PTR)

  ! get the compilation log
  ierr2 = clGetProgramBuildInfo(prog, device(dinl%devid), CL_PROGRAM_BUILD_LOG, sizeof(source), C_LOC(source), cnum)
  ! if(cnum > 1) call Message(trim(source(1:cnum))//'test',frm='(A)')
  call CLerror_check('InnerProdGPU:clBuildProgram', ierr)
  call CLerror_check('InnerProdGPU:clGetProgramBuildInfo', ierr2)

  ! finally get the kernel and release the program
  kernelname = 'InnerProd'
  ckernelname = kernelname
  ckernelname(10:10) = C_NULL_CHAR
  kernel = clCreateKernel(prog, C_LOC(ckernelname), ierr)
  call CLerror_check('InnerProdGPU:clCreateKernel', ierr)

  ierr = clReleaseProgram(prog)
  call CLerror_check('InnerProdGPU:clReleaseProgram', ierr)

  ! the remainder is done in the InnerProdGPU routine
  !=========================================

  !=========================================
  ! ALLOCATION AND INITIALIZATION OF ARRAYS
  !=========================================
  call Message('--> Allocating various arrays for indexing')

  allocate(expt(Ne*correctsize),stat=istat)
  if (istat .ne. 0) stop 'Could not allocate array for experimental patterns'
  expt = 0.0

  allocate(dict1(Nd*correctsize),dict2(Nd*correctsize),dicttranspose(Nd*correctsize),stat=istat)
  if (istat .ne. 0) stop 'Could not allocate array for dictionary patterns'
  dict1 = 0.0
  dict2 = 0.0
  dict => dict1
  dicttranspose = 0.0

  allocate(results(Ne*Nd),stat=istat)
  if (istat .ne. 0) stop 'Could not allocate array for results'
  results = 0.0

  allocate(mask(binx,biny),masklin(L),stat=istat)
  if (istat .ne. 0) stop 'Could not allocate arrays for masks'
  mask = 1.0
  masklin = 0.0

  allocate(imageexpt(L),imageexptflt(correctsize),imagedictflt(correctsize),imagedictfltflip(correctsize),stat=istat)
  allocate(tmpimageexpt(correctsize),stat=istat)
  if (istat .ne. 0) stop 'Could not allocate array for reading experimental image patterns'
  imageexpt = 0.0
  imageexptflt = 0.0

  allocate(meandict(correctsize),meanexpt(correctsize),imagedict(correctsize),stat=istat)
  if (istat .ne. 0) stop 'Could not allocate array for mean dictionary and experimental patterns'
  meandict = 0.0
  meanexpt = 0.0

  ! allocate(EBSDpattern(dinl%numsx,dinl%numsy),binned(binx,biny),stat=istat)
  allocate(EBSDpattern(binx,biny),binned(binx,biny),stat=istat)
  if (istat .ne. 0) stop 'Could not allocate array for EBSD pattern'
  EBSDpattern = 0.0
  binned = 0.0

  allocate(resultarray(1:Nd),stat=istat)
  if (istat .ne. 0) stop 'could not allocate result arrays'

  resultarray = 0.0

  allocate(indexarray(1:Nd),stat=istat)
  if (istat .ne. 0) stop 'could not allocate index arrays'

  indexarray = 0

  allocate(indexlist(1:Nd*(ceiling(float(FZcnt)/float(Nd)))),stat=istat)
  if (istat .ne. 0) stop 'could not allocate indexlist arrays'

  indexlist = 0

  do ii = 1,Nd*ceiling(float(FZcnt)/float(Nd))
      indexlist(ii) = ii
  end do

  allocate(resultmain(nnk,Ne*ceiling(float(totnumexpt)/float(Ne))),stat=istat)
  if (istat .ne. 0) stop 'could not allocate main result array'

  resultmain = -2.0

  allocate(indexmain(nnk,Ne*ceiling(float(totnumexpt)/float(Ne))),stat=istat)
  if (istat .ne. 0) stop 'could not allocate main index array'

  indexmain = 0

  allocate(resulttmp(2*nnk,Ne*ceiling(float(totnumexpt)/float(Ne))),stat=istat)
  if (istat .ne. 0) stop 'could not allocate temporary result array'

  resulttmp = -2.0

  allocate(indextmp(2*nnk,Ne*ceiling(float(totnumexpt)/float(Ne))),stat=istat)
  if (istat .ne. 0) stop 'could not allocate temporary index array'

  indextmp = 0

  allocate(eulerarray(1:3,Nd*ceiling(float(FZcnt)/float(Nd))),stat=istat)
  if (istat .ne. 0) stop 'could not allocate euler array'
  eulerarray = 0.0
  if (trim(dinl%indexingmode).eq.'static') then
      eulerarray(1:3,1:FZcnt) = eulerarray2(1:3,1:FZcnt)
      deallocate(eulerarray2)
  end if

  allocate(exptIQ(totnumexpt), exptCI(totnumexpt), exptFit(totnumexpt), stat=istat)
  if (istat .ne. 0) stop 'could not allocate exptIQ array'

  allocate(rdata(binx,biny),fdata(binx,biny),stat=istat)
  if (istat .ne. 0) stop 'could not allocate arrays for Hi-Pass filter'
  rdata = 0.D0
  fdata = 0.D0

  ! also, allocate the arrays used to create the average dot product map; this will require
  ! reading the actual EBSD HDF5 file to figure out how many rows and columns there
  ! are in the region of interest.  For now we get those from the nml until we actually
  ! implement the HDF5 reading bit
  ! this portion of code was first tested in IDL.
  allocate(EBSDpatterninteger(binx,biny))
  EBSDpatterninteger = 0
  allocate(EBSDpatternad(binx,biny),EBSDpatternintd(binx,biny))
  EBSDpatternad = 0.0
  EBSDpatternintd = 0.0

  !=====================================================
  ! determine loop variables to avoid having to duplicate
  ! large sections of mostly identical code
  !=====================================================
  ratio = float(FZcnt)/float(Nd)
  cratio = ceiling(ratio)
  fratio = floor(ratio)

  ratioE = float(totnumexpt)/float(Ne)
  cratioE = ceiling(ratioE)
  fratioE = floor(ratioE)

  allocate(ppend(cratio),ppendE(cratioE))
  ppend = (/ (Nd, i=1,cratio) /)
  if (fratio.lt.cratio) then
    ppend(cratio) = MODULO(FZcnt,Nd)
  end if

  ppendE = (/ (Ne, i=1,cratioE) /)
  if (fratioE.lt.cratioE) then
    ppendE(cratioE) = MODULO(totnumexpt,Ne)
  end if

  !=====================================================
  ! define the circular mask if necessary and convert to 1D vector
  !=====================================================

  if (trim(dinl%maskfile).ne.'undefined') then
  ! read the mask from file; the mask can be defined by a 2D array of 0 and 1 values
  ! that is stored in row form as strings, e.g.
  !    0000001110000000
  !    0000011111000000
  ! ... etc
  !
      f_exists = .FALSE.
      fname = trim(EMsoft_getEMdatapathname())//trim(dinl%maskfile)
      fname = EMsoft_toNativePath(fname)
      inquire(file=trim(fname), exist=f_exists)
      if (f_exists.eqv..TRUE.) then
        mask = 0.0
        open(unit=dataunit,file=trim(fname),status='old',form='formatted')
        do jj=biny,1,-1
          read(dataunit,"(A)") charline
          do ii=1,binx
            if (charline(ii:ii).eq.'1') mask(ii,jj) = 1.0
          end do
        end do
        close(unit=dataunit,status='keep')
      else
        call FatalError('EBSDDISubroutine','maskfile '//trim(fname)//' does not exist')
      end if
  else
      if (dinl%maskpattern.eq.'y') then
        do ii = 1,biny
            do jj = 1,binx
                if((ii-biny/2)**2 + (jj-binx/2)**2 .ge. dinl%maskradius**2) then
                    mask(jj,ii) = 0.0
                end if
            end do
        end do
      end if
  end if

  ! convert the mask to a linear (1D) array
  do ii = 1,biny
      do jj = 1,binx
          masklin((ii-1)*binx+jj) = mask(jj,ii)
      end do
  end do

  !=====================================================
  ! Preprocess all the experimental patterns and store
  ! them in a temporary file as vectors; also, create
  ! an average dot product map to be stored in the h5ebsd output file
  !=====================================================
  call h5open_EMsoft(hdferr)
  call PreProcessPatterns(dinl%nthreads, .FALSE., dinl, binx, biny, masklin, correctsize, totnumexpt, exptIQ=exptIQ)
  call h5close_EMsoft(hdferr)

  !=====================================================
  call Message(' -> computing Average Dot Product map (ADP)')
  call Message(' ')

  ! re-open the temporary file
  fname = trim(EMsoft_getEMtmppathname())//trim(dinl%tmpfile)
  fname = EMsoft_toNativePath(fname)

  open(unit=itmpexpt,file=trim(fname),&
       status='old',form='unformatted',access='direct',recl=recordsize_correct,iostat=ierr)

  ! use the getADPmap routine in the filters module
  if (ROIselected.eqv..TRUE.) then
    allocate(dpmap(dinl%ROI(3)*dinl%ROI(4)))
    call getADPmap(itmpexpt, dinl%ROI(3)*dinl%ROI(4), L, dinl%ROI(3), dinl%ROI(4), dpmap)
  else
    allocate(dpmap(totnumexpt))
    call getADPmap(itmpexpt, totnumexpt, L, dinl%ipf_wd, dinl%ipf_ht, dpmap)
  end if

  ! we will leave the itmpexpt file open, since we'll be reading from it again...

  !=====================================================
  ! MAIN COMPUTATIONAL LOOP
  !
  ! Some explanation is necessary here... the bulk of this code is
  ! executed in OpenMP multithreaded mode, with nthreads threads.
  ! Thread 0 has a special role described below; threads 1 ... nthreads-1
  ! share the computation of the dictionary patterns, and wait for
  ! thread 0 to finish, if necessary.
  !
  ! Thread 0 takes the dictionary patterns computed by the other threads
  ! in the previous step in the dictionaryloop and sends them to the GPU,
  ! along with as many chunks of experimental data are to be handled (experimentalloop
  ! inside the thread 0 portion of the code); the experimental patterns
  ! are then read from the temporary file (unit itmpexpt).  Once all dot
  ! products have been computed by the GPU, thread 0 will rank them largest
  ! to smallest and keep only the top nnk values along with their indices
  ! into the array of Euler angle triplets.  If the other threads are still
  ! computing dictionary patterns, thread 0 will join them; otherwise
  ! thread 0 will immediately take the next batch of dictionary patterns
  ! and start all over.
  !
  ! The trick is for the user to determine the array chunk sizes so that
  ! threads 1 ... nthreads-1 do not have to wait for thread 0 to finish;
  ! this requires a bit of experimenting and observing the load on all the
  ! system cores.  The load should always be approximately 100% x nthreads-1
  ! for an efficient execution.  The appropriate number of threads will depend
  ! on how powerful the GPU card is...
  !=====================================================

  call Time_tick(tickstart)
  call Time_tick(tickstart2)

  if (trim(dinl%indexingmode).eq.'dynamic') then
      call OMP_SET_NUM_THREADS(dinl%nthreads)
      io_int(1) = dinl%nthreads
  else
      call OMP_SET_NUM_THREADS(2)
      io_int(1) = 2
  end if
  call WriteValue(' -> Number of threads set to ',io_int,1,"(I3)")

  ! define the jpar array of integer parameters
  jpar(1) = dinl%binning
  jpar(2) = dinl%numsx
  jpar(3) = dinl%numsy
  jpar(4) = mpnl%npx
  jpar(5) = npy
  jpar(6) = EBSDMCdata%numEbins
  jpar(7) = dinl%nE

  call timestamp()

  dictionaryloop: do ii = 1,cratio+1
      results = 0.0

  ! if ii is odd, then we use dict1 for the dictionary computation, and dict2 for the GPU
  ! (assuming ii>1); when ii is even we switch the two pointers
      if (mod(ii,2).eq.1) then
        dict => dict1
        dict1 = 0.0
        T0dict => dict2   ! these are the patterns to be sent to the GPU
        if (verbose.eqv..TRUE.) call WriteValue('','dict => dict1; T0dict => dict2')
      else
        dict => dict2
        dict2 = 0.0
        T0dict => dict1   ! these are the patterns to be sent to the GPU
        if (verbose.eqv..TRUE.) call WriteValue('','dict => dict2; T0dict => dict1')
      end if

      if (verbose.eqv..TRUE.) then
        io_int(1) = ii
        call WriteValue('Dictionaryloop index = ',io_int,1)
      end if

  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(TID,iii,jj,ll,mm,pp,ierr,io_int, tock, ttime) &
  !$OMP& PRIVATE(binned, ma, mi, EBSDpatternintd, EBSDpatterninteger, EBSDpatternad, quat, imagedictflt,imagedictfltflip)

          TID = OMP_GET_THREAD_NUM()

        if ((ii.eq.1).and.(TID.eq.0)) write(*,*) ' actual number of OpenMP threads  = ',OMP_GET_NUM_THREADS()

  ! the master thread should be the one working on the GPU computation
  !$OMP MASTER
      if (ii.gt.1) then
        iii = ii-1        ! the index ii is already one ahead, since the GPU thread lags one cycle behind the others...
        if (verbose.eqv..TRUE.) then
          if (associated(T0dict,dict1)) then
            write(*,"('   GPU thread is working on dict1')")
          else
            write(*,"('   GPU thread is working on dict2')")
          end if
        end if

        dicttranspose = 0.0

        do ll = 1,correctsize
          do mm = 1,Nd
              dicttranspose((ll-1)*Nd+mm) = T0dict((mm-1)*correctsize+ll)
          end do
        end do

        ierr = clEnqueueWriteBuffer(command_queue, cl_dict, CL_TRUE, 0_8, size_in_bytes_dict, C_LOC(dicttranspose(1)), &
                                    0, C_NULL_PTR, C_NULL_PTR)
        call CLerror_check('EBSDDISubroutine:clEnqueueWriteBuffer', ierr)

        mvres = 0.0

        experimentalloop: do jj = 1,cratioE

          expt = 0.0

          do pp = 1,ppendE(jj)   ! Ne or MODULO(totnumexpt,Ne)
            read(itmpexpt,rec=(jj-1)*Ne+pp) tmpimageexpt
            expt((pp-1)*correctsize+1:pp*correctsize) = tmpimageexpt
          end do

          ierr = clEnqueueWriteBuffer(command_queue, cl_expt, CL_TRUE, 0_8, size_in_bytes_expt, C_LOC(expt(1)), &
                                      0, C_NULL_PTR, C_NULL_PTR)
          call CLerror_check('EBSDDISubroutine:clEnqueueWriteBuffer', ierr)

          call InnerProdGPU(cl_expt,cl_dict,Ne,Nd,correctsize,results,numd,dinl%devid,kernel,context,command_queue)

          dp =  maxval(results)
          if (dp.gt.mvres) mvres = dp

  ! this might be simplified later for the remainder of the patterns
          do qq = 1,ppendE(jj)
              jjj = (jj-1)*Ne+qq
              resultarray(1:Nd) = results((qq-1)*Nd+1:qq*Nd)
              indexarray(1:Nd) = indexlist((iii-1)*Nd+1:iii*Nd)

              call SSORT(resultarray,indexarray,Nd,-2)
              resulttmp(nnk+1:2*nnk,jjj) = resultarray(1:nnk)
              indextmp(nnk+1:2*nnk,jjj) = indexarray(1:nnk)

              call SSORT(resulttmp(:,jjj),indextmp(:,jjj),2*nnk,-2)

              resultmain(1:nnk,jjj) = resulttmp(1:nnk,jjj)
              indexmain(1:nnk,jjj) = indextmp(1:nnk,jjj)
          end do
        end do experimentalloop

        io_real(1) = mvres
        io_real(2) = float(iii)/float(cratio)*100.0
        call WriteValue('',io_real,2,"(' max. dot product = ',F10.6,';',F6.1,'% complete')")

        if (mod(iii,10) .eq. 0) then
  ! do a remaining time estimate
  ! and print information
          if (iii.eq.10) then
              tock = Time_tock(tickstart)
              ttime = float(tock) * float(cratio) / float(iii)
              tstop = ttime
              io_int(1:4) = (/iii,cratio, int(ttime/3600.0), int(mod(ttime,3600.0)/60.0)/)
              call WriteValue('',io_int,4,"(' -> Completed cycle ',I5,' out of ',I5,'; est. total time ', &
                             I4,' hrs',I3,' min')")
          else
              ttime = tstop * float(cratio-iii) / float(cratio)
              io_int(1:4) = (/iii,cratio, int(ttime/3600.0), int(mod(ttime,3600.0)/60.0)/)
              call WriteValue('',io_int,4,"(' -> Completed cycle ',I5,' out of ',I5,'; est. remaining time ', &
                             I4,' hrs',I3,' min')")
          end if
        end if
      else
         if (verbose.eqv..TRUE.) call WriteValue('','        GPU thread is idling')
      end if  ! ii.gt.1

  !$OMP END MASTER

  ! here we carry out the dictionary pattern computation, unless we are in the ii=cratio+1 step
      if (ii.lt.cratio+1) then
       if (verbose.eqv..TRUE.) then
         if (associated(dict,dict1)) then
           write(*,"('    Thread ',I2,' is working on dict1')") TID
         else
           write(*,"('    Thread ',I2,' is working on dict2')") TID
         end if
       end if

  if (trim(dinl%indexingmode).eq.'dynamic') then
  !$OMP DO SCHEDULE(DYNAMIC)

       do pp = 1,ppend(ii)  !Nd or MODULO(FZcnt,Nd)
         binned = 0.0
         quat = ro2qu(FZarray(1:4,(ii-1)*Nd+pp))

         call CalcEBSDPatternSingleFull(jpar,quat,accum_e_MC,mLPNH,mLPSH,EBSDdetector%rgx,&
                                        EBSDdetector%rgy,EBSDdetector%rgz,binned,Emin,Emax,mask,prefactor)

         if (dinl%scalingmode .eq. 'gam') then
           binned = binned**dinl%gammavalue
         end if

  ! hi pass filtering
  !      rdata = dble(binned)
  !      fdata = HiPassFilter(rdata,dims,w)
  !      binned = sngl(fdata)

  ! adaptive histogram equalization
         ma = maxval(binned)
         mi = minval(binned)

         EBSDpatternintd = ((binned - mi)/ (ma-mi))
         EBSDpatterninteger = nint(EBSDpatternintd*255.0)
         EBSDpatternad =  adhisteq(dinl%nregions,binx,biny,EBSDpatterninteger)
         binned = float(EBSDpatternad)

         imagedictflt = 0.0
         imagedictfltflip = 0.0
         do ll = 1,biny
           do mm = 1,binx
             imagedictflt((ll-1)*binx+mm) = binned(mm,ll)
           end do
         end do

  ! normalize and apply circular mask
         imagedictflt(1:L) = imagedictflt(1:L) * masklin(1:L)
         vlen = vecnorm(imagedictflt(1:correctsize))
         if (vlen.ne.0.0) then
           imagedictflt(1:correctsize) = imagedictflt(1:correctsize)/vlen
         else
           imagedictflt(1:correctsize) = 0.0
         end if

         dict((pp-1)*correctsize+1:pp*correctsize) = imagedictflt(1:correctsize)

         eulerarray(1:3,(ii-1)*Nd+pp) = 180.0/cPi*ro2eu(FZarray(1:4,(ii-1)*Nd+pp))
       end do
  !$OMP END DO

  else  ! we are doing static indexing, so only 2 threads in total

  ! get a set of patterns from the precomputed dictionary file...
  ! we'll use a hyperslab to read a block of preprocessed patterns from file

     if (TID .ne. 0) then
  ! read data from the hyperslab
       dataset = SC_EBSDpatterns
       dims2 = (/ correctsize, ppend(ii) /)
       offset2 = (/ 0, (ii-1)*Nd /)

       if(allocated(EBSDdictpatflt)) deallocate(EBSDdictpatflt)
       EBSDdictpatflt = HDF_readHyperslabFloatArray2D(dataset, offset2, dims2, HDF_head)

       do pp = 1,ppend(ii)  !Nd or MODULO(FZcnt,Nd)
         dict((pp-1)*correctsize+1:pp*correctsize) = EBSDdictpatflt(1:correctsize,pp)
       end do
     end if

  end if

       if (verbose.eqv..TRUE.) then
         io_int(1) = TID
         call WriteValue('',io_int,1,"('       Thread ',I2,' is done')")
       end if
     else
       if (verbose.eqv..TRUE.) then
         io_int(1) = TID
         call WriteValue('',io_int,1,"('       Thread ',I2,' idling')")
       end if
     end if

  ! and we end the parallel section here (all threads will synchronize).
  !$OMP END PARALLEL

  end do dictionaryloop

  if (dinl%keeptmpfile.eq.'n') then
      close(itmpexpt,status='delete')
  else
      close(itmpexpt,status='keep')
  end if

  ! release the OpenCL kernel
  ierr = clReleaseKernel(kernel)
  call CLerror_check('InnerProdGPU:clReleaseKernel', ierr)

  if (trim(dinl%indexingmode).eq.'static') then
  ! close file and nullify pointer
      call HDF_pop(HDF_head,.TRUE.)
      call h5close_EMsoft(hdferr)
  end if

  ! perform some timing stuff
  tstop = Time_tock(tickstart2)
  io_real(1) = tstop
  call WriteValue('Indexing duration (system_clock, s)                : ',io_real,1,"(/,F14.3)")
  io_real(1) = float(totnumexpt)*float(FZcnt) / tstop
  call WriteValue('Number of pattern comparisons per second           : ',io_real,1,"(/,F14.3)")
  io_real(1) = float(totnumexpt) / tstop
  call WriteValue('Number of experimental patterns indexed per second : ',io_real,1,"(/,F14.3,/)")

  ! ===================
  ! MAIN OUTPUT SECTION
  ! ===================

  ! fill the ipar array with integer parameters that are needed to write the h5ebsd file
  ! (anything other than what is already in the dinl structure)
  ipar = 0
  ipar(1) = nnk
  ipar(2) = Ne*ceiling(float(totnumexpt)/float(Ne))
  ipar(3) = totnumexpt
  ipar(4) = Nd*ceiling(float(FZcnt)/float(Nd))
  ipar(5) = FZcnt
  ipar(6) = pgnum
  if (ROIselected.eqv..TRUE.) then
    ipar(7) = dinl%ROI(3)
    ipar(8) = dinl%ROI(4)
  else
    ipar(7) = dinl%ipf_wd
    ipar(8) = dinl%ipf_ht
  end if

  allocate(OSMmap(jjend, iiiend))

  ! Initialize FORTRAN interface.
  call h5open_EMsoft(hdferr)

  if (dinl%datafile.ne.'undefined') then
    vendor = 'TSL'
    call h5ebsd_writeFile(vendor, dinl, mcnl%xtalname, dstr, tstrb, ipar, resultmain, exptIQ, indexmain, eulerarray, &
                          dpmap, progname, nmldeffile, OSMmap)
    call Message('Data stored in h5ebsd file : '//trim(dinl%datafile))
  end if

  if (dinl%ctffile.ne.'undefined') then
    call ctfebsd_writeFile(dinl,mcnl%xtalname,ipar,indexmain,eulerarray,resultmain, OSMmap, exptIQ)
    call Message('Data stored in ctf file : '//trim(dinl%ctffile))
  end if

  if (dinl%angfile.ne.'undefined') then
      call angebsd_writeFile(dinl,mcnl%xtalname,ipar,indexmain,eulerarray,resultmain,exptIQ)
      call Message('Data stored in ang file : '//trim(dinl%angfile))
  end if

  ! close the fortran HDF5 interface
  call h5close_EMsoft(hdferr)

  ! if requested, we notify the user that this program has completed its run
  if (trim(EMsoft_getNotify()).ne.'Off') then
    if (trim(dinl%Notify).eq.'On') then
      NumLines = 3
      allocate(MessageLines(NumLines))

      call hostnm(c)

      MessageLines(1) = 'EMEBSDDI program has ended successfully'
      MessageLines(2) = 'Indexed data stored in '//trim(dinl%datafile)
      write (exectime,"(F15.0)") tstop
      MessageLines(3) = 'Total execution time [s]: '//trim(exectime)
      TitleMessage = 'EMsoft on '//trim(c)
      i = PostMessage(MessageLines, NumLines, TitleMessage)
    end if
  end if

  end subroutine EBSDDISubroutine

!--------------------------------------------------------------------------
!
! SUBROUTINE:InnerProdGPU
!
!> @author Saransh Singh, Carnegie Mellon University
!
!> @brief Perform the inner product computations for the dictionary approach
!
!> @param expt vector with list of observed patterns
!> @param dict vector with list of calculated patterns
!> @param Ne number of patterns in the expt vector
!> @param Nd number of patterns in the dict vector
!> @param L size of one single pattern
!> @param result result of the matrix multiplication
!> @param kernel opencl kernel pointer
!> @param context opencl context type
!> @param command_queue opencl command queue
!
!> @date 12/09/14  SS 1.0 original
!> @date 27/01/15  SS 1.1 modified to call the subroutine from mastersubroutine
!> @date 02/24/16 MDG 1.2 converted OpenCL calls to clfortran from fortrancl
!> @date 03/03/16 MDG 1.3 added C_NULL_CHAR to kernelname
!> @date 06/07/17 MDG 1.4 removed progoptions from Build Program call; caused some issues on Linux in Release mode
!> @date 11/13/17 MDG 2.0 moved several OpenCL init statements to main calling program
!--------------------------------------------------------------------------
recursive subroutine InnerProdGPU(cl_expt,cl_dict,Ne,Nd,correctsize,results,numd,selnumd,kernel,context,command_queue)
!DEC$ ATTRIBUTES DLLEXPORT :: InnerProdGPU

use local
use clfortran
use CLsupport
use ISO_C_BINDING
use io
use error

IMPLICIT NONE

integer(kind=4),INTENT(IN)                          :: Ne
integer(kind=4),INTENT(IN)                          :: Nd
real(kind=4),INTENT(OUT),target                     :: results(Ne*Nd)
integer(c_intptr_t),target,INTENT(INOUT)            :: cl_expt
!f2py intent(in,out) ::  cl_expt
integer(c_intptr_t),target,INTENT(INOUT)            :: cl_dict
!f2py intent(in,out) ::  cl_dict

integer(kind=4),INTENT(IN)                          :: correctsize
integer(kind=irg),INTENT(IN)                        :: numd, selnumd
integer(c_intptr_t),target,INTENT(INOUT)            :: context
!f2py intent(in,out) ::  context
integer(c_intptr_t),target,INTENT(INOUT)            :: kernel
!f2py intent(in,out) ::  kernel
integer(c_intptr_t),target,INTENT(INOUT)            :: command_queue
!f2py intent(in,out) ::  command_queue

integer(c_int32_t)                                  :: ierr, ierr2, pcnt
integer(c_intptr_t),target                          :: cl_result

real(kind=4)                                        :: dicttranspose(Nd*correctsize)
integer(kind=4),parameter                           :: iunit = 40
character(fnlen)                                    :: info ! info about the GPU
integer(kind=8),target                              :: globalsize(2),localsize(2)
integer(kind=4)                                     :: num,istat,i,j,ii,jj,kk, io_int(1)
integer(kind=4),target                              :: Wexp,Wdict
integer(kind=8)                                     :: size_in_bytes_expt,size_in_bytes_dict,size_in_bytes_result
integer(kind=irg)                                   :: irec

size_in_bytes_result = Ne*Nd*sizeof(results(1))
Wexp = correctsize
Wdict = Nd
localsize = (/16,16/)
globalsize = (/Ne,Nd/)

!=====================
! INITIALIZATION [mostly performed in the calling program]
!=====================

! create buffer
cl_result = clCreateBuffer(context, CL_MEM_READ_WRITE, size_in_bytes_result, C_NULL_PTR, ierr)
call CLerror_check('InnerProdGPU:clCreateBuffer', ierr)

! ----

! set kernel arguments
ierr =  clSetKernelArg(kernel, 0, sizeof(cl_expt), C_LOC(cl_expt))
call CLerror_check('InnerProdGPU:clSetKernelArg:cl_expt', ierr)

ierr = clSetKernelArg(kernel, 1, sizeof(cl_dict), C_LOC(cl_dict))
call CLerror_check('InnerProdGPU:clSetKernelArg:cl_dict', ierr)

ierr = clSetKernelArg(kernel, 2, sizeof(Wexp), C_LOC(Wexp))
call CLerror_check('InnerProdGPU:clSetKernelArg:Wexp', ierr)

ierr = clSetKernelArg(kernel, 3, sizeof(Wdict), C_LOC(Wdict))
call CLerror_check('InnerProdGPU:clSetKernelArg:Wdict', ierr)

ierr = clSetKernelArg(kernel, 4, sizeof(cl_result), C_LOC(cl_result))
call CLerror_check('InnerProdGPU:clSetKernelArg:cl_result', ierr)

!execute the kernel
ierr = clEnqueueNDRangeKernel(command_queue, kernel, 2, C_NULL_PTR, C_LOC(globalsize), C_LOC(localsize), &
                              0, C_NULL_PTR, C_NULL_PTR)
call CLerror_check('InnerProdGPU:clEnqueueNDRangeKernel', ierr)

! wait for the commands to finish
ierr = clFinish(command_queue)
call CLerror_check('InnerProdGPU:clFinish', ierr)

! read the resulting vector from device memory
ierr = clEnqueueReadBuffer(command_queue,cl_result,CL_TRUE,0_8,size_in_bytes_result,C_LOC(results(1)),0,C_NULL_PTR,C_NULL_PTR)
call CLerror_check('InnerProdGPU:clEnqueueReadBuffer', ierr)

ierr = clReleaseMemObject(cl_result)
call CLerror_check('InnerProdGPU:clReleaseMemObject:cl_result', ierr)

! ---
end subroutine InnerProdGPU
!--------------------------------------------------------------------------

end module

max_stars_repo_name: hwagit/EMsoft

max_stars_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense

max_issues_repo_name: hwagit/EMsoft

max_issues_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense

max_forks_repo_name: hwagit/EMsoft

max_forks_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense


text:
C    Copyright (c) 2014, Sandia Corporation.
C    Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
C    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 Sandia Corporation 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

CC* FILE: [.QMESH]INNERH.FOR
CC* MODIFIED BY: TED BLACKER
CC* MODIFICATION DATE: 7/6/90
CC* MODIFICATION: COMPLETED HEADER INFORMATION
C
      SUBROUTINE INNERH (MXND, NXH, NUID, LXK, KXL, NXL, LXN, KKK, LLL,
     &   NNN, NNNOLD, NH, ISTART, IAVAIL, NAVAIL, NOROOM, ERR)
C***********************************************************************
C
C  SUBROUTINE INNERH  =  INSERT A ROW OF ELEMENTS AROUND A HOLE
C
C***********************************************************************
C
      DIMENSION NUID(MXND), LXK(4, MXND)
      DIMENSION KXL(2, 3*MXND), NXL(2, 3*MXND), LXN(4, MXND)
      DIMENSION NXH(MXND)
C
      LOGICAL ERR, NOROOM
C
      ERR = .FALSE.
      IF ((KKK + NH .GT. MXND) .OR. (LLL + 2*NH .GT. 3*MXND) .OR.
     &   (NNN + NH .GT. MXND)) THEN
         NOROOM = .TRUE.
         GO TO 200
      ELSE
         NOROOM = .FALSE.
      END IF
C
C  GENERATE LINES
C
      LLLOLD = LLL
      I1 = ISTART
      DO 100 I = 1, NH
         LLL = LLL + 1
         NXL(1, LLL) = NXH(I1)
         CALL ADDLXN(MXND, LXN, NUID, NAVAIL, IAVAIL, NXL(1, LLL),
     &      LLL, NNN, ERR, NOROOM)
         IF (NOROOM .OR. ERR) GO TO 200
         NXL(2, LLL) = NNNOLD + I
         I1 = I1 + 1
         IF (I1 .GT. NH) I1 = 1
  100 CONTINUE
C
      DO 110 I = 1, NH - 1
         LLL = LLL + 1
         NXL(1, LLL) = NNNOLD + I
         NXL(2, LLL) = NNNOLD + I + 1
  110 CONTINUE
      LLL = LLL + 1
      NXL(1, LLL) = NNN
      NXL(2, LLL) = NNNOLD + 1
C
C  MAKE SPACE IN LXN LIST FOR NEW NODES
C
      IF (IAVAIL .LT. NNNOLD) THEN
         JJ = IAVAIL
         MAXLNK = 0
  120    CONTINUE
         IF (LXN(4, JJ) .LE. NNNOLD) THEN
            JJ = LXN(4, JJ)
            MAXLNK = MAXLNK + 1
            IF (MAXLNK .GT. NNNOLD) THEN
               CALL MESAGE ('INNERH - LXN LINKED LIST ERROR')
               ERR = .TRUE.
               GO TO 200
            ELSE
               GO TO 120
            END IF
         END IF
         LXN(4, JJ) = NNNOLD + NH + 1
         NAVAIL = NAVAIL - NH
      ELSE IF (IAVAIL .GT. NNNOLD) THEN
         DO 130 J = 1, IAVAIL - 1
            IF (LXN(4, J) .LT. 0 .AND. ABS(LXN(4, J)) .GT. NNNOLD)
     &         LXN(4, J) = LXN(4, J) - NH
  130    CONTINUE
         JJ = NNNOLD + NH
         DO 150 J = NNNOLD + 1, IAVAIL - 1
            JJ = JJ + 1
            DO 140 K = 1, 4
               LXN(K, JJ) = LXN(K, J)
  140       CONTINUE
  150    CONTINUE
         IAVAIL = JJ + 1
         NAVAIL = NAVAIL - NH
      ELSE
         IAVAIL = NNN + 1
         NAVAIL = MXND - NNN
      END IF
C
C  MARK NODES ON HOLE BOUNDARY
C
      LXN(1, NNNOLD + 1) = LLLOLD + 1
      LXN(2, NNNOLD + 1) = -(LLLOLD + 2*NH)
      LXN(3, NNNOLD + 1) = LLLOLD + NH + 1
      LXN(4, NNNOLD + 1) = 0
      NNN = NNNOLD + 1
      DO 160 I = 2, NH
         NNN = NNN + 1
         LXN(1, NNN) = LLLOLD + I
         LXN(2, NNN) = -(LLLOLD + NH + I - 1)
         LXN(3, NNN) = LLLOLD + NH + I
         LXN(4, NNN) = 0
  160 CONTINUE
C
C  GENERATE ELEMENTS
C
      DO 180 I = LLLOLD + 1, LLL
         DO 170 J = 1, 2
            KXL(J, I) = 0
  170    CONTINUE
  180 CONTINUE
C
      I1 = ISTART
      I2 = I1 + 1
      IF (I2 .GT. NH) I2 = 1
      DO 190 I = 1, NH - 1
         CALL FNDLIN_FQ (MXND, LXN, NXH(I1), NXH(I2), LINE, ERR)
         IF (ERR) GO TO 200
         KKK = KKK + 1
         LXK(1, KKK) = LINE
         LXK(2, KKK) = LLLOLD + I
         LXK(3, KKK) = LLLOLD + I + 1
         LXK(4, KKK) = LLLOLD + I + NH
         I1 = I1 + 1
         IF (I1 .GT. NH) I1 = 1
         I2 = I1 + 1
         IF (I2 .GT. NH) I2 = 1
C
         IF (KXL(1, LINE) .EQ. 0) THEN
            KXL(1, LINE) = KKK
         ELSE IF (KXL(2, LINE) .EQ. 0) THEN
            KXL(2, LINE) = KKK
         ELSE
            CALL MESAGE ('KXL TABLE FULL')
         END IF
         KXL(2, LLLOLD + I) = KKK
         KXL(1, LLLOLD + I + 1) = KKK
         KXL(1, LLLOLD + I + NH) = KKK
  190 CONTINUE
      CALL FNDLIN_FQ (MXND, LXN, NXH(I1), NXH(I2), LINE, ERR)
      IF (ERR) GO TO 200
      KKK = KKK + 1
      LXK(1, KKK) = LINE
      LXK(2, KKK) = LLLOLD + NH
      LXK(3, KKK) = LLLOLD + 1
      LXK(4, KKK) = LLLOLD + 2*NH
C
      IF (KXL(1, LINE) .EQ. 0) THEN
         KXL(1, LINE) = KKK
      ELSE IF (KXL(2, LINE) .EQ. 0) THEN
         KXL(2, LINE) = KKK
      ELSE
         CALL MESAGE ('KXL TABLE FULL')
      END IF
      KXL(2, LLLOLD + NH) = KKK
      KXL(1, LLLOLD + 1) = KKK
      KXL(1, LLLOLD + 2*NH) = KKK
C
  200 CONTINUE
      RETURN
      END

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
c Obsolete
c
      subroutine tqrad(np,x,px,y,py,z,g,dv,sx,sy,sz,l,al,ak,
     1                 dx,dy,theta,cost,sint,radlvl,
     1                 f1in,f2in,f1out,f2out,mfring)
      use ffs_flag
      use tmacro
      implicit real*8(a-h,o-z)
      integer*4 np,l
      real*8 x(np),px(np),y(np),py(np),z(np),dv(np),g(np),
     $     sx(np),sy(np),sz(np)
      LOGICAL GAUSSR
      real*8 byx,bxa,bya,dprad,dpradx,dprady
      real*8 f1in,f2in,f1out,f2out
      include 'inc/TENT.inc'
      call ttfrin(np,x,px,y,py,z,g,4,ak,al,0.d0)
      al2=al*.5d0
      do 2010 i=1,np
        a=px(i)**2+py(i)**2
        dpz=a*(-.5d0-a*(.125d0+a*.0625d0))
        dpz=(dpz**2-a)/(2.d0+2.d0*dpz)
        dpz=(dpz**2-a)/(2.d0+2.d0*dpz)
        r=-dpz/(1.d0+dpz)*al2
        x(i)=x(i)+px(i)*r
        y(i)=y(i)+py(i)*r
        z(i)=z(i)-(3.d0+dpz)*a/2.d0/(2.d0+dpz)*r
2010  continue
      if(mfring .eq. 1 .or. mfring .eq. 3)then
        do 2110 i=1,np
c          p=(1.d0+g(i))**2
          p=(1.d0+g(i))
          a=f1in/p
          ea=exp(a)
          b=f2in/p
          pxf=px(i)/ea
          pyf=py(i)*ea
          z(i)=z(i)-(a*x(i)+b*(1.d0+.5d0*a)*pxf)*px(i)
     $             +(a*y(i)+b*(1.d0-.5d0*a)*pyf)*py(i)
          x(i)=ea*x(i)+b*px(i)
          y(i)=y(i)/ea-b*py(i)
          px(i)=pxf
          py(i)=pyf
2110    continue
      endif
      GAUSSR=RADLVL .EQ. 2.D0
      akk0=ak/al
      sqrtk=sqrt(abs(akk0))
      byx=akk0*brhoz
      b1=(byx)**2
      ur=urad*p0**3
      an=anrad*p0
      sp=0.d0
      do 100 i=1,np
        alsum=0.d0
        al0=al
c        p=g(i)*(2.d0+g(i))
        p=g(i)
        pr=1.d0+p
c        sqrtp=1.d0+g(i)
        sqrtp=sqrt(1.d0+g(i))
 110    brad0=x(i)**2+y(i)**2
        brad1=x(i)*px(i)+y(i)*py(i)
        brad2=(px(i)**2+py(i)**2-akk0*(x(i)-y(i))*(x(i)+y(i)))/3.d0
        brad=b1*(brad0+al0*(brad1+al0*brad2))
        if(brad .ne. 0.d0)then
          rho1=brhoz/sqrt(brad)
          alx=min(al0,.09d0*rho1/an)
111       brad=b1*(brad0+alx*(brad1+alx*brad2))
          rho=brhoz/sqrt(brad)
          alr=alx*(1.d0+(px(i)**2+py(i)**2)*.5d0)
          prob=alr*an/rho
          if(prob .gt. .1d0)then
            alx=alx*(.09d0/prob)
            go to 111
          endif
          if(tran() .gt. 1.d0-prob)then
            bxa=byx*(y(i)*(1.d0+akk0*alx**2/6.d0)+py(i)*alx*.5d0)
            bya=byx*(x(i)*(1.d0-akk0*alx**2/6.d0)+px(i)*alx*.5d0)
            call tsynchrad(pr,alr,bxa,bya,
     $           dprad,dpradx,dprady,
     $           i,l,alsum,0.d0,0.d0,theta,
     $           x(i),y(i),px(i),py(i))
            p=max(-.999d0,p-dprad*.5d0)
            px(i)=px(i)-dpradx*.5d0
            py(i)=py(i)-dprady*.5d0
            pr=1.d0+p
            sp=sp-dprad
            sqrtp=sqrt(pr)
          else
            dprad=0.d0
          endif
        else
          alx=al0
          dprad=0.d0
        endif
        alsum=alsum+alx
        if(ak .gt. 0.d0)then
          akk=sqrtk/sqrtp
          phi=akk*alx*.5d0
          s=phi**2
          t=1.d0-s/(3.d0-s/(5.d0-s/(7.d0-s/9.d0)))
          th=1.d0+s/(3.d0+s/(5.d0+s/(7.d0+s/9.d0)))
          u=t**2+phi**2
          a11=(t-phi)*(t+phi)/u
          a12=2.d0*t*phi/u
          u=(th-phi)*(th+phi)
          b11=(th**2+phi**2)/u
          b12=2.d0*th*phi/u
          a21=-a12*akk
          a12=a12/akk
          b21=b12*akk
          b12=b12/akk
c    Here I ignore the change of dv(i)
          xi=x(i)
          yi=y(i)
          x(i)=a11*xi+a12*px(i)
          y(i)=b11*yi+b12*py(i)
          z(i)=z(i)-dv(i)*alx-
     1       (((px(i)**2+(xi*akk)**2)*(alx+a11*a12)+
     1         (py(i)-yi*akk)*(py(i)+yi*akk)*(alx+b11*b12))*.5d0+
     1         xi*x(i)*a21+yi*y(i)*b21)*.5d0
          px(i)=a21*xi+a11*px(i)
          py(i)=b21*yi+b11*py(i)
        elseif(ak .lt. 0.d0)then
          akk=sqrtk/sqrtp
          phi=akk*alx
          phi=akk*alx*.5d0
          s=phi**2
          t=1.d0-s/(3.d0-s/(5.d0-s/(7.d0-s/9.d0)))
          th=1.d0+s/(3.d0+s/(5.d0+s/(7.d0+s/9.d0)))
          u=t**2+phi**2
          b11=(t-phi)*(t+phi)/u
          b12=2.d0*t*phi/u
          u=(th-phi)*(th+phi)
          a11=(th**2+phi**2)/u
          a12=2.d0*th*phi/u
          b21=-b12*akk
          b12=b12/akk
          a21=a12*akk
          a12=a12/akk
          xi=x(i)
          yi=y(i)
          x(i)=xi*a11+px(i)*a12
          y(i)=yi*b11+py(i)*b12
          z(i)=z(i)-dv(i)*alx-
     1       (((px(i)-xi*akk)*(px(i)+xi*akk)*(alx+a11*a12)+
     1         (py(i)**2+(yi*akk)**2)*(alx+b11*b12))*.5d0+
     1         xi*x(i)*a21+yi*y(i)*b21)*.5d0
          px(i)=a21*xi+a11*px(i)
          py(i)=b21*yi+b11*py(i)
        endif
        if(dprad .ne. 0.d0)then
          p=max(-.999d0,p-dprad*.5d0)
          px(i)=px(i)-dpradx*.5d0
          py(i)=py(i)-dprady*.5d0
          pr=1.d0+p
          sqrtp=sqrt(pr)
        endif
        al0=al0-alx
        if(al0 .gt. 0.d0)then
          go to 110
        endif
        g(i)=p
100   continue
      if(radcod)then
        sp=0.d0
      else
        sp=sp/np
      endif
      do 310 i=1,np
        dp=max(-.99d0,g(i)-sp)
        pr=1.d0+dp
        r=(1.d0+g(i))/pr
c       px(i)=px(i)*r
c       py(i)=py(i)*r
c        g(i)=dp/(1.d0+sqrt(pr))
        g(i)=dp
        p1=pr*p0
        h=sqrt(1.d0+p1**2)
        dv(i)=-dp*(1.d0+pr)/h/(h+pr*h0)+dvfs
310   continue
      if(mfring .eq. 2 .or. mfring .eq. 3)then
        do 2120 i=1,np
c          p=(1.d0+g(i))**2
          p=(1.d0+g(i))
          a=f1out/p
          ea=exp(a)
          b=f2out/p
          pxf=px(i)/ea
          pyf=py(i)*ea
          z(i)=z(i)-(a*x(i)+b*(1.d0+.5d0*a)*pxf)*px(i)
     $             +(a*y(i)+b*(1.d0-.5d0*a)*pyf)*py(i)
          x(i)=ea*x(i)+b*px(i)
          y(i)=y(i)/ea-b*py(i)
          px(i)=pxf
          py(i)=pyf
2120    continue
      endif
      do 2020 i=1,np
        a=px(i)**2+py(i)**2
        dpz=a*(-.5d0-a*(.125d0+a*.0625d0))
        dpz=(dpz**2-a)/(2.d0+2.d0*dpz)
        dpz=(dpz**2-a)/(2.d0+2.d0*dpz)
        r=-dpz/(1.d0+dpz)*al2
        x(i)=x(i)+px(i)*r
        y(i)=y(i)+py(i)*r
        z(i)=z(i)-(3.d0+dpz)*a/2.d0/(2.d0+dpz)*r
2020  continue
      call ttfrin(np,x,px,y,py,z,g,4,-ak,al,0.d0)
      include 'inc/TEXIT.inc'
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


text:
In WikiWikiWeb other Wikis, like C2, insults are not appreciated.  Here, feel free  if that place down the street gave you food poisoning, put how you felt and exactly how long you wanted to kill the cook as you hunkered down over the toilet. (Make sure you understand the process of food poisoning before leveling any charges however.) Profanity is discouraged, but there are plenty of insults that can be made without them.

Theres a difference between insults that are merely insults and a heated response.  I dont like unnecessarily derisive and insulting language, because it is usually written in haste.  What might be considered an insult to someone may end up useful to someone else.  Being libel careful with your words is important.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE FAILRS (FTHR,ULTSTN,STRESL,FINDEX)
C
C     THIS ROUTINE COMPUTES THE FAILURE INDEX OF A LAYER IN A LAMINATED
C     COMPOSITE ELEMENT USING ONE OF THE FOLLOWING FIVE FAILURE THEORIES
C     CURRENTLY AVAILABLE
C        1.   HILL
C        2.   HOFFMAN
C        3.   TSAI-WU
C        4.   MAX STRESS
C        5.   MAX STRAIN
C
C
C     DEFINITIONS
C
C     XT = ULTIMATE UNIAXIAL TENSILE STRENGTH IN THE FIBER DIRECTION
C     XC = ULTIMATE UNIAXIAL COMPRESSIVE STRENGTH IN THE FIBER DIRECTION
C     YT = ULTIMATE UNIAXIAL TENSILE STRENGTH PERPENDICULAR TO THE FIBER
C          DIRECTION
C     YC = ULTIMATE UNIAXIAL COMPRESSIVE STRENGTH PERPENDICULAR TO THE
C          FIBER DIRECTION
C     S  = ULTIMATE PLANAR SHEAR STRENGTH UNDER PURE SHEAR LOADING
C
C     SIMILARILY FOR THE ULTIMATE STRAINS
C
C
      INTEGER FTHR
      REAL    ULTSTN(6),STRESL(3)
C
C
C     CHECK FOR ZERO STRENGTH VALUES
C
      DO 50 I = 1,5
   50 IF (ULTSTN(I) .EQ. 0.0) GO TO 700
C
C     ULTIMATE STRENGTH VALUES
C
      XT   = ULTSTN(1)
      XC   = ULTSTN(2)
      YT   = ULTSTN(3)
      YC   = ULTSTN(4)
      S    = ULTSTN(5)
      F12  = ULTSTN(6)
C
C     LAYER STRESSES
C
      SIG1 = STRESL(1)
      SIG2 = STRESL(2)
      TAU12= STRESL(3)
C
C     LAYER STRAINS
C
      EPS1 = STRESL(1)
      EPS2 = STRESL(2)
      GAMA = STRESL(3)
C
      GO TO (100,200,300,400,500), FTHR
C
C     HILL FAILURE THEORY
C     -------------------
C
  100 X = XT
      IF (SIG1 .LT. 0.0) X = XC
C
      Y = YT
      IF (SIG2 .LT. 0.0) Y = YC
C
      XX = XT
      IF (SIG1*SIG2 .LT. 0.0) XX = XC
C
      FINDEX = (SIG1*SIG1)/(X*X)
      FINDEX = FINDEX + (SIG2 * SIG2)/(Y * Y)
      FINDEX = FINDEX - (SIG1 * SIG2)/(XX*XX)
      FINDEX = FINDEX + (TAU12*TAU12)/(S * S)
      GO TO 600
C
C
C     HOFFMAN FAILURE THEORY
C     ----------------------
C
  200 FINDEX = (1.0/XT-1.0/XC)*SIG1
      FINDEX = FINDEX + (1.0/YT-1.0/YC)*SIG2
      FINDEX = FINDEX + (SIG1 * SIG1)/(XT*XC)
      FINDEX = FINDEX + (SIG2 * SIG2)/(YT*YC)
      FINDEX = FINDEX + (TAU12*TAU12)/(S * S)
      FINDEX = FINDEX + (SIG1 * SIG2)/(XT*XC)
      GO TO 600
C
C
C     TSAI-WU FAILURE THEORY
C     ----------------------
C
C     CHECK STABILITY CRITERIA FOR THE INTERACTION TERM F12
C
  300 IF (F12 .EQ. 0.0) GO TO 350
C
      CRIT = (1.0/(XT*XC))*(1.0/(YT*YC)) - F12*F12
      IF (CRIT .GT. 0.0) GO TO 350
C
C     IF STABILITY CRITERIA IS VIOLATED THEN SET THE F12 THE INTERACTION
C     TERM TO ZERO
C
      F12 = 0.0
C
C
  350 FINDEX = (1.0/XT-1.0/XC)*SIG1
      FINDEX = FINDEX + (1.0/YT-1.0/YC)*SIG2
      FINDEX = FINDEX + (SIG1 * SIG1)/(XT*XC)
      FINDEX = FINDEX + (SIG2 * SIG2)/(YT*YC)
      FINDEX = FINDEX + (TAU12*TAU12)/(S * S)
      IF (F12 .EQ. 0.0) GO TO 600
      FINDEX = FINDEX + 2.0*F12*SIG1*SIG2
      GO TO 600
C
C
C     MAX STRESS FAILURE THEORY
C     -------------------------
C
  400 FI1 = SIG1/XT
      IF (SIG1 .LT. 0.0) FI1 = ABS(SIG1/XC)
C
      FI2 = SIG2/YT
      IF (SIG2 .LT. 0.0) FI2 = ABS(SIG2/YC)
C
      FI12 = ABS(TAU12)/S
C
      FINDEX = FI1
      IF (FI2  .GT. FINDEX) FINDEX = FI2
      IF (FI12 .GT. FINDEX) FINDEX = FI12
      GO TO 600
C
C
C     MAX STRAIN FAILURE THEORY
C     -------------------------
C
  500 FI1 = EPS1/XT
      IF (EPS1 .LT. 0.0) FI1 = ABS(EPS1/XC)
C
      FI2 = EPS2/YT
      IF (EPS2 .LT. 0.0) FI2 = ABS(EPS2/YC)
C
      FI12 = ABS(GAMA)/S
C
      FINDEX = FI1
      IF (FI2  .GT. FINDEX) FINDEX = FI2
      IF (FI12 .GT. FINDEX) FINDEX = FI12
C
  600 CONTINUE
C
      RETURN
C
C
C     NON-FATAL ERROR
C
  700 FINDEX = 0.0
      RETURN
      END


! ###############################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! 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 and documentation.

! 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.
! _______________________________________________________________________________________________________

! End MIT license text.

   MODULE TRANSFORM_NODE_FORCES_Interface

   INTERFACE

      SUBROUTINE TRANSFORM_NODE_FORCES ( COORD_SYS )

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, MELGP, NCORD
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  TRANSFORM_NODE_FORCES_BEGEND
      USE CONSTANTS_1, ONLY           :  ZERO
      USE MODEL_STUF, ONLY            :  CAN_ELEM_TYPE_OFFSET, GRID, CORD, BGRID, ELDOF, ELGP, OFFDIS, OFFSET, PEB, PEG, PEL, TE,  &
                                         TYPE

      IMPLICIT NONE

      CHARACTER(LEN=*), INTENT(IN)    :: COORD_SYS         ! 'B" for basic, 'G' for global

      INTEGER(LONG), PARAMETER        :: NCOLA     = 3     ! An input to subr MATMULT_FFF called herein
      INTEGER(LONG), PARAMETER        :: NCOLB     = 1     ! An input to subr MATMULT_FFF called herein
      INTEGER(LONG), PARAMETER        :: NROWA     = 3     ! An input to subr MATMULT_FFF called herein
      INTEGER(LONG), PARAMETER        :: NROW      = 3     ! An input to subr MATPUT, MATGET called herein
      INTEGER(LONG), PARAMETER        :: NCOL      = 1     ! An input to subr MATPUT, MATGET called herein
      INTEGER(LONG), PARAMETER        :: PCOL      = 1     ! An input to subr MATPUT, MATGET called herein
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = TRANSFORM_NODE_FORCES_BEGEND

      REAL(DOUBLE)                    :: THETAD,PHID       ! Returns from subr GEN_T0L (not used here)

      END SUBROUTINE TRANSFORM_NODE_FORCES

   END INTERFACE

   END MODULE TRANSFORM_NODE_FORCES_Interface