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

      SUBROUTINE BD_CTETRA ( CARD, LARGE_FLD_INP, NUM_GRD )

! Processes CTETRA Bulk Data Cards
!  1) Sets ETYPE for this element type
!  2) Calls subr ELEPRO to read element ID, property ID and connection data into array EDAT

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, JCARD_LEN, FATAL_ERR, NCTETRA4, NCTETRA10, NEDAT, NELE
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  BD_CTETRA_BEGEND
      USE MODEL_STUF, ONLY            :  ETYPE

      USE BD_CTETRA_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'BD_CTETRA'
      CHARACTER(LEN=*), INTENT(INOUT) :: CARD              ! A Bulk Data card
      CHARACTER(LEN=*), INTENT(IN)    :: LARGE_FLD_INP     ! If 'Y', CARD is large field format
      CHARACTER(LEN(CARD))            :: CHILD             ! "Child" card read in subr NEXTC, called herein
      CHARACTER(LEN=JCARD_LEN)        :: JCARD(10)         ! The 10 fields of characters making up CARD
      CHARACTER(LEN(JCARD))           :: ID                ! Character value of element ID (field 2 of parent card)
      CHARACTER(LEN(JCARD))           :: JCARD_EDAT(10)    ! JCARD values sent to subr ELEPRO
      CHARACTER(LEN(JCARD))           :: NAME              ! JCARD(1) from parent entry

      INTEGER(LONG), INTENT(OUT)      :: NUM_GRD           ! Number of GRID's + SPOINT's for the elem
      INTEGER(LONG)                   :: ICONT     = 0     ! Indicator of whether a cont card exists. Output from subr NEXTC
      INTEGER(LONG)                   :: IERR      = 0     ! Error indicator returned from subr NEXTC called herein
      INTEGER(LONG)                   :: I                 ! DO loop index
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = BD_CTETRA_BEGEND

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
! CTETRA element Bulk Data Card routine

!   FIELD   ITEM                   ARRAY ELEMENT
!   -----   ------------   ---------------------------------

! Mandatory parent card

!    1      Element type   ETYPE(nele) = 'TETRA4' or 'TETRA10'
!    2      Element ID     EDAT(nedat+1)
!    3      Property ID    EDAT(nedat+1)
!    4-9    Grids 1-6      EDAT(nedat+3) thru EDAT(nedat+8)

! Possible 2nd card if this is a 10 node TETRA:

!    2-5    Grids 7-10     EDAT(nedat+9) thru EDAT(nedat+16) if this is a TETRA 10 node element

! Make JCARD from CARD

      CALL MKJCARD ( SUBR_NAME, CARD, JCARD )
      NAME = JCARD(1)
      ID   = JCARD(2)

! Set JCARD_EDAT to JCARD

      DO I=1,10
         JCARD_EDAT(I) = JCARD(I)
      ENDDO

! Read and check data

      IF (JCARD(8)(1:) == ' ') THEN                        ! If field 8 is blank then this must be a TETRA 4
         CALL ELEPRO ( 'Y', JCARD_EDAT, 6, 6, 'Y', 'Y', 'Y', 'Y', 'Y', 'Y', 'N', 'N' )
         ETYPE(NELE) = 'TETRA4  '
         NCTETRA4 = NCTETRA4 + 1
         NUM_GRD = 4
         CALL BD_IMBEDDED_BLANK ( JCARD,2,3,4,5,6,7,0,0 )
         CALL CARD_FLDS_NOT_BLANK ( JCARD,0,0,0,0,0,0,8,9 )
      ELSE                                                 ! This must be a TETRA 10 if field 8 is not blank
         CALL ELEPRO ( 'Y', JCARD_EDAT, 8, 8, 'Y', 'Y', 'Y', 'Y', 'Y', 'Y', 'Y', 'Y' )
         ETYPE(NELE) = 'TETRA10 '
         NCTETRA10 = NCTETRA10 + 1
         NUM_GRD = 10
         CALL BD_IMBEDDED_BLANK ( JCARD,2,3,4,5,6,7,8,9 )
      ENDIF
      CALL CRDERR ( CARD )

! 2nd card - required if TETRA 10:

      IF (ETYPE(NELE) == 'TETRA10 ') THEN

         IF (LARGE_FLD_INP == 'N') THEN
            CALL NEXTC  ( CARD, ICONT, IERR )
         ELSE
            CALL NEXTC2 ( CARD, ICONT, IERR, CHILD )
            CARD = CHILD
         ENDIF
         CALL MKJCARD ( SUBR_NAME, CARD, JCARD )
         IF (ICONT == 1) THEN

            DO I=1,10
               JCARD_EDAT(I) = JCARD(I)
            ENDDO

            IF (ETYPE(NELE) == 'TETRA10  ') THEN
               CALL ELEPRO ( 'N', JCARD_EDAT, 4, 4, 'Y', 'Y', 'Y', 'Y', 'N', 'N', 'N', 'N' )
               CALL BD_IMBEDDED_BLANK ( JCARD,2,3,4,5,6,7,8,9 )  ! Make sure that there are no imbedded blanks in fields 2-9
            ENDIF

         ELSE

            FATAL_ERR = FATAL_ERR + 1
            WRITE(ERR,1136) NAME, ID
            WRITE(F06,1136) NAME, ID

         ENDIF

      ENDIF

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
 1136 FORMAT(' *ERROR  1136: REQUIRED CONTINUATION FOR ',A,' ID = ',A,' MISSING')

! **********************************************************************************************************************************

      END SUBROUTINE BD_CTETRA


text:
      subroutine ampmidWZ(p1,p2,p3,p4,p5,p6,p7,p8,za,zb,amp)
      implicit none
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      include 'cmplxmass.f'
      include 'masses.f'
      include 'runstring.f'
      include 'ewcharge.f'
      include 'zcouple.f'
      include 'WWbits.f'
      integer jdu1,jdu2,i1,i2,i3,i4,
     & p1,p2,p3,p4,p5,p6,p7,p8
      double complex zab2,zba2,amp(2,2,2),sqwmass,rxw,cosW,
     & propw34,propz56,propz28,propw17,
     & propw3456,proph1347,gamZ56ne(2,2),
     & gamZ56ee(2,2),gamZ28nq(2,2,2),gamZ28eq(2,2,2),gamZ28qW(2,2),
     & gamZ56eW(2),gamZ28eqZ(2,2,2),gamZ28nqZ(2,2,2)
      double precision qn,t3,t4,s34,s56,s17,s28,s137,s147,s167,
     & s238,s248,s258,s268,s456,s345,s356,s3456,s1347,
     & twop17Dp3456,twop34Dp3456,q3,q4,l3,l4
      parameter(qn=0d0)
C-----Begin statement functions
      zab2(i1,i2,i3,i4)=za(i1,i2)*zb(i2,i4)+za(i1,i3)*zb(i3,i4)
      zba2(i1,i2,i3,i4)=zb(i1,i2)*za(i2,i4)+zb(i1,i3)*za(i3,i4)
      t3(i1,i2,i3)=s(i1,i2)+s(i2,i3)+s(i3,i1)
      t4(i1,i2,i3,i4)=s(i1,i2)+s(i1,i3)+s(i1,i4)
     &               +s(i2,i3)+s(i2,i4)+s(i3,i4)
C-----end statement functions

c--- special fix for Madgraph check
      if (index(runstring,'mad') .gt. 0) then
        sqwmass=dcmplx(wmass**2,0d0)
      else
        sqwmass=cwmass2
      endif

      rxw=sqrt((cone-cxw)/cxw)
      cosW=sqrt((cone-cxw))
      s3456=t4(p3,p4,p5,p6)
      s1347=t4(p1,p3,p4,p7)
      s137=t3(p1,p3,p7)
      s147=t3(p1,p4,p7)
      s167=t3(p1,p6,p7)
      s238=t3(p2,p3,p8)
      s248=t3(p2,p4,p8)
      s258=t3(p2,p5,p8)
      s268=t3(p2,p6,p8)
      s345=t3(p3,p4,p5)
      s356=t3(p3,p5,p6)
      s456=t3(p4,p5,p6)
      s34=s(p3,p4)
      s56=s(p5,p6)
      s17=s(p1,p7)
      s28=s(p2,p8)

      twop17Dp3456=s28-s17-s3456
      twop34Dp3456=s34-s56+s3456

      proph1347=dcmplx(s1347-hmass**2,hmass*hwidth)
      propw34=s34-dcmplx(wmass**2,-wmass*wwidth)
      propz56=s56-dcmplx(zmass**2,-zmass*zwidth)
      propw3456=s3456-dcmplx(wmass**2,-wmass*wwidth)
      propw17=s17-dcmplx(wmass**2,-wmass*wwidth)
      propz28=s28-dcmplx(zmass**2,-zmass*zwidth)

      q3=qn
      l3=ln
      q4=qe
      l4=le
      do jdu1=1,2
      gamZ28nq(jdu1,1,1)=q3*Q(jdu1)/s28+l3*L(jdu1)/propZ28
      gamZ28nq(jdu1,1,2)=q3*Q(jdu1)/s28+l3*R(jdu1)/propZ28

      gamZ28eq(jdu1,1,1)=q4*Q(jdu1)/s28+l4*L(jdu1)/propZ28
      gamZ28eq(jdu1,1,2)=q4*Q(jdu1)/s28+l4*R(jdu1)/propZ28

      gamZ28nqZ(jdu1,1,1)=qn*Q(jdu1)/s28+ln*L(jdu1)/propZ28
      gamZ28nqZ(jdu1,1,2)=qn*Q(jdu1)/s28+ln*R(jdu1)/propZ28

      gamZ28eqZ(jdu1,1,1)=qe*Q(jdu1)/s28+le*L(jdu1)/propZ28
      gamZ28eqZ(jdu1,1,2)=qe*Q(jdu1)/s28+le*R(jdu1)/propZ28

      gamZ28qW(jdu1,1)=Q(jdu1)/s28+rxw*L(jdu1)/propZ28
      gamZ28qW(jdu1,2)=Q(jdu1)/s28+rxw*R(jdu1)/propZ28

      gamZ56eW(1)=qe/s56+rxw*le/propZ56
      gamZ56eW(2)=qe/s56+rxw*re/propZ56
      enddo

      gamZ56ne(1,1)=q3*qe/s56+l3*le/propZ56
      gamZ56ne(1,2)=q3*qe/s56+l3*re/propZ56

      gamZ56ee(1,1)=q4*qe/s56+l4*le/propZ56
      gamZ56ee(1,2)=q4*qe/s56+l4*re/propZ56

      do jdu2=1,2
      amp(jdu2,1,1)= + L(jdu2)*cxw**(-2)*Hbit*cosW**(-2)*propw17**(-1)*
     & propw34**(-1)*propZ28**(-1) * (  - 2.D0*za(p3,p7)*za(p5,p8)*zb(
     &    p1,p4)*zb(p2,p6)*le*sqwmass*proph1347**(-1)*propZ56**(-1) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28qW(jdu2,1)*cwmass2**(-1)*
     & cxw**(-2)*Bbit*propw17**(-1)*propw34**(-1)*propw3456**(-1)*
     & s345**(-1) * ( za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p3,p4)*zab2(p3,
     &    p4,p5,p6)*s3456 + za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p3,p4)*
     &    zab2(p3,p4,p5,p6)*twop17Dp3456 - za(p3,p5)*za(p7,p8)*zb(p1,p2
     &    )*zb(p4,p5)*zab2(p5,p3,p4,p6)*s3456 - za(p3,p5)*za(p7,p8)*zb(
     &    p1,p2)*zb(p4,p5)*zab2(p5,p3,p4,p6)*twop17Dp3456 )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28qW(jdu2,1)*cxw**(-2)*Bbit*
     & propw17**(-1)*propw34**(-1)*propw3456**(-1)*s345**(-1) * ( 2.D0*
     &    za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p1,p6)*zba2(p4,p3,p5,p1) -
     &    za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p3,p4)*zab2(p3,p4,p5,p6) +
     &    za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p4,p5)*zab2(p5,p3,p4,p6) - 2.
     &    D0*za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p6,p7)*zba2(p4,p3,p5,p7)
     &     + 2.D0*za(p3,p5)*zb(p1,p6)*zab2(p8,p1,p7,p2)*zba2(p4,p3,p5,
     &    p7) - 2.D0*za(p3,p5)*zb(p2,p6)*zab2(p7,p2,p8,p1)*zba2(p4,p3,
     &    p5,p8) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28qW(jdu2,1)*gamZ56eW(1)*
     & cwmass2**(-1)*cxw**(-1)*Bbit*propw17**(-1)*propw34**(-1)*
     & propw3456**(-1) * (  - 2.D0*za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p4,
     &    p6)*s3456**2 + 2.D0*za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p4,p6)*
     &    twop34Dp3456*s3456 - 2.D0*za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p4
     &    ,p6)*twop17Dp3456*s3456 + 2.D0*za(p3,p5)*za(p7,p8)*zb(p1,p2)*
     &    zb(p4,p6)*twop17Dp3456*twop34Dp3456 )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28qW(jdu2,1)*gamZ56eW(1)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw34**(-1)*propw3456**(-1) * (
     &     - 4.D0*za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p4,p6)*s1347 + 2.D0*
     &    za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p4,p6)*s3456 + 4.D0*za(p3,p5
     &    )*za(p7,p8)*zb(p1,p2)*zb(p4,p6)*s34 + 4.D0*za(p3,p5)*za(p7,p8
     &    )*zb(p1,p2)*zb(p4,p6)*s17 - 2.D0*za(p3,p5)*za(p7,p8)*zb(p1,p2
     &    )*zb(p4,p6)*twop34Dp3456 + 2.D0*za(p3,p5)*za(p7,p8)*zb(p1,p2)
     &    *zb(p4,p6)*twop17Dp3456 + 4.D0*za(p3,p5)*zb(p4,p6)*zab2(p7,p2
     &    ,p8,p1)*zab2(p8,p3,p4,p2) - 4.D0*za(p3,p5)*zb(p4,p6)*zab2(p7,
     &    p3,p4,p1)*zab2(p8,p1,p7,p2) - 4.D0*za(p3,p7)*zb(p1,p4)*zab2(
     &    p5,p3,p4,p6)*zab2(p8,p1,p7,p2) + 4.D0*za(p3,p8)*zb(p2,p4)*
     &    zab2(p5,p3,p4,p6)*zab2(p7,p2,p8,p1) + 4.D0*za(p5,p7)*zb(p1,p6
     &    )*zab2(p3,p5,p6,p4)*zab2(p8,p1,p7,p2) - 4.D0*za(p5,p8)*zb(p2,
     &    p6)*zab2(p3,p5,p6,p4)*zab2(p7,p2,p8,p1) - 4.D0*za(p7,p8)*zb(
     &    p1,p2)*zab2(p3,p1,p7,p4)*zab2(p5,p3,p4,p6) + 4.D0*za(p7,p8)*
     &    zb(p1,p2)*zab2(p3,p5,p6,p4)*zab2(p5,p1,p7,p6) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28qW(jdu2,1)*gamZ56eW(1)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw34**(-1) * ( 2.D0*za(p3,p5)*
     &    za(p7,p8)*zb(p1,p2)*zb(p4,p6) - 4.D0*za(p3,p7)*za(p5,p8)*zb(
     &    p1,p4)*zb(p2,p6) + 2.D0*za(p3,p8)*za(p5,p7)*zb(p1,p6)*zb(p2,
     &    p4) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28qW(jdu2,1)*gamZ56ne(1,1)*
     & cwmass2**(-1)*cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*
     & s356**(-1) * (  - 2.D0*za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p3,p6)*
     &    zab2(p3,p5,p6,p4)*s3456 - 2.D0*za(p3,p5)*za(p7,p8)*zb(p1,p2)*
     &    zb(p3,p6)*zab2(p3,p5,p6,p4)*twop17Dp3456 - 2.D0*za(p3,p5)*za(
     &    p7,p8)*zb(p1,p2)*zb(p5,p6)*zab2(p5,p3,p6,p4)*s3456 - 2.D0*za(
     &    p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p5,p6)*zab2(p5,p3,p6,p4)*
     &    twop17Dp3456 )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28qW(jdu2,1)*gamZ56ne(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*s356**(-1) * (  - 4.
     &    D0*za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p1,p4)*zba2(p6,p3,p5,p1)
     &     + 2.D0*za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p3,p6)*zab2(p3,p5,p6
     &    ,p4) + 4.D0*za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p4,p7)*zba2(p6,
     &    p3,p5,p7) + 2.D0*za(p3,p5)*za(p7,p8)*zb(p1,p2)*zb(p5,p6)*
     &    zab2(p5,p3,p6,p4) - 4.D0*za(p3,p5)*zb(p1,p4)*zab2(p8,p1,p7,p2
     &    )*zba2(p6,p3,p5,p7) + 4.D0*za(p3,p5)*zb(p2,p4)*zab2(p7,p2,p8,
     &    p1)*zba2(p6,p3,p5,p8) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28qW(jdu2,1)*gamZ56ee(1,1)*
     & cwmass2**(-1)*cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*
     & s456**(-1) * ( 2.D0*za(p4,p5)*za(p7,p8)*zb(p1,p2)*zb(p4,p6)*
     &    zab2(p3,p5,p6,p4)*s3456 + 2.D0*za(p4,p5)*za(p7,p8)*zb(p1,p2)*
     &    zb(p4,p6)*zab2(p3,p5,p6,p4)*twop17Dp3456 - 2.D0*za(p5,p6)*za(
     &    p7,p8)*zb(p1,p2)*zb(p4,p6)*zab2(p3,p4,p5,p6)*s3456 - 2.D0*za(
     &    p5,p6)*za(p7,p8)*zb(p1,p2)*zb(p4,p6)*zab2(p3,p4,p5,p6)*
     &    twop17Dp3456 )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28qW(jdu2,1)*gamZ56ee(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*s456**(-1) * ( 4.D0
     &    *za(p1,p3)*za(p7,p8)*zb(p1,p2)*zb(p4,p6)*zba2(p1,p4,p6,p5) -
     &    4.D0*za(p3,p7)*za(p7,p8)*zb(p1,p2)*zb(p4,p6)*zba2(p7,p4,p6,p5
     &    ) - 4.D0*za(p3,p7)*zb(p4,p6)*zab2(p8,p1,p7,p2)*zba2(p1,p4,p6,
     &    p5) + 4.D0*za(p3,p8)*zb(p4,p6)*zab2(p7,p2,p8,p1)*zba2(p2,p4,
     &    p6,p5) - 2.D0*za(p4,p5)*za(p7,p8)*zb(p1,p2)*zb(p4,p6)*zab2(p3
     &    ,p5,p6,p4) + 2.D0*za(p5,p6)*za(p7,p8)*zb(p1,p2)*zb(p4,p6)*
     &    zab2(p3,p4,p5,p6) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28eq(jdu2,1,1)*gamZ56ne(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*s356**(-1) * (  - 4.D0*za(p3,p5)*
     &    zb(p2,p4)*zba2(p1,p2,p4,p8)*zba2(p6,p3,p5,p7)*s248**(-1) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28eq(jdu2,1,1)*gamZ56ee(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*s456**(-1) * ( 4.D0*za(p3,p7)*zb(p4
     &    ,p6)*zba2(p1,p3,p7,p8)*zba2(p2,p4,p6,p5)*s137**(-1) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28eq(jdu2,1,1)*gamZ56ee(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1) * (  - 4.D0*za(p3,p7)*zb(p2,p4)*
     &    zba2(p1,p3,p7,p5)*zba2(p6,p2,p4,p8)*s137**(-1)*s248**(-1) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28nq(jdu2,1,1)*gamZ56ne(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*s147**(-1)*s356**(-1) * (  - 4.D0*
     &    za(p3,p5)*zb(p1,p4)*zba2(p2,p1,p4,p7)*zba2(p6,p3,p5,p8) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28nq(jdu2,1,1)*gamZ56ne(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*s147**(-1) * (  - 4.D0*za(p3,p8)*
     &    zb(p1,p4)*zba2(p2,p3,p8,p5)*zba2(p6,p1,p4,p7)*s238**(-1) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28nq(jdu2,1,1)*gamZ56ee(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*s456**(-1) * ( 4.D0*za(p3,p8)*zb(p4
     &    ,p6)*zba2(p1,p4,p6,p5)*zba2(p2,p3,p8,p7)*s238**(-1) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28eqZ(jdu2,1,1)*cxw**(-2)*
     & Bbit*propw17**(-1)*propw34**(-1)*s258**(-1) * (  - 2.D0*za(p5,p8
     &    )*zb(p1,p6)*zba2(p2,p5,p8,p3)*zba2(p4,p1,p6,p7)*s167**(-1) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28eqZ(jdu2,1,1)*cxw**(-2)*
     & Bbit*propw17**(-1)*propw34**(-1)*s268**(-1)*s345**(-1) * ( 2.D0*
     &    za(p3,p5)*zb(p2,p6)*zba2(p1,p2,p6,p8)*zba2(p4,p3,p5,p7) )
      amp(jdu2,1,1) = amp(jdu2,1,1) + gamZ28nqZ(jdu2,1,1)*cxw**(-2)*
     & Bbit*propw17**(-1)*propw34**(-1)*s345**(-1) * ( 2.D0*za(p3,p5)*
     &    zb(p1,p6)*zba2(p2,p1,p6,p7)*zba2(p4,p3,p5,p8)*s167**(-1) )
      amp(jdu2,1,2)= + L(jdu2)*cxw**(-2)*Hbit*cosW**(-2)*propw17**(-1)*
     & propw34**(-1)*propZ28**(-1) * (  - 2.D0*za(p3,p7)*za(p6,p8)*zb(
     &    p1,p4)*zb(p2,p5)*re*sqwmass*proph1347**(-1)*propZ56**(-1) )
      amp(jdu2,1,2) = amp(jdu2,1,2) + gamZ28qW(jdu2,1)*gamZ56eW(2)*
     & cwmass2**(-1)*cxw**(-1)*Bbit*propw17**(-1)*propw34**(-1)*
     & propw3456**(-1) * (  - 2.D0*za(p3,p6)*za(p7,p8)*zb(p1,p2)*zb(p4,
     &    p5)*s3456**2 + 2.D0*za(p3,p6)*za(p7,p8)*zb(p1,p2)*zb(p4,p5)*
     &    twop34Dp3456*s3456 - 2.D0*za(p3,p6)*za(p7,p8)*zb(p1,p2)*zb(p4
     &    ,p5)*twop17Dp3456*s3456 + 2.D0*za(p3,p6)*za(p7,p8)*zb(p1,p2)*
     &    zb(p4,p5)*twop17Dp3456*twop34Dp3456 )
      amp(jdu2,1,2) = amp(jdu2,1,2) + gamZ28qW(jdu2,1)*gamZ56eW(2)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw34**(-1)*propw3456**(-1) * (
     &     - 4.D0*za(p3,p6)*za(p7,p8)*zb(p1,p2)*zb(p4,p5)*s1347 + 2.D0*
     &    za(p3,p6)*za(p7,p8)*zb(p1,p2)*zb(p4,p5)*s3456 + 4.D0*za(p3,p6
     &    )*za(p7,p8)*zb(p1,p2)*zb(p4,p5)*s34 + 4.D0*za(p3,p6)*za(p7,p8
     &    )*zb(p1,p2)*zb(p4,p5)*s17 - 2.D0*za(p3,p6)*za(p7,p8)*zb(p1,p2
     &    )*zb(p4,p5)*twop34Dp3456 + 2.D0*za(p3,p6)*za(p7,p8)*zb(p1,p2)
     &    *zb(p4,p5)*twop17Dp3456 + 4.D0*za(p3,p6)*zb(p4,p5)*zab2(p7,p2
     &    ,p8,p1)*zab2(p8,p3,p4,p2) - 4.D0*za(p3,p6)*zb(p4,p5)*zab2(p7,
     &    p3,p4,p1)*zab2(p8,p1,p7,p2) - 4.D0*za(p3,p7)*zb(p1,p4)*zab2(
     &    p8,p1,p7,p2)*zba2(p5,p3,p4,p6) + 4.D0*za(p3,p8)*zb(p2,p4)*
     &    zab2(p7,p2,p8,p1)*zba2(p5,p3,p4,p6) + 4.D0*za(p6,p7)*zb(p1,p5
     &    )*zab2(p3,p5,p6,p4)*zab2(p8,p1,p7,p2) - 4.D0*za(p6,p8)*zb(p2,
     &    p5)*zab2(p3,p5,p6,p4)*zab2(p7,p2,p8,p1) - 4.D0*za(p7,p8)*zb(
     &    p1,p2)*zab2(p3,p1,p7,p4)*zba2(p5,p3,p4,p6) + 4.D0*za(p7,p8)*
     &    zb(p1,p2)*zab2(p3,p5,p6,p4)*zba2(p5,p1,p7,p6) )
      amp(jdu2,1,2) = amp(jdu2,1,2) + gamZ28qW(jdu2,1)*gamZ56eW(2)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw34**(-1) * ( 2.D0*za(p3,p6)*
     &    za(p7,p8)*zb(p1,p2)*zb(p4,p5) - 4.D0*za(p3,p7)*za(p6,p8)*zb(
     &    p1,p4)*zb(p2,p5) + 2.D0*za(p3,p8)*za(p6,p7)*zb(p1,p5)*zb(p2,
     &    p4) )
      amp(jdu2,1,2) = amp(jdu2,1,2) + gamZ28qW(jdu2,1)*gamZ56ne(1,2)*
     & cwmass2**(-1)*cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*
     & s356**(-1) * (  - 2.D0*za(p3,p6)*za(p7,p8)*zb(p1,p2)*zb(p3,p5)*
     &    zab2(p3,p5,p6,p4)*s3456 - 2.D0*za(p3,p6)*za(p7,p8)*zb(p1,p2)*
     &    zb(p3,p5)*zab2(p3,p5,p6,p4)*twop17Dp3456 + 2.D0*za(p3,p6)*za(
     &    p7,p8)*zb(p1,p2)*zb(p5,p6)*zab2(p6,p3,p5,p4)*s3456 + 2.D0*za(
     &    p3,p6)*za(p7,p8)*zb(p1,p2)*zb(p5,p6)*zab2(p6,p3,p5,p4)*
     &    twop17Dp3456 )
      amp(jdu2,1,2) = amp(jdu2,1,2) + gamZ28qW(jdu2,1)*gamZ56ne(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*s356**(-1) * (  - 4.
     &    D0*za(p3,p6)*za(p7,p8)*zb(p1,p2)*zb(p1,p4)*zba2(p5,p3,p6,p1)
     &     + 2.D0*za(p3,p6)*za(p7,p8)*zb(p1,p2)*zb(p3,p5)*zab2(p3,p5,p6
     &    ,p4) + 4.D0*za(p3,p6)*za(p7,p8)*zb(p1,p2)*zb(p4,p7)*zba2(p5,
     &    p3,p6,p7) - 2.D0*za(p3,p6)*za(p7,p8)*zb(p1,p2)*zb(p5,p6)*
     &    zab2(p6,p3,p5,p4) - 4.D0*za(p3,p6)*zb(p1,p4)*zab2(p8,p1,p7,p2
     &    )*zba2(p5,p3,p6,p7) + 4.D0*za(p3,p6)*zb(p2,p4)*zab2(p7,p2,p8,
     &    p1)*zba2(p5,p3,p6,p8) )
      amp(jdu2,1,2) = amp(jdu2,1,2) + gamZ28qW(jdu2,1)*gamZ56ee(1,2)*
     & cwmass2**(-1)*cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*
     & s456**(-1) * ( 2.D0*za(p4,p6)*za(p7,p8)*zb(p1,p2)*zb(p4,p5)*
     &    zab2(p3,p5,p6,p4)*s3456 + 2.D0*za(p4,p6)*za(p7,p8)*zb(p1,p2)*
     &    zb(p4,p5)*zab2(p3,p5,p6,p4)*twop17Dp3456 + 2.D0*za(p5,p6)*za(
     &    p7,p8)*zb(p1,p2)*zb(p4,p5)*zab2(p3,p4,p6,p5)*s3456 + 2.D0*za(
     &    p5,p6)*za(p7,p8)*zb(p1,p2)*zb(p4,p5)*zab2(p3,p4,p6,p5)*
     &    twop17Dp3456 )
      amp(jdu2,1,2) = amp(jdu2,1,2) + gamZ28qW(jdu2,1)*gamZ56ee(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*s456**(-1) * ( 4.D0
     &    *za(p1,p3)*za(p7,p8)*zb(p1,p2)*zb(p4,p5)*zba2(p1,p4,p5,p6) -
     &    4.D0*za(p3,p7)*za(p7,p8)*zb(p1,p2)*zb(p4,p5)*zba2(p7,p4,p5,p6
     &    ) - 4.D0*za(p3,p7)*zb(p4,p5)*zab2(p8,p1,p7,p2)*zba2(p1,p4,p5,
     &    p6) + 4.D0*za(p3,p8)*zb(p4,p5)*zab2(p7,p2,p8,p1)*zba2(p2,p4,
     &    p5,p6) - 2.D0*za(p4,p6)*za(p7,p8)*zb(p1,p2)*zb(p4,p5)*zab2(p3
     &    ,p5,p6,p4) - 2.D0*za(p5,p6)*za(p7,p8)*zb(p1,p2)*zb(p4,p5)*
     &    zab2(p3,p4,p6,p5) )
      amp(jdu2,1,2) = amp(jdu2,1,2) + gamZ28eq(jdu2,1,1)*gamZ56ne(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*s356**(-1) * (  - 4.D0*za(p3,p6)*
     &    zb(p2,p4)*zba2(p1,p2,p4,p8)*zba2(p5,p3,p6,p7)*s248**(-1) )
      amp(jdu2,1,2) = amp(jdu2,1,2) + gamZ28eq(jdu2,1,1)*gamZ56ee(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*s456**(-1) * ( 4.D0*za(p3,p7)*zb(p4
     &    ,p5)*zba2(p1,p3,p7,p8)*zba2(p2,p4,p5,p6)*s137**(-1) )
      amp(jdu2,1,2) = amp(jdu2,1,2) + gamZ28eq(jdu2,1,1)*gamZ56ee(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1) * (  - 4.D0*za(p3,p7)*zb(p2,p4)*
     &    zba2(p1,p3,p7,p6)*zba2(p5,p2,p4,p8)*s137**(-1)*s248**(-1) )
      amp(jdu2,1,2) = amp(jdu2,1,2) + gamZ28nq(jdu2,1,1)*gamZ56ne(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*s147**(-1)*s356**(-1) * (  - 4.D0*
     &    za(p3,p6)*zb(p1,p4)*zba2(p2,p1,p4,p7)*zba2(p5,p3,p6,p8) )
      amp(jdu2,1,2) = amp(jdu2,1,2) + gamZ28nq(jdu2,1,1)*gamZ56ne(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*s147**(-1) * (  - 4.D0*za(p3,p8)*
     &    zb(p1,p4)*zba2(p2,p3,p8,p6)*zba2(p5,p1,p4,p7)*s238**(-1) )
      amp(jdu2,1,2) = amp(jdu2,1,2) + gamZ28nq(jdu2,1,1)*gamZ56ee(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*s456**(-1) * ( 4.D0*za(p3,p8)*zb(p4
     &    ,p5)*zba2(p1,p4,p5,p6)*zba2(p2,p3,p8,p7)*s238**(-1) )
      amp(jdu2,2,1)= + cxw**(-2)*Hbit*cosW**(-2)*propw17**(-1)*
     & propw34**(-1)*propZ28**(-1) * (  - 2.D0*R(jdu2)*za(p2,p5)*za(p3,
     &    p7)*zb(p1,p4)*zb(p6,p8)*le*sqwmass*proph1347**(-1)*
     &    propZ56**(-1) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28qW(jdu2,2)*cwmass2**(-1)*
     & cxw**(-2)*Bbit*propw17**(-1)*propw34**(-1)*propw3456**(-1)*
     & s345**(-1) * (  - za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p3,p4)*zab2(
     &    p3,p4,p5,p6)*s3456 - za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p3,p4)*
     &    zab2(p3,p4,p5,p6)*twop17Dp3456 + za(p2,p7)*za(p3,p5)*zb(p1,p8
     &    )*zb(p4,p5)*zab2(p5,p3,p4,p6)*s3456 + za(p2,p7)*za(p3,p5)*zb(
     &    p1,p8)*zb(p4,p5)*zab2(p5,p3,p4,p6)*twop17Dp3456 )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28qW(jdu2,2)*cxw**(-2)*Bbit*
     & propw17**(-1)*propw34**(-1)*propw3456**(-1)*s345**(-1) * (  - 2.D
     &    0*za(p2,p7)*za(p3,p5)*zb(p1,p6)*zb(p1,p8)*zba2(p4,p3,p5,p1)
     &     + za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p3,p4)*zab2(p3,p4,p5,p6)
     &     - za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p4,p5)*zab2(p5,p3,p4,p6)
     &     + 2.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p6,p7)*zba2(p4,p3,p5
     &    ,p7) + 2.D0*za(p3,p5)*zb(p1,p6)*zba2(p4,p3,p5,p7)*zba2(p8,p1,
     &    p7,p2) + 2.D0*za(p3,p5)*zb(p6,p8)*zab2(p7,p2,p8,p1)*zba2(p4,
     &    p3,p5,p2) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28qW(jdu2,2)*gamZ56eW(1)*
     & cwmass2**(-1)*cxw**(-1)*Bbit*propw17**(-1)*propw34**(-1)*
     & propw3456**(-1) * ( 2.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p4,p6)
     &    *s3456**2 - 2.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p4,p6)*
     &    twop34Dp3456*s3456 + 2.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p4
     &    ,p6)*twop17Dp3456*s3456 - 2.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8)*
     &    zb(p4,p6)*twop17Dp3456*twop34Dp3456 )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28qW(jdu2,2)*gamZ56eW(1)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw34**(-1)*propw3456**(-1) * ( 4.
     &    D0*za(p2,p3)*zb(p4,p8)*zab2(p5,p3,p4,p6)*zab2(p7,p2,p8,p1) -
     &    4.D0*za(p2,p5)*zb(p6,p8)*zab2(p3,p5,p6,p4)*zab2(p7,p2,p8,p1)
     &     + 4.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p4,p6)*s1347 - 2.D0*
     &    za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p4,p6)*s3456 - 4.D0*za(p2,p7
     &    )*za(p3,p5)*zb(p1,p8)*zb(p4,p6)*s34 - 4.D0*za(p2,p7)*za(p3,p5
     &    )*zb(p1,p8)*zb(p4,p6)*s17 + 2.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8
     &    )*zb(p4,p6)*twop34Dp3456 - 2.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8)
     &    *zb(p4,p6)*twop17Dp3456 + 4.D0*za(p2,p7)*zb(p1,p8)*zab2(p3,p1
     &    ,p7,p4)*zab2(p5,p3,p4,p6) - 4.D0*za(p2,p7)*zb(p1,p8)*zab2(p3,
     &    p5,p6,p4)*zab2(p5,p1,p7,p6) + 4.D0*za(p3,p5)*zb(p4,p6)*zab2(
     &    p7,p2,p8,p1)*zba2(p8,p3,p4,p2) - 4.D0*za(p3,p5)*zb(p4,p6)*
     &    zab2(p7,p3,p4,p1)*zba2(p8,p1,p7,p2) - 4.D0*za(p3,p7)*zb(p1,p4
     &    )*zab2(p5,p3,p4,p6)*zba2(p8,p1,p7,p2) + 4.D0*za(p5,p7)*zb(p1,
     &    p6)*zab2(p3,p5,p6,p4)*zba2(p8,p1,p7,p2) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28qW(jdu2,2)*gamZ56eW(1)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw34**(-1) * ( 2.D0*za(p2,p3)*
     &    za(p5,p7)*zb(p1,p6)*zb(p4,p8) - 4.D0*za(p2,p5)*za(p3,p7)*zb(
     &    p1,p4)*zb(p6,p8) - 2.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p4,
     &    p6) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28qW(jdu2,2)*gamZ56ne(1,1)*
     & cwmass2**(-1)*cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*
     & s356**(-1) * ( 2.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p3,p6)*
     &    zab2(p3,p5,p6,p4)*s3456 + 2.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8)*
     &    zb(p3,p6)*zab2(p3,p5,p6,p4)*twop17Dp3456 + 2.D0*za(p2,p7)*za(
     &    p3,p5)*zb(p1,p8)*zb(p5,p6)*zab2(p5,p3,p6,p4)*s3456 + 2.D0*za(
     &    p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p5,p6)*zab2(p5,p3,p6,p4)*
     &    twop17Dp3456 )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28qW(jdu2,2)*gamZ56ne(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*s356**(-1) * ( 4.D0
     &    *za(p2,p7)*za(p3,p5)*zb(p1,p4)*zb(p1,p8)*zba2(p6,p3,p5,p1) -
     &    2.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p3,p6)*zab2(p3,p5,p6,p4
     &    ) - 4.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p4,p7)*zba2(p6,p3,
     &    p5,p7) - 2.D0*za(p2,p7)*za(p3,p5)*zb(p1,p8)*zb(p5,p6)*zab2(p5
     &    ,p3,p6,p4) - 4.D0*za(p3,p5)*zb(p1,p4)*zba2(p6,p3,p5,p7)*zba2(
     &    p8,p1,p7,p2) - 4.D0*za(p3,p5)*zb(p4,p8)*zab2(p7,p2,p8,p1)*
     &    zba2(p6,p3,p5,p2) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28qW(jdu2,2)*gamZ56ee(1,1)*
     & cwmass2**(-1)*cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*
     & s456**(-1) * (  - 2.D0*za(p2,p7)*za(p4,p5)*zb(p1,p8)*zb(p4,p6)*
     &    zab2(p3,p5,p6,p4)*s3456 - 2.D0*za(p2,p7)*za(p4,p5)*zb(p1,p8)*
     &    zb(p4,p6)*zab2(p3,p5,p6,p4)*twop17Dp3456 + 2.D0*za(p2,p7)*za(
     &    p5,p6)*zb(p1,p8)*zb(p4,p6)*zab2(p3,p4,p5,p6)*s3456 + 2.D0*za(
     &    p2,p7)*za(p5,p6)*zb(p1,p8)*zb(p4,p6)*zab2(p3,p4,p5,p6)*
     &    twop17Dp3456 )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28qW(jdu2,2)*gamZ56ee(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*s456**(-1) * (  - 4.
     &    D0*za(p1,p3)*za(p2,p7)*zb(p1,p8)*zb(p4,p6)*zba2(p1,p4,p6,p5)
     &     - 4.D0*za(p2,p3)*zb(p4,p6)*zab2(p7,p2,p8,p1)*zba2(p8,p4,p6,
     &    p5) + 4.D0*za(p2,p7)*za(p3,p7)*zb(p1,p8)*zb(p4,p6)*zba2(p7,p4
     &    ,p6,p5) + 2.D0*za(p2,p7)*za(p4,p5)*zb(p1,p8)*zb(p4,p6)*zab2(
     &    p3,p5,p6,p4) - 2.D0*za(p2,p7)*za(p5,p6)*zb(p1,p8)*zb(p4,p6)*
     &    zab2(p3,p4,p5,p6) - 4.D0*za(p3,p7)*zb(p4,p6)*zba2(p1,p4,p6,p5
     &    )*zba2(p8,p1,p7,p2) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28eq(jdu2,1,2)*gamZ56ne(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*s356**(-1) * ( 4.D0*za(p3,p5)*zb(p4
     &    ,p8)*zba2(p1,p4,p8,p2)*zba2(p6,p3,p5,p7)*s248**(-1) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28eq(jdu2,1,2)*gamZ56ee(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*s456**(-1) * ( 4.D0*za(p3,p7)*zb(p4
     &    ,p6)*zba2(p1,p3,p7,p2)*zba2(p8,p4,p6,p5)*s137**(-1) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28eq(jdu2,1,2)*gamZ56ee(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1) * ( 4.D0*za(p3,p7)*zb(p4,p8)*zba2(
     &    p1,p3,p7,p5)*zba2(p6,p4,p8,p2)*s137**(-1)*s248**(-1) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28nq(jdu2,1,2)*gamZ56ne(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*s147**(-1)*s356**(-1) * (  - 4.D0*
     &    za(p3,p5)*zb(p1,p4)*zba2(p6,p3,p5,p2)*zba2(p8,p1,p4,p7) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28nq(jdu2,1,2)*gamZ56ne(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*s147**(-1) * ( 4.D0*za(p2,p3)*zb(p1
     &    ,p4)*zba2(p6,p1,p4,p7)*zba2(p8,p2,p3,p5)*s238**(-1) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28nq(jdu2,1,2)*gamZ56ee(1,1)*
     & cxw**(-1)*Bbit*propw17**(-1)*s456**(-1) * (  - 4.D0*za(p2,p3)*
     &    zb(p4,p6)*zba2(p1,p4,p6,p5)*zba2(p8,p2,p3,p7)*s238**(-1) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28eqZ(jdu2,1,2)*cxw**(-2)*
     & Bbit*propw17**(-1)*propw34**(-1)*s258**(-1) * ( 2.D0*za(p2,p5)*
     &    zb(p1,p6)*zba2(p4,p1,p6,p7)*zba2(p8,p2,p5,p3)*s167**(-1) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28eqZ(jdu2,1,2)*cxw**(-2)*
     & Bbit*propw17**(-1)*propw34**(-1)*s268**(-1)*s345**(-1) * (  - 2.D
     &    0*za(p3,p5)*zb(p6,p8)*zba2(p1,p6,p8,p2)*zba2(p4,p3,p5,p7) )
      amp(jdu2,2,1) = amp(jdu2,2,1) + gamZ28nqZ(jdu2,1,2)*cxw**(-2)*
     & Bbit*propw17**(-1)*propw34**(-1)*s345**(-1) * ( 2.D0*za(p3,p5)*
     &    zb(p1,p6)*zba2(p4,p3,p5,p2)*zba2(p8,p1,p6,p7)*s167**(-1) )
      amp(jdu2,2,2)= + cxw**(-2)*Hbit*cosW**(-2)*propw17**(-1)*
     & propw34**(-1)*propZ28**(-1) * (  - 2.D0*R(jdu2)*za(p2,p6)*za(p3,
     &    p7)*zb(p1,p4)*zb(p5,p8)*re*sqwmass*proph1347**(-1)*
     &    propZ56**(-1) )
      amp(jdu2,2,2) = amp(jdu2,2,2) + gamZ28qW(jdu2,2)*gamZ56eW(2)*
     & cwmass2**(-1)*cxw**(-1)*Bbit*propw17**(-1)*propw34**(-1)*
     & propw3456**(-1) * ( 2.D0*za(p2,p7)*za(p3,p6)*zb(p1,p8)*zb(p4,p5)
     &    *s3456**2 - 2.D0*za(p2,p7)*za(p3,p6)*zb(p1,p8)*zb(p4,p5)*
     &    twop34Dp3456*s3456 + 2.D0*za(p2,p7)*za(p3,p6)*zb(p1,p8)*zb(p4
     &    ,p5)*twop17Dp3456*s3456 - 2.D0*za(p2,p7)*za(p3,p6)*zb(p1,p8)*
     &    zb(p4,p5)*twop17Dp3456*twop34Dp3456 )
      amp(jdu2,2,2) = amp(jdu2,2,2) + gamZ28qW(jdu2,2)*gamZ56eW(2)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw34**(-1)*propw3456**(-1) * ( 4.
     &    D0*za(p2,p3)*zb(p4,p8)*zab2(p7,p2,p8,p1)*zba2(p5,p3,p4,p6) -
     &    4.D0*za(p2,p6)*zb(p5,p8)*zab2(p3,p5,p6,p4)*zab2(p7,p2,p8,p1)
     &     + 4.D0*za(p2,p7)*za(p3,p6)*zb(p1,p8)*zb(p4,p5)*s1347 - 2.D0*
     &    za(p2,p7)*za(p3,p6)*zb(p1,p8)*zb(p4,p5)*s3456 - 4.D0*za(p2,p7
     &    )*za(p3,p6)*zb(p1,p8)*zb(p4,p5)*s34 - 4.D0*za(p2,p7)*za(p3,p6
     &    )*zb(p1,p8)*zb(p4,p5)*s17 + 2.D0*za(p2,p7)*za(p3,p6)*zb(p1,p8
     &    )*zb(p4,p5)*twop34Dp3456 - 2.D0*za(p2,p7)*za(p3,p6)*zb(p1,p8)
     &    *zb(p4,p5)*twop17Dp3456 + 4.D0*za(p2,p7)*zb(p1,p8)*zab2(p3,p1
     &    ,p7,p4)*zba2(p5,p3,p4,p6) - 4.D0*za(p2,p7)*zb(p1,p8)*zab2(p3,
     &    p5,p6,p4)*zba2(p5,p1,p7,p6) + 4.D0*za(p3,p6)*zb(p4,p5)*zab2(
     &    p7,p2,p8,p1)*zba2(p8,p3,p4,p2) - 4.D0*za(p3,p6)*zb(p4,p5)*
     &    zab2(p7,p3,p4,p1)*zba2(p8,p1,p7,p2) - 4.D0*za(p3,p7)*zb(p1,p4
     &    )*zba2(p5,p3,p4,p6)*zba2(p8,p1,p7,p2) + 4.D0*za(p6,p7)*zb(p1,
     &    p5)*zab2(p3,p5,p6,p4)*zba2(p8,p1,p7,p2) )
      amp(jdu2,2,2) = amp(jdu2,2,2) + gamZ28qW(jdu2,2)*gamZ56eW(2)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw34**(-1) * ( 2.D0*za(p2,p3)*
     &    za(p6,p7)*zb(p1,p5)*zb(p4,p8) - 4.D0*za(p2,p6)*za(p3,p7)*zb(
     &    p1,p4)*zb(p5,p8) - 2.D0*za(p2,p7)*za(p3,p6)*zb(p1,p8)*zb(p4,
     &    p5) )
      amp(jdu2,2,2) = amp(jdu2,2,2) + gamZ28qW(jdu2,2)*gamZ56ne(1,2)*
     & cwmass2**(-1)*cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*
     & s356**(-1) * ( 2.D0*za(p2,p7)*za(p3,p6)*zb(p1,p8)*zb(p3,p5)*
     &    zab2(p3,p5,p6,p4)*s3456 + 2.D0*za(p2,p7)*za(p3,p6)*zb(p1,p8)*
     &    zb(p3,p5)*zab2(p3,p5,p6,p4)*twop17Dp3456 - 2.D0*za(p2,p7)*za(
     &    p3,p6)*zb(p1,p8)*zb(p5,p6)*zab2(p6,p3,p5,p4)*s3456 - 2.D0*za(
     &    p2,p7)*za(p3,p6)*zb(p1,p8)*zb(p5,p6)*zab2(p6,p3,p5,p4)*
     &    twop17Dp3456 )
      amp(jdu2,2,2) = amp(jdu2,2,2) + gamZ28qW(jdu2,2)*gamZ56ne(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*s356**(-1) * ( 4.D0
     &    *za(p2,p7)*za(p3,p6)*zb(p1,p4)*zb(p1,p8)*zba2(p5,p3,p6,p1) -
     &    2.D0*za(p2,p7)*za(p3,p6)*zb(p1,p8)*zb(p3,p5)*zab2(p3,p5,p6,p4
     &    ) - 4.D0*za(p2,p7)*za(p3,p6)*zb(p1,p8)*zb(p4,p7)*zba2(p5,p3,
     &    p6,p7) + 2.D0*za(p2,p7)*za(p3,p6)*zb(p1,p8)*zb(p5,p6)*zab2(p6
     &    ,p3,p5,p4) - 4.D0*za(p3,p6)*zb(p1,p4)*zba2(p5,p3,p6,p7)*zba2(
     &    p8,p1,p7,p2) - 4.D0*za(p3,p6)*zb(p4,p8)*zab2(p7,p2,p8,p1)*
     &    zba2(p5,p3,p6,p2) )
      amp(jdu2,2,2) = amp(jdu2,2,2) + gamZ28qW(jdu2,2)*gamZ56ee(1,2)*
     & cwmass2**(-1)*cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*
     & s456**(-1) * (  - 2.D0*za(p2,p7)*za(p4,p6)*zb(p1,p8)*zb(p4,p5)*
     &    zab2(p3,p5,p6,p4)*s3456 - 2.D0*za(p2,p7)*za(p4,p6)*zb(p1,p8)*
     &    zb(p4,p5)*zab2(p3,p5,p6,p4)*twop17Dp3456 - 2.D0*za(p2,p7)*za(
     &    p5,p6)*zb(p1,p8)*zb(p4,p5)*zab2(p3,p4,p6,p5)*s3456 - 2.D0*za(
     &    p2,p7)*za(p5,p6)*zb(p1,p8)*zb(p4,p5)*zab2(p3,p4,p6,p5)*
     &    twop17Dp3456 )
      amp(jdu2,2,2) = amp(jdu2,2,2) + gamZ28qW(jdu2,2)*gamZ56ee(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*propw3456**(-1)*s456**(-1) * (  - 4.
     &    D0*za(p1,p3)*za(p2,p7)*zb(p1,p8)*zb(p4,p5)*zba2(p1,p4,p5,p6)
     &     - 4.D0*za(p2,p3)*zb(p4,p5)*zab2(p7,p2,p8,p1)*zba2(p8,p4,p5,
     &    p6) + 4.D0*za(p2,p7)*za(p3,p7)*zb(p1,p8)*zb(p4,p5)*zba2(p7,p4
     &    ,p5,p6) + 2.D0*za(p2,p7)*za(p4,p6)*zb(p1,p8)*zb(p4,p5)*zab2(
     &    p3,p5,p6,p4) + 2.D0*za(p2,p7)*za(p5,p6)*zb(p1,p8)*zb(p4,p5)*
     &    zab2(p3,p4,p6,p5) - 4.D0*za(p3,p7)*zb(p4,p5)*zba2(p1,p4,p5,p6
     &    )*zba2(p8,p1,p7,p2) )
      amp(jdu2,2,2) = amp(jdu2,2,2) + gamZ28eq(jdu2,1,2)*gamZ56ne(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*s356**(-1) * ( 4.D0*za(p3,p6)*zb(p4
     &    ,p8)*zba2(p1,p4,p8,p2)*zba2(p5,p3,p6,p7)*s248**(-1) )
      amp(jdu2,2,2) = amp(jdu2,2,2) + gamZ28eq(jdu2,1,2)*gamZ56ee(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*s456**(-1) * ( 4.D0*za(p3,p7)*zb(p4
     &    ,p5)*zba2(p1,p3,p7,p2)*zba2(p8,p4,p5,p6)*s137**(-1) )
      amp(jdu2,2,2) = amp(jdu2,2,2) + gamZ28eq(jdu2,1,2)*gamZ56ee(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1) * ( 4.D0*za(p3,p7)*zb(p4,p8)*zba2(
     &    p1,p3,p7,p6)*zba2(p5,p4,p8,p2)*s137**(-1)*s248**(-1) )
      amp(jdu2,2,2) = amp(jdu2,2,2) + gamZ28nq(jdu2,1,2)*gamZ56ne(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*s147**(-1)*s356**(-1) * (  - 4.D0*
     &    za(p3,p6)*zb(p1,p4)*zba2(p5,p3,p6,p2)*zba2(p8,p1,p4,p7) )
      amp(jdu2,2,2) = amp(jdu2,2,2) + gamZ28nq(jdu2,1,2)*gamZ56ne(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*s147**(-1) * ( 4.D0*za(p2,p3)*zb(p1
     &    ,p4)*zba2(p5,p1,p4,p7)*zba2(p8,p2,p3,p6)*s238**(-1) )
      amp(jdu2,2,2) = amp(jdu2,2,2) + gamZ28nq(jdu2,1,2)*gamZ56ee(1,2)*
     & cxw**(-1)*Bbit*propw17**(-1)*s456**(-1) * (  - 4.D0*za(p2,p3)*
     &    zb(p4,p5)*zba2(p1,p4,p5,p6)*zba2(p8,p2,p3,p7)*s238**(-1) )
      enddo

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 02.03.03/06                           *
C  *    TEST TITLE : Modelling transformation and clipping *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

C clipping indicator
C                noclip      clip
      INTEGER    PNCLIP,     PCLIP
      PARAMETER (PNCLIP = 0, PCLIP = 1)

C modelling clipping operator
      INTEGER    PMCREP,     PMCINT
      PARAMETER (PMCREP = 1, PMCINT = 2)

C composition type
C                preconcatenate  postconcatenate  replace
      INTEGER    PCPRE,          PCPOST,          PCREPL
      PARAMETER (PCPRE = 0,      PCPOST = 1,      PCREPL = 2)

      INTEGER     IPTHSZ
      PARAMETER  (IPTHSZ = 20)

      INTEGER     M1LEN,M2LEN, PMLEN, STPATH(2,5), IDUM1, IDUM(2)
      INTEGER     FPTHSZ,FPATH(2,IPTHSZ), EXLEN,EXPATH(2,IPTHSZ)
      INTEGER     FALEN

      REAL        M1(4,4), M2(4,4), PMLOCX(5), PMLOCY(5), PMLOCZ(5)
      REAL        FALOCX(10), FALOCY(10), FALOCZ(10), TX,TY,TZ
      REAL        TPLOCX(10), TPLOCY(10), TPLOCZ(10), RDUM
      REAL        MCOMP(4,4), EXLOCX,EXLOCY,EXLOCZ, SDIST, PTREGD
      REAL        SRPX,SRPY,SRPZ, HALFSP(6,10), NVECX,NVECY,NVECZ
      REAL        CLPTX,CLPTY,CLPTZ, CLVCX,CLVCY,CLVCZ
      REAL        DVECX,DVECY,DVECZ, NVECL, VECL

      LOGICAL     IAREQL, SPHDIS

      CALL INITGL ('02.03.03/06')

C open PHIGS
      CALL XPOPPH (ERRFIL, MEMUN)

      IF (.NOT. SPHDIS(123)) THEN
         CALL INMSG ('Skipping test because implementation '     //
     1               'apparently uses non-Euclidean metric for ' //
     2               'search distance.')
         GOTO 666
      ENDIF

C  m1 = random transformation matrix
C     =   1.43  -2.41  3.25 -1.15
C        -3.25  -3.13  2.54 -1.35
C        -2.43  -1.54  2.31  1.42
C         3.51  -2.32  2.41  3.25

      CALL SETRVS (' 1.43,  -3.25, -2.43,  3.51,  ' //
     1             '-2.41,  -3.13, -1.54, -2.32,  ' //
     2             ' 3.25,   2.54,  2.31,  2.41,  ' //
     3             '-1.15,  -1.35,  1.42,  3.25', M1, M1LEN)

C  m2 = random transformation matrix
C     =   -1.23   1.23  -1.02   9.30
C          1.28  -9.34   0.93  -4.89
C          2.03   9.85   4.01  -9.23
C          4.18  -0.12  -9.34   8.01

      CALL SETRVS ( '-1.23,  1.28,  2.03,  4.18, ' //
     1              ' 1.23, -9.34,  9.85, -0.12, ' //
     2              '-1.02,  0.93,  4.01, -9.34, ' //
     3              ' 9.30, -4.89, -9.23,  8.01', M2, M2LEN)

C  *** *** *** ***   Modelling transformations   *** *** *** ***
C
C  set up structure #101:
C  1. set local transformation = m1
C  2. execute structure 102
      CALL POPST  (101)
      CALL PSLMT3 (M1, PCREPL)
      CALL PEXST  (102)
      CALL PCLST

C  pmloc = polymarker location = -6.347, 4.637, 1.082
      PMLOCX(1) = -6.347
      PMLOCY(1) =  4.637
      PMLOCZ(1) =  1.082

C  set up structure #102:
C  1. set global transformation = m2
C  2. polymarker at pmloc
      CALL POPST  (102)
      CALL PSGMT3 (M2)
      CALL PPM3   (1, PMLOCX, PMLOCY, PMLOCZ)
      CALL PCLST

C  exloc = expected WC location = m2 X m1 X pmloc
      CALL ECOM3 (M2,M1, MCOMP)
      CALL ETP3  (PMLOCX(1),PMLOCY(1),PMLOCZ(1), MCOMP,
     1            EXLOCX,EXLOCY,EXLOCZ)

      SDIST = 0.1
      SRPX = EXLOCX + SDIST
      SRPY = EXLOCY
      SRPZ = EXLOCZ

      CALL SETMSG ('2 3 4 5 6 8 10', 'ISS should detect a primitive ' //
     1             'according to its location in world coordinates, ' //
     2             'computed with the current composite modelling '   //
     3             'transformation, including the effect of '         //
     4             'inherited transformations.')

C  ISS twice with SRP, varying search distance
C    starting path:  101,1, 102,1
C    search ceiling: 2
C    expected path:  101,1, 102,2

      CALL ISSAB (SRPX,SRPY,SRPZ, SDIST, '101,1, 102,1', PNCLIP, 2,
     1            '101,1, 102,2')

C  exloc = expected WC location = m2 X pmloc
      CALL ETP3  (PMLOCX(1),PMLOCY(1),PMLOCZ(1), M2,
     1            EXLOCX,EXLOCY,EXLOCZ)

      SDIST = 0.1
      SRPX = EXLOCX + SDIST
      SRPY = EXLOCY
      SRPZ = EXLOCZ

      CALL SETMSG ('2 3 4 5 6 10', 'ISS should detect a primitive '   //
     1             'according to its location in world coordinates, ' //
     2             'computed with the current composite modelling '   //
     3             'transformation when there are no inherited '      //
     4             'transformations.')

C  ISS twice with SRP, varying search distance
C    starting path:  102,1
C    search ceiling: 1
C    expected path:  102,2

      CALL ISSAB (SRPX,SRPY,SRPZ, SDIST, '102,1', PNCLIP, 1, '102,2')

C  *** *** *** ***   Modelling clipping   *** *** *** ***
C
C  clear CSS
      CALL PDAS

C  set up structure #101:
C    1. set modelling clipping CLIP
C    2. set modelling clipping volume with:
C          point         =  5, 5, 5
C          normal vector = -1,-1,-1
C    3. set local transformation = shift by +3,+3,+3
C    4. polymarker at MC: 4,4,4 (WC will be 7,7,7)
C    5. polymarker at MC: 1,1,1 (WC will be 4,4,4)

      CALL POPST  (101)
      CALL PSMCLI (PCLIP)
      CALL SETRVS ('5,5,5, -1,-1,-1', HALFSP, IDUM1)
      CALL PSMCV3 (PMCREP, 1, HALFSP)
      CALL ETR3   (3.0,3.0,3.0, M2)
      CALL PSLMT3 (M2, PCREPL)
      CALL SETRVS ('4,1', PMLOCX, PMLEN)
      CALL SETRVS ('4,1', PMLOCY, PMLEN)
      CALL SETRVS ('4,1', PMLOCZ, PMLEN)
      CALL PPM3   (1, PMLOCX(1), PMLOCY(1), PMLOCZ(1))
      CALL PPM3   (1, PMLOCX(2), PMLOCY(2), PMLOCZ(2))

      CALL SETMSG ('2 3 4 5 6 10', 'If the modelling clip flag is '   //
     1             'NOCLIP, ISS should detect a primitive whether '   //
     2             'or not its location in world coordinates is '     //
     3             'outside the current modelling clipping volume, '  //
     4             'even if the current modelling clipping '          //
     5             'indicator is CLIP.')

C  ISS twice with
C    modelling clip flag: NOCLIP
C    SRP:                 5.5,5.5,5.5
C    search distance:     1.02*1.5
C
C  pass/fail depending on
C    (found at 101,4 with starting path = 101,0  and
C     found at 101,5 with starting path = 101,4)

      CALL SETVS ('101,0', STPATH, IDUM1)
      CALL SETVS ('101,4', EXPATH, EXLEN)
      CALL PISS3 (5.5,5.5,5.5, 1.02*1.5, 1, STPATH, PNCLIP, 1,
     1            0,IDUM,IDUM,IDUM,IDUM, 0,IDUM,IDUM,IDUM,IDUM,
     2            IPTHSZ, ERRIND, FPTHSZ, FPATH)

      IF (ERRIND .EQ. 0 .AND. IAREQL(2*FPTHSZ,FPATH, EXLEN,EXPATH)) THEN
C        OK so far
      ELSE
         CALL FAIL
         GOTO 300
      ENDIF

      CALL SETVS ('101,4', STPATH, IDUM1)
      CALL SETVS ('101,5', EXPATH, EXLEN)
      CALL PISS3 (5.5,5.5,5.5, 1.02*1.5, 1, STPATH, PNCLIP, 1,
     1            0,IDUM,IDUM,IDUM,IDUM, 0,IDUM,IDUM,IDUM,IDUM,
     2            IPTHSZ, ERRIND, FPTHSZ, FPATH)

      CALL IFPF (ERRIND.EQ.0 .AND. IAREQL(2*FPTHSZ,FPATH, EXLEN,EXPATH))

300   CONTINUE

C  change element #1 of structure 101 to:
C    1. set modelling clipping NOCLIP
      CALL PSEP (1)
      CALL PDEL
      CALL PSMCLI (PNCLIP)
      CALL PCLST

      CALL SETMSG ('2 3 4 5 6 10', 'If the modelling clip flag is '   //
     1             'CLIP, ISS should detect only those primitives '   //
     2             'whose location in world coordinates is inside '   //
     3             'the current modelling clipping volume, even if '  //
     4             'the current modelling clipping indicator is '     //
     5             'NOCLIP.')

C  ISS with
C    modelling clip flag: CLIP
C    SRP:                 5.5,5.5,5.5
C    search distance:     1.02*1.5
C
C  pass/fail depending on
C    (found at 101,5 with starting path = 101,0)

      CALL SETVS ('101,0', STPATH, IDUM1)
      CALL SETVS ('101,5', EXPATH, EXLEN)
      CALL PISS3 (5.5,5.5,5.5, 1.02*1.5, 1, STPATH, PCLIP, 1,
     1            0,IDUM,IDUM,IDUM,IDUM, 0,IDUM,IDUM,IDUM,IDUM,
     2            IPTHSZ, ERRIND, FPTHSZ, FPATH)

      CALL IFPF (ERRIND.EQ.0 .AND. IAREQL(2*FPTHSZ,FPATH, EXLEN,EXPATH))

C  *** *** *** ***   Complex primitive   *** *** *** ***
C
C  clear CSS
      CALL PDAS

C  set up vertices for complex fill area:
C
C     1----------2    p1: 0, 1
C     |          |    p2: 3, 1
C     |  5---6   |    p3: 3,-2
C     |  |   |   |    p4: 1,-2
C     |  |   |   |    p5: 1, 0
C     8--+---7   |    p6: 2, 0
C        |       |    p7: 2,-1
C        4-------3    p8: 0,-1

      CALL SETRVS (' 0, 3, 3, 1, 1, 2, 2, 0', FALOCX, FALEN)
      CALL SETRVS (' 1, 1,-2,-2, 0, 0,-1,-1', FALOCY, FALEN)
      CALL SETRVS (' 0, 0, 0, 0, 0, 0, 0, 0', FALOCZ, FALEN)

C  tp1-tp8 = WC vertices = m1 X p1-p8
      CALL ARRTP3 (FALEN, FALOCX, FALOCY, FALOCZ, M1,
     1                    TPLOCX, TPLOCY, TPLOCZ)

C  nvec = unit vector normal to WC plane of fill area
      CALL PTSPL (FALEN, TPLOCX,TPLOCY,TPLOCZ, TX,TY,TZ, RDUM)
      CALL VEC1  (TX,TY,TZ, NVECX,NVECY,NVECZ)

C  set up structure #101:
C  1. set global transformation = m1
C  2. fill area at p1-p8
      CALL POPST  (101)
      CALL PSGMT3 (M1)
      CALL PFA    (FALEN, FALOCX,FALOCY)
      CALL PCLST

C  SRP   = point over hole in fill area = average of pt5,pt6,pt7 + 0.02*nvec
      SRPX = (TPLOCX(5)+TPLOCX(6)+TPLOCX(7)) / 3 + 0.02 * NVECX
      SRPY = (TPLOCY(5)+TPLOCY(6)+TPLOCY(7)) / 3 + 0.02 * NVECY
      SRPZ = (TPLOCZ(5)+TPLOCZ(6)+TPLOCZ(7)) / 3 + 0.02 * NVECZ
C  Note: the particular values chosen for FALOC preserve the "hole"
C  at p5-7, but since M1 is a distorting transformation, this will
C  not always be the case with other values for FALOC.

C  sdist = distance from SRP to fill area, as computed by PTREGD
      SDIST = PTREGD (SRPX,SRPY,SRPZ, FALEN, TPLOCX,TPLOCY,TPLOCZ)

      CALL SETMSG ('2 3 4 5 6 10 13', 'ISS should detect a complex ' //
     1             'primitive according to its location in world '   //
     2             'coordinates.')

C  ISS twice with SRP, NOCLIP, varying search distance
C    starting path:  101,0
C    search ceiling: 1
C    expected path:  101,2

      CALL ISSAB (SRPX,SRPY,SRPZ, SDIST, '101,0', PNCLIP, 1, '101,2')

C  *** *** *** ***   Complex clipping   *** *** *** ***
C
C  set up vertices for fill area:
C
C        1----------2     p1: 3,2
C       /           /     p2: 6,2
C      /           /      p3: 4,-1
C     4           /       p4: 1,1
C      \         /
C       \       /
C        \     /
C         \   /
C          \ /
C           3

      CALL SETRVS ('3,6, 4,1', FALOCX, FALEN)
      CALL SETRVS ('2,2,-1,1', FALOCY, FALEN)
      CALL SETRVS ('0,0, 0,0', FALOCZ, FALEN)

C  tp1-tp4 = WC vertices = m1 X p1-p4
      CALL ARRTP3 (FALEN, FALOCX, FALOCY, FALOCZ, M1,
     1                    TPLOCX, TPLOCY, TPLOCZ)

C  clpt  = WC clip point  = point on tp2:tp4 line segment closest to tp1
C    (tp1 normal projection onto line segment)
      CALL PTLNDS ( TPLOCX(1),TPLOCY(1),TPLOCZ(1),
     1  TPLOCX(2)-TPLOCX(4), TPLOCY(2)-TPLOCY(4), TPLOCZ(2)-TPLOCZ(4),
     2  TPLOCX(4),TPLOCY(4),TPLOCZ(4), CLPTX,CLPTY,CLPTZ, RDUM)

C  clvec = WC clip vector = vector from clpt to tp1 (in WC plane of
C     fill area, normal to line between tp2 and tp4, towards tp1).
      CLVCX = TPLOCX(1) - CLPTX
      CLVCY = TPLOCY(1) - CLPTY
      CLVCZ = TPLOCZ(1) - CLPTZ

C  nvec  = unit vector normal to WC fill area plane
      CALL PTSPL (FALEN, TPLOCX,TPLOCY,TPLOCZ, TX,TY,TZ, RDUM)
      CALL VEC1  (TX,TY,TZ, NVECX,NVECY,NVECZ)

C  set up structure #102:
C  1. set modelling clipping volume with:
C        point         = clpt
C        normal vector = clvec
C  2. set global transformation = m1
C  3. fill area at p1-p4
      CALL POPST  (102)

      HALFSP (1,1) = CLPTX
      HALFSP (2,1) = CLPTY
      HALFSP (3,1) = CLPTZ
      HALFSP (4,1) = CLVCX
      HALFSP (5,1) = CLVCY
      HALFSP (6,1) = CLVCZ
      CALL PSMCV3 (PMCREP, 1, HALFSP)

      CALL PSGMT3 (M1)
      CALL PFA    (FALEN, FALOCX,FALOCY)
      CALL PCLST

C  dvec = vector away from clipping volume and up from fill area plane
      NVECL = VECL(CLVCX, CLVCY, CLVCZ) / 10
      DVECX = NVECL*NVECX - CLVCX/2
      DVECY = NVECL*NVECY - CLVCY/2
      DVECZ = NVECL*NVECZ - CLVCZ/2

C  SRP = clpt + dvec (outside clipping volume)
      SRPX = CLPTX + DVECX
      SRPY = CLPTY + DVECY
      SRPZ = CLPTZ + DVECZ

C  sdist = length of dvec
      SDIST = VECL (DVECX,DVECY,DVECZ)
C  Note: this is not a general-purpose formula - only works when
C  clip point is inside transformed quadrilateral.  Since M1 is a
C  distorting transformation, this will not always be the case,
C  although it is for the particular values chosen for FALOC.

      CALL SETMSG ('2 3 4 5 6 10 13', 'If the modelling clip flag '   //
     1             'is CLIP, ISS should detect only those portions '  //
     2             'of a complex primitive whose location in world '  //
     3             'coordinates is inside the current modelling '     //
     4             'clipping volume.')

C  ISS twice with SRP, CLIP, varying search distance
C    starting path:  102,0
C    search ceiling: 1
C    expected path:  102,3
      CALL ISSAB (SRPX,SRPY,SRPZ, SDIST, '102,0', PCLIP, 1, '102,3')

666   CONTINUE
      CALL ENDIT
      END


text:
!*******************************************************************
!                                                                  *
      SUBROUTINE W1JJG(K1,QM1,QM2,IK,BK,ID,BD,WW)
!                                                                  *
!   Written by  G. Gaigalas                                        *
!   Transform to fortran 90/95 by G. Gaigalas       December 2012  *
!   The last modification made by G. Gaigalas       October  2017  *
!                                                                  *
!*******************************************************************
!
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY: DOUBLE
      USE CONS_C,          ONLY: ZERO, TENTH, HALF, TWO, EPS
      USE ribojj_C
      USE ribojj9_C
      USE ribojj11_C
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE mes_I
      USE rumtjj_I
      USE ixjtik_I
      USE ittk_I
      USE a1jj_I
      USE sixj_I
      USE izas1_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER,      INTENT(IN)               :: K1
      INTEGER,      INTENT(IN), DIMENSION(7) :: IK, ID
      REAL(DOUBLE), INTENT(IN)               :: QM1, QM2
      REAL(DOUBLE), INTENT(IN), DIMENSION(3) :: BK, BD
      REAL(DOUBLE), INTENT(OUT)              :: WW
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: KK1,IQ,IQM,IE,IT,ITK,ITD,ITP,ITP1,ITG,ITG1,IBTT
      INTEGER,      DIMENSION(7) :: IBT
      REAL(DOUBLE)               :: ENQP, D1, SI1, W
      REAL(DOUBLE), DIMENSION(3) :: BT
!-----------------------------------------------
      WW=ZERO
      IF(IZAS1(ID(7),BD(3),IK(7),BK(3)) == 0)RETURN
      ENQP=ZERO
      KK1=K1*2
      IQ=QM2*TWO+QM2*TENTH
      IF(ID(3) > 37) RETURN
      IF(ITTK(ID(6),IK(6),KK1) == 0)RETURN
      ITK=IK(1)
      ITD=ID(1)
      IF(ID(3) == 9) THEN
        IF(ID(4) > 2) CALL MES(1)
        IF(IK(4) > 2) CALL MES(1)
        ITK=ITK-300
        ITD=ITD-300
        ITP1=IMPNJJ9(ITK)
        ITP=IMPNJJ9(ITD)
        IF(ITP1 /= ITP)RETURN
        ITG1=IMGNJJ9(ITK)
        ITG=IMGNJJ9(ITD)
       ELSE
        IF(ID(4) > 2) CALL MES(1)
        IF(IK(4) > 2) CALL MES(1)
        ITP1=IMPNJJ11(ITK)
        ITP=IMPNJJ11(ITD)
        IF(ITP1 /= ITP)RETURN
        ITG1=IMGNJJ11(ITK)
        ITG=IMGNJJ11(ITD)
      ENDIF
      IF(ITG1 /= ITG)RETURN
      IBT(2)=ID(2)
      IBT(3)=ID(3)
      IBT(4)=ID(4)+IQ
      BT(3)=BD(3)+HALF*DBLE(IQ)
      IQM=TWO*DABS(BT(3))+TENTH
      DO IT=ITP,ITG
        CALL RUMTJJ(IT,IBT(3),IBT(7),IBTT,IBT(6))
        IF(IQM <= IBT(7)) THEN
          IF(IXJTIK(IK(3),IK(3),KK1,ID(6),IK(6),IBT(6)) /= 0) THEN
            IBT(1)=IT
            BT(2)=DBLE(IBT(6))/TWO
            BT(1)=DBLE(IBT(7))/TWO
            CALL A1JJ(IK,BK,IBT,BT,QM1,D1)
            IF(DABS(D1) > EPS) THEN
              CALL A1JJ(IBT,BT,ID,BD,QM2,W)
              IF(DABS(W) > EPS) THEN
                D1=D1*W
                CALL SIXJ(IK(3),IK(3),KK1,ID(6),IK(6),IBT(6),0,SI1)
                ENQP=ENQP+D1*SI1
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      END DO
      WW=ENQP*DSQRT(DBLE((KK1+1)))
      IE=KK1+IK(6)+ID(6)
      IF(((IE/4)*4) /= IE)WW=-WW
      RETURN
      END SUBROUTINE W1JJG

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text: The Moss Room is one of the Meeting Rooms conference rooms in the Memorial Union. Look for on on the west side of the third floor. It is named for Ralph H. Moss.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*> \brief \b SSYTRS_AA_2STAGE
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download SSYTRS_AA_2STAGE + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*      SUBROUTINE SSYTRS_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB, IPIV,
*                                   IPIV2, B, LDB, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          UPLO
*       INTEGER            N, NRHS, LDA, LTB, LDB, INFO
*       ..
*       .. Array Arguments ..
*       INTEGER            IPIV( * ), IPIV2( * )
*       REAL               A( LDA, * ), TB( * ), B( LDB, * )
*       ..
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SSYTRS_AA_2STAGE solves a system of linear equations A*X = B with a real
*> symmetric matrix A using the factorization A = U*T*U**T or
*> A = L*T*L**T computed by SSYTRF_AA_2STAGE.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*>          Specifies whether the details of the factorization are stored
*>          as an upper or lower triangular matrix.
*>          = 'U':  Upper triangular, form is A = U*T*U**T;
*>          = 'L':  Lower triangular, form is A = L*T*L**T.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in] NRHS
*> \verbatim
*>          NRHS is INTEGER
*>          The number of right hand sides, i.e., the number of columns
*>          of the matrix B.  NRHS >= 0.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*>          A is REAL array, dimension (LDA,N)
*>          Details of factors computed by SSYTRF_AA_2STAGE.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,N).
*> \endverbatim
*>
*> \param[out] TB
*> \verbatim
*>          TB is REAL array, dimension (LTB)
*>          Details of factors computed by SSYTRF_AA_2STAGE.
*> \endverbatim
*>
*> \param[in] LTB
*> \verbatim
*>          The size of the array TB. LTB >= 4*N.
*> \endverbatim
*>
*> \param[in] IPIV
*> \verbatim
*>          IPIV is INTEGER array, dimension (N)
*>          Details of the interchanges as computed by
*>          SSYTRF_AA_2STAGE.
*> \endverbatim
*>
*> \param[in] IPIV2
*> \verbatim
*>          IPIV2 is INTEGER array, dimension (N)
*>          Details of the interchanges as computed by
*>          SSYTRF_AA_2STAGE.
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*>          B is REAL array, dimension (LDB,NRHS)
*>          On entry, the right hand side matrix B.
*>          On exit, the solution matrix X.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the array B.  LDB >= max(1,N).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit
*>          < 0:  if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2017
*
*> \ingroup realSYcomputational
*
*  =====================================================================
      SUBROUTINE SSYTRS_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB,
     $                             IPIV, IPIV2, B, LDB, INFO )
*
*  -- LAPACK computational routine (version 3.8.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2017
*
      IMPLICIT NONE
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            N, NRHS, LDA, LTB, LDB, INFO
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * ), IPIV2( * )
      REAL               A( LDA, * ), TB( * ), B( LDB, * )
*     ..
*
*  =====================================================================
*
      REAL               ONE
      PARAMETER          ( ONE = 1.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            LDTB, NB
      LOGICAL            UPPER
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. External Subroutines ..
      EXTERNAL           SGBTRS, SLASWP, STRSM, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX
*     ..
*     .. Executable Statements ..
*
      INFO = 0
      UPPER = LSAME( UPLO, 'U' )
      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LTB.LT.( 4*N ) ) THEN
         INFO = -7
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -11
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SSYTRS_AA_2STAGE', -INFO )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( N.EQ.0 .OR. NRHS.EQ.0 )
     $   RETURN
*
*     Read NB and compute LDTB
*
      NB = INT( TB( 1 ) )
      LDTB = LTB/N
*
      IF( UPPER ) THEN
*
*        Solve A*X = B, where A = U*T*U**T.
*
         IF( N.GT.NB ) THEN
*
*           Pivot, P**T * B
*
            CALL SLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 )
*
*           Compute (U**T \P**T * B) -> B    [ (U**T \P**T * B) ]
*
            CALL STRSM( 'L', 'U', 'T', 'U', N-NB, NRHS, ONE, A(1, NB+1),
     $                 LDA, B(NB+1, 1), LDB)
*
         END IF
*
*        Compute T \ B -> B   [ T \ (U**T \P**T * B) ]
*
         CALL SGBTRS( 'N', N, NB, NB, NRHS, TB, LDTB, IPIV2, B, LDB,
     $               INFO)
         IF( N.GT.NB ) THEN
*
*           Compute (U \ B) -> B   [ U \ (T \ (U**T \P**T * B) ) ]
*
            CALL STRSM( 'L', 'U', 'N', 'U', N-NB, NRHS, ONE, A(1, NB+1),
     $                  LDA, B(NB+1, 1), LDB)
*
*           Pivot, P * B  [ P * (U \ (T \ (U**T \P**T * B) )) ]
*
            CALL SLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 )
*
         END IF
*
      ELSE
*
*        Solve A*X = B, where A = L*T*L**T.
*
         IF( N.GT.NB ) THEN
*
*           Pivot, P**T * B
*
            CALL SLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 )
*
*           Compute (L \P**T * B) -> B    [ (L \P**T * B) ]
*
            CALL STRSM( 'L', 'L', 'N', 'U', N-NB, NRHS, ONE, A(NB+1, 1),
     $                 LDA, B(NB+1, 1), LDB)
*
         END IF
*
*        Compute T \ B -> B   [ T \ (L \P**T * B) ]
*
         CALL SGBTRS( 'N', N, NB, NB, NRHS, TB, LDTB, IPIV2, B, LDB,
     $               INFO)
         IF( N.GT.NB ) THEN
*
*           Compute (L**T \ B) -> B   [ L**T \ (T \ (L \P**T * B) ) ]
*
            CALL STRSM( 'L', 'L', 'T', 'U', N-NB, NRHS, ONE, A(NB+1, 1),
     $                  LDA, B(NB+1, 1), LDB)
*
*           Pivot, P * B  [ P * (L**T \ (T \ (L \P**T * B) )) ]
*
            CALL SLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 )
*
         END IF
      END IF
*
      RETURN
*
*     End of SSYTRS_AA_2STAGE
*
      END

max_stars_repo_name: kortschak/OpenBLAS

max_issues_repo_name: kortschak/OpenBLAS

max_forks_repo_name: kortschak/OpenBLAS


text:
C
C
C  Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.
C
C  This program is free software; you can redistribute it and/or modify it
C  under the terms of version 2.1 of the GNU Lesser General Public License
C  as published by the Free Software Foundation.
C
C  This program is distributed in the hope that it would be useful, but
C  WITHOUT ANY WARRANTY; without even the implied warranty of
C  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
C
C  Further, this software is distributed without any warranty that it is
C  free of the rightful claim of any third person regarding infringement
C  or the like.  Any license provided herein, whether implied or
C  otherwise, applies only to this software file.  Patent licenses, if
C  any, provided herein do not apply to combinations of this program with
C  other software, or any other product whatsoever.
C
C  You should have received a copy of the GNU Lesser General Public
C  License along with this program; if not, write the Free Software
C  Foundation, Inc., 59 Temple Place - Suite 330, Boston MA 02111-1307,
C  USA.
C
C  Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
C  Mountain View, CA 94043, or:
C
C  http://www.sgi.com
C
C  For further information regarding this notice, see:
C
C  http://oss.sgi.com/projects/GenInfo/NoticeExplan
C
C

!DIR$ ID "@(#) libu/multi/posix/mactlstf.f	92.0	10/08/98 14:57:41"

!DIR$ FREE

!========================================
!     _MACT_DO_TSKLIST
!
!     This routine reports the status of all active macrotasks.  It
!     must be compiled by the CraySoft Fortran-90 compiler, so that
!     the output it produces is properly synchronized with any other
!     output produced by the program.
!
!     This routine is derived from the original TSKLIST from the PVP
!     environment.
!
!     This source file must be compiled with the f90 -eQ switch, to
!     allow identifiers with leading underbars to be acceptable.  This
!     means that the external symbol for _MACT_DO_TSKLIST, in C terms,
!     is "_MACT_DO_TSKLIST".
!
      SUBROUTINE _MACT_DO_TSKLIST(NUM_TASKS,                      &
                                  MTIDS, STATES, STIDS, UTIDS,    &
                                  WAIT_OBJS, WAIT_COUNTS          &
                                 )
      INTEGER               NUM_TASKS
      INTEGER               MTIDS      (NUM_TASKS)
      INTEGER               STATES     (NUM_TASKS)
      INTEGER               STIDS      (NUM_TASKS)
      INTEGER               UTIDS      (NUM_TASKS)
      INTEGER               WAIT_OBJS  (NUM_TASKS)
      INTEGER               WAIT_COUNTS(NUM_TASKS)

      INTEGER, PARAMETER :: TS_NONEXISTENT = 0
      INTEGER, PARAMETER :: TS_READY       = 1
      INTEGER, PARAMETER :: TS_RUNNING     = 2
      INTEGER, PARAMETER :: TS_WAITLOCK    = 3
      INTEGER, PARAMETER :: TS_WAITEV      = 4
      INTEGER, PARAMETER :: TS_WAITBAR     = 5
      INTEGER, PARAMETER :: TS_WAITTASK    = 6
      INTEGER, PARAMETER :: TS_COMPLETE    = 7

      INTEGER               OUTFILE

 9000 FORMAT(//' STATUS OF EXISTING TASKS',                          &
             //' Internal  User defined task val',38X,'# tasks'      &
              /' task ID   Hexadecimal     ASCII    Task status',    &
               23X,'waiting'                                         &
              /' -------   ---------------------    -----------',    &
               23X,'-------'                                         &
            )
 9001 FORMAT(1X,I6,4X,Z8,9X,A4,4X,'ready to run',              23X,'      ')
 9002 FORMAT(1X,I6,4X,Z8,9X,A4,4X,'running on logical CPU ',I5, 6X,'      ')
 9003 FORMAT(1X,I6,4X,Z8,9X,A4,4X,'waiting for lock at ',Z8,    7X,I6)
 9004 FORMAT(1X,I6,4X,Z8,9X,A4,4X,'waiting for event at ',Z8,   6X,I6)
 9005 FORMAT(1X,I6,4X,Z8,9X,A4,4X,'waiting for barrier at ',Z8, 4X,I6)
 9006 FORMAT(1X,I6,4X,Z8,9X,A4,4X,'waiting for task ',I6,      12X,I6)
 9007 FORMAT(1X,I6,4X,Z8,9X,A4,4X,'complete',                  27X,'      ')
 9008 FORMAT(1X,I6,4X,Z8,9X,A4,4X,'unknown',                   28X,'      ')

      OUTFILE = 6

!
!     Write information about each task.
!
      WRITE(OUTFILE,9000)
      DO I = 1, NUM_TASKS
        SELECT CASE (STATES(I))
          CASE (TS_READY)
               WRITE(OUTFILE, 9001) MTIDS(I), UTIDS(I), FIX_UTID(UTIDS(I))

          CASE (TS_RUNNING)
               WRITE(OUTFILE, 9002) MTIDS(I), UTIDS(I), FIX_UTID(UTIDS(I)),   &
                                    WAIT_OBJS(I)

          CASE (TS_WAITLOCK)
               WRITE(OUTFILE, 9003) MTIDS(I), UTIDS(I), FIX_UTID(UTIDS(I)),   &
                                    WAIT_OBJS(I), WAIT_COUNTS(I)

          CASE (TS_WAITEV)
               WRITE(OUTFILE, 9004) MTIDS(I), UTIDS(I), FIX_UTID(UTIDS(I)),   &
                                    WAIT_OBJS(I), WAIT_COUNTS(I)

          CASE (TS_WAITBAR)
               WRITE(OUTFILE, 9005) MTIDS(I), UTIDS(I), FIX_UTID(UTIDS(I)),   &
                                    WAIT_OBJS(I), WAIT_COUNTS(I)

          CASE (TS_WAITTASK)
               WRITE(OUTFILE, 9006) MTIDS(I), UTIDS(I), FIX_UTID(UTIDS(I)),   &
                                    WAIT_OBJS(I), WAIT_COUNTS(I)

          CASE (TS_COMPLETE)
               WRITE(OUTFILE, 9007) MTIDS(I), UTIDS(I), FIX_UTID(UTIDS(I))

          CASE DEFAULT
               WRITE(OUTFILE, 9008) MTIDS(I), UTIDS(I), FIX_UTID(UTIDS(I))
        END SELECT
      ENDDO

      RETURN

!     ==============================================================
      CONTAINS

        CHARACTER(4) FUNCTION FIX_UTID(UTID)
        INTEGER      UTID
        CHARACTER(4) A_UTID
        INTEGER      I

        INTRINSIC    TRANSFER
        INTRINSIC    ICHAR

        A_UTID = TRANSFER(UTID, A_UTID)
        DO I = 1, 4
          IF (ICHAR(A_UTID(I:I)) .LE. 31) THEN
            A_UTID(I:I) = '.'
          ENDIF
        ENDDO

        FIX_UTID = A_UTID

        END FUNCTION FIX_UTID

      END SUBROUTINE _MACT_DO_TSKLIST

max_stars_repo_name: sharugupta/OpenUH

max_issues_repo_name: sharugupta/OpenUH

max_forks_repo_name: sharugupta/OpenUH


text:
      integer function hepchg(id)

C...Purpose: to give three times the charge for a particle/parton.

C      ID     = particle ID
C      hepchg = particle charge times 3

      integer id,kqa,kq1,kq2,kq3,kqj,irt,kqx

      integer ichg(109)
      data ichg/-1,2,-1,2,-1,2,-1,2,2*0,
     1 -3,0,-3,0,-3,0,-3,3*0,
     2 3*0,3,6*0,
     3 3*0,3,0,0,3,3*0,
     4 10*0,
     5 0,6,3,6,6*0,
     6 49*0/
      save ichg

C...Initial values. Simple case of direct readout.
      hepchg=0

      kqa=iabs(id)
      kqn=mod(kqa/1000000000,10)
      kqx=mod(kqa/1000000,10)
      kq3=mod(kqa/1000,10)
      kq2=mod(kqa/100,10)
      kq1=mod(kqa/10,10)
      kqj=mod(kqa,10)
      irt = mod(kqa,10000)

C...illegal or ion
      if(kqa.eq.0.or.kqa.ge.10000000) then
C...set ion charge to zero - not enough information
        if(kqn.eq.1) hepchg = 0
C... direct translation
      elseif(kqa.le.100) then
        hepchg = ichg(kqa)
C... KS and KL (and undefined)
      elseif(kqj.eq.0) then
        hepchg = 0
C... direct translation
      elseif(kqx.gt.0 .and. irt.le.100) then
        hepchg = ichg(irt)
C... not everything maps, therefore deal with special cases
        if(kqa.eq.1000017 .or. kqa.eq.1000018) hepchg = 0
        if(kqa.eq.1000034 .or. kqa.eq.1000052) hepchg = 0
        if(kqa.eq.1000053 .or. kqa.eq.1000054) hepchg = 0
        if(kqa.eq.9900061 .or. kqa.eq.9900062) hepchg = 6
C... pentaquarks
      elseif(kqa.eq.9221132) then
        hepchg = 3
      elseif(kqa.eq.9331122) then
        hepchg = -6

C...Construction from quark content for heavy meson, diquark, baryon.
      elseif(kq3.eq.0) then
C...Mesons.
        hepchg = ichg(kq2) - ichg(kq1)
C...Strange or beauty mesons.
        if(kq2.eq.3 .or. kq2.eq.5) hepchg = ichg(kq1) - ichg(kq2)
      elseif(kq1.eq.0) then
C...Diquarks.
        hepchg = ichg(kq3) + ichg(kq2)
      else
C...Baryons
        hepchg = ichg(kq3) + ichg(kq2) + ichg(kq1)
      endif

C... fix sign of charge
      if(id.lt.0 .and. hepchg.ne.0) hepchg = -hepchg

      return
      end


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   EVB-QMDFF - RPMD molecular dynamics and rate constant calculations on
!               black-box generated potential energy surfaces
!
!   Copyright (c) 2021 by Julien Steffen (steffen@pctc.uni-kiel.de)
!                         Stefan Grimme (grimme@thch.uni-bonn.de) (QMDFF code)
!
!   Permission is hereby granted, free of charge, to any person obtaining a
!   copy of this software and associated documentation files (the "Software"),
!   to deal in the Software without restriction, including without limitation
!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!   and/or sell copies of the Software, and to permit persons to whom the
!   Software is furnished to do so, subject to the following conditions:
!
!   The above copyright notice and this permission notice shall be included in
!   all copies or substantial portions of the Software.
!
!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
!   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
!   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!   DEALINGS IN THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!     subroutine rdsolvff_three: read in informations for QMDFF3
!
!     part of QMDFF
!
subroutine rdsolvff_three(n,xyz,at,q,imass,d,scalehb,scalexb,fname)
use qmdff
implicit none
integer::n,at(n),imass(n)
real(kind=8)::xyz(3,n),q(n),d
real(kind=8)::scalehb(94),scalexb(94)
character(len=*)::fname

integer::i,j,nn,nterm,idum
real(kind=8)::qdum,r,xx(10),dum1,dum2
character(len=80)::atmp

open(unit=1,file=fname)

read (1,'(a)') atmp
call readl(atmp,xx,nn)
n=idint(xx(1))
qdum=xx(2)
d=xx(3)

if (nn.eq.3.and.d.lt.0) then
   d=xx(3)
   write(*,*) 'taking solvent density from file ',d
end if
read (1,'(a)') atmp

do i=1,n
   read (1,*)at(i),xyz(1:3,i),q(i),imass(i)
end do
read (1,*) nbond_three,nangl_three,ntors_three,nhb_three,nnci_three
if (nbond_three.gt.ndim.or.nangl_three.gt.ndim.or.ntors_three.gt.ndim) &
    &  stop 'too many FF terms'
do i=1,nbond_three
   read (1,*)bond_three(1:2,i),vbond_three(1:3,i)
end do
do i=1,nangl_two
   read (1,*)angl_three(1,i),angl_three(2,i),angl_three(3,i), &
     &     vangl_three(1:2,i)
enddo
do i=1,ntors_two
   read (1,*)tors_three(1:6,i),vtors_three(1:2,i), &
       &   (vtors_three(3*(j-1)+3,i), &
       &    vtors_three(3*(j-1)+4,i), &
       &    vtors_three(3*(j-1)+5,i),j=1,tors_three(5,i))
end do
do i=1,ntors_three
   do j=1,tors_three(5,i)
      vtors_three(3*(j-1)+4,i)=vtors_three(3*(j-1)+4,i)*pi
   end do
end do
if(nhb_three.gt.0) read(1,'(6(3i4,2x))')hb_three(1:3,1:nhb_three)
do i=1,nhb_three
   if (at(hb_three(3,i)).eq.1) then
!
!     subroutine hbpara is totally FF-unspecific!
!
      call hbpara(10.0d0,5.0d0,q(hb_three(1,i)),dum1)
      vhb_three(1,i)=dum1*scalehb(at(hb_three(1,i)))
      call hbpara(10.0d0,5.0d0,q(hb_three(2,i)),dum1)
      vhb_three(2,i)=dum1*scalehb(at(hb_three(2,i)))
   else
      call hbpara(-6.5d0,1.0d0,q(hb_three(3,i)),dum1)
      vhb_three(1,i)=scalexb(at(hb_three(3,i)))*dum1
   end if
end do
read(1,'(8(3i4,2x))')nci_three(1:3,1:nnci_three)
close(1)

q=q*qdum

return
end subroutine rdsolvff_three

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense


text:
c------------------------------------------------------------------------------
c Copyright (c) 2012, liisa holm
c All rights reserved.
c
c Redistribution and use in source and binary forms, with or without
c modification, are permitted provided that the following conditions are met:
c
c 1. Redistributions of source code must retain the above copyright notice, this
c    list of conditions and the following disclaimer.
c 2. Redistributions in binary form must reproduce the above copyright notice,
c    this list of conditions and the following disclaimer in the documentation
c    and/or other materials provided with the distribution.
c
c THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
c ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
c WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
c DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
c ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
c (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
c LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
c ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
c (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
c SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
c------------------------------------------------------------------------------
c       in: uniprot.ISA [8 bytes], uniprot.pin
c       out: uniprot.SAP [4 bytes] uniprot.SRES [2 bytes]
c
c       USAGE: kb7_31 uniprot.pin uniprot.ISA uniprot.SAP uniprot.sres
c
c       algo: SAP[ISA[seq]]=prot(seq)
c
c       gfortran -O3 sap8.f -o sap8
c
        program kb7_31
        implicit none
        INTEGER, PARAMETER :: long = SELECTED_INT_KIND(18)
        INTEGER, PARAMETER :: giga=1000000000 ! 1 G
        INTEGER(long), PARAMETER :: MEMORY=500_long*giga ! 500 G
        integer(long) :: seq,naa,offset,ptr,k,i,chunksize
        integer bufsize,ibuf,irec,seqlen,jrec,nprot
        integer, dimension(:), allocatable :: sap,sres
        parameter(bufsize=250000000) ! 250 M = 2 Gbytes
        character*256 sname,fn_pin,fn_ISA, fn_SAP,fn_sres
        integer(long) :: isabuffer(bufsize),lex_from,lex_to
        integer iprot,ires,ok,outbuf,outbuffer(bufsize),divisor
        integer*2 outbuffer_sres(bufsize)

        ! mandatory arguments
        if(iargc().lt.1) then
          write(*,*) 'sap8 generates protein labels for suffix array'
          write(*,*) ' '
          write(*,*) 'usage: sap8 database-name '
          stop
        end if
        call getarg(1,sname)
        fn_pin=trim(sname)//'.pin'
        fn_ISA=trim(sname)//'.ISA'
        fn_SAP=trim(sname)//'.SAP'
        fn_sres=trim(sname)//'.SRES'
        write(*,*) '# inputs: ',fn_pin,fn_ISA
        write(*,*) '# outputs: ',fn_SAP,fn_sres
        naa=0

        ! open input ISA
        call openfile(10,fn_ISA,8*bufsize)

        ! open output SAP
        call openfile(20,fn_SAP,4*bufsize)
        ! open output sres
        call openfile(40,fn_sres,2*bufsize)

        ! open input pin
        open(30,file=fn_pin,status='old')

        ! read uniprot.pin, fill SAP[ISA[seq]]=prot(seq)
        ! get naa
        nprot=0
10      read(30,*,end=19) ptr,seqlen
        nprot=nprot+1
        goto 10
19      continue
        naa=ptr+seqlen+1
        close(30)

        ! split database to fit MEMORY
        ! sap(chunksize)+isabuffer(bufsize)+outbuffer(bufsize).lt.MEMORY
        chunksize=naa
        divisor=1
        do while(chunksize*8_long+bufsize*14_long.gt.MEMORY)
                divisor=divisor+1
                chunksize=1+naa/divisor
        end do
        write(*,*) '# chunksize,naa,chunks',chunksize,naa,divisor
        if(chunksize*divisor.lt.naa) stop 'chunksize bug'

        ! dynamic memory allocation
        allocate(sap(chunksize),stat=ok)
        if(ok /= 0) stop 'sap memory allocation failed'
        allocate(sres(chunksize),stat=ok)
        if(ok /= 0) stop 'sres memory allocation failed'

        ! scan chunks
        lex_from=0
        lex_to=chunksize
        irec=0 ! output
        outbuf=0
	offset=0
        do while(lex_from.lt.naa)
          write(*,*) '#processing chunk',lex_from,lex_to,naa
          ! rewind
          open(30,file=fn_pin,status='old')
          jrec=0 ! input
          ibuf=bufsize
          iprot=0
          do iprot=1,nprot
20          read(30,*,end=29) ptr,seqlen
            seq=ptr
            !write(*,*) '#processing',iprot,ptr,seqlen
            do ires=1,seqlen+1 ! aaseq plus newline
                seq=seq+1
                ibuf=ibuf+1
                if(ibuf.gt.bufsize) then
                        jrec=jrec+1
                        k=min(bufsize,naa-offset)
	  write(*,*) '#to read',jrec,offset,k,offset+k,naa,iprot
                        read(10,rec=jrec,err=29) (isabuffer(i),i=1,k)
	                offset=offset+k
                        ibuf=1
                end if
                i=isabuffer(ibuf)-lex_from
                if(i.gt.0.and.i.le.chunksize) then
                        sap(i)=iprot
                        sres(i)=ires
          !if(iprot.eq.1) write(*,*) '#sap',iprot,ires,seq,i,i+lex_from
                end if
            end do
            if(mod(iprot,100000).eq.0)
     $          write(*,*) iprot,naa,ptr,seqlen,seq,float(seq)/naa
          end do
29        continue
          close(30)

          ! buffered output SAP; flush full records
          do i=1,chunksize
                if(outbuf.ge.bufsize) then ! flush
                  irec=irec+1
                  write(20,rec=irec) (outbuffer(k),k=1,bufsize)
                  write(40,rec=irec) (outbuffer_sres(k),k=1,bufsize)
                  outbuf=0
                end if
                outbuf=outbuf+1
                outbuffer(outbuf)=sap(i)
                outbuffer_sres(outbuf)=sres(i)
          end do
          lex_from=lex_from+chunksize
          lex_to=lex_to+chunksize
        end do

        deallocate(sap)

        ! output last record
        if(outbuf.gt.0) then
                irec=irec+1
                write(20,rec=irec) (outbuffer(k),k=1,outbuf)
                write(40,rec=irec) (outbuffer_sres(k),k=1,outbuf)
        end if

        close(10)
        close(20)
        close(40)

        end program kb7_31

c
c====================================================================
c
        subroutine openfile(unit10,fn,recordlength)
c
c       open binary file for direct access
c
        implicit none
        integer unit10,recordlength
        character*(*) fn

        open(unit=unit10,
     $          file=fn,
     $          form='unformatted',
     $          access='direct',
     $          recl=recordlength)

        return
        end
c
c======================================================================
c

max_stars_repo_name: nestorzaburannyi/annotate

max_issues_repo_name: nestorzaburannyi/annotate

max_forks_repo_name: nestorzaburannyi/annotate


text:
! reduction on a vector
      program redonvect
      integer n
      parameter(n=10)
      real*8 A(n,n), S(n)
      integer i,j

!hpf$ processors P(2,2)
!hpf$ distribute A onto P

!hpf$ independent(j,i)
      do j=1, n
         do i=1, n
            A(i,j) = n-i/(j+1)
         enddo
      enddo

!hpf$ independent
      do i=1, n
         S(i) = 0.0
      enddo

!hpf$ independent(j,i), reduction(S)
      do j=1, n
         do i=1, n
            S(i) = S(i) + A(i,j)
         enddo
      enddo

      print *, (S(i),i=1,n)

      end


text:
      subroutine gafscl(cdrsd,cscale,ndrsmx,ndrsrd,nst,nstmax)
c
c     This subroutine calculates the cscale array of affinity scaling
c     factors. Affinity scaling is used to help make decisions on
c     which phases are the best choices to precipitate when there
c     are multiple supersaturations.
c
c     This subroutine is called by:
c
c       EQ6/eq6.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       cdrsd  = array of reaction coefficients ('d' set)
c
c     Principal output:
c
c       cscale = array of affinity scaling factors
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer ndrsmx,nstmax
c
      integer ndrsrd(2,nstmax)
      integer nst
c
      real*8 cdrsd(ndrsmx),cscale(nstmax)
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer n,ns,nr1,nr2
c
      real*8 cx
c
c-----------------------------------------------------------------------
c
      do ns = 1,nst
        cx = 0.
        nr1 = ndrsrd(1,ns)
        nr2 = ndrsrd(2,ns)
        do n = nr1 + 1,nr2
          cx = cx + abs(cdrsd(n))
        enddo
        if (cx .le. 0.) cx = 1.0
        cscale(ns) = cx
      enddo
c
      end


text:
C (C) Copyright 1996-2016 ECMWF.
C
C This software is licensed under the terms of the Apache Licence Version 2.0
C which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
C In applying this licence, ECMWF does not waive the privileges and immunities
C granted to it by virtue of its status as an intergovernmental organisation nor
C does it submit to any jurisdiction.
C

      program  magics

C define the data for the 2 graphs
      real ymin(3), ymax(3)
      character*22 x(3)

C define the text
      character*20 lines(1)

      data ymin /12.,17.,15./
      data ymax /22.,21.,17./

      data x /"2011-03-02 00:00:00",
     +        "2011-03-02 12:00:00",
     +        "2011-03-03 00:00:00"/

      data lines /"My other Graph"/

C open magics
      call popen

C Setting of the output file name
      call psetc('output_name', 'date_graph')

C define the cartesian projection
      call psetc("subpage_map_projection", "cartesian")
	  ??? Position the subpage ???

C define horizontal axis
      call psetc("axis_orientation","horizontal")
      call psetc("axis_type","???")
      call psetc("axis_grid", "on")
	  call psetc("axis_grid_colour", "grey")
	  call pseti("axis_grid_thickness", 1)
      call p???("axis_date_min_value", "2011-03-01 12:00:00")
      call p???("???date_max???", "2011-03-03 12:00:00")
	  call psetc("axis_grid_line_style", "dot")
      call p???

C define horizontal axis
      call psetc("axis_orientation", "???")
        ???
      call psetc("axis_grid", "on")
      call psetc("axis_grid_colour", "grey")
      call pseti("axis_grid_thickness", 1)
      call psetc("axis_grid_line_style", "dot")
      call psetr("??min??", 10.)
      call psetr("???max??", 30.)
      call p???

C Define the blue gurve with the min values
      call psetc("graph_line_colour", "???")
      call pseti("graph_line_thickness", 8)
      call psetc("graph_symbol", "on")
      call psetc("legend_user_text",
     + "<font ???? > Min </font>")
      call pseti("graph_symbol_marker_index", 1)
      call psetr("graph_symbol_height", 0.5)
      call psetc("graph_symbol_colour", "black")
      call p???("graph_curve_date_x_values", x, 3)
	  call p???("graph_curve_y_values", ymin, 3)
      call pgraph

C Define the blue gurve with the max values
      call psetc("graph_line_colour", "red")
      call pseti("graph_line_thickness", 8)
      call psetc("graph_symbol", "on")
      call psetc("legend_user_text",
     + "<font colour='red'> Max </font>")
      call pseti("graph_symbol_marker_index", 1)
      call psetr("graph_symbol_height", 0.5)
      call psetc("graph_symbol_colour", "black")
      call pset1c("???date_x???", x, ???)
	  call pset1r("???y???", ymax, ???)
      call p???

C Define the text
      call pset1c("text_lines", lines, 1)
      call psetc("text_html", "true")
      call psetc("text_colour", "black")
      call psetr("text_font_size", 0.6)
      call psetc("text_mode", "positional")
      call psetr("text_box_x_position", 1.5)
      call psetr("text_box_y_position", 15.5)
      call psetr("text_box_x_length", 20.)
      call psetr("text_box_y_length", 2.5)
      call psetc("text_border", "off")
      call psetc("text_justification", "left")
      call p????

      call psetc("legend_text_colour", "black")

      call p???

      end


text:
program prototype
use, intrinsic :: iso_fortran_env, only : dp => REAL64, error_unit
use :: mpi
implicit none
integer, parameter :: MEM_ALLOC_ERR = 1
integer, parameter :: default_n = 5, nr_dims = 2, nr_iters = 500
real(kind=dp), parameter :: tol = 1.0e-4_dp
integer :: ierr, iter_nr, n, comm_rank, comm_size
real(Kind=dp), dimension(:, :), allocatable :: temp
character(len=1024) :: buffer
integer, dimension(nr_dims) :: dims
logical, dimension(nr_dims) :: periodic
logical :: reorder
logical, parameter :: is_verbose = .false.
integer :: cart_comm, row_type, col_type

! setup MPI and get size and initial rank
call MPI_Init(ierr)
call MPI_Comm_size(MPI_COMM_WORLD, comm_size, ierr)
call MPI_Comm_rank(MPI_COMM_WORLD, comm_rank, ierr)

! setup Carthesian communicator
dims = 0
periodic = .true.
reorder = .true.
call MPI_Dims_create(comm_size, nr_dims, dims, ierr)
call MPI_Cart_create(MPI_COMM_WORLD, nr_dims, dims, periodic, reorder, &
                     cart_comm, ierr)
call MPI_Comm_rank(cart_comm, comm_rank, ierr)

! rank 0 handles command line arguments, and broadcasts them to other
! ranks
if (comm_rank == 0) then
    if (command_argument_count() == 1) then
        call get_command_argument(1, buffer)
        read (buffer, '(I5)') n
    else
        n = default_n
    end if
end if
call MPI_Bcast(n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)

! create and register types
call create_types(row_type, col_type, n)

! create and initialize local arrays for temperature and old temperature
call allocate_matrix(temp, n)
call init_matrix(temp, cart_comm)

! show matrix
if (is_verbose) call print_Matrix(temp, cart_comm)

! allow diffusion until changes are smaller than tolerance, or the maximum
! number of iterations is exceeded
call halo_exchange(temp, cart_comm, row_type, col_type)

! show matrix
if (is_verbose) call print_Matrix(temp, cart_comm)

! take down MPI
call MPI_Finalize(ierr)

contains

! routine to allocate an nxn matrix
    subroutine allocate_matrix(matrix, n)
        implicit none
        real(kind=dp), dimension(:, :), allocatable, intent(inout) :: matrix
        integer, intent(in) :: n
        allocate(matrix(n, n), stat=ierr)
        if (ierr /= 0) then
            write (error_unit, '(A, I0)') &
                'can not allocate temp in ', comm_rank
            call MPI_Abort(MPI_COMM_WORLD, MEM_ALLOC_ERR, ierr)
        end if
    end subroutine allocate_matrix

! create MPI types for halo exchange
    subroutine create_types(row_type, col_type, n)
        implicit none
        integer, intent(inout) :: row_type, col_type
        integer, intent(in) :: n
        integer :: ierr
        call MPI_Type_vector(n, 1, n, MPI_DOUBLE_PRECISION, row_type, ierr)
        call MPI_Type_commit(row_type, ierr)
        call MPI_type_contiguous(n, MPI_DOUBLE_PRECISION, col_type, ierr)
        call MPI_Type_commit(col_type, ierr)
     end subroutine create_types

! initialize a matrix, mostly to zero, some elements randomly to a
! non-zero value
    subroutine init_matrix(matrix, comm)
        implicit none
        real(kind=dp), dimension(:, :), intent(inout) :: matrix
        integer, intent(in) :: comm
        integer :: i, j, n, ierr, offset, stride
        integer, dimension(nr_dims) :: comm_coords, comm_dims
        logical, dimension(nr_dims) :: periodic
        n = size(matrix, 1)
        call MPI_Cart_get(comm, nr_dims, comm_dims, periodic, &
                          comm_coords, ierr)
        stride = comm_dims(2)*n
        offset = comm_coords(1)*n + stride*comm_coords(2)*n
        do j = 1, size(matrix, 2)
            do i = 1, size(matrix, 1)
                matrix(i, j) = offset + j + (i - 1)*stride
            end do
        end do
    end subroutine init_matrix

    subroutine halo_exchange(matrix, comm, row_type, col_type)
        implicit none
        real(kind=dp), dimension(:, :), intent(inout) :: matrix
        integer, intent(in) :: comm, row_type, col_type
        integer, parameter :: tag = 17
        integer :: up, right, down, left, comm_rank, ierr, nr_rows, nr_cols
        nr_rows = size(matrix, 1)
        nr_cols = size(matrix, 2)
        call MPI_Comm_rank(comm, comm_rank, ierr)
        call MPI_Cart_shift(comm, 0, 1, left, right, ierr)
        call MPI_Cart_shift(comm, 1, 1, up, down, ierr)
        call MPI_Sendrecv(matrix(2, 1), 1, row_type, up, tag, &
                          matrix(1, 1), 1, row_type, down, tag, &
                          comm, MPI_STATUS_IGNORE, ierr)
        call MPI_Sendrecv(matrix(nr_rows - 1, 1), 1, row_type, down, tag, &
                          matrix(nr_rows, 1), 1, row_type, up, tag, &
                          comm, MPI_STATUS_IGNORE, ierr)
        call MPI_Sendrecv(matrix(1, 2), 1, col_type, left, tag, &
                          matrix(1, 1), 1, col_type, right, tag, &
                          comm, MPI_STATUS_IGNORE, ierr)
        call MPI_Sendrecv(matrix(1, nr_cols - 1), 1, col_type, right, tag, &
                          matrix(1, nr_cols), 1, col_type, left, tag, &
                          comm, MPI_STATUS_IGNORE, ierr)
    end subroutine halo_exchange

    subroutine print_matrix(matrix, comm)
        implicit none
        real(kind=dp), dimension(:, :), intent(in) :: matrix
        integer, intent(in) :: comm
        integer :: i, comm_rank, comm_size, ierr, rank
        call MPI_Comm_size(comm, comm_size, ierr)
        call MPI_Comm_rank(comm, comm_rank, ierr)
        do rank = 0, comm_size - 1
            if (comm_rank == rank) then
                print '(A, I0)', 'rank ', comm_rank
                do i = 1, size(matrix, 1)
                    print *, matrix(i, :)
                end do
            end if
            call MPI_Barrier(comm, ierr)
        end do
    end subroutine print_Matrix

end program prototype

max_stars_repo_name: Gjacquenot/training-material

max_issues_repo_name: Gjacquenot/training-material

max_forks_repo_name: Gjacquenot/training-material