text:
subroutine frequencyCalculation2(timesRepeated, generationNum, APlayerPayoff, BPlayerPayoff)
    use global
    implicit none
    integer,intent(in)::timesRepeated, generationNum
    real,intent(in)::APlayerPayoff(numPlayerInGame),BPlayerPayoff(numPlayerInGame)
    !runs (generations) in one repeat
    integer::run
    ! number of runs (repetitions)
    integer::indexRepeated
    !Output
    character(4)::timeRepeatFile
    character(5)::strEta
    !!!
    integer::numA
    !!! the distribution of the number of A players
    integer(kind=8)::numADist(Ntotal+1)
    !!!
    real::payoffUpdate
    !!
    integer, allocatable::nodeRandUpdate(:)
    integer::numNodeRandExceedNtotal, indexNodeRandExceedNtotal
    integer::iNodeRandExceed, nodeRandReDraw(1)
    !!
    integer::nodeUpdateNow
    !!
    real, allocatable::randNumUseUpdate(:)
    integer::numGenerationLoop, indexGenerationLoop
    !!!
    type(structuredPopulations)::networkMain(Ntotal)
    !!!
    !
    numGenerationLoop = 200

    !Testing
    !write(*,*) "For Finite Ring!"
    write(*,*) "For Well-mixed! Always Focal Game!"
    write(*,*) "Always 200 Generation Loops, Sample after 100 Loops"

    !!! Read networks
    !!! For structured populations, uncomment the next line
    !call readNetwork(networkMain)

    call RNSET(0)

    allocate(nodeRandUpdate(generationNum))
    allocate(randNumUseUpdate(generationNum))

    open(unit = 121,file = "numA_Strategy_Dist")

    timeRepeatedOneRealization:do indexRepeated=1, timesRepeated
        write(timeRepeatFile,"(I4)")  indexRepeated
        write(*,*) "-------------------"
        write(*,*) "Repeat = "//timeRepeatFile

        numADist = 0

        !-------Initialize strategies-------!
        call Inital_strategies_not_equal()
        !
        numA = sum(strategyPGG)
        !
        wholeGenerationLoop: do indexGenerationLoop = 1, numGenerationLoop
            !reset runs
            run = 1

            !! ------------------Speed up the simulation--------------------------!!!!
            nodeRandUpdate = -1.0
            !! All the random chosen node
            call RNUND(Ntotal, nodeRandUpdate)

            !!! Check the random chosen nodes
            if(any(nodeRandUpdate > Ntotal))then
                !write(*,*) "Updated Nodes exceed Population size -- ReDraw!!"
                numNodeRandExceedNtotal = count(nodeRandUpdate > Ntotal)
                !write(*,*) "numNodeRandExceedNtotal", numNodeRandExceedNtotal
                !!! redraw the nodes that exceed population size
                do iNodeRandExceed = 1, numNodeRandExceedNtotal
                    indexNodeRandExceedNtotal = maxloc(nodeRandUpdate, dim=1)
                    call RNUND(Ntotal, nodeRandReDraw)
                    nodeRandUpdate(indexNodeRandExceedNtotal) = nodeRandReDraw(1)
                end do
            end if
            if(any(nodeRandUpdate > Ntotal))then
                write(*,*) "Updated Nodes Still exceed Population size!!"
                pause
            end if

            !! All the random numbers used in the updating
            randNumUseUpdate = -1.0
            call RNUN(randNumUseUpdate)
            !! --------------------------------------------------------!!!!

            !!!
            generationLoop:do while(run <= generationNum)
                !!! ---------------------------- !!!
                !!! ------Game and Update------- !!!
                nodeUpdateNow = nodeRandUpdate(run)

                !!! For Regular graph, uncomment this line
                !call PGG_AroundUpdateNode_RegularGraph(APlayerPayoff,BPlayerPayoff,nodeUpdateNow,payoffUpdate, networkMain)

                !!! For well-mixed populations, uncomment this line
                call PGG_AroundUpdateNode_WellMixed(APlayerPayoff,BPlayerPayoff,nodeUpdateNow,payoffUpdate)

                !!! For both well-mixed & Regular graphs
                call AspirationUpdate(nodeUpdateNow, payoffUpdate, randNumUseUpdate(run))
                !!------------------------------!!
                numA = sum(strategyPGG)
                !!

                !!! Sample after a relaxtion time of 100 * generationNum time steps
                if(indexGenerationLoop > 100)then
                    numADist(numA + 1) = numADist(numA + 1) + 1
                end if

                run = run + 1
            end do generationLoop

        end do wholeGenerationLoop

        write(121,"(101(1X,I8))") numADist

    end do timeRepeatedOneRealization

    deallocate(nodeRandUpdate)
    deallocate(randNumUseUpdate)

    close(121)

    return
end subroutine frequencyCalculation2

max_stars_repo_name: leizhougetbetter/IAD

max_issues_repo_name: leizhougetbetter/IAD

max_forks_repo_name: leizhougetbetter/IAD


text:
!--- This is a driver function for photon isolation, the relevant photon isolation
!--- critera are selected and applied

      logical function iso(p,phot_id,isub,nd)
      implicit none
      include 'constants.f'
      include 'frag.f'
      include 'useet.f'
      include 'phot_dip.f'
      include 'z_dip.f'
      include 'lastphot.f'
      double precision p(mxpart,4)
      integer isub,nd,imode
      integer phot_id ! refers to which photon we are isolating
      logical first
      logical photo_iso_phys,photo_iso_z!,photo_iso,iso_old
      data first/.true./
      save first
!$omp threadprivate(first)

c--- imode: set the behaviour of the isolation requirement (default: imode=0)
c---   imode=0 : scaling cut for epsilon_h < 1, fixed cut for epsilon_h > 1
c---   imode=1 : scaling cut for all epsilon_h
c---   imode=2 : fixed cut for all epsilon_h
      imode=0

      iso=.true.

! Check if no isolation required
      if ((abs(cone_ang) .lt. 1d-4).or.(abs(epsilon_h) .lt. 1d-4)) then
        if (first) then
        write(6,*)'****************************************************'
        write(6,*)'*                                                  *'
        write(6,*)'*         No photon isolation cuts applied         *'
        write(6,*)'*                                                  *'
        write(6,*)'****************************************************'
        first=.false.
      endif
        return
      endif

c--- for imode=0, decide which cut to use depending on epsilon_h
!--   epsilon_h < 1 : epsilon_h corresponds to a pt fraction i.e. pt(had) < epsilon_h pt_gamma
!--   epsilon_h > 1 : treat it as an E_t max in cone   i.e. pt(had) < epsilon_h
       if (imode .eq. 0) then
         if (epsilon_h .lt. 0.9999d0) then
         imode=1
         else
         imode=2
         endif
       endif

!===== NEW ISOLATION : WORK OUT WHICH ROUTINE TO CALL

!========== are we doing fragmentation and is this the
!========== photon produced by fragmentation?
       if((fragint_mode) .and. (phot_id .eq. lastphot)) then
          iso=photo_iso_z(p,phot_id,z_frag,imode,isub)
!========== are we doing a photon dipole and is this the
!========== photon produced by fragmentation?
       elseif(phot_dip(nd) .and. (phot_id .eq. lastphot)) then
          iso=photo_iso_z(p,phot_id,z_dip(nd),imode,isub)
!========== we are in default mode, BORN, virt, real, non-photon dipole
       else
          iso=photo_iso_phys(p,phot_id,imode,isub)
       endif

!===== OLD ISOLATION
c       iso_old=photo_iso(p,isub,phot_dip(nd),phot_id,nd,imode)
c       if (iso .neqv. iso_old) then
c         write(6,*)'WARNING:',iso,'(new) vs',iso_old,' (old) for nd=',nd
c       else
c         write(6,*) 'OKAY'
c       endif

c--- write out isolation parameters
       if    (first .and. (imode .eq. 1)) then
        write(6,*)'************** Photons Isolated     ****************'
        write(6,*)'*                                                  *'
        write(6,99)'*    E_t(had) in cone',cone_ang,' < ',epsilon_h,
     &   ' E_t(phot)     *'
        write(6,*)'*                                                  *'
        write(6,*)'****************************************************'
        first=.false.
      elseif (first .and. (imode .eq. 2)) then
        write(6,*)'************** Photons Isolated     ****************'
        write(6,*)'*                                                  *'
        write(6,96)'* E_t (had) in cone',cone_ang,' < ',epsilon_h,
     &   'GeV    *'
        write(6,*)'*                                                  *'
        write(6,*)'****************************************************'
        first=.false.
      endif

!      endif

      return

 99   format(1x,a21,f6.2,a3,f6.2,a16)
 96   format(1x,a19,f6.2,a4,f6.2,a17)
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
	SUBROUTINE IM_ISUB  ( garea, iret )
C************************************************************************
C* IM_ISUB								*
C*									*
C* This subroutine subsets the current projection area for an image.	*
C* 									*
C*									*
C* IM_ISUB  ( GAREA, IRET )						*
C*									*
C* Input parameters:							*
C*	GAREA		CHAR*		Graphics area			*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -5 = invalid projection	*
C*					 -6 = invalid graphics area	*
C**									*
C* Log:									*
C* J. Cowie/COMET	 1/95	Original				*
C* S. Jacobs/NMC	 9/95	Fixed a typo yout(3) to yout(1) and	*
C*				  yout(4) to yout(2)			*
C* G. Krueger/EAI	 6/96	Add default projection			*
C* J. Cowie/COMET	 1/97	Changed IMGDEF common variable names	*
C* S. Jacobs/NCEP	11/98	Added check for MER proj to use GSMMAP	*
C* E. Safford/GSC	05/00	Removed check for MER proj		*
C* S. Jacobs/NCEP	 6/00	Increased imgnam from 80 to 256 chars	*
C* S. Jacobs/NCEP	 6/00	Changed to check projection set to SAT	*
C* S. Guan/NCEP          1/16   Modified for NETCDF4 Himawari data      *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'IMGDEF.CMN'
	INCLUDE		'AREAFL.CMN'
        INCLUDE         'SATDEF.CMN'
C*
	CHARACTER*(*)	garea
C*
	CHARACTER	nvtype*8, imgnam*256, loc*132, cproj*4,
     +			cdproj*30, goe4*4
	REAL		gltln(4), xin(2), xout(2), yin(2), yout(2)
	REAL		centrd (2)
        EQUIVALENCE     ( ianav(1),  goe4 )
C------------------------------------------------------------------------
C
	iret = 0
C
C*	Exit if graphics area is DSET
C
	CALL ST_LCUC ( garea, loc, ier)
	IF  ( loc .eq. 'DSET' ) RETURN
C
C*	Convert garea into lat/lon bounds
C
	CALL LC_GARE ( garea, gltln, cdproj, centrd, ier )
	IF ( ier .ne. 0 ) THEN
	    iret = -6
	    RETURN
	END IF
C
C*	Find view coords for these bounds
C
	xin (1) = gltln (1)
	xin (2) = gltln (3)
	yin (1) = gltln (2)
	yin (2) = gltln (4)
C
C       For NETCDF4 Himawari data
        IF ( goe4 .eq. 'GOE4' ) THEN
           yin (1) = (-1.0) *gltln (4)
           yin (2) = (-1.0) *gltln (2)
        END IF
C
	CALL GTRANS ( 'M', 'P', 2, xin, yin, xout, yout, ier )
