text:
c Interface to LHAPDF, v4.1 onwards. The inclusion of LHAGLUE allows the
c use of PDFLIB format. This file has been derived from mcatnlo_mlmtopdf.f.
c Calls to newmode and pdftomlm are trivial here, but are kept for easier
c future changes
c
c Modified on 16/1/2008: implements an alternative strategy to get Lambda_QCD;
c does not use LHAPDF common blocks, but calls alpha_S(M_Z) and computes
c Lambda_QCD. The old strategy is still available (see the parameter ilam)
c
      subroutine mlmpdf(mode,ih,q2,x,fx,nf)
      implicit real * 8 (a-h,o-z)
      real * 4 q2,x,fx(-nf:nf)
      common/trans/nptype,ngroup,nset
c x*pdf(x) in PDG format
      dimension fxp(-6:6)
c used by pdfset
      dimension val(20)
c used by pdfset
      character * 20 parm(20)
c used by structp
      parameter (zero=0.d0)
      parameter (izero=0)
      logical ini
      data ini/.true./,imode/0/
c
      if(x.le.0.or.x.ge.1) then
        do j=-nf,nf
          fx(j) = 0
        enddo
        return
      endif
      if(ih.ne.5)then
c incoming particle is not an electron: use PDFLIB
        if(ini.or.imode.ne.mode) then
          ini = .false.
          imode = mode
c pass from our conventions to PDFLIB conventions for parameter settings.
c See subroutine newmode for comments on the conventions adopted. Proper
c settings of /LHAPDFC/ and /LHACONTROL/ must have been adopted before
c runtime (and before the call to PDFSET)
          call newmode(ih,mode)
          parm(1) = 'DEFAULT'
          val (1) =  nset
          call pdfset(parm,val)
        endif
        xd = dble(x)
        qd = dble(sqrt(q2))
c first convert structm format into pftopdg format (fxp), in order to use
c what already done in mcatnlo_mlmtopdf.f where pftopdg is available. For
c photons, use structp and consider only on-shell photons
        if(abs(ih).le.4)then
          if(abs(ih).le.3)then
            call structm(xd,qd,upv,dnv,usea,dsea,str,chm,bot,top,glu)
          elseif(ih.eq.4)then
            q2d = dble(q2)
            call structp(xd,q2d,zero,izero,
     #                   upv,dnv,usea,dsea,str,chm,bot,top,glu)
          else
            write(*,*)'Unknow particles in LHAPDF'
            stop
          endif
c If a nucleon, define fxp assuming LHAPDF returns proton densities, i.e.
c up=upv+usea, upbar=usea, dn=dnv+dsea, dnbar=dsea. If a pion or photon,
c use the SAME conventions (which are wrong) according to PDFLIB manual.
c Correct below in order to obtain the physical densities from which
c to calculate the luminosities
          fxp( 0)=glu
          fxp( 1)=dnv+dsea
          fxp( 2)=upv+usea
          fxp( 3)=str
          fxp( 4)=chm
          fxp( 5)=bot
          fxp( 6)=top
          fxp(-1)=dsea
          fxp(-2)=usea
          fxp(-3)=str
          fxp(-4)=chm
          fxp(-5)=bot
          fxp(-6)=top
        else
          write(*,*)'Unknow particles in LHAPDF'
          stop
        endif
c in our conventions 1=up, 2=down; PDFLIB 1=down, 2=up. With
c f(1)<-->f(2) we mean also f(-1)<-->f(-2)
c in the following lines, deals with particles only (no antiparticles)
c proton(ih=1) ==> f(1)<-->f(2)
c neutron(ih=2) ==> no action (f(1)<-->f(2) for PDFLIB convention and
c    f(1)<-->f(2) for isospin symmetry (u_proton=d_neutron....)
c pion+(ih=3) ==> f(2)<-->f(-2), since PDFLIB has d=u=q_v+q_sea,
c    ubar=dbar=q_sea
c photon(ih=4) ==> f(-1)<-->f(-2) and f(i)=f(-i)/2, i=1,2 since PDFLIB
c    has f(i)=2*f(-i), and f(1)<-->f(2)
c Notice that in the jet package pions and neutrons are not used. If
c selected, they are rejected by the routine pdfpar. This routine
c is therefore a completely general interface with PDFLIB
        if(abs(ih).eq.1) then
          tmp     = fxp(1)
          fxp(1)  = fxp(2)
          fxp(2)  = tmp
          tmp     = fxp(-1)
          fxp(-1) = fxp(-2)
          fxp(-2) = tmp
        elseif(abs(ih).eq.2) then
          continue
        elseif(abs(ih).eq.3) then
          tmp = fxp(-2)
          fxp(-2) = fxp(2)
          fxp(2)  = tmp
        elseif(abs(ih).eq.4) then
          tmp     = fxp(-1)
          fxp(-1) = fxp(-2)
          fxp(-2) = tmp
          fxp(1)=fxp(-1)
          fxp(2)=fxp(-2)
c this is (p+n)/2
        elseif(ih.eq.0) then
          va  = (fxp(1)+fxp(2))/2
          sea = (fxp(-1)+fxp(-2))/2
          fxp(1)  = va
          fxp(2)  = va
          fxp(-1) = sea
          fxp(-2) = sea
        else
          write(*,*) ' ih was', ih,' instead of 0, +-1, +-2, +-3, or 4'
          stop
        endif
c for particles, ich=1, for antiparticles, ich=-1
        if(ih.lt.0) then
          ich = -1
        else
          ich = 1
        endif
c divide by x and exchange q with qbar in the case of antiparticles
        do j=-nf,nf
          fx(j) = sngl(fxp(j*ich)/xd)
        enddo
      else
c incoming "hadron" is an electron
        if(ini.or.imode.ne.mode) then
          ini = .false.
          imode = mode
          iset=mode-50
        endif
        if(iset.eq.1) then
          call elpdf_lac1(q2,x,fx,nf)
        elseif(iset.eq.2) then
          call elpdf_grv(q2,x,fx,nf)
        elseif(iset.eq.3) then
          call elpdf_user(q2,x,fx,nf)
        else
          write(*,*)'Electron set non implemented'
          stop
        endif
      endif
      return
      end

      subroutine pdfpar(mode,ih,xlam,scheme,iret)
      implicit real * 8 (a-h,o-z)
      parameter (xmz=91.1876d0)
      parameter (ifive=5)
c If ilam=0, obtains Lambda_5 from LHAPDF common blocks
c If ilam=1, computes Lambda_5 using as(MZ) as computed by LHAPDF
      parameter (ilam=1)
      common/trans/nptype,ngroup,nset
      common/w50512/qcdl4,qcdl5
      dimension val(20)
      character * 20 parm(20)
      character * 2 scheme
      logical ini
      data ini/.true./
c iret#0 when problem occur
      iret = 0
      if(ini)then
        write(*,*)'This is an interface to LHAPDF'
        ini = .false.
      endif
c modify the following if non-nucleon PDFs will be included in LHAPDF
      if(abs(ih).gt.5)then
        write(*,*) ' Hadron tpye ',ih,' not implemented'
        iret=1
        return
      endif
c fake values. If kept, the main program crashes
      scheme='XX'
      xlam=0.0
c incoming particle is not an electron: use LHAPDF
      if(ih.ne.5)then
c the scheme of the LHAPDF set is not given in any common block.
c Set it by hand in the main program
        scheme = '**'
        call newmode(ih,mode)
        parm(1) = 'DEFAULT'
        val (1) =  nset
c set the parameters
        call pdfset(parm,val)
        if(ilam.eq.0)then
c Lambda_QCD_5, as given by LHAPDF common blocks; this may be wrong,
c so be careful
          xlam = qcdl5
        elseif(ilam.eq.1)then
          asmzlha=alphasPDF(xmz)
          xlam=fk88xlambd(asmzlha,xmz,ifive)
        else
          write(*,*)'Error in pdfpar: unknown ilam',ilam
          stop
        endif
      else
c incoming particle is an electron
        if(mode.eq.51)then
          scheme='MS'
          xlam=.130
        elseif(mode.eq.52)then
          scheme='DG'
          xlam=.130
        elseif(mode.eq.53)then
          scheme='**'
          xlam=.001
        else
          write(*,*)'Electron set not implemented'
          iret=1
        endif
      endif
      return
      end

      subroutine getlamlha(xlam,xlamlha)
c This is the part of pdfpar which returns Lambda_QCD as used in MC@NLO
c code (xlam) using the information provided by LHAPDF. A separate
c routine is necessary since pdfpar is not called by MC@NLO codes
c linked with LHAPDF
      implicit none
      real*8 xlam,xlamlha,xmz,asmzlha,alphasPDF,fk88xlambd
      parameter (xmz=91.1876d0)
      integer ifive,ilam
      parameter (ifive=5)
c If ilam=0, obtains Lambda_5 from LHAPDF common blocks
c If ilam=1, computes Lambda_5 using as(MZ) as computed by LHAPDF
      parameter (ilam=1)
c
      if(ilam.eq.0.and.xlamlha.le.0.d0)then
        write(*,*)'Error in getlamlha:',xlamlha
        write(*,*)' This function must be called after PDFSET'
        stop
      endif
c
      if(ilam.eq.0)then
c Lambda_QCD_5, as given by LHAPDF common blocks; this may be wrong,
c so be careful
        xlam = xlamlha
      elseif(ilam.eq.1)then
        asmzlha=alphasPDF(xmz)
        xlam=fk88xlambd(asmzlha,xmz,ifive)
      else
        write(*,*)'Error in pdfpar: unknown ilam',ilam
        stop
      endif
c
      return
      end

      subroutine prntsf
c     prints details of the structure function sets
c
      write(*,100)
     #  '  Refer to LHAPDF on-line manual to see'
     # ,'  PDFLIB-compatible labeling conventions'
 100  format(1x,a,100(/,1x,a))
      return
      end

      subroutine newmode(ih,mode)
c Trivial except for ih in the case of LHAPDF, since PDFLIB-type set
c numbers are specified without using a group number; this format is
c therefore compatible with mlm's one. For particle types (LHAPDF v4.1
c doesn't seem to have a particle type classification -- retain here
c that of PDFLIB, to allow more flexibility for the future).
c
c                        MC@NLO               PDFLIB
c
c  nucleons           -2,-1,0,1,2                1
c  pions                  -3,3                   2
c  photons                  4                    3
c
      implicit real * 8 (a-h,o-z)
      common/trans/nptype,ngroup,nset
c
      if(abs(ih).le.2)then
        nptype=1
      elseif(abs(ih).eq.3)then
        nptype=2
      elseif(ih.eq.4)then
        nptype=3
      else
        write(6,*)'Hadron type not implemented in PDFLIB/LHAPDF'
        stop
      endif
      ngroup=-1
      nset=mode
      return
      end

      subroutine pdftomlm(ipdfih,ipdfgroup,ipdfndns,ihmlm,ndnsmlm)
c Performs the inverse operation of newmode. From LHAPDFv4.0, proton,
c pions, and photons are included; the common block w50511 is not
c filled by LHAGLUE, and thus the information on particle type is not
c available in versions up to v4.1; keep the following for more flexibility
      implicit real * 8 (a-h,o-z)
c
      if(ipdfih.eq.1)then
        ihmlm=1
      elseif(ipdfih.eq.2)then
        ihmlm=3
      elseif(ipdfih.eq.3)then
        ihmlm=4
      else
        write(*,*)'Wrong hadron type in LHAPDF:',ipdfih
        stop
      endif
      ndnsmlm=ipdfndns
      return
      end

      subroutine setlhacblk(strin)
c Must be called BEFORE SETPAR. It sets the variables that control LHAPDF,
c stored in the common block /LHACONTROL/, taken from LHAGLUE. /LHAPDFC/
c can also be set here if needed (in this version, use the default and
c define logical link at runtime instead)
      implicit none
      character*70 strin,strout
      character*20 lhaparm(20)
      double precision lhavalue(20)
      common/lhacontrol/lhaparm,lhavalue
c
c change the following to collect statistics for PDFs outside range
      lhaparm(16)='NOSTAT'
c 'EXTRAPOLATE' will extrapolate PDFs outside their range
      call fk88low_to_upp(strin,strout)
      lhaparm(18)=strout
c set to 'SILENT' to suppress output completely
      lhaparm(19)='LOWKEY'
c
      return
      end

C------- ALPHA QCD -------------------------------------
c Program to calculate alfa strong with nf flavours,
c as a function of lambda with 5 flavors.
c The value of alfa is matched at the thresholds q = mq.
c When invoked with nf < 0 it chooses nf as the number of
c flavors with mass less then q.
c
      function alfas(q2,xlam,inf)
      implicit real * 8 (a-h,o-z)
      data olam/0.d0/,pi/3.14159d0/
      data xmb/5.d0/,xmc/1.5d0/
      if(xlam.ne.olam) then
        olam = xlam
        b5  = (33-2*5)/pi/12
        bp5 = (153 - 19*5) / pi / 2 / (33 - 2*5)
        b4  = (33-2*4)/pi/12
        bp4 = (153 - 19*4) / pi / 2 / (33 - 2*4)
        b3  = (33-2*3)/pi/12
        bp3 = (153 - 19*3) / pi / 2 / (33 - 2*3)
        xlc = 2 * log(xmc/xlam)
        xlb = 2 * log(xmb/xlam)
        xllc = log(xlc)
        xllb = log(xlb)
        c45  =  1/( 1/(b5 * xlb) - xllb*bp5/(b5 * xlb)**2 )
     #        - 1/( 1/(b4 * xlb) - xllb*bp4/(b4 * xlb)**2 )
        c35  =  1/( 1/(b4 * xlc) - xllc*bp4/(b4 * xlc)**2 )
     #        - 1/( 1/(b3 * xlc) - xllc*bp3/(b3 * xlc)**2 ) + c45
      endif
      q   = sqrt(q2)
      xlq = 2 * log( q/xlam )
      xllq = log( xlq )
      nf = inf
      if( nf .lt. 0) then
        if( q .gt. xmb ) then
          nf = 5
        elseif( q .gt. xmc ) then
          nf = 4
        else
          nf = 3
        endif
      endif
      if    ( nf .eq. 5 ) then
        alfas = 1/(b5 * xlq) -  bp5/(b5 * xlq)**2 * xllq
      elseif( nf .eq. 4 ) then
        alfas = 1/( 1/(1/(b4 * xlq) - bp4/(b4 * xlq)**2 * xllq) + c45 )
      elseif( nf .eq. 3 ) then
        alfas = 1/( 1/(1/(b3 * xlq) - bp3/(b3 * xlq)**2 * xllq) + c35 )
      else
        print *,'error in alfa: unimplemented # of light flavours',nf
        stop
      endif
      return
      end

      function fk88xlambd(as,q,nf)
      implicit real*8(a-h,o-z)
      data pi/3.1415926535897931d0/
      b  = (33-2*nf)/pi/12
      bp = (153 - 19*nf) / pi / 2 / (33 - 2*nf)
      t  = 1/b/as
    1 xlt = log(t)
      ot = t
c-----------------------------------------------------------
c Value and Derivative of alfa with respect to t
      as0  = 1/b/t - bp*xlt/(b*t)**2
      as1  = - 1/b/t**2 -bp/b**2*(1-2*xlt)/t**3
      t  = (as-as0)/as1 + t
      if(abs(ot-t)/ot.gt.1.d-5)goto 1
      fk88xlambd = q/exp(t/2)
      return
      end

max_stars_repo_name: JaySandesara/madminer

max_issues_repo_name: JaySandesara/madminer

max_forks_repo_name: JaySandesara/madminer


!*==sgttrf.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b SGTTRF
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download SGTTRF + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE SGTTRF( N, DL, D, DU, DU2, IPIV, INFO )
!
!       .. Scalar Arguments ..
!       INTEGER            INFO, N
!       ..
!       .. Array Arguments ..
!       INTEGER            IPIV( * )
!       REAL               D( * ), DL( * ), DU( * ), DU2( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> SGTTRF computes an LU factorization of a real tridiagonal matrix A
!> using elimination with partial pivoting and row interchanges.
!>
!> The factorization has the form
!>    A = L * U
!> where L is a product of permutation and unit lower bidiagonal
!> matrices and U is upper triangular with nonzeros in only the main
!> diagonal and first two superdiagonals.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the matrix A.
!> \endverbatim
!>
!> \param[in,out] DL
!> \verbatim
!>          DL is REAL array, dimension (N-1)
!>          On entry, DL must contain the (n-1) sub-diagonal elements of
!>          A.
!>
!>          On exit, DL is overwritten by the (n-1) multipliers that
!>          define the matrix L from the LU factorization of A.
!> \endverbatim
!>
!> \param[in,out] D
!> \verbatim
!>          D is REAL array, dimension (N)
!>          On entry, D must contain the diagonal elements of A.
!>
!>          On exit, D is overwritten by the n diagonal elements of the
!>          upper triangular matrix U from the LU factorization of A.
!> \endverbatim
!>
!> \param[in,out] DU
!> \verbatim
!>          DU is REAL array, dimension (N-1)
!>          On entry, DU must contain the (n-1) super-diagonal elements
!>          of A.
!>
!>          On exit, DU is overwritten by the (n-1) elements of the first
!>          super-diagonal of U.
!> \endverbatim
!>
!> \param[out] DU2
!> \verbatim
!>          DU2 is REAL array, dimension (N-2)
!>          On exit, DU2 is overwritten by the (n-2) elements of the
!>          second super-diagonal of U.
!> \endverbatim
!>
!> \param[out] IPIV
!> \verbatim
!>          IPIV is INTEGER array, dimension (N)
!>          The pivot indices; for 1 <= i <= n, row i of the matrix was
!>          interchanged with row IPIV(i).  IPIV(i) will always be either
!>          i or i+1; IPIV(i) = i indicates a row interchange was not
!>          required.
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -k, the k-th argument had an illegal value
!>          > 0:  if INFO = k, U(k,k) is exactly zero. The factorization
!>                has been completed, but the factor U is exactly
!>                singular, and division by zero will occur if it is used
!>                to solve a system of equations.
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup realGTcomputational
!
!  =====================================================================
      SUBROUTINE SGTTRF(N,Dl,D,Du,Du2,Ipiv,Info)
      IMPLICIT NONE
!*--SGTTRF128
!
!  -- LAPACK computational routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      INTEGER Info , N
!     ..
!     .. Array Arguments ..
      INTEGER Ipiv(*)
      REAL D(*) , Dl(*) , Du(*) , Du2(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      REAL ZERO
      PARAMETER (ZERO=0.0E+0)
!     ..
!     .. Local Scalars ..
      INTEGER i
      REAL fact , temp
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS
!     ..
!     .. External Subroutines ..
      EXTERNAL XERBLA
!     ..
!     .. Executable Statements ..
!
      Info = 0
      IF ( N<0 ) THEN
         Info = -1
         CALL XERBLA('SGTTRF',-Info)
         RETURN
      ENDIF
!
!     Quick return if possible
!
      IF ( N==0 ) RETURN
!
!     Initialize IPIV(i) = i and DU2(I) = 0
!
      DO i = 1 , N
         Ipiv(i) = i
      ENDDO
      DO i = 1 , N - 2
         Du2(i) = ZERO
      ENDDO
!
      DO i = 1 , N - 2
         IF ( ABS(D(i))<ABS(Dl(i)) ) THEN
!
!           Interchange rows I and I+1, eliminate DL(I)
!
            fact = D(i)/Dl(i)
            D(i) = Dl(i)
            Dl(i) = fact
            temp = Du(i)
            Du(i) = D(i+1)
            D(i+1) = temp - fact*D(i+1)
            Du2(i) = Du(i+1)
            Du(i+1) = -fact*Du(i+1)
            Ipiv(i) = i + 1
!
!           No row interchange required, eliminate DL(I)
!
         ELSEIF ( D(i)/=ZERO ) THEN
            fact = Dl(i)/D(i)
            Dl(i) = fact
            D(i+1) = D(i+1) - fact*Du(i)
         ENDIF
      ENDDO
      IF ( N>1 ) THEN
         i = N - 1
         IF ( ABS(D(i))<ABS(Dl(i)) ) THEN
            fact = D(i)/Dl(i)
            D(i) = Dl(i)
            Dl(i) = fact
            temp = Du(i)
            Du(i) = D(i+1)
            D(i+1) = temp - fact*D(i+1)
            Ipiv(i) = i + 1
         ELSEIF ( D(i)/=ZERO ) THEN
            fact = Dl(i)/D(i)
            Dl(i) = fact
            D(i+1) = D(i+1) - fact*Du(i)
         ENDIF
      ENDIF
!
!     Check for a zero on the diagonal of U.
!
      DO i = 1 , N
         IF ( D(i)==ZERO ) THEN
            Info = i
            EXIT
         ENDIF
      ENDDO
!
!
!     End of SGTTRF
!
      END SUBROUTINE SGTTRF