C
C*	Check for bad garea plot bounds
C
	IF ( ( ier .ne. 0 ) .or. ( xout (1) .eq. RMISSD ) .or.
     +	     ( xout (2) .eq. RMISSD ) .or.
     +       ( yout (1) .eq. RMISSD ) .or.
     +       ( yout (2) .eq. RMISSD ) ) THEN
	    iret = -6
	    RETURN
	END IF
C
C*	Get plotting bounds
C
	CALL GQBND ( 'P', xl, yb, xr, yt, ier )
C
C*	Use full image size from common
C
	imxsiz = imnpix
	imysiz = imnlin
C
C*	Map plot coords to image coords
C
	imleft = NINT (imxsiz - (xr - xout (1)) /
     +                (xr - xl) * (imxsiz - 1))
	imtop = NINT (1 - (yt - yout (2)) / (yt - yb) * (1 - imysiz))
	imrght  = NINT (imxsiz - (xr - xout (2)) /
     +                (xr - xl) * (imxsiz - 1))
	imbot  = NINT (1 - (yt - yout (1)) / (yt - yb) * (1 - imysiz))
C
C*	Find out the map projection
C
 	CALL GQMPRJ ( cproj, angle1, angle2, angle3,
     +	   	      rlat1, rlon1, rlat2, rlon2, ier )
C
C*	If no map projection is set, it might be MCIDAS GOES or GVAR.
C
	IF ( cproj .eq. 'SAT' ) THEN
C
	    CALL GQSATN ( nvtype, imgnam, ier )
C
	    IF ( nvtype(1:2) .eq. 'MC') THEN
C
C              For NETCDF4 Himawari data
               IF ( goe4 .eq. 'GOE4' ) THEN
                  itmp = imleft
                  imleft = imrght
                  imrght = itmp
	          CALL GSATMG4 ( imgnam, iadir, ianav, imleft, imbot,
     +              imrght, imtop, ier )
               ELSE
                  CALL GSATMG ( imgnam, iadir, ianav, imleft, imtop,
     +                          imrght, imbot, ier )
               END IF

	    END IF
C
	ELSE
C
	    CALL GSMPRJ ( cproj, angle1, angle2, angle3,
     +			      gltln(1), gltln(2),
     +			      gltln(3), gltln(4), ier )
	END IF
C
C*	Check for projection error
C
	IF ( ier .ne. 0 ) THEN
	    iret = -5
	    RETURN
	END IF
C
C*	Check for image out of garea
C
	IF ( (imleft .gt. imxsiz) .or. (imrght .lt. 1) .or.
     +	     (imtop  .gt. imysiz) .or. (imbot  .lt. 1) ) THEN
	    CALL ER_WMSG ( 'IM', 2, ' ', ier )
	END IF

	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
*DECK DBVSUP
      SUBROUTINE DBVSUP (Y, NROWY, NCOMP, XPTS, NXPTS, A, NROWA, ALPHA,
     +   NIC, B, NROWB, BETA, NFC, IGOFX, RE, AE, IFLAG, WORK, NDW,
     +   IWORK, NDIW, NEQIVP)
C***BEGIN PROLOGUE  DBVSUP
C***PURPOSE  Solve a linear two-point boundary value problem using
C            superposition coupled with an orthonormalization procedure
C            and a variable-step integration scheme.
C***LIBRARY   SLATEC
C***CATEGORY  I1B1
C***TYPE      DOUBLE PRECISION (BVSUP-S, DBVSUP-D)
C***KEYWORDS  ORTHONORMALIZATION, SHOOTING,
C             TWO-POINT BOUNDARY VALUE PROBLEM
C***AUTHOR  Scott, M. R., (SNLA)
C           Watts, H. A., (SNLA)
C***DESCRIPTION
C
C **********************************************************************
C
C     Subroutine DBVSUP solves a linear two-point boundary-value problem
C     of the form
C                        DY/DX = MATRIX(X,U)*Y(X) + G(X,U)
C                A*Y(XINITIAL) = ALPHA ,  B*Y(XFINAL) = BETA
C
C     coupled with the solution of the initial value problem
C
C                        DU/DX = F(X,U)
C                      U(XINITIAL) = ETA
C
C **********************************************************************
C     ABSTRACT
C        The method of solution uses superposition coupled with an
C     orthonormalization procedure and a variable-step integration
C     scheme.  Each time the superposition solutions start to
C     lose their numerical linear independence, the vectors are
C     reorthonormalized before integration proceeds.  The underlying
C     principle of the algorithm is then to piece together the
C     intermediate (orthogonalized) solutions, defined on the various
C     subintervals, to obtain the desired solutions.
C
C **********************************************************************
C     INPUT to DBVSUP
C **********************************************************************
C
C     NROWY = actual row dimension of Y in calling program.
C             NROWY must be .GE. NCOMP
C
C     NCOMP = number of components per solution vector.
C             NCOMP is equal to number of original differential
C             equations.  NCOMP = NIC + NFC.
C
C     XPTS = desired output points for solution. They must be monotonic.
C            XINITIAL = XPTS(1)
C            XFINAL = XPTS(NXPTS)
C
C     NXPTS = number of output points.
C
C     A(NROWA,NCOMP) = boundary condition matrix at XINITIAL
C                      must be contained in (NIC,NCOMP) sub-matrix.
C
C     NROWA = actual row dimension of A in calling program,
C             NROWA must be .GE. NIC.
C
C     ALPHA(NIC+NEQIVP) = boundary conditions at XINITIAL.
C                         If NEQIVP .GT. 0 (see below), the boundary
C                         conditions at XINITIAL for the initial value
C                         equations must be stored starting in
C                         position (NIC + 1) of ALPHA.
C                         Thus,  ALPHA(NIC+K) = ETA(K).
C
C     NIC = number of boundary conditions at XINITIAL.
C
C     B(NROWB,NCOMP) = boundary condition matrix at XFINAL.
C                      Must be contained in (NFC,NCOMP) sub-matrix.
C
C     NROWB = actual row dimension of B in calling program,
C             NROWB must be .GE. NFC.
C
C     BETA(NFC) = boundary conditions at XFINAL.
C
C     NFC = number of boundary conditions at XFINAL.
C
C     IGOFX =0 -- The inhomogeneous term G(X) is identically zero.
C           =1 -- The inhomogeneous term G(X) is not identically zero.
C                 (if IGOFX=1, then Subroutine DGVEC (or DUVEC) must be
C                  supplied).
C
C     RE = relative error tolerance used by the integrator.
C          (see one of the integrators)
C
C     AE = absolute error tolerance used by the integrator.
C          (see one of the integrators)
C **NOTE-  RE and AE should not both be zero.
C
C     IFLAG = a status parameter used principally for output.
C             However, for efficient solution of problems which
C             are originally defined as COMPLEX*16 valued (but
C             converted to double precision systems to use this code),
C             the user must set IFLAG=13 on input. See the comment
C             below for more information on solving such problems.
C
C     WORK(NDW) = floating point array used for internal storage.
C
C     NDW = actual dimension of work array allocated by user.
C           An estimate for NDW can be computed from the following
C            NDW = 130 + NCOMP**2 * (6 + NXPTS/2 + expected number of
C                                           orthonormalizations/8)
C           For the disk or tape storage mode,
C            NDW = 6 * NCOMP**2 + 10 * NCOMP + 130
C  However, when the ADAMS integrator is to be used, the estimates are
C            NDW = 130 + NCOMP**2 * (13 + NXPTS/2 + expected number of
C                                           orthonormalizations/8)
C    and     NDW = 13 * NCOMP**2 + 22 * NCOMP + 130   , respectively.
C
C     IWORK(NDIW) = integer array used for internal storage.
C
C     NDIW = actual dimension of IWORK array allocated by user.
C            An estimate for NDIW can be computed from the following
C            NDIW = 68 + NCOMP * (1 + expected number of
C                                            orthonormalizations)
C **NOTE --  the amount of storage required is problem dependent and may
C            be difficult to predict in advance.  Experience has shown
C            that for most problems 20 or fewer orthonormalizations
C            should suffice. If the problem cannot be completed with the
C            allotted storage, then a message will be printed which
C            estimates the amount of storage necessary. In any case, the
C            user can examine the IWORK array for the actual storage
C            requirements, as described in the output information below.
C
C     NEQIVP = number of auxiliary initial value equations being added
C              to the boundary value problem.
C **NOTE -- Occasionally the coefficients  matrix  and/or  G  may be
C           functions which depend on the independent variable  X  and
C           on  U, the solution of an auxiliary initial value problem.
C           In order to avoid the difficulties associated with
C           interpolation, the auxiliary equations may be solved
C           simultaneously with the given boundary value problem.
C           This initial value problem may be linear or nonlinear.
C                 See SAND77-1328 for an example.
C
C
C     The user must supply subroutines DFMAT, DGVEC, DUIVP and DUVEC,
C     when needed (they must be so named), to evaluate the derivatives
C     as follows
C
C        A. DFMAT must be supplied.
C
C              SUBROUTINE DFMAT(X,Y,YP)
C              X = independent variable (input to DFMAT)
C              Y = dependent variable vector (input to DFMAT)
C              YP = DY/DX = derivative vector (output from DFMAT)
C
C            Compute the derivatives for the homogeneous problem
C              YP(I) = DY(I)/DX = MATRIX(X) * Y(I)  , I = 1,...,NCOMP
C
C            When (NEQIVP .GT. 0) and  matrix  is dependent on  U  as
C            well as on  X, the following common statement must be
C            included in DFMAT
C                    COMMON /DMLIVP/ NOFST
C            for convenience, the  U  vector is stored at the bottom
C            of the  Y  array.  Thus, during any call to DFMAT,
C            U(I) is referenced by  Y(NOFST + I).
C
C
C            Subroutine DBVDER calls DFMAT NFC times to evaluate the
C            homogeneous equations and, if necessary, it calls DFMAT
C            once in evaluating the particular solution. since X remains
C            unchanged in this sequence of calls it is possible to
C            realize considerable computational savings for complicated
C            and expensive evaluations of the matrix entries. To do this
C            the user merely passes a variable, say XS, via common where
C            XS is defined in the main program to be any value except
C            the initial X. Then the non-constant elements of matrix(x)
C            appearing in the differential equations need only be
C            computed if X is unequal to XS, whereupon XS is reset to X.
C
C
C        B. If  NEQIVP .GT. 0 ,  DUIVP must also be supplied.
C
C              SUBROUTINE DUIVP(X,U,UP)
C              X = independent variable (input to DUIVP)
C              U = dependent variable vector (input to DUIVP)
C              UP = DU/DX = derivative vector (output from DUIVP)
C
C            Compute the derivatives for the auxiliary initial value eqs
C              UP(I) = DU(I)/DX, I = 1,...,NEQIVP.
C
C            Subroutine DBVDER calls DUIVP once to evaluate the
C            derivatives for the auxiliary initial value equations.
C
C
C        C. If  NEQIVP = 0  and  IGOFX = 1 ,  DGVEC must be supplied.
C
C              SUBROUTINE DGVEC(X,G)
C              X = independent variable (input to DGVEC)
C              G = vector of inhomogeneous terms G(X) (output from
C              DGVEC)
C
C            Compute the inhomogeneous terms G(X)
C                G(I) = G(X) values for I = 1,...,NCOMP.
C
C            Subroutine DBVDER calls DGVEC in evaluating the particular
C            solution provided G(X) is not identically zero. Thus, when
C            IGOFX=0, the user need not write a DGVEC subroutine. Also,
C            the user does not have to bother with the computational
C            savings scheme for DGVEC as this is automatically achieved
C            via the DBVDER subroutine.
C
C
C        D. If  NEQIVP .GT. 0  and  IGOFX = 1 ,  DUVEC must be supplied.
C
C             SUBROUTINE DUVEC(X,U,G)
C             X = independent variable (input to DUVEC)
C             U = dependent variable vector from the auxiliary initial
C                 value problem    (input to DUVEC)
C             G = array of inhomogeneous terms G(X,U)(output from DUVEC)
C
C            Compute the inhomogeneous terms G(X,U)
C                G(I) = G(X,U) values for I = 1,...,NCOMP.
C
C            Subroutine DBVDER calls DUVEC in evaluating the particular
C            solution provided G(X,U) is not identically zero.  Thus,
C            when IGOFX=0, the user need not write a DUVEC subroutine.
C
C
C
C     The following is optional input to DBVSUP to give user more
C     flexibility in use of code.  See SAND75-0198, SAND77-1328,
C     SAND77-1690, SAND78-0522, and SAND78-1501 for more information.
C
C ****CAUTION -- The user must zero out IWORK(1),...,IWORK(15)
C                prior to calling DBVSUP. These locations define
C                optional input and must be zero unless set to special
C                values by the user as described below.
C
C     IWORK(1) -- number of orthonormalization points.
C                 A value need be set only if IWORK(11) = 1
C
C     IWORK(9) -- integrator and orthonormalization parameter
C                 (default value is 1)
C                 1 = RUNGE-KUTTA-FEHLBERG code using GRAM-SCHMIDT test.
C                 2 = ADAMS code using GRAM-SCHMIDT test.
C
C     IWORK(11) -- orthonormalization points parameter
C                  (default value is 0)
C                  0 - orthonormalization points not pre-assigned.
C                  1 - orthonormalization points pre-assigned in
C                      the first IWORK(1) positions of work.
C
C     IWORK(12) -- storage parameter
C                  (default value is 0)
C                  0 - all storage in core.
C                  LUN - homogeneous and inhomogeneous solutions at
C                      output points and orthonormalization information
C                      are stored on disk.  The logical unit number to
C                      be used for disk I/O (NTAPE) is set to IWORK(12).
C
C     WORK(1),... -- pre-assigned orthonormalization points, stored
C                    monotonically, corresponding to the direction
C                    of integration.
C
C
C
C                 ******************************************************
C                 *** COMPLEX*16 VALUED PROBLEM ***
C                 ******************************************************
C **NOTE***
C       Suppose the original boundary value problem is NC equations
C     of the form
C                   DW/DX = MAT(X,U)*W(X) + H(X,U)
C                 R*W(XINITIAL)=GAMMA , S*W(XFINAL)=DELTA
C     where all variables are COMPLEX*16 valued. The DBVSUP code can be
C     used by converting to a double precision system of size 2*NC. To
C     solve the larger dimensioned problem efficiently, the user must
C     initialize IFLAG=13 on input and order the vector components
C     according to Y(1)=DOUBLE PRECISION(W(1)),...,Y(NC)=DOUBLE
C     PRECISION(W(NC)),Y(NC+1)=IMAG(W(1)),...., Y(2*NC)=IMAG(W(NC)).
C     Then define
C                        ...............................................
C                        . DOUBLE PRECISION(MAT)    -IMAG(MAT) .
C            MATRIX  =   .                         .
C                        . IMAG(MAT)     DOUBLE PRECISION(MAT) .
C                        ...............................................
C
C     The matrices A,B and vectors G,ALPHA,BETA must be defined
C     similarly. Further details can be found in SAND78-1501.
C
C
C **********************************************************************
C     OUTPUT from DBVSUP
C **********************************************************************
C
C     Y(NROWY,NXPTS) = solution at specified output points.
C
C     IFLAG Output Values
C            =-5 algorithm ,for obtaining starting vectors for the
C                special COMPLEX*16 problem structure, was unable to
C                obtain the initial vectors satisfying the necessary
C                independence criteria.
C            =-4 rank of boundary condition matrix A is less than NIC,
C                as determined by DLSSUD.
C            =-2 invalid input parameters.
C            =-1 insufficient number of storage locations allocated for
C                WORK or IWORK.
C
C            =0 indicates successful solution.
C
C            =1 a computed solution is returned but uniqueness of the
C               solution of the boundary-value problem is questionable.
C               For an eigenvalue problem, this should be treated as a
C               successful execution since this is the expected mode
C               of return.
C            =2 a computed solution is returned but the existence of the
C               solution to the boundary-value problem is questionable.
C            =3 a nontrivial solution approximation is returned although
C               the boundary condition matrix B*Y(XFINAL) is found to be
C               nonsingular (to the desired accuracy level) while the
C               right hand side vector is zero. To eliminate this type
C               of return, the accuracy of the eigenvalue parameter
C               must be improved.
C            ***NOTE-We attempt to diagnose the correct problem behavior
C               and report possible difficulties by the appropriate
C               error flag.  However, the user should probably resolve
C               the problem using smaller error tolerances and/or
C               perturbations in the boundary conditions or other
C               parameters. This will often reveal the correct
C               interpretation for the problem posed.
C
C            =13 maximum number of orthonormalizations attained before
C                reaching XFINAL.
C            =20-flag from integrator (DDERKF or DDEABM) values can
C                range from 21 to 25.
C            =30 solution vectors form a dependent set.
C
C     WORK(1),...,WORK(IWORK(1)) = orthonormalization points
C                                  determined by DBVPOR.
C
C     IWORK(1) = number of orthonormalizations performed by DBVPOR.
C
C     IWORK(2) = maximum number of orthonormalizations allowed as
C                calculated from storage allocated by user.
C
C     IWORK(3),IWORK(4),IWORK(5),IWORK(6)   give information about
C                actual storage requirements for WORK and IWORK
C                arrays.  In particular,
C                       required storage for  work array is
C        IWORK(3) + IWORK(4)*(expected number of orthonormalizations)
C
C                       required storage for IWORK array is
C        IWORK(5) + IWORK(6)*(expected number of orthonormalizations)
C
C     IWORK(8) = final value of exponent parameter used in tolerance
C                test for orthonormalization.
C
C     IWORK(16) = number of independent vectors returned from DMGSBV.
C                It is only of interest when IFLAG=30 is obtained.
C
C     IWORK(17) = numerically estimated rank of the boundary
C                 condition matrix defined from B*Y(XFINAL)
C
C **********************************************************************
C
C     Necessary machine constants are defined in the Function
C     Routine D1MACH. The user must make sure that the values
C     set in D1MACH are relevant to the computer being used.
C
C **********************************************************************
C **********************************************************************
C
C***REFERENCES  M. R. Scott and H. A. Watts, SUPORT - a computer code
C                 for two-point boundary-value problems via
C                 orthonormalization, SIAM Journal of Numerical
C                 Analysis 14, (1977), pp. 40-70.
C               B. L. Darlow, M. R. Scott and H. A. Watts, Modifications
C                 of SUPORT, a linear boundary value problem solver
C                 Part I - pre-assigning orthonormalization points,
C                 auxiliary initial value problem, disk or tape storage,
C                 Report SAND77-1328, Sandia Laboratories, Albuquerque,
C                 New Mexico, 1977.
C               B. L. Darlow, M. R. Scott and H. A. Watts, Modifications
C                 of SUPORT, a linear boundary value problem solver
C                 Part II - inclusion of an Adams integrator, Report
C                 SAND77-1690, Sandia Laboratories, Albuquerque,
C                 New Mexico, 1977.
C               M. E. Lord and H. A. Watts, Modifications of SUPORT,
C                 a linear boundary value problem solver Part III -
C                 orthonormalization improvements, Report SAND78-0522,
C                 Sandia Laboratories, Albuquerque, New Mexico, 1978.
C               H. A. Watts, M. R. Scott and M. E. Lord, Computational
C                 solution of complex*16 valued boundary problems,
C                 Report SAND78-1501, Sandia Laboratories,
C                 Albuquerque, New Mexico, 1978.
C***ROUTINES CALLED  DEXBVP, DMACON, XERMSG
C***COMMON BLOCKS    DML15T, DML17B, DML18J, DML5MC, DML8SZ
C***REVISION HISTORY  (YYMMDD)
C   750601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   890921  Realigned order of variables in certain COMMON blocks.
C           (WRB)
C   890921  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900510  Convert XERRWV calls to XERMSG calls, remove some extraneous
C           comments.  (RWC)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  DBVSUP
C **********************************************************************
C
      INTEGER ICOCO, IFLAG, IGOFX, IGOFXD, INDPVT, INFO, INHOMO, INTEG,
     1     IS, ISTKOP, IVP, IWORK(*), J, K, K1, K10, K11, K2,
     2     K3, K4, K5, K6, K7, K8, K9, KKKCOE, KKKCOF, KKKG, KKKINT,
     3     KKKS, KKKSTO, KKKSUD, KKKSVC, KKKU, KKKV, KKKWS, KKKYHP,
     4     KKKZPW, KNSWOT, KOP, KPTS, L1, L2, LLLCOF, LLLINT, LLLIP,
     5     LLLIWS, LLLSUD, LLLSVC, LOTJP, LPAR, MNSWOT,
     6     MXNON, MXNONI, MXNONR, NCOMP, NCOMPD, NDEQ, NDISK, NDIW,
     7     NDW, NEEDIW, NEEDW, NEQ, NEQIVD, NEQIVP, NFC, NFCC,
     8     NFCD, NIC, NICD, NITEMP, NON, NOPG, NPS, NROWA, NROWB,
     9     NROWY, NRTEMP, NSWOT, NTAPE, NTP, NUMORT, NXPTS, NXPTSD,
     1     NXPTSM
      DOUBLE PRECISION A(NROWA,*), AE, AED, ALPHA(*),
     1     B(NROWB,*), BETA(*), C, EPS, FOURU, PWCND, PX, RE,
     2     RED, SQOVFL, SRU, TND, TOL, TWOU, URO, WORK(NDW), X, XBEG,
     3     XEND, XOP, XOT, XPTS(*), XSAV, Y(NROWY,*)
      CHARACTER*8 XERN1, XERN2, XERN3, XERN4
C
C     ******************************************************************
C         THE COMMON BLOCK BELOW IS USED TO COMMUNICATE WITH SUBROUTINE
C         DBVDER.  THE USER SHOULD NOT ALTER OR USE THIS COMMON BLOCK IN
C         THE CALLING PROGRAM.
C
      COMMON /DML8SZ/ C,XSAV,IGOFXD,INHOMO,IVP,NCOMPD,NFCD
C
C     ******************************************************************
C         THESE COMMON BLOCKS AID IN REDUCING THE NUMBER OF SUBROUTINE
C         ARGUMENTS PREVALENT IN THIS MODULAR STRUCTURE
C
      COMMON /DML18J/ AED,RED,TOL,NXPTSD,NICD,NOPG,MXNON,NDISK,NTAPE,
     1                NEQ,INDPVT,INTEG,NPS,NTP,NEQIVD,NUMORT,NFCC,
     2                ICOCO
      COMMON /DML17B/ KKKZPW,NEEDW,NEEDIW,K1,K2,K3,K4,K5,K6,K7,K8,K9,
     1                K10,K11,L1,L2,KKKINT,LLLINT
C
C     ******************************************************************
C         THIS COMMON BLOCK IS USED IN SUBROUTINES DBVSUP,DBVPOR,DRKFAB,
C         DREORT, AND DSTWAY. IT CONTAINS INFORMATION NECESSARY
C         FOR THE ORTHONORMALIZATION TESTING PROCEDURE AND A BACKUP
C         RESTARTING CAPABILITY.
C
      COMMON /DML15T/ PX,PWCND,TND,X,XBEG,XEND,XOT,XOP,INFO(15),ISTKOP,
     1                KNSWOT,KOP,LOTJP,MNSWOT,NSWOT
C
C     ******************************************************************
C         THIS COMMON BLOCK CONTAINS THE MACHINE DEPENDENT PARAMETERS
C         USED BY THE CODE
C
      COMMON /DML5MC/ URO,SRU,EPS,SQOVFL,TWOU,FOURU,LPAR
C
C      *****************************************************************
C          SET UP MACHINE DEPENDENT CONSTANTS.
C
C***FIRST EXECUTABLE STATEMENT  DBVSUP
                        CALL DMACON
C
C                       ************************************************
C                           TEST FOR INVALID INPUT
C
                        IF (NROWY .LT. NCOMP) GO TO 80
                        IF (NCOMP .NE. NIC + NFC) GO TO 80
                        IF (NXPTS .LT. 2) GO TO 80
                        IF (NIC .LE. 0) GO TO 80
                        IF (NROWA .LT. NIC) GO TO 80
                        IF (NFC .LE. 0) GO TO 80
                        IF (NROWB .LT. NFC) GO TO 80
                        IF (IGOFX .LT. 0 .OR. IGOFX .GT. 1) GO TO 80
                        IF (RE .LT. 0.0D0) GO TO 80
                        IF (AE .LT. 0.0D0) GO TO 80
                        IF (RE .EQ. 0.0D0 .AND. AE .EQ. 0.0D0) GO TO 80
C                          BEGIN BLOCK PERMITTING ...EXITS TO 70
                              IS = 1
                              IF (XPTS(NXPTS) .LT. XPTS(1)) IS = 2
                              NXPTSM = NXPTS - 1
                              DO 30 K = 1, NXPTSM
                                 IF (IS .EQ. 2) GO TO 10
C                          .........EXIT
                                    IF (XPTS(K+1) .LE. XPTS(K)) GO TO 70
                                 GO TO 20
   10                            CONTINUE
C                          .........EXIT
                                    IF (XPTS(K) .LE. XPTS(K+1)) GO TO 70
   20                            CONTINUE
   30                         CONTINUE
C
C                             ******************************************
C                                 CHECK FOR DISK STORAGE
C
                              KPTS = NXPTS
                              NDISK = 0
                              IF (IWORK(12) .EQ. 0) GO TO 40
                                 NTAPE = IWORK(12)
                                 KPTS = 1
                                 NDISK = 1
   40                         CONTINUE
C
C                             ******************************************
C                                 SET INTEG PARAMETER ACCORDING TO
C                                 CHOICE OF INTEGRATOR.
C
                              INTEG = 1
                              IF (IWORK(9) .EQ. 2) INTEG = 2
C
C                             ******************************************
C                                 COMPUTE INHOMO
C
C                 ............EXIT
                              IF (IGOFX .EQ. 1) GO TO 100
                              DO 50 J = 1, NIC
C                 ...............EXIT
                                 IF (ALPHA(J) .NE. 0.0D0) GO TO 100
   50                         CONTINUE
                              DO 60 J = 1, NFC
C                    ............EXIT
                                 IF (BETA(J) .NE. 0.0D0) GO TO 90
   60                         CONTINUE
                              INHOMO = 3
C              ...............EXIT
                              GO TO 110
   70                      CONTINUE
   80                   CONTINUE
                        IFLAG = -2
C     ..................EXIT
                        GO TO 220
   90                CONTINUE
                     INHOMO = 2
C              ......EXIT
                     GO TO 110
  100             CONTINUE
                  INHOMO = 1
  110          CONTINUE
C
C              *********************************************************
C                  TO TAKE ADVANTAGE OF THE SPECIAL STRUCTURE WHEN
C                  SOLVING A COMPLEX*16 VALUED PROBLEM,WE INTRODUCE
C                  NFCC=NFC WHILE CHANGING THE INTERNAL VALUE OF NFC
C
               NFCC = NFC
               IF (IFLAG .EQ. 13) NFC = NFC/2
C
C              *********************************************************
C                  DETERMINE NECESSARY STORAGE REQUIREMENTS
C
C              FOR BASIC ARRAYS IN DBVPOR
               KKKYHP = NCOMP*(NFC + 1) + NEQIVP
               KKKU = NCOMP*NFC*KPTS
               KKKV = NCOMP*KPTS
               KKKCOE = NFCC
               KKKS = NFC + 1
               KKKSTO = NCOMP*(NFC + 1) + NEQIVP + 1
               KKKG = NCOMP
C
C              FOR ORTHONORMALIZATION RELATED MATTERS
               NTP = (NFCC*(NFCC + 1))/2
               KKKZPW = 1 + NTP + NFCC
               LLLIP = NFCC
C
C              FOR ADDITIONAL REQUIRED WORK SPACE
C                (DLSSUD)
               KKKSUD = 4*NIC + (NROWA + 1)*NCOMP
               LLLSUD = NIC
C              (DVECS)
               KKKSVC = 1 + 4*NFCC + 2*NFCC**2
               LLLSVC = 2*NFCC
C
               NDEQ = NCOMP*NFC + NEQIVP
               IF (INHOMO .EQ. 1) NDEQ = NDEQ + NCOMP
               GO TO (120,130), INTEG
C              (DDERKF)
  120          CONTINUE
                  KKKINT = 33 + 7*NDEQ
                  LLLINT = 34
               GO TO 140
C              (DDEABM)
  130          CONTINUE
                  KKKINT = 130 + 21*NDEQ
                  LLLINT = 51
  140          CONTINUE
C
C              (COEF)
               KKKCOF = 5*NFCC + NFCC**2
               LLLCOF = 3 + NFCC
C
               KKKWS = MAX(KKKSUD,KKKSVC,KKKINT,KKKCOF)
               LLLIWS = MAX(LLLSUD,LLLSVC,LLLINT,LLLCOF)
C
               NEEDW = KKKYHP + KKKU + KKKV + KKKCOE + KKKS + KKKSTO
     1                 + KKKG + KKKZPW + KKKWS
               NEEDIW = 17 + LLLIP + LLLIWS
C              *********************************************************
C                  COMPUTE THE NUMBER OF POSSIBLE ORTHONORMALIZATIONS
C                  WITH THE ALLOTTED STORAGE
C
               IWORK(3) = NEEDW
               IWORK(4) = KKKZPW
               IWORK(5) = NEEDIW
               IWORK(6) = LLLIP
               NRTEMP = NDW - NEEDW
               NITEMP = NDIW - NEEDIW
C           ...EXIT
               IF (NRTEMP .LT. 0) GO TO 180
C           ...EXIT
               IF (NITEMP .LT. 0) GO TO 180
C
               IF (NDISK .EQ. 0) GO TO 150
                  NON = 0
                  MXNON = NRTEMP
               GO TO 160
  150          CONTINUE
C
                  MXNONR = NRTEMP/KKKZPW
                  MXNONI = NITEMP/LLLIP
                  MXNON = MIN(MXNONR,MXNONI)
                  NON = MXNON
  160          CONTINUE
C
               IWORK(2) = MXNON
C
C              *********************************************************
C                  CHECK FOR PRE-ASSIGNED ORTHONORMALIZATION POINTS
C
               NOPG = 0
C        ......EXIT
               IF (IWORK(11) .NE. 1) GO TO 210
               IF (MXNON .LT. IWORK(1)) GO TO 170
                  NOPG = 1
                  MXNON = IWORK(1)
                  WORK(MXNON+1) = 2.0D0*XPTS(NXPTS) - XPTS(1)
C        .........EXIT
                  GO TO 210
  170          CONTINUE
  180       CONTINUE
C
            IFLAG = -1
      IF (NDISK .NE. 1) THEN
         WRITE (XERN1, '(I8)') NEEDW
         WRITE (XERN2, '(I8)') KKKZPW
         WRITE (XERN3, '(I8)') NEEDIW
         WRITE (XERN4, '(I8)') LLLIP
         CALL XERMSG ('SLATEC', 'DBVSUP',
     *      'REQUIRED STORAGE FOR WORK ARRAY IS '  // XERN1 // ' + ' //
     *      XERN2 // '*(EXPECTED NUMBER OF ORTHONORMALIZATIONS) $$'  //
     *      'REQUIRED STORAGE FOR IWORK ARRAY IS ' // XERN3 // ' + ' //
     *      XERN4 // '*(EXPECTED NUMBER OF ORTHONORMALIZATIONS)', 1, 0)
      ELSE
         WRITE (XERN1, '(I8)') NEEDW
         WRITE (XERN2, '(I8)') NEEDIW
         CALL XERMSG ('SLATEC', 'DBVSUP',
     *      'REQUIRED STORAGE FOR WORK ARRAY IS '  // XERN1 //
     *      ' + NUMBER OF ORTHONOMALIZATIONS. $$'  //
     *      'REQUIRED STORAGE FOR IWORK ARRAY IS ' // XERN2, 1, 0)
      ENDIF
      RETURN
C
C        ***************************************************************
C            ALLOCATE STORAGE FROM WORK AND IWORK ARRAYS
C
C         (Z)
  210    K1 = 1 + (MXNON + 1)
C        (P)
         K2 = K1 + NTP*(NON + 1)
C        (W)
         K3 = K2 + NFCC*(NON + 1)
C        (YHP)
         K4 = K3 + KKKYHP
C        (U)
         K5 = K4 + KKKU
C        (V)
         K6 = K5 + KKKV
C        (COEF)
         K7 = K6 + KKKCOE
C        (S)
         K8 = K7 + KKKS
C        (STOWA)
         K9 = K8 + KKKSTO
C        (G)
         K10 = K9 + KKKG
         K11 = K10 + KKKWS
C                  REQUIRED ADDITIONAL DOUBLE PRECISION WORK SPACE
C                  STARTS AT WORK(K10) AND EXTENDS TO WORK(K11-1)
C
C           FIRST 17 LOCATIONS OF IWORK ARE USED FOR OPTIONAL
C           INPUT AND OUTPUT ITEMS
C        (IP)
         L1 = 18 + NFCC*(NON + 1)
         L2 = L1 + LLLIWS
C                   REQUIRED INTEGER WORK SPACE STARTS AT IWORK(L1)
C                   AND EXTENDS TO IWORK(L2-1)
C
C        ***************************************************************
C            SET INDICATOR FOR NORMALIZATION OF PARTICULAR SOLUTION
C
         NPS = 0
         IF (IWORK(10) .EQ. 1) NPS = 1
C
C        ***************************************************************
C            SET PIVOTING PARAMETER
C
         INDPVT = 0
         IF (IWORK(15) .EQ. 1) INDPVT = 1
C
C        ***************************************************************
C            SET OTHER COMMON BLOCK PARAMETERS
C
         NFCD = NFC
         NCOMPD = NCOMP
         IGOFXD = IGOFX
         NXPTSD = NXPTS
         NICD = NIC
         RED = RE
         AED = AE
         NEQIVD = NEQIVP
         MNSWOT = 20
         IF (IWORK(13) .EQ. -1) MNSWOT = MAX(1,IWORK(14))
         XBEG = XPTS(1)
         XEND = XPTS(NXPTS)
         XSAV = XEND
         ICOCO = 1
         IF (INHOMO .EQ. 3 .AND. NOPG .EQ. 1) WORK(MXNON+1) = XEND
C
C        ***************************************************************
C
         CALL DEXBVP(Y,NROWY,XPTS,A,NROWA,ALPHA,B,NROWB,BETA,IFLAG,WORK,
     1               IWORK)
         NFC = NFCC
         IWORK(17) = IWORK(L1)
  220 CONTINUE
      RETURN
      END


text:
SUBROUTINE c2num
!======================================
! 2nd Neighbour Number
!======================================
! Calculates the second nearest neighbour number for each atom in the cluster.
! Calculates the derivatives.
! Use the c2num light version disregard the derivative loop and improve performance.
!======================================
   USE PARACLUSTER
   USE CLUSTER
   USE POTENTIAL
   USE ENFORCE
   USE DISTANCE
   USE META

   IMPLICIT NONE

   integer               :: i, j, mminusn_pwr
   double precision      :: ratio_ton, ratio_tom, rnp, rmp, num
   double precision      :: ratio, comfactor
   real(8), dimension(3) :: rij0 ! 3 possible cases: atoms A-A, B-B, A-B
   real(8), dimension(3) :: dij0 ! 3 possible cases: atoms A-A, B-B, A-B
   real(8), DIMENSION(nsiz,nsiz) :: dij, df_dx, df_dy, df_dz

   if (ipas.eq.1) write(*,*) '> c2num: first calculation of window function 1 (C2N)'

   !------------------------------------------------------------------------
   mminusn_pwr = m2n_pwr - n2n_pwr

   ! Evaluating the values of rij0 parameter in dependence of atom pair kind
   rij0 = dist * r2n * dsqrt(2.d0) ! rij0 and dist are 3 element vectors

   ! Evaluating the ref distance for the couple of atoms (THIS IS DIFFERENT)
   dij0 = dist * d2n * dsqrt(2.d0)
   !------------------------------------------------------------------------

!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(i, j), SCHEDULE(STATIC)
   DO i = 1, natom
      DO j = 1, natom
         ! Calculating the distance between atom i and atom j
         dij(j,i) = pair_dist (j,i) - dij0(pairkindmat(j,i))
      END DO
   END DO
!$OMP END DO

!$OMP DO PRIVATE(i), SCHEDULE(STATIC)
   ! Elements on the diagonal
   DO i = 1, natom
      s2n(i,i) = 0.d0
   END DO
!$OMP END DO

!$OMP DO PRIVATE(i, j, ratio, ratio_ton, ratio_tom, rnp, rmp, num), SCHEDULE(DYNAMIC,20)
   DO i = 1, natom-1
      DO j = i+1, natom
         ! Calculating the analytic C2N
         ratio = dij(j,i)/rij0(pairkindmat(j,i))
         ! I need to rise to (X_pwr-1) for the derivative
         ratio_ton = ratio**(n2n_pwr-1)
         ratio_tom = ratio_ton*(ratio**(mminusn_pwr))

         ! Calculating some pieces of derivative
         rnp = n2n_pwr * ratio_ton
         rmp = m2n_pwr * ratio_tom

         ! Now ratio_ton becomes really ratio to n_pwr
         ratio_ton = ratio_ton * ratio

         ! Another piece of derivative
         num = rmp - rnp - mminusn_pwr*ratio_ton*ratio_tom
         ! same as:
         ! num = rmp - rnp + ( n_pwr - m_pwr )*ratio**(n_pwr + m_pwr - 1)

         ! Now ratio_tom becomes really ratio to m_pwr
         ratio_tom = ratio_tom * ratio

         ! The sum of these s(i,j) is the analytic coord no.
         s2n(j,i) = (1.d0-ratio_ton)/(1.d0-ratio_tom)  ! this is f(rij)
         s2n(i,j) = s2n(j,i)

         ! Derivative of s
         dS(j,i) = num/( rij0(pairkindmat(j,i))* (1.d0 - ratio_tom)**2)
         dS(i,j) = dS(j,i)
      END DO ! on j
   END DO ! on i
!$OMP END DO

!$OMP DO PRIVATE(i), SCHEDULE(STATIC)
   ! LP: Calculating the sum of s2n(i,j), which is C2Natom(i), global.
   ! This is the 2nd neigh num of any single atom i
   ! No problem because s2n(i,i) is always zero
   DO i = 1, natom
      C2Natom(i) = SUM(s2n(:,i))
   END DO
!$OMP END DO

!$OMP DO PRIVATE(i,j,comfactor), SCHEDULE(STATIC)
   !-------------------------------------------------------------------------
   ! Calculating the sum on j of derivatives of s
   ! with respect to coordinates of every atoms
   DO i = 1, natom
      DO j = 1, natom
         IF (j .NE.i ) THEN
            ! Calculating the derivative of coordination number with
            ! respect to distance between atoms i and j, multiplied for
            ! the derivative of distance between atoms i and j
            ! with respect to position of atom i
            comfactor  = dS(j,i)/ pair_dist(j,i)    ! Common factor for the following...
            df_dx(j,i) = comfactor* xyz_dist(1,j,i) ! order of indexes depends on the definition of xyz_dist
            df_dy(j,i) = comfactor* xyz_dist(2,j,i)
            df_dz(j,i) = comfactor* xyz_dist(3,j,i)
         ELSE
            df_dx(j,i) = 0.d0
            df_dy(j,i) = 0.d0
            df_dz(j,i) = 0.d0
         END IF
      END DO
   END DO
!$OMP END DO

!$OMP SECTIONS PRIVATE(i)
!$OMP SECTION
   DO i = 1, natom
      dS2n_dx(i) = 2.d0 * SUM(df_dx(1:natom, i))
   END DO
!$OMP SECTION
   DO i = 1, natom
      dS2n_dy(i) = 2.d0 * SUM(df_dy(1:natom, i))
   END DO
!$OMP SECTION
   DO i = 1, natom
      dS2n_dz(i) = 2.d0 * SUM(df_dz(1:natom, i))
   END DO
!$OMP END SECTIONS
!$OMP END PARALLEL
   !---------------------------------------------------------------------------
END SUBROUTINE c2num

max_stars_repo_name: kcl-tscm/LoDiS

max_issues_repo_name: kcl-tscm/LoDiS

max_forks_repo_name: kcl-tscm/LoDiS


text:
Hoagland Hall Annex, or just The Annex is where the graduate students in Atmospheric Science live used to live, though theyve now moved upstairs.  Its essentially a garage converted into offices, which explains why they moved out.  The Biometeorology Group definitely had the best room, as its the only one where the air conditioner really works in the Summer.

Rumor has it that the Art Department will be taking over the space, which is way across campus from most of its offices.

It also has a workshop associated with UC Agriculture and Natural Resources DANR on the west side.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      program user_error03

C     Check that PIPS survives the user error in the call graph because
C     B is called with only one argument and that PIPS can then compute the
C     effects for B.

      call B(x)

      end

      subroutine B(x, y)

      print *, x+y

      end


text: #redirect Police Department Trading Cards


text:
      PROGRAM LOOPWHILE
      INTEGER I
C     FORTRAN 66 does not have IF block.
      I = 1024
 10   CONTINUE
      IF (I .LE. 0) GOTO 20
      WRITE (*,*) I
      I = I / 2
      GOTO 10
 20   CONTINUE
      STOP
      END

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
	SUBROUTINE SNMCLR  ( ncol, prmlst, colors, icolor, iret )
C************************************************************************
C* SNMCLR								*
C*									*
C* This subroutine translates the colors.				*
C*									*
C* SNMCLR  ( NCOL, PRMLST, COLORS, ICOLOR, IRET )			*
C*									*
C* Input parameters:							*
C*	NCOL		INTEGER		Number of colors		*
C*	PRMLST(*)	CHAR*		Parameter names			*
C*	COLORS		CHAR*		COLORS input			*
C*									*
C* Output parameters:							*
C*	ICOLOR (*)	INTEGER		Colors array			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal			*
C**									*
C* Log:									*
C* I. Graffman/RDS	 7/87	GEMPAK4					*
C* M. desJardins/GSFC	 6/88	Reorganized				*
C* S. Schotz/GSC	 8/90	Reorganized as SFMCLR			*
C************************************************************************
	CHARACTER*(*)	colors, prmlst (*)
	INTEGER		icolor (*)
C*
	INTEGER		icbuff (11)
C------------------------------------------------------------------------
	CALL IN_COLR  ( colors, ncol, icbuff, iret )
C
C*	Turn off colors for 'BLNK' and 'SPAC'.
C
	ic  = 0
	DO  i = 1, ncol
            IF  ( (prmlst (i) .eq. 'BLNK') .or.
     1            (prmlst(i) .eq. 'SPAC') ) THEN
                icolor (i) = 0
              ELSE
		ic = ic + 1
		icolor (i) = icbuff (ic)
	    END IF
	END DO
C*
	RETURN
	END


text:
+incdir+${TUE_HOME}/src
${TUE_HOME}/src/tue_pkg.sv

max_stars_repo_name: little-blue/Surelog

max_issues_repo_name: little-blue/Surelog

max_forks_repo_name: little-blue/Surelog


text:
C    @(#)kpovlt.f	20.3 2/13/96
      function kpovlt(i,j)

C *** This function determines the sort key for overloaded transformer *
C *** array OLTR.                                                      *

      include 'ipfinc/parametr.inc'

      include 'ipfinc/anlys.inc'
      include 'ipfinc/arcntl.inc'
      include 'ipfinc/area.inc'
      include 'ipfinc/branch.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/busanl.inc'

      common /sortsw/ sortsw, vltsrt(MAXBUS)
      integer sortsw, vltsrt

      character own1 * 3, own2 * 3

      j1 = koltr(1,i)
      j2 = koltr(1,j)

      if (sortsw .eq. 1) then  !          Sort by ownership-bus
         nbr = iabs (brnch_ptr(j1))
         write (own1,100) kbrnch(3,nbr)
  100    format (a3)
         nbr = iabs (brnch_ptr(j2))
         write (own2,100) kbrnch(3,nbr)
         kpovlt = kompr (own1,own2,kpovlt)
         if (kpovlt .eq. 0) then
            kpovlt = inp2alf(kx(j1)) - inp2alf(kx(j2))
         endif
         if (kpovlt .eq. 0) then
            kpovlt = inp2alf(ky(j1)) - inp2alf(ky(j2))
         endif
         if (kpovlt .eq. 0) then
            kpovlt = kompr (brid(j1), brid(j2), junk)
         endif
         if (kpovlt .eq. 0) then
            kpovlt = isign (brsect(j1), brnch_ptr(j1))
     &             - isign (brsect(j2), brnch_ptr(j2))
         endif

      else if (sortsw .eq. 2) then   !    Sort by bus
         kpovlt = inp2alf(kx(j1)) - inp2alf(kx(j2))
         if (kpovlt .eq. 0) then
            kpovlt = inp2alf(ky(j1)) - inp2alf(ky(j2))
         endif
         if (kpovlt .eq. 0) then
            kpovlt = kompr (brid(j1), brid(j2), junk)
         endif
         if (kpovlt .eq. 0) then
            kpovlt = isign (brsect(j1), brnch_ptr(j1))
     &             - isign (brsect(j2), brnch_ptr(j2))
         endif

      else if (sortsw .eq. 3) then   !    Sort by zone-bus
         k1 = kx(j1)
         k2 = kx(j2)
         kpovlt = kompr (zone(k1),zone(k2),kpovlt)
         if (kpovlt .eq. 0) then
            m1 = ky(j1)
            m2 = ky(j2)
            kpovlt = kompr (zone(m1),zone(m2),kpovlt)
            if (kpovlt .eq. 0) then
               kpovlt = inp2alf(kx(j1)) - inp2alf(kx(j2))
            endif
            if (kpovlt .eq. 0) then
               kpovlt = inp2alf(ky(j1)) - inp2alf(ky(j2))
            endif
            if (kpovlt .eq. 0) then
               kpovlt = kompr (brid(j1), brid(j2), junk)
            endif
            if (kpovlt .eq. 0) then
               kpovlt = isign (brsect(j1), brnch_ptr(j1))
     &                - isign (brsect(j2), brnch_ptr(j2))
            endif
         endif
      else if (sortsw .eq. 4) then   !    Sort by area-bus
         k1 = kx(j1)
         k2 = kx(j2)
         l1 = jarzn(k1)
         l2 = jarzn(k2)
         kpovlt = kompr (arcnam(l1),arcnam(l2),kpovlt)
         if (kpovlt .eq. 0) then
            m1 = ky(j1)
            m2 = ky(j2)
            l1 = jarzn(k1)
            l2 = jarzn(k2)
            kpovlt = kompr (arcnam(l1),arcnam(l2),kpovlt)
            if (kpovlt .eq. 0) then
               kpovlt = inp2alf(kx(j1)) - inp2alf(kx(j2))
            endif
            if (kpovlt .eq. 0) then
               kpovlt = inp2alf(ky(j1)) - inp2alf(ky(j2))
            endif
            if (kpovlt .eq. 0) then
               kpovlt = kompr (brid(j1), brid(j2), junk)
            endif
            if (kpovlt .eq. 0) then
               kpovlt = isign (brsect(j1), brnch_ptr(j1))
     &                - isign (brsect(j2), brnch_ptr(j2))
            endif
         endif
      endif
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
      SUBROUTINE PZPOCON( UPLO, N, A, IA, JA, DESCA, ANORM, RCOND, WORK,
     $                    LWORK, RWORK, LRWORK, INFO )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 25, 2001
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            IA, INFO, JA, LRWORK, LWORK, N
      DOUBLE PRECISION   ANORM, RCOND
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * )
      DOUBLE PRECISION   RWORK( * )
      COMPLEX*16         A( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  PZPOCON estimates the reciprocal of the condition number (in the
*  1-norm) of a complex Hermitian positive definite distributed matrix
*  using the Cholesky factorization A = U**H*U or A = L*L**H computed by
*  PZPOTRF.
*
*  An estimate is obtained for norm(inv(A(IA:IA+N-1,JA:JA+N-1))), and
*  the reciprocal of the condition number is computed as
*             RCOND = 1 / ( norm( A(IA:IA+N-1,JA:JA+N-1)      ) *
*                           norm( inv(A(IA:IA+N-1,JA:JA+N-1)) ) ).
*
*  Notes
*  =====
*
*  Each global data object is described by an associated description
*  vector.  This vector stores the information required to establish
*  the mapping between an object element and its corresponding process
*  and memory location.
*
*  Let A be a generic term for any 2D block cyclicly distributed array.
*  Such a global array has an associated description vector DESCA.
*  In the following comments, the character _ should be read as
*  "of the global array".
*
*  NOTATION        STORED IN      EXPLANATION
*  --------------- -------------- --------------------------------------
*  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
*                                 DTYPE_A = 1.
*  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
*                                 the BLACS process grid A is distribu-
*                                 ted over. The context itself is glo-
*                                 bal, but the handle (the integer
*                                 value) may vary.
*  M_A    (global) DESCA( M_ )    The number of rows in the global
*                                 array A.
*  N_A    (global) DESCA( N_ )    The number of columns in the global
*                                 array A.
*  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
*                                 the rows of the array.
*  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
*                                 the columns of the array.
*  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
*                                 row of the array A is distributed.
*  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
*                                 first column of the array A is
*                                 distributed.
*  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
*
*  Let K be the number of rows or columns of a distributed matrix,
*  and assume that its process grid has dimension p x q.
*  LOCr( K ) denotes the number of elements of K that a process
*  would receive if K were distributed over the p processes of its
*  process column.
*  Similarly, LOCc( K ) denotes the number of elements of K that a
*  process would receive if K were distributed over the q processes of
*  its process row.
*  The values of LOCr() and LOCc() may be determined via a call to the
*  ScaLAPACK tool function, NUMROC:
*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
*  An upper bound for these quantities may be computed by:
*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
*  Arguments
*  =========
*
*  UPLO    (global input) CHARACTER
*          Specifies whether the factor stored in
*          A(IA:IA+N-1,JA:JA+N-1) is upper or lower triangular.
*          = 'U':  Upper triangular
*          = 'L':  Lower triangular
*
*  N       (global input) INTEGER
*          The order of the distributed matrix A(IA:IA+N-1,JA:JA+N-1).
*          N >= 0.
*
*  A       (local input) COMPLEX*16 pointer into the local memory to
*          an array of dimension ( LLD_A, LOCc(JA+N-1) ). On entry, this
*          array contains the local pieces of the factors L or U from
*          the Cholesky factorization A(IA:IA+N-1,JA:JA+N-1) = U'*U or
*          L*L', as computed by PZPOTRF.
*
*  IA      (global input) INTEGER
*          The row index in the global array A indicating the first
*          row of sub( A ).
*
*  JA      (global input) INTEGER
*          The column index in the global array A indicating the
*          first column of sub( A ).
*
*  DESCA   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix A.
*
*  ANORM   (global input) DOUBLE PRECISION
*          The 1-norm (or infinity-norm) of the hermitian distributed
*          matrix A(IA:IA+N-1,JA:JA+N-1).
*
*  RCOND   (global output) DOUBLE PRECISION
*          The reciprocal of the condition number of the distributed
*          matrix A(IA:IA+N-1,JA:JA+N-1), computed as
*             RCOND = 1 / ( norm( A(IA:IA+N-1,JA:JA+N-1)      ) *
*                           norm( inv(A(IA:IA+N-1,JA:JA+N-1)) ) ).
*
*  WORK    (local workspace/local output) COMPLEX*16 array,
*                                                   dimension (LWORK)
*          On exit, WORK(1) returns the minimal and optimal LWORK.
*
*  LWORK   (local or global input) INTEGER
*          The dimension of the array WORK.
*          LWORK is local input and must be at least
*          LWORK >= 2*LOCr(N+MOD(IA-1,MB_A)) +
*          MAX( 2, MAX(NB_A*MAX(1,CEIL(P-1,Q)),LOCc(N+MOD(JA-1,NB_A)) +
*          NB_A*MAX(1,CEIL(Q-1,P))) ).
*
*          If LWORK = -1, then LWORK is global input and a workspace
*          query is assumed; the routine only calculates the minimum
*          and optimal size for all work arrays. Each of these
*          values is returned in the first entry of the corresponding
*          work array, and no error message is issued by PXERBLA.
*
*  RWORK   (local workspace/local output) DOUBLE PRECISION array,
*                                                  dimension (LRWORK)
*          On exit, RWORK(1) returns the minimal and optimal LRWORK.
*
*  LRWORK  (local or global input) INTEGER
*          The dimension of the array RWORK.
*          LRWORK is local input and must be at least
*          LRWORK >= 2*LOCc(N+MOD(JA-1,NB_A)).
*
*          If LRWORK = -1, then LRWORK is global input and a workspace
*          query is assumed; the routine only calculates the minimum
*          and optimal size for all work arrays. Each of these
*          values is returned in the first entry of the corresponding
*          work array, and no error message is issued by PXERBLA.
*
*
*  INFO    (global output) INTEGER
*          = 0:  successful exit
*          < 0:  If the i-th argument is an array and the j-entry had
*                an illegal value, then INFO = -(i*100+j), if the i-th
*                argument is a scalar and had an illegal value, then
*                INFO = -i.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $                   LLD_, MB_, M_, NB_, N_, RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
     $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
      DOUBLE PRECISION   ONE, ZERO
      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY, UPPER
      CHARACTER          CBTOP, COLCTOP, NORMIN, ROWCTOP
      INTEGER            IACOL, IAROW, ICOFF, ICTXT, IIA, IPNL, IPNU,
     $                   IPV, IPW, IPX, IROFF, IV, IX, IXX, JJA, JV,
     $                   JX, KASE, LRWMIN, LWMIN, MYCOL, MYROW, NP,
     $                   NPCOL, NPROW, NPMOD, NQ, NQMOD
      DOUBLE PRECISION   AINVNM, SCALE, SL, SU, SMLNUM
      COMPLEX*16         WMAX, ZDUM
*     ..
*     .. Local Arrays ..
      INTEGER            DESCV( DLEN_ ), DESCX( DLEN_ ), IDUM1( 3 ),
     $                   IDUM2( 3 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, CHK1MAT, DESCSET, INFOG2L,
     $                   PB_TOPGET, PB_TOPSET, PXERBLA, PCHK1MAT,
     $                   PZAMAX, PZLATRS, PZLACON, PZDRSCL,
     $                   ZGEBR2D, ZGEBS2D
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ICEIL, INDXG2P, NUMROC
      DOUBLE PRECISION   PDLAMCH
      EXTERNAL           ICEIL, INDXG2P, LSAME, NUMROC, PDLAMCH
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, DBLE, DIMAG, ICHAR, MAX, MOD
*     ..
*     .. Statement Functions ..
      DOUBLE PRECISION   CABS1
*     ..
*     .. Statement Function definitions ..
      CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
*     ..
*     .. Executable Statements ..
*
*     Get grid parameters
*
      ICTXT = DESCA( CTXT_ )
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
*     Test the input parameters
*
      INFO = 0
      IF( NPROW.EQ.-1 ) THEN
         INFO = -(600+CTXT_)
      ELSE
         CALL CHK1MAT( N, 2, N, 2, IA, JA, DESCA, 6, INFO )
         IF( INFO.EQ.0 ) THEN
            UPPER = LSAME( UPLO, 'U' )
            IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ),
     $                       NPROW )
            IACOL = INDXG2P( JA, DESCA( NB_ ), MYCOL, DESCA( CSRC_ ),
     $                       NPCOL )
            NPMOD = NUMROC( N + MOD( IA-1, DESCA( MB_ ) ), DESCA( MB_ ),
     $                      MYROW, IAROW, NPROW )
            NQMOD = NUMROC( N + MOD( JA-1, DESCA( NB_ ) ), DESCA( NB_ ),
     $                      MYCOL, IACOL, NPCOL )
            LWMIN = 2*NPMOD +
     $              MAX( 2, MAX( DESCA( NB_ )*
     $                   MAX( 1, ICEIL( NPROW-1, NPCOL ) ), NQMOD +
     $                   DESCA( NB_ )*
     $                   MAX( 1, ICEIL( NPCOL-1, NPROW ) ) ) )
            WORK( 1 ) = DBLE( LWMIN )
            LRWMIN = 2*NQMOD
            RWORK( 1 ) = DBLE( LRWMIN )
            LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 )
*
            IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
               INFO = -1
            ELSE IF( ANORM.LT.ZERO ) THEN
               INFO = -7
            ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
               INFO = -10
            ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
               INFO = -12
            END IF
         END IF
*
         IF( UPPER ) THEN
            IDUM1( 1 ) = ICHAR( 'U' )
         ELSE
            IDUM1( 1 ) = ICHAR( 'L' )
         END IF
         IDUM2( 1 ) = 1
         IF( LWORK.EQ.-1 ) THEN
            IDUM1( 2 ) = -1
         ELSE
            IDUM1( 2 ) = 1
         END IF
         IDUM2( 2 ) = 10
         IF( LRWORK.EQ.-1 ) THEN
            IDUM1( 3 ) = -1
         ELSE
            IDUM1( 3 ) = 1
         END IF
         IDUM2( 3 ) = 12
         CALL PCHK1MAT( N, 2, N, 2, IA, JA, DESCA, 6, 3, IDUM1, IDUM2,
     $                  INFO )
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL PXERBLA( ICTXT, 'PZPOCON', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
*
*     Quick return if possible
*
      RCOND = ZERO
      IF( N.EQ.0 ) THEN
         RCOND = ONE
         RETURN
      ELSE IF( ANORM.EQ.ZERO ) THEN
         RETURN
      ELSE IF( N.EQ.1 ) THEN
         RCOND = ONE
         RETURN
      END IF
*
      CALL PB_TOPGET( ICTXT, 'Combine', 'Columnwise', COLCTOP )
      CALL PB_TOPGET( ICTXT, 'Combine', 'Rowwise',    ROWCTOP )
      CALL PB_TOPSET( ICTXT, 'Combine', 'Columnwise', '1-tree' )
      CALL PB_TOPSET( ICTXT, 'Combine', 'Rowwise',    '1-tree' )
*
      SMLNUM = PDLAMCH( ICTXT, 'Safe minimum' )
      IROFF = MOD( IA-1, DESCA( MB_ ) )
      ICOFF = MOD( JA-1, DESCA( NB_ ) )
      CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, IIA, JJA,
     $              IAROW, IACOL )
      NP = NUMROC( N+IROFF, DESCA( MB_ ), MYROW, IAROW, NPROW )
      NQ = NUMROC( N+ICOFF, DESCA( NB_ ), MYCOL, IACOL, NPCOL )
      IV = IROFF + 1
      IX = IV
      JV = ICOFF + 1
      JX = JV
*
      IPX = 1
      IPV = IPX + NP
      IPW = IPV + NP
      IPNL = 1
      IPNU = IPNL + NQ
*
      CALL DESCSET( DESCV, N+IROFF, 1, DESCA( MB_ ), 1, IAROW, MYCOL,
     $              ICTXT, MAX( 1, NP ) )
      CALL DESCSET( DESCX, N+IROFF, 1, DESCA( MB_ ), 1, IAROW, MYCOL,
     $              ICTXT, MAX( 1, NP ) )
*
*     Estimate the 1-norm (or I-norm) of inv(A).
*
      AINVNM = ZERO
      KASE   = 0
      NORMIN = 'N'
*
   10 CONTINUE
      CALL PZLACON( N, WORK( IPV ), IV, JV, DESCV, WORK( IPX ), IX, JX,
     $              DESCX, AINVNM, KASE )
      IF( KASE.NE.0 ) THEN
         IF( UPPER ) THEN
*
*           Multiply by inv(U').
*
            DESCX( CSRC_ ) = IACOL
            CALL PZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
     $                    NORMIN, N, A, IA, JA, DESCA, WORK( IPX ),
     $                    IX, JX, DESCX, SL, RWORK( IPNL ),
     $                    WORK( IPW ) )
            DESCX( CSRC_ ) = MYCOL
            NORMIN = 'Y'
*
*           Multiply by inv(U).
*
            DESCX( CSRC_ ) = IACOL
            CALL PZLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN,
     $                    N, A, IA, JA, DESCA, WORK( IPX ), IX, JX,
     $                    DESCX, SU, RWORK( IPNU ), WORK( IPW ) )
            DESCX( CSRC_ ) = MYCOL
         ELSE
*
*           Multiply by inv(L).
*
            DESCX( CSRC_ ) = IACOL
            CALL PZLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN,
     $                    N, A, IA, JA, DESCA, WORK( IPX ), IX, JX,
     $                    DESCX, SL, RWORK( IPNL ), WORK( IPW ) )
            DESCX( CSRC_ ) = MYCOL
            NORMIN = 'Y'
*
*           Multiply by inv(L').
*
            DESCX( CSRC_ ) = IACOL
            CALL PZLATRS( 'Lower', 'Conjugate transpose', 'Non-unit',
     $                    NORMIN, N, A, IA, JA, DESCA, WORK( IPX ),
     $                    IX, JX, DESCX, SU, RWORK( IPNU ),
     $                    WORK( IPW ) )
            DESCX( CSRC_ ) = MYCOL
         END IF
*
*        Multiply by 1/SCALE if doing so will not cause overflow.
*
         SCALE = SL*SU
         IF( SCALE.NE.ONE ) THEN
            CALL PZAMAX( N, WMAX, IXX, WORK( IPX ), IX, JX, DESCX, 1 )
            IF( DESCX( M_ ).EQ.1 .AND. N.EQ.1 ) THEN
               CALL PB_TOPGET( ICTXT, 'Broadcast', 'Columnwise', CBTOP )
               IF( MYROW.EQ.IAROW ) THEN
                  CALL ZGEBS2D( ICTXT, 'Column', CBTOP, 1, 1, WMAX, 1 )
               ELSE
                  CALL ZGEBR2D( ICTXT, 'Column', CBTOP, 1, 1, WMAX, 1,
     $                          IAROW, MYCOL )
               END IF
            END IF
            IF( SCALE.LT.CABS1( WMAX )*SMLNUM .OR. SCALE.EQ.ZERO )
     $         GO TO 20
            CALL PZDRSCL( N, SCALE, WORK( IPX ), IX, JX, DESCX, 1 )
         END IF
         GO TO 10
      END IF
*
*     Compute the estimate of the reciprocal condition number.
*
      IF( AINVNM.NE.ZERO )
     $   RCOND = ( ONE / AINVNM ) / ANORM
*
   20 CONTINUE
*
      CALL PB_TOPSET( ICTXT, 'Combine', 'Columnwise', COLCTOP )
      CALL PB_TOPSET( ICTXT, 'Combine', 'Rowwise',    ROWCTOP )
*
      RETURN
*
*     End of PZPOCON
*
      END

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


text:
module usub_utils_mod
implicit none

    private

    public  :: get_fid
    public  :: check_iostat
    public  :: write_node_info

    contains

    ! Use Abaqus' utility routine getoutdir to determine the full path to the base_name that resides
    ! in the current working directory. Open that file and return the file identifier.
    function get_fid(base_name, the_action) result(file_id)
    use abaqus_utils_mod
    implicit none
        character(len=*), intent(in)            :: base_name    ! Base name for file
        character(len=*), intent(in), optional  :: the_action   ! action for open, default: 'read'
        integer                                 :: file_id      ! File identifier for file to read

        integer                                 :: io_status    ! Used to check that file opens
                                                                ! sucessfully
        integer                                 :: cwd_length   ! Length of current path (used by
                                                                ! getoutdir)
        character(len=256)                      :: filename     ! Filename (full path)
        character(len=20)                       :: int_action   ! Internal action

        if (present(the_action)) then
            int_action = the_action
        else
            int_action = 'read'
        endif

        call getoutdir(filename, cwd_length)

        if ((len(trim(filename)) + len(trim(base_name)) + 1) > len(filename)) then
            write(*,*) 'Current path = "'//trim(filename)//'" is too long'
            call xit()
        endif

        write(filename, *) trim(filename)//'/'//trim(base_name)

        open(newunit=file_id, file=trim(filename), iostat=io_status, action=int_action)
        call check_iostat(io_status, 'Error opening "'//base_name//'"')

    end function get_fid

    subroutine check_iostat(io_status, error_message)
    use abaqus_utils_mod
    implicit none
        integer, intent(in)             :: io_status        ! value returned by io_stat
        character(len=*), intent(in)    :: error_message    ! error message, if io_stat /= 0

        if (io_status /= 0) then
            write(*,*) error_message
            call xit()
        endif

    end subroutine

    subroutine write_node_info(node_label, node_coords, node_jdof, kstep, kinc)
    implicit none
        integer, intent(in)                 :: node_label
        double precision, intent(in)        :: node_coords(3)
        integer, optional, intent(in)       :: node_jdof
        integer, optional, intent(in)       :: kstep
        integer, optional, intent(in)       :: kinc

        write(*,"(A,I0)") 'Node label = ', node_label
        write(*,"(A,F10.5,F10.5,F10.5)") 'Node coord = ', node_coords(1), node_coords(2), node_coords(3)
        if (present(node_jdof)) then
            write(*,"(A,I0)") 'Node jdof  = ', node_jdof
        endif
        if (present(kstep)) then
            write(*,"(A,I0)") 'Step nr.   = ', kstep
        endif
        if (present(kinc)) then
            write(*,"(A,I0)") 'Increment  = ', kinc
        endif

    end subroutine

end module usub_utils_mod

max_stars_repo_name: KnutAM/AbaqusRolloverSimulation

max_issues_repo_name: KnutAM/AbaqusRolloverSimulation

max_forks_repo_name: KnutAM/AbaqusRolloverSimulation


text:
      SUBROUTINE CMIWRT (ICODE,NAME1,NAME2,LOC,NW,A,IZ)
C
C     THIS SUBROUTINE WRITES FORMATTED SOF ITEMS.
C     ICODE = 1 FOR EQSS    ICODE = 2 FOR BGSS    ICODE = 3 FOR CSTM
C     ICODE = 4 FOR PLTS    ICODE = 5 FOR LODS    ICODE = 7 FOR LOAP
C     NAME1 IS PSEUDOSTRUCTURE NAME, NAME2 IS COMPONENT NAME
C
      EXTERNAL        ANDF
      INTEGER         OUTT,ANDF
      DIMENSION       NAME1(2),NAME2(2),A(1),IZ(1),IBITS(32),IPL(6),
     1                IH1(96),IH2(96),IH3(96),IH4(96),IH5(96),IH6(96)
      COMMON /SYSTEM/ XXX,OUTT,JUNK1(6),NLPP,JUNK2(2),NLINE
      COMMON /OUTPUT/ ITITL(96),IHEAD(96)
      DATA IH1 / 9*4H    ,4H EQS,4HS IT,4HEM F,4HOR S,4HUBST,4HRUCT,
     1           4HURE ,2*4H    ,4H COM,4HPONE,4HNT  ,11*4H    ,4HGRID,
     2           4H POI,4HNT  ,4H INT,4HERNA,4HL   ,4H  CO,4HMPON,
     3           4HENT ,2*4H    ,4H GRI,4HD PO,4HINT ,4H  IN,4HTERN,
     4           4HAL  ,4H   C,4HOMPO,4HNENT,2*4H    ,4H  GR,4HID P,
     5           4HOINT,4H   I,4HNTER,4HNAL ,4H    ,4HCOMP,4HONEN,
     6           4HT   ,4H    ,4HID  ,4H    ,4H POI,4HNT I,4HD   ,
     7           4H    ,4H DOF,4*4H    ,4H ID ,4H    ,4H  PO,4HINT ,
     8           4HID  ,4H    ,4H  DO,4HF   ,3*4H    ,4H  ID,4H    ,
     9           4H   P,4HOINT,4H ID ,4H    ,4H   D,4HOF  ,4H     /
      DATA IH2 / 11*4H    ,4HBGSS,4H ITE,4HM FO,4HR SU,4HBSTR,4HUCTU,
     1           4HRE  ,21*4H    ,4HINTE,4HRNAL,4H    ,4H CST,4HM ID,
     2           4*4H    ,4H  C ,4HO O ,4HR D ,4HI N ,4HA T ,4HE S ,
     3           17*4H    ,4HPOIN,4HT ID,4H    ,4H   N,4HO.  ,3*4H    ,
     4           4HX1  ,3*4H    ,4HX2  ,3*4H    ,4HX3  ,8*4H     /
      DATA IH3 / 12*4H    ,4HCSTM,4H ITE,4HM FO,4HR SU,4HBSTR,
     1           4HUCTU,4HRE  ,13*4H    ,2*4H    ,4H CST,4HM   ,4HTYPE,
     2           2*4H    ,4HC O ,4HO R ,4HD I ,4HN A ,4HT E ,4HS   ,
     3           4HO F ,4H  O ,4HR I ,4HG I ,4HN   ,3*4H    ,4H   T,
     4           4H R A,4H N S,4H F O,4H R M,4H A T,4H I O,4H N  ,
     5           5*4H    ,4H  ID,5*4H    ,4HX1  ,3*4H    ,4HX2  ,
     6           3*4H    ,4HX3  ,6*4H    ,4H   M,4H A T,4H R I,4H X  ,
     7           5*4H      /
      DATA IH4 / 12*4H    ,4HPLTS,4H ITE,4HM FO,4HR SU,4HBSTR,
     1           4HUCTU,4HRE  ,13*4H    ,2*4H    ,4HCOMP,4HONEN,4HT   ,
     2           4H    ,4H C O,4H O R,4H D I,4H N A,4H T E,4H S  ,
     3           4H O F,4H   O,4HR I ,4HG I ,4HN   ,3*4H    ,4H   T,
     4           4H R A,4H N S,4H F O,4H R M,4H A T,4H I O,4H N  ,
     5           6*4H    ,4H  NA,4HME  ,3*4H    ,4H X1 ,3*4H    ,
     6           4H X2 ,3*4H    ,4H X3 ,6*4H    ,4H   M,4H A T,4H R I,
     7           4H X  ,6*4H      /
      DATA IH5 / 12*4H    ,4HLODS,4H ITE,4HM FO,4HR SU,4HBSTR,
     1           4HUCTU,4HRE  ,18*4H    ,4H COM,4HPONE,4HNT  ,4H  NU,
     2           4HMBER,4H OF ,21*4H    ,5*4H    ,4H   N,4HAME ,
     3           4H    ,4H  LO,4HAD S,4HETS ,4H  L ,4HO A ,4HD   ,
     4           4HS E ,4HT   ,4HI D ,4HE N ,4HT I ,4HF I ,4HC A ,
     5           4HT I ,4HO N ,4H  N ,4HU M ,4HB E ,4HR S ,5*4H      /
      DATA IH6 / 9*4H         ,4HEQSS,4H ITE,4HM - ,4HSCAL,4HAR I,
     1           4HNDEX,4H LIS,4HT FO,4HR SU,4HBSTR,4HUCTU,4HRE  ,
     2           11*4H        ,4H INT,4HERNA,4HL   ,4H INT,4HERNA,
     3           4HL   ,4H  CO,4HMPON,4HENT ,2*4H         ,4H  IN,
     4           4HTERN,4HAL  ,4H  IN,4HTERN,4HAL  ,4H   C,4HOMPO,
     5           4HNENT,2*4H         ,4H   I,4HNTER,4HNAL ,4H   I,
     6           4HNTER,4HNAL ,4H    ,4HCOMP,4HONEN,4HT   ,4H POI,
     7           4HNT I,4HD   ,4H  SI,4HL ID,2*4H         ,4H DOF,
     8           3*4H         ,4H  PO,4HINT ,4HID  ,4H   S,4HIL I,
     9           4HD   ,4H    ,4H  DO,4HF   ,2*4H         ,4H   P,
     A           4HOINT,4H ID ,4H    ,4HSIL ,4HID  ,4H    ,4H   D,
     B           4HOF  ,4H     /
      DATA LOAP/ 4HLOAP/
C
      IST  = LOC
      IFIN = LOC + NW - 1
      GO TO (1,2,3,4,5,6,5,8), ICODE
C
C     EQSS ITEM
C
    1 DO 100 I = 1,96
  100 IHEAD(I) = IH1(I)
C
C     INSERT NAMES INTO HEADING
C
      IHEAD(17) = NAME1(1)
      IHEAD(18) = NAME1(2)
      IHEAD(22) = NAME2(1)
      IHEAD(23) = NAME2(2)
      CALL PAGE
      IF (NW .NE. 0) GO TO 140
      WRITE (OUTT,1009)
      GO TO 700
C
  140 DO 101 I = IST,IFIN,9
      NLINE = NLINE + 1
      IF (NLINE .LE. NLPP) GO TO 150
      CALL PAGE
      NLINE = NLINE + 1
  150 CONTINUE
      ICOMP = ANDF(IZ(I+2),63)
      CALL BITPAT (ICOMP,IBITS(1))
      I2 = 3
      IF (I+5 .GT. IFIN) GO TO 151
      ICOMP = ANDF(IZ(I+5),63)
      CALL BITPAT (ICOMP,IBITS(4))
      I2 = 6
      IF (I+8 .GT. IFIN) GO TO 151
      ICOMP = ANDF(IZ(I+8),63)
      CALL BITPAT (ICOMP,IBITS(7))
      I2 = 9
  151 CONTINUE
      WRITE (OUTT,1000) (IZ(I+J-1),IZ(I+J),IBITS(J),IBITS(J+1),J=1,I2,3)
  101 CONTINUE
      GO TO 700
C
C     EQSS - SCALER INDEX LIST
C
    8 DO 600 I = 1,96
  600 IHEAD(I) = IH6(I)
      IHEAD(22) = NAME1(1)
      IHEAD(23) = NAME1(2)
      CALL PAGE
C
      IP = 0
      DO 603 I = IST,IFIN,6
      NLINE = NLINE + 1
      IF (NLINE .LE. NLPP) GO TO 601
      CALL PAGE
      NLINE = NLINE + 1
  601 CONTINUE
      KCODE = IZ(I+1)
      CALL BITPAT (KCODE,IBITS(1))
      I2 = 2
      IPL(1) = IP + 1
      IF (I+3 .GT. IFIN) GO TO 602
      KCODE = IZ(I+3)
      CALL BITPAT (KCODE,IBITS(3))
      I2 = 4
      IPL(3) = IP + 2
      IF (I+5 .GT. IFIN) GO TO 602
      KCODE = IZ(I+5)
      CALL BITPAT (KCODE,IBITS(5))
      I2 = 6
      IPL(5) = IP + 3
  602 CONTINUE
      WRITE (OUTT,1000) (IPL(J),IZ(I+J-1),IBITS(J),IBITS(J+1),J=1,I2,2)
      IP = IP + 3
  603 CONTINUE
      GO TO 700
C
C     BGSS ITEM
C
    2 DO 200 I = 1,96
  200 IHEAD(I)  = IH2(I)
      IHEAD(20) = NAME1(1)
      IHEAD(21) = NAME1(2)
      CALL PAGE
      J = 0
      DO 201 I = IST,IFIN,4
      J = J + 1
      NLINE = NLINE + 1
      IF (NLINE .LE. NLPP) GO TO 250
      CALL PAGE
      NLINE = NLINE + 1
  250 CONTINUE
      WRITE (OUTT,1001) J,IZ(I),A(I+1),A(I+2),A(I+3)
  201 CONTINUE
      GO TO 700
C
C     CSTM ITEM
C
    3 DO 300 I = 1,96
  300 IHEAD(I)  = IH3(I)
      IHEAD(20) = NAME1(1)
      IHEAD(21) = NAME1(2)
      CALL PAGE
      DO 301 I = IST,IFIN,14
      NLINE = NLINE + 4
      IF (NLINE .LE. NLPP) GO TO 350
      CALL PAGE
      NLINE = NLINE + 4
  350 CONTINUE
      I1 = I + 2
      I2 = I + 13
      WRITE (OUTT,1002) IZ(I),IZ(I+1),(A(KK),KK= I1,I2)
  301 CONTINUE
      GO TO 700
C
    4 DO 400 I = 1,96
C
C     PLTS ITEM
C
  400 IHEAD(I)  = IH4(I)
      IHEAD(20) = NAME1(1)
      IHEAD(21) = NAME1(2)
      CALL PAGE
      DO 401 I = IST,IFIN,14
      NLINE = NLINE + 4
      IF (NLINE .LE. NLPP) GO TO 450
      CALL PAGE
      NLINE = NLINE + 4
  450 CONTINUE
      I1 = I + 2
      I2 = I + 13
      WRITE (OUTT,1004) IZ(I),IZ(I+1),(A(J),J=I1,I2)
  401 CONTINUE
      GO TO 700
C
C     LODS AND LOAP ITEMS
C
    5 DO 500 I = 1,96
  500 IHEAD(I)  = IH5(I)
      IHEAD(20) = NAME1(1)
      IHEAD(21) = NAME1(2)
      IF (ICODE .EQ. 7) IHEAD(13) = LOAP
      CALL PAGE
    6 IF (NW.EQ.0 .OR. NW.EQ.1) GO TO 520
      NL = NW/5 + 3
      NLINE = NLINE + NL
      IF (NLINE .LE. NLPP) GO TO 550
      CALL PAGE
      NLINE = NLINE + NL
  550 CONTINUE
      IST1 = IST + 1
      WRITE (OUTT,1006) NAME2(1),NAME2(2),IZ(IST),(IZ(J),J=IST1,IFIN)
      GO TO 700
C
  520 NLINE = NLINE + 2
      IF (NLINE .LE. NLPP) GO TO 560
      CALL PAGE
      NLINE = NLINE + 2
  560 CONTINUE
      WRITE (OUTT,1008) NAME2(1),NAME2(2)
  700 RETURN
C
 1000 FORMAT (6X,I8,4X,I8,6X,A4,A2,2(13X,I8,4X,I8,6X,A4,A2))
 1001 FORMAT (33X,I8,4X,I8,3X,3(3X,E13.6))
 1002 FORMAT (/10X,I8,3X,I4,3X,3(3X,E13.6),4X,3(3X,E13.6),
     1        /80X,3(3X,E13.6), /80X,3(3X,E13.6))
 1004 FORMAT (/14X,2A4,3X,3(3X,E13.6),4X,3(3X,E13.6)
     1        /77X,3(3X,E13.6), /77X,3(3X,E13.6))
 1006 FORMAT (/26X,2A4,3X,I8,5X,6(2X,I8)/(50X,2X,I8,2X,I8,2X,I8,2X,I8,
     1        2X,I8,2X,I8,/))
 1008 FORMAT (/26X,2A4,17X,32HNO LOAD SETS FOR THIS COMPONENT. )
 1009 FORMAT (/30X,64HALL DEGREES OF FREEDOM FOR THIS COMPONENT HAVE BEE
     1N REDUCED OUT. )
      END