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

!>    @brief commonly used global variables and constants
module mod_genrl
      double precision tslocal
!

      !> @brief Memory of data in 8Bytes (std. double precision size).
      !> @details
      !> Memory of data in 8Bytes (std. double precision size). \n
      double precision memory

!     file output compression tool
      integer compress_out

      !> @brief Number of cells in x-direction.
      !> @details
      !> Number of cells in x-direction. \n
      !> Read in under `# grid`, first entry.
      integer :: i0

      !> @brief Number of cells in y-direction.
      !> @details
      !> Number of cells in y-direction. \n
      !> Read in under `# grid`, second entry.
      integer :: j0

      !> @brief Number of cells in z-direction.
      !> @details
      !> Number of cells in z-direction. \n
      !> Read in under `# grid`, third entry.
      integer :: k0

      !> @brief Output orientation of monitoring output.
      !> @details
      !> Output orientation of monitoring output. \n \n
      !> `out_orientation` values: \n
      !>   0: normal (pv column wise), each time step has its own block
      !>      (mp row wise) \n
      !>   1: transposed (pv row wise), each time step has its own
      !>      block (mp column wise) \n
      !>   2: normal (pv column wise), new file for each time step (mp
      !>      row wise) \n
      !>   3: transposed (pv row wise), new file for each time step
      !>      (mp column wise) \n
      !>   4: normal (pv column wise), new file for each monitor point
      !>      (time row wise) \n
      integer :: out_orientation

      !> @brief Runmode integer.
      !> @details
      !> Runmode integer. \n
      !> Read in under `# runmode`. \n\n
      !>
      !> 0: forward modeling \n
      !> 1: forward modeling and data fit \n
      !> 2: inverse modeling or extra steady state \n
      !> 3: inverse modeling \n
      integer :: runmode

      !> @brief Maximum number of nonlinear iterations.
      !> @details
      !> Maximum number of nonlinear iterations. \n
      !> Important for the Picard iterations.\n\n
      !>
      !> Read under `# nlsolve`, first entry.
      integer :: maxiter_nl

      integer nladapt

      !> @brief Switch for checking non-linear iteration convergence.
      !> @details
      !> Switch for checking non-linear iteration convergence. \n
      !> - 0: check convergence (default).\n
      !> - other: do not check convergence.\n
      !>
      !> Read under `# nlsolve`, second line, only entry.
      integer nlconverge
!
      logical head_active, temp_active
      logical trac_active,chem_active,trans_active, pres_active
      logical vtk_out,tec_out,hdf_out,txt_out

      !> @brief Switch informing if an initial flow transformation is needed.
      !> @details
      !>  Switch informing if an initial flow transformation is needed. \n
      !> This variable is .false. if the not-computed flow-variable
      !> (f.e. pres in head-based computation) is given in the input
      !> file. Otherwise `head2pres` or `pres2head` will be called
      !> inside /forward/forward_init.f90.
      logical is_init_flow_trafo_needed

      !> @brief Switch for disabling multiple outputs.
      !> @details
      !> Switch for disabling multiple outputs. \n
      !> Disables output in write_logs, write_data, write_outt,
      !> write_tecdiff, write_monitor.
      logical write_disable
      logical write_smonitor, write_param, write_eoutt

      !> @brief Switch for disabling iteration output.
      !> @details
      !> Switch for disabling iteration output. \n
      !> Disables output from forward iteration.
      logical write_iter_disable

      !> @brief Switch for reading external input files.
      !> @details
      !> Switch for reading external input files. \n
      !> If and only if read_external_input is .true., the input file
      !> is searched for external input files names in the subroutine
      !> read_control.
      logical read_external_input

      !> @brief Hashtag character.
      !> @details
      !> Hashtag character. \n
      !> Used as input keyword signifier.
      character (len=1), parameter :: key_char = '#'

      !> @brief Number of samples.
      !> @details
      !> Number of samples. \n
      !> For example, the number of samples in a Monte Carlo
      !> simulation. Also used in OpenMP parallelization.
      integer nsmpl

      ! Useful constants
      ! ----------------

      !> @brief Very big number.
      !> @details
      !> Very big number. \n\n
      double precision, parameter :: big = 1.0d30

      !> @brief Very small number.
      !> @details
      !> Very small number. \n\n
      double precision, parameter :: small = 1.0d-300

      !> @brief Circle number pi.
      !> @details
      !> Circle number pi. \n\n
      double precision, parameter :: pi = 3.141592653589793d0

      !> @brief Flag for log permeability output.
      !> @details
      !> Flag for log permeability output. \n\n
      logical, parameter :: klogflag = .true.

      !> @brief Number of possible old saved arrays.
      !> @details
      !> Number of possible old saved arrays. \n\n
      !>
      !> max number of "cgen_*"
      integer, parameter :: ncgen = 5

      !> @brief Index constant for old saved arrays.
      !> @details
      !> Index constant for old saved arrays. \n\n
      !>
      !> index constants for HEAD, TEMP and PRES
      !> forward iteration [x_i-1,j];
      !> do not change it without a change of any
      !> "headold" into "headold(1,cgen_fw)" ccc
      integer, parameter :: cgen_fw = 1

      !> @brief Index constant for old saved arrays.
      !> @details
      !> Index constant for old saved arrays. \n\n
      !>
      !> Used around nonlinear forward iteration in `forward_iter`.
      !>
      !> index constants for HEAD, TEMP and PRES
      !> forward iteration [x_i-1,j];
      !> do not change it without a change of any
      !> "headold" into "headold(1,cgen_fw)" ccc\n\n
      !>
      !>       time iteration
      integer, parameter :: cgen_time = 2

      !> @brief Index constant for old saved arrays.
      !> @details
      !> Index constant for old saved arrays. \n\n
      !>
      !> index constants for HEAD, TEMP and PRES
      !> forward iteration [x_i-1,j];
      !> do not change it without a change of any
      !> "headold" into "headold(1,cgen_fw)" ccc\n\n
      !>
      !>       nonlinear-forward iteration (use of ad)
      integer, parameter :: cgen_fwad = 3
!
!       optimisation/realisation iteration [x_*,j]
!         index constants for g_HEAD, g_TEMP and g_PRES
        integer, parameter :: cgen_opti = 4
!       special DD use [f(x_*,j)]
        integer, parameter :: cgen_dd = 5

      !       index of the <master> sample
      integer, parameter :: idx_master=1

        !> @brief Doubling threshold for variable time step size
        !> @details
        !> Doubling threshold for variable time step size. \n\n
        !>
        !> Once delt_count reaches delt_double, the time step size is
        !> doubled. \n\n
        !>
        !> Convergence flag, doubling value for variable time step
        !> size
        integer delt_double

        !> @brief Variable step size nonlinear iteration counter.
        !> @details
        !> Variable step size nonlinear iteration counter. \n\n
        !>
        !> This is a helper-counter that together with the current
        !> number of iterations (in `forward_picard`) decides whether
        !> the doubling counter of the variable time step size utitily
        !> should be moved towards or away from doubling the time step
        !> size.
        integer :: iter_nlold

        !> @brief Flag for Variable step size
        !> @details
        !> Flag for Variable step size. \n\n
        !>
        !> True if variable step size is used.
        logical :: delt_vary

        !> @brief Initial time step size for variable step size
        !> @details
        !> Initial time step size for variable step size. \n\n
        !>
        !> First input under `# variable step size`.
        double precision :: delt_start

        !> @brief Minimum time step size for variable step size
        !> @details
        !> Minimum time step size for variable step size. \n\n
        !>
        !> Second input under `# variable step size`.
        double precision :: delt_min

        !> @brief Maximum time step size for variable step size
        !> @details
        !> Maximum time step size for variable step size. \n\n
        !>
        !> Third input under `# variable step size`.
        double precision :: delt_max

end module mod_genrl


text:
module m3
end module m3

module m1
    use m2
    use m3
!$claw
end module m1

max_stars_repo_name: clementval/claw-compiler

max_issues_repo_name: skosukhin/claw-compiler

max_forks_repo_name: skosukhin/claw-compiler


text:
*DECK TINVIT
      SUBROUTINE TINVIT (NM, N, D, E, E2, M, W, IND, Z, IERR, RV1, RV2,
     +   RV3, RV4, RV6)
C***BEGIN PROLOGUE  TINVIT
C***PURPOSE  Compute the eigenvectors of symmetric tridiagonal matrix
C            corresponding to specified eigenvalues, using inverse
C            iteration.
C***LIBRARY   SLATEC (EISPACK)
C***CATEGORY  D4C3
C***TYPE      SINGLE PRECISION (TINVIT-S)
C***KEYWORDS  EIGENVECTORS, EISPACK
C***AUTHOR  Smith, B. T., et al.
C***DESCRIPTION
C
C     This subroutine is a translation of the inverse iteration tech-
C     nique in the ALGOL procedure TRISTURM by Peters and Wilkinson.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
C
C     This subroutine finds those eigenvectors of a TRIDIAGONAL
C     SYMMETRIC matrix corresponding to specified eigenvalues,
C     using inverse iteration.
C
C     On Input
C
C        NM must be set to the row dimension of the two-dimensional
C          array parameter, Z, as declared in the calling program
C          dimension statement.  NM is an INTEGER variable.
C
C        N is the order of the matrix.  N is an INTEGER variable.
C          N must be less than or equal to NM.
C
C        D contains the diagonal elements of the symmetric tridiagonal
C          matrix.  D is a one-dimensional REAL array, dimensioned D(N).
C
C        E contains the subdiagonal elements of the symmetric
C          tridiagonal matrix in its last N-1 positions.  E(1) is
C          arbitrary.  E is a one-dimensional REAL array, dimensioned
C          E(N).
C
C        E2 contains the squares of the corresponding elements of E,
C          with zeros corresponding to negligible elements of E.
C          E(I) is considered negligible if it is not larger than
C          the product of the relative machine precision and the sum
C          of the magnitudes of D(I) and D(I-1).  E2(1) must contain
C          0.0e0 if the eigenvalues are in ascending order, or 2.0e0
C          if the eigenvalues are in descending order.  If  BISECT,
C          TRIDIB, or  IMTQLV  has been used to find the eigenvalues,
C          their output E2 array is exactly what is expected here.
C          E2 is a one-dimensional REAL array, dimensioned E2(N).
C
C        M is the number of specified eigenvalues for which eigenvectors
C          are to be determined.  M is an INTEGER variable.
C
C        W contains the M eigenvalues in ascending or descending order.
C          W is a one-dimensional REAL array, dimensioned W(M).
C
C        IND contains in its first M positions the submatrix indices
C          associated with the corresponding eigenvalues in W --
C          1 for eigenvalues belonging to the first submatrix from
C          the top, 2 for those belonging to the second submatrix, etc.
C          If  BISECT  or  TRIDIB  has been used to determine the
C          eigenvalues, their output IND array is suitable for input
C          to TINVIT.  IND is a one-dimensional INTEGER array,
C          dimensioned IND(M).
C
C     On Output
C
C       ** All input arrays are unaltered.**
C
C        Z contains the associated set of orthonormal eigenvectors.
C          Any vector which fails to converge is set to zero.
C          Z is a two-dimensional REAL array, dimensioned Z(NM,M).
C
C        IERR is an INTEGER flag set to
C          Zero       for normal return,
C          -J         if the eigenvector corresponding to the J-th
C                     eigenvalue fails to converge in 5 iterations.
C
C        RV1, RV2 and RV3 are one-dimensional REAL arrays used for
C          temporary storage.  They are used to store the main diagonal
C          and the two adjacent diagonals of the triangular matrix
C          produced in the inverse iteration process.  RV1, RV2 and
C          RV3 are dimensioned RV1(N), RV2(N) and RV3(N).
C
C        RV4 and RV6 are one-dimensional REAL arrays used for temporary
C          storage.  RV4 holds the multipliers of the Gaussian
C          elimination process.  RV6 holds the approximate eigenvectors
C          in this process.  RV4 and RV6 are dimensioned RV4(N) and
C          RV6(N).
C
C     Questions and comments should be directed to B. S. Garbow,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C     ------------------------------------------------------------------
C
C***REFERENCES  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
C                 Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
C                 system Routines - EISPACK Guide, Springer-Verlag,
C                 1976.
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   760101  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  TINVIT
C
      INTEGER I,J,M,N,P,Q,R,S,II,IP,JJ,NM,ITS,TAG,IERR,GROUP
      INTEGER IND(*)
      REAL D(*),E(*),E2(*),W(*),Z(NM,*)
      REAL RV1(*),RV2(*),RV3(*),RV4(*),RV6(*)
      REAL U,V,UK,XU,X0,X1,EPS2,EPS3,EPS4,NORM,ORDER
C
C***FIRST EXECUTABLE STATEMENT  TINVIT
      IERR = 0
      IF (M .EQ. 0) GO TO 1001
      TAG = 0
      ORDER = 1.0E0 - E2(1)
      Q = 0
C     .......... ESTABLISH AND PROCESS NEXT SUBMATRIX ..........
  100 P = Q + 1
C
      DO 120 Q = P, N
         IF (Q .EQ. N) GO TO 140
         IF (E2(Q+1) .EQ. 0.0E0) GO TO 140
  120 CONTINUE
C     .......... FIND VECTORS BY INVERSE ITERATION ..........
  140 TAG = TAG + 1
      S = 0
C
      DO 920 R = 1, M
         IF (IND(R) .NE. TAG) GO TO 920
         ITS = 1
         X1 = W(R)
         IF (S .NE. 0) GO TO 510
C     .......... CHECK FOR ISOLATED ROOT ..........
         XU = 1.0E0
         IF (P .NE. Q) GO TO 490
         RV6(P) = 1.0E0
         GO TO 870
  490    NORM = ABS(D(P))
         IP = P + 1
C
         DO 500 I = IP, Q
  500    NORM = MAX(NORM, ABS(D(I)) + ABS(E(I)))
C     .......... EPS2 IS THE CRITERION FOR GROUPING,
C                EPS3 REPLACES ZERO PIVOTS AND EQUAL
C                ROOTS ARE MODIFIED BY EPS3,
C                EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW ..........
         EPS2 = 1.0E-3 * NORM
         EPS3 = NORM
  502    EPS3 = 0.5E0*EPS3
         IF (NORM + EPS3 .GT. NORM) GO TO 502
         UK = SQRT(REAL(Q-P+5))
         EPS3 = UK * EPS3
         EPS4 = UK * EPS3
         UK = EPS4 / UK
         S = P
  505    GROUP = 0
         GO TO 520
C     .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
  510    IF (ABS(X1-X0) .GE. EPS2) GO TO 505
         GROUP = GROUP + 1
         IF (ORDER * (X1 - X0) .LE. 0.0E0) X1 = X0 + ORDER * EPS3
C     .......... ELIMINATION WITH INTERCHANGES AND
C                INITIALIZATION OF VECTOR ..........
  520    V = 0.0E0
C
         DO 580 I = P, Q
            RV6(I) = UK
            IF (I .EQ. P) GO TO 560
            IF (ABS(E(I)) .LT. ABS(U)) GO TO 540
C     .......... WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
C                E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY ..........
            XU = U / E(I)
            RV4(I) = XU
            RV1(I-1) = E(I)
            RV2(I-1) = D(I) - X1
            RV3(I-1) = 0.0E0
            IF (I .NE. Q) RV3(I-1) = E(I+1)
            U = V - XU * RV2(I-1)
            V = -XU * RV3(I-1)
            GO TO 580
  540       XU = E(I) / U
            RV4(I) = XU
            RV1(I-1) = U
            RV2(I-1) = V
            RV3(I-1) = 0.0E0
  560       U = D(I) - X1 - XU * V
            IF (I .NE. Q) V = E(I+1)
  580    CONTINUE
C
         IF (U .EQ. 0.0E0) U = EPS3
         RV1(Q) = U
         RV2(Q) = 0.0E0
         RV3(Q) = 0.0E0
C     .......... BACK SUBSTITUTION
C                FOR I=Q STEP -1 UNTIL P DO -- ..........
  600    DO 620 II = P, Q
            I = P + Q - II
            RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
            V = U
            U = RV6(I)
  620    CONTINUE
C     .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS
C                MEMBERS OF GROUP ..........
         IF (GROUP .EQ. 0) GO TO 700
         J = R
C
         DO 680 JJ = 1, GROUP
  630       J = J - 1
            IF (IND(J) .NE. TAG) GO TO 630
            XU = 0.0E0
C
            DO 640 I = P, Q
  640       XU = XU + RV6(I) * Z(I,J)
C
            DO 660 I = P, Q
  660       RV6(I) = RV6(I) - XU * Z(I,J)
C
  680    CONTINUE
C
  700    NORM = 0.0E0
C
         DO 720 I = P, Q
  720    NORM = NORM + ABS(RV6(I))
C
         IF (NORM .GE. 1.0E0) GO TO 840
C     .......... FORWARD SUBSTITUTION ..........
         IF (ITS .EQ. 5) GO TO 830
         IF (NORM .NE. 0.0E0) GO TO 740
         RV6(S) = EPS4
         S = S + 1
         IF (S .GT. Q) S = P
         GO TO 780
  740    XU = EPS4 / NORM
C
         DO 760 I = P, Q
  760    RV6(I) = RV6(I) * XU
C     .......... ELIMINATION OPERATIONS ON NEXT VECTOR
C                ITERATE ..........
  780    DO 820 I = IP, Q
            U = RV6(I)
C     .......... IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
C                WAS PERFORMED EARLIER IN THE
C                TRIANGULARIZATION PROCESS ..........
            IF (RV1(I-1) .NE. E(I)) GO TO 800
            U = RV6(I-1)
            RV6(I-1) = RV6(I)
  800       RV6(I) = U - RV4(I) * RV6(I-1)
  820    CONTINUE
C
         ITS = ITS + 1
         GO TO 600
C     .......... SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
  830    IERR = -R
         XU = 0.0E0
         GO TO 870
C     .......... NORMALIZE SO THAT SUM OF SQUARES IS
C                1 AND EXPAND TO FULL ORDER ..........
  840    U = 0.0E0
C
         DO 860 I = P, Q
  860    U = U + RV6(I)**2
C
         XU = 1.0E0 / SQRT(U)
C
  870    DO 880 I = 1, N
  880    Z(I,R) = 0.0E0
C
         DO 900 I = P, Q
  900    Z(I,R) = RV6(I) * XU
C
         X0 = X1
  920 CONTINUE
C
      IF (Q .LT. N) GO TO 100
 1001 RETURN
      END

ext: f


text:
      integer function f()
      f=77
      end function f

      program t
      implicit none
      integer::f
      write (*,'(i0)') f()
      end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
!
!    version 1.0
!    built Fri Apr 16 18:03:05 CDT 2004
!
!  by B.R. Fabijonas
!     Department of Mathematics
!     Southern Methodist University
!     bfabi@smu.edu
!
!  see the file airy_README for an explanation
!***
!************************************************************************
!***
      module airy_functions_real_single
        implicit none
        private
        public :: airy_ai, airy_bi, airy_ai_zero, airy_bi_zero
        public :: airy_info, airy_aux, airy_aux_info
        integer, parameter :: prd = kind(0.0e0)
        include 'airy_head'
!*
!  interfaces
!*
        interface airy_ai
          module procedure airy_air
        end interface
        interface airy_bi
          module procedure airy_bir
        end interface
        interface airy_ai_zero
  	  module procedure ai_zeror
          module procedure ai_zerorv
        end interface
        interface airy_bi_zero
  	  module procedure bi_zeror
          module procedure bi_zerorv
        end interface
        interface airy_info
          module procedure airy_paramsr
        end interface
        interface airy_aux
  	  module procedure airy_auxr
        end interface
        interface airy_aux_info
          module procedure airy_paramsr_aux
        end interface
        contains
         include 'airy_real'
         include 'airy_parameters'
      end module airy_functions_real_single
!***
!************************************************************************
!***
      module airy_functions_complex_single
        implicit none
        private
        public :: airy_ai, airy_ai_zero, airy_bi_zero
        integer, parameter :: prd = kind(0.0e0)
        include 'airy_head'
!*
!  interfaces
!*
        interface airy_ai
          module procedure airy_aic
          module procedure airy_ai_rayc
        end interface
        interface airy_ai_zero
   	  module procedure ai_zeroc
          module procedure ai_zerocv
        end interface
        interface airy_bi_zero
   	  module procedure bi_zeroc
          module procedure bi_zerocv
        end interface
        contains
         include 'airy_complex'
         include 'airy_parameters'
      end module airy_functions_complex_single
!***
!************************************************************************
!***
      module airy_functions_real_double
        implicit none
        private
        public :: airy_ai, airy_bi, airy_ai_zero, airy_bi_zero
        public :: airy_info, airy_aux, airy_aux_info
        integer, parameter :: prd = kind(0.0d0)
        include 'airy_head'
!*
!  interfaces
!*
        interface airy_ai
          module procedure airy_air
        end interface
        interface airy_bi
          module procedure airy_bir
        end interface
        interface airy_ai_zero
  	  module procedure ai_zeror
          module procedure ai_zerorv
        end interface
        interface airy_bi_zero
  	  module procedure bi_zeror
          module procedure bi_zerorv
        end interface
        interface airy_info
          module procedure airy_paramsr
        end interface
        interface airy_aux
  	  module procedure airy_auxr
        end interface
        interface airy_aux_info
          module procedure airy_paramsr_aux
        end interface
        contains
         include 'airy_real'
         include 'airy_parameters'
      end module airy_functions_real_double
!***
!************************************************************************
!***
      module airy_functions_complex_double
        implicit none
        private
        public :: airy_ai, airy_ai_zero, airy_bi_zero
        integer, parameter :: prd = kind(0.0d0)
        include 'airy_head'
!*
!  interfaces
!*
        interface airy_ai
          module procedure airy_aic
          module procedure airy_ai_rayc
        end interface
        interface airy_ai_zero
   	  module procedure ai_zeroc
          module procedure ai_zerocv
        end interface
        interface airy_bi_zero
   	  module procedure bi_zeroc
          module procedure bi_zerocv
        end interface
        contains
         include 'airy_complex'
         include 'airy_parameters'
      end module airy_functions_complex_double
!***
!************************************************************************
!***
!      module airy_functions_real_quad
!        implicit none
!        private
!        public :: airy_ai, airy_bi, airy_ai_zero, airy_bi_zero
!        public :: airy_info, airy_aux, airy_aux_info
!        integer, parameter :: prd = kind(0.0q0)
!        include 'airy_head'
!!*
!!  interfaces
!!*
!        interface airy_ai
!          module procedure airy_air
!        end interface
!        interface airy_bi
!          module procedure airy_bir
!        end interface
!        interface airy_ai_zero
!  	  module procedure ai_zeror
!          module procedure ai_zerorv
!        end interface
!        interface airy_bi_zero
!  	  module procedure bi_zeror
!          module procedure bi_zerorv
!        end interface
!        interface airy_info
!          module procedure airy_paramsr
!        end interface
!        interface airy_aux
!  	  module procedure airy_auxr
!        end interface
!        interface airy_aux_info
!          module procedure airy_paramsr_aux
!        end interface
!        contains
!         include 'airy_real'
!         include 'airy_parameters'
!      end module airy_functions_real_quad
!!***
!!************************************************************************
!!***
!      module airy_functions_complex_quad
!        implicit none
!        private
!        public :: airy_ai, airy_ai_zero, airy_bi_zero
!        integer, parameter :: prd = kind(0.0q0)
!        include 'airy_head'
!!*
!!  interfaces
!!*
!        interface airy_ai
!          module procedure airy_aic
!          module procedure airy_ai_rayc
!        end interface
!        interface airy_ai_zero
!   	  module procedure ai_zeroc
!          module procedure ai_zerocv
!        end interface
!        interface airy_bi_zero
!   	  module procedure bi_zeroc
!          module procedure bi_zerocv
!        end interface
!        contains
!         include 'airy_complex'
!         include 'airy_parameters'
!      end module airy_functions_complex_quad
!!***
!!************************************************************************
!!***

max_stars_repo_name: IRFM/ALOHA

max_issues_repo_name: IRFM/ALOHA

max_forks_repo_name: IRFM/ALOHA


text:
* ============================================================================ *
      SUBROUTINE DSM(M,N,NPAIRS,INDROW,INDCOL,NGRP,MAXGRP,MINGRP,
     *               INFO,IPNTR,JPNTR,IWA,LIWA)
      INTEGER M,N,NPAIRS,MAXGRP,MINGRP,INFO,LIWA
      INTEGER INDROW(NPAIRS),INDCOL(NPAIRS),NGRP(N),
     *        IPNTR(M+1),JPNTR(N+1),IWA(LIWA)
C     **********
C
C     SUBROUTINE DSM
C
C     THE PURPOSE OF DSM IS TO DETERMINE AN OPTIMAL OR NEAR-
C     OPTIMAL CONSISTENT PARTITION OF THE COLUMNS OF A SPARSE
C     M BY N MATRIX A.
C
C     THE SPARSITY PATTERN OF THE MATRIX A IS SPECIFIED BY
C     THE ARRAYS INDROW AND INDCOL. ON INPUT THE INDICES
C     FOR THE NON-ZERO ELEMENTS OF A ARE
C
C           INDROW(K),INDCOL(K), K = 1,2,...,NPAIRS.
C
C     THE (INDROW,INDCOL) PAIRS MAY BE SPECIFIED IN ANY ORDER.
C     DUPLICATE INPUT PAIRS ARE PERMITTED, BUT THE SUBROUTINE
C     ELIMINATES THEM.
C
C     THE SUBROUTINE PARTITIONS THE COLUMNS OF A INTO GROUPS
C     SUCH THAT COLUMNS IN THE SAME GROUP DO NOT HAVE A
C     NON-ZERO IN THE SAME ROW POSITION. A PARTITION OF THE
C     COLUMNS OF A WITH THIS PROPERTY IS CONSISTENT WITH THE
C     DIRECT DETERMINATION OF A.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE DSM(M,N,NPAIRS,INDROW,INDCOL,NGRP,MAXGRP,MINGRP,
C                      INFO,IPNTR,JPNTR,IWA,LIWA)
C
C     WHERE
C
C       M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF ROWS OF A.
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       NPAIRS IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE
C         NUMBER OF (INDROW,INDCOL) PAIRS USED TO DESCRIBE THE
C         SPARSITY PATTERN OF A.
C
C       INDROW IS AN INTEGER ARRAY OF LENGTH NPAIRS. ON INPUT INDROW
C         MUST CONTAIN THE ROW INDICES OF THE NON-ZERO ELEMENTS OF A.
C         ON OUTPUT INDROW IS PERMUTED SO THAT THE CORRESPONDING
C         COLUMN INDICES ARE IN NON-DECREASING ORDER. THE COLUMN
C         INDICES CAN BE RECOVERED FROM THE ARRAY JPNTR.
C
C       INDCOL IS AN INTEGER ARRAY OF LENGTH NPAIRS. ON INPUT INDCOL
C         MUST CONTAIN THE COLUMN INDICES OF THE NON-ZERO ELEMENTS OF
C         A. ON OUTPUT INDCOL IS PERMUTED SO THAT THE CORRESPONDING
C         ROW INDICES ARE IN NON-DECREASING ORDER. THE ROW INDICES
C         CAN BE RECOVERED FROM THE ARRAY IPNTR.
C
C       NGRP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE PARTITION OF THE COLUMNS OF A. COLUMN JCOL BELONGS
C         TO GROUP NGRP(JCOL).
C
C       MAXGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES THE
C         NUMBER OF GROUPS IN THE PARTITION OF THE COLUMNS OF A.
C
C       MINGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES A LOWER
C         BOUND FOR THE NUMBER OF GROUPS IN ANY CONSISTENT PARTITION
C         OF THE COLUMNS OF A.
C
C       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS. FOR
C         NORMAL TERMINATION INFO = 1. IF M, N, OR NPAIRS IS NOT
C         POSITIVE OR LIWA IS LESS THAN MAX(M,6*N), THEN INFO = 0.
C         IF THE K-TH ELEMENT OF INDROW IS NOT AN INTEGER BETWEEN
C         1 AND M OR THE K-TH ELEMENT OF INDCOL IS NOT AN INTEGER
C         BETWEEN 1 AND N, THEN INFO = -K.
C
C       IPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH M + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL.
C         THE COLUMN INDICES FOR ROW I ARE
C
C               INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1.
C
C         NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       JPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW.
C         THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       IWA IS AN INTEGER WORK ARRAY OF LENGTH LIWA.
C
C       LIWA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN
C         MAX(M,6*N).
C
C     SUBPROGRAMS CALLED
C
C       MINPACK-SUPPLIED ... DEGR,IDO,NUMSRT,SEQ,SETR,SLO,SRTDAT
C
C       FORTRAN-SUPPLIED ... MAX
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER I,IR,J,JP,K,MAXCLQ,NNZ,NUMGRP
C
C     CHECK THE INPUT DATA.
C
      INFO = 0
      IF (M .LT. 1 .OR. N .LT. 1 .OR. NPAIRS .LT. 1 .OR.
     *    LIWA .LT. MAX(M,6*N)) RETURN
      DO 10 K = 1, NPAIRS
         INFO = -K
         IF (INDROW(K) .LT. 1 .OR. INDROW(K) .GT. M .OR.
     *       INDCOL(K) .LT. 1 .OR. INDCOL(K) .GT. N) RETURN
   10    CONTINUE
      INFO = 1
C
C     SORT THE DATA STRUCTURE BY COLUMNS.
C
      CALL SRTDAT(N,NPAIRS,INDROW,INDCOL,JPNTR,IWA)
C
C     COMPRESS THE DATA AND DETERMINE THE NUMBER OF
C     NON-ZERO ELEMENTS OF A.
C
      DO 20 I = 1, M
         IWA(I) = 0
   20    CONTINUE
      NNZ = 1
      DO 40 J = 1, N
         K = NNZ
         DO 30 JP = JPNTR(J), JPNTR(J+1)-1
            IR = INDROW(JP)
            IF (IWA(IR) .NE. J) THEN
               INDROW(NNZ) = IR
               NNZ = NNZ + 1
               IWA(IR) = J
               END IF
   30       CONTINUE
         JPNTR(J) = K
   40    CONTINUE
      JPNTR(N+1) = NNZ
C
C     EXTEND THE DATA STRUCTURE TO ROWS.
C
      CALL SETR(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA)
C
C     DETERMINE A LOWER BOUND FOR THE NUMBER OF GROUPS.
C
      MINGRP = 0
      DO 50 I = 1, M
         MINGRP = MAX(MINGRP,IPNTR(I+1)-IPNTR(I))
   50    CONTINUE
C
C     DETERMINE THE DEGREE SEQUENCE FOR THE INTERSECTION
C     GRAPH OF THE COLUMNS OF A.
C
      CALL DEGR(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(5*N+1),IWA(N+1))
C
C     COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF A
C     WITH THE SMALLEST-LAST (SL) ORDERING.
C
      CALL SLO(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(5*N+1),IWA(4*N+1),
     *         MAXCLQ,IWA(1),IWA(N+1),IWA(2*N+1),IWA(3*N+1))
      CALL SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),NGRP,MAXGRP,
     *         IWA(N+1))
      MINGRP = MAX(MINGRP,MAXCLQ)
C
C     EXIT IF THE SMALLEST-LAST ORDERING IS OPTIMAL.
C
      IF (MAXGRP .EQ. MINGRP) RETURN
C
C     COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF A
C     WITH THE INCIDENCE-DEGREE (ID) ORDERING.
C
      CALL IDO(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA(5*N+1),IWA(4*N+1),
     *         MAXCLQ,IWA(1),IWA(N+1),IWA(2*N+1),IWA(3*N+1))
      CALL SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),IWA(1),NUMGRP,
     *         IWA(N+1))
      MINGRP = MAX(MINGRP,MAXCLQ)
C
C     RETAIN THE BETTER OF THE TWO ORDERINGS SO FAR.
C
      IF (NUMGRP .LT. MAXGRP) THEN
         MAXGRP = NUMGRP
         DO 60 J = 1, N
            NGRP(J) = IWA(J)
   60       CONTINUE
C
C        EXIT IF THE INCIDENCE-DEGREE ORDERING IS OPTIMAL.
C
         IF (MAXGRP .EQ. MINGRP) RETURN
         END IF
C
C     COLOR THE INTERSECTION GRAPH OF THE COLUMNS OF A
C     WITH THE LARGEST-FIRST (LF) ORDERING.
C
      CALL NUMSRT(N,N-1,IWA(5*N+1),-1,IWA(4*N+1),IWA(2*N+1),IWA(N+1))
      CALL SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,IWA(4*N+1),IWA(1),NUMGRP,
     *         IWA(N+1))
C
C     RETAIN THE BEST OF THE THREE ORDERINGS AND EXIT.
C
      IF (NUMGRP .LT. MAXGRP) THEN
         MAXGRP = NUMGRP
         DO 70 J = 1, N
            NGRP(J) = IWA(J)
   70       CONTINUE
         END IF
      RETURN
C
C     LAST CARD OF SUBROUTINE DSM.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE DEGR(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,IWA)
      INTEGER N
      INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(*),NDEG(N),IWA(N)
C     **********
C
C     SUBROUTINE DEGR
C
C     GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A,
C     THIS SUBROUTINE DETERMINES THE DEGREE SEQUENCE FOR
C     THE INTERSECTION GRAPH OF THE COLUMNS OF A.
C
C     IN GRAPH-THEORY TERMINOLOGY, THE INTERSECTION GRAPH OF
C     THE COLUMNS OF A IS THE LOOPLESS GRAPH G WITH VERTICES
C     A(J), J = 1,2,...,N WHERE A(J) IS THE J-TH COLUMN OF A
C     AND WITH EDGE (A(I),A(J)) IF AND ONLY IF COLUMNS I AND J
C     HAVE A NON-ZERO IN THE SAME ROW POSITION.
C
C     NOTE THAT THE VALUE OF M IS NOT NEEDED BY DEGR AND IS
C     THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE DEGR(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,IWA)
C
C     WHERE
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW
C         INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW.
C         THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE
C         COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL.
C         THE COLUMN INDICES FOR ROW I ARE
C
C               INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1.
C
C         NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       NDEG IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH
C         SPECIFIES THE DEGREE SEQUENCE. THE DEGREE OF THE
C         J-TH COLUMN OF A IS NDEG(J).
C
C       IWA IS AN INTEGER WORK ARRAY OF LENGTH N.
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER IC,IP,IR,JCOL,JP
C
C     INITIALIZATION BLOCK.
C
      DO 10 JP = 1, N
         NDEG(JP) = 0
         IWA(JP) = 0
   10    CONTINUE
C
C     COMPUTE THE DEGREE SEQUENCE BY DETERMINING THE CONTRIBUTIONS
C     TO THE DEGREES FROM THE CURRENT(JCOL) COLUMN AND FURTHER
C     COLUMNS WHICH HAVE NOT YET BEEN CONSIDERED.
C
      DO 40 JCOL = 2, N
         IWA(JCOL) = N
C
C        DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND
C        TO NON-ZEROES IN THE MATRIX.
C
         DO 30 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1
            IR = INDROW(JP)
C
C           FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC)
C           WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX.
C
            DO 20 IP = IPNTR(IR), IPNTR(IR+1)-1
               IC = INDCOL(IP)
C
C              ARRAY IWA MARKS COLUMNS WHICH HAVE CONTRIBUTED TO
C              THE DEGREE COUNT OF COLUMN JCOL. UPDATE THE DEGREE
C              COUNTS OF THESE COLUMNS AS WELL AS COLUMN JCOL.
C
               IF (IWA(IC) .LT. JCOL) THEN
                  IWA(IC) = JCOL
                  NDEG(IC) = NDEG(IC) + 1
                  NDEG(JCOL) = NDEG(JCOL) + 1
                  END IF
   20          CONTINUE
   30       CONTINUE
   40    CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE DEGR.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE IDO(M,N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST,
     *               MAXCLQ,IWA1,IWA2,IWA3,IWA4)
      INTEGER M,N,MAXCLQ
      INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(M+1),NDEG(N),
     *        LIST(N),IWA1(0:N-1),IWA2(N),IWA3(N),IWA4(N)
C     **********
C
C     SUBROUTINE IDO
C
C     GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, THIS
C     SUBROUTINE DETERMINES AN INCIDENCE-DEGREE ORDERING OF THE
C     COLUMNS OF A.
C
C     THE INCIDENCE-DEGREE ORDERING IS DEFINED FOR THE LOOPLESS
C     GRAPH G WITH VERTICES A(J), J = 1,2,...,N WHERE A(J) IS THE
C     J-TH COLUMN OF A AND WITH EDGE (A(I),A(J)) IF AND ONLY IF
C     COLUMNS I AND J HAVE A NON-ZERO IN THE SAME ROW POSITION.
C
C     THE INCIDENCE-DEGREE ORDERING IS DETERMINED RECURSIVELY BY
C     LETTING LIST(K), K = 1,...,N BE A COLUMN WITH MAXIMAL
C     INCIDENCE TO THE SUBGRAPH SPANNED BY THE ORDERED COLUMNS.
C     AMONG ALL THE COLUMNS OF MAXIMAL INCIDENCE, IDO CHOOSES A
C     COLUMN OF MAXIMAL DEGREE.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE IDO(M,N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST,
C                      MAXCLQ,IWA1,IWA2,IWA3,IWA4)
C
C     WHERE
C
C       M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF ROWS OF A.
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW
C         INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW.
C         THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE
C         COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL.
C         THE COLUMN INDICES FOR ROW I ARE
C
C               INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1.
C
C         NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       NDEG IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE DEGREE SEQUENCE. THE DEGREE OF THE J-TH COLUMN
C         OF A IS NDEG(J).
C
C       LIST IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE INCIDENCE-DEGREE ORDERING OF THE COLUMNS OF A. THE J-TH
C         COLUMN IN THIS ORDER IS LIST(J).
C
C       MAXCLQ IS AN INTEGER OUTPUT VARIABLE SET TO THE SIZE
C         OF THE LARGEST CLIQUE FOUND DURING THE ORDERING.
C
C       IWA1,IWA2,IWA3, AND IWA4 ARE INTEGER WORK ARRAYS OF LENGTH N.
C
C     SUBPROGRAMS CALLED
C
C       MINPACK-SUPPLIED ... NUMSRT
C
C       FORTRAN-SUPPLIED ... MAX
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER IC,IP,IR,JCOL,JP,
     *        MAXINC,MAXLST,NCOMP,NUMINC,NUMLST,NUMORD,NUMWGT
C
C     SORT THE DEGREE SEQUENCE.
C
      CALL NUMSRT(N,N-1,NDEG,-1,IWA4,IWA2,IWA3)
C
C     INITIALIZATION BLOCK.
C
C     CREATE A DOUBLY-LINKED LIST TO ACCESS THE INCIDENCES OF THE
C     COLUMNS. THE POINTERS FOR THE LINKED LIST ARE AS FOLLOWS.
C
C     EACH UN-ORDERED COLUMN IC IS IN A LIST (THE INCIDENCE LIST)
C     OF COLUMNS WITH THE SAME INCIDENCE.
C
C     IWA1(NUMINC) IS THE FIRST COLUMN IN THE NUMINC LIST
C     UNLESS IWA1(NUMINC) = 0. IN THIS CASE THERE ARE
C     NO COLUMNS IN THE NUMINC LIST.
C
C     IWA2(IC) IS THE COLUMN BEFORE IC IN THE INCIDENCE LIST
C     UNLESS IWA2(IC) = 0. IN THIS CASE IC IS THE FIRST
C     COLUMN IN THIS INCIDENCE LIST.
C
C     IWA3(IC) IS THE COLUMN AFTER IC IN THE INCIDENCE LIST
C     UNLESS IWA3(IC) = 0. IN THIS CASE IC IS THE LAST
C     COLUMN IN THIS INCIDENCE LIST.
C
C     IF IC IS AN UN-ORDERED COLUMN, THEN LIST(IC) IS THE
C     INCIDENCE OF IC TO THE GRAPH INDUCED BY THE ORDERED
C     COLUMNS. IF JCOL IS AN ORDERED COLUMN, THEN LIST(JCOL)
C     IS THE INCIDENCE-DEGREE ORDER OF COLUMN JCOL.
C
      MAXINC = 0
      DO 10 JP = N, 1, -1
         IC = IWA4(JP)
         IWA1(N-JP) = 0
         IWA2(IC) = 0
         IWA3(IC) = IWA1(0)
         IF (IWA1(0) .GT. 0) IWA2(IWA1(0)) = IC
         IWA1(0) = IC
         IWA4(JP) = 0
         LIST(JP) = 0
   10    CONTINUE
C
C     DETERMINE THE MAXIMAL SEARCH LENGTH FOR THE LIST
C     OF COLUMNS OF MAXIMAL INCIDENCE.
C
      MAXLST = 0
      DO 20 IR = 1, M
         MAXLST = MAXLST + (IPNTR(IR+1) - IPNTR(IR))**2
   20    CONTINUE
      MAXLST = MAXLST/N
      MAXCLQ = 0
      NUMORD = 1
C
C     BEGINNING OF ITERATION LOOP.
C
   30 CONTINUE
C
C        UPDATE THE SIZE OF THE LARGEST CLIQUE
C        FOUND DURING THE ORDERING.
C
         IF (MAXINC .EQ. 0) NCOMP = 0
         NCOMP = NCOMP + 1
         IF (MAXINC + 1 .EQ. NCOMP) MAXCLQ = MAX(MAXCLQ,NCOMP)
C
C        CHOOSE A COLUMN JCOL OF MAXIMAL DEGREE AMONG THE
C        COLUMNS OF MAXIMAL INCIDENCE MAXINC.
C
   40    CONTINUE
            JP = IWA1(MAXINC)
            IF (JP .GT. 0) GO TO 50
            MAXINC = MAXINC - 1
            GO TO 40
   50    CONTINUE
         NUMWGT = -1
         DO 60 NUMLST = 1, MAXLST
            IF (NDEG(JP) .GT. NUMWGT) THEN
               NUMWGT = NDEG(JP)
               JCOL = JP
               END IF
            JP = IWA3(JP)
            IF (JP .LE. 0) GO TO 70
   60       CONTINUE
   70    CONTINUE
         LIST(JCOL) = NUMORD
         NUMORD = NUMORD + 1
C
C        TERMINATION TEST.
C
         IF (NUMORD .GT. N) GO TO 100
C
C        DELETE COLUMN JCOL FROM THE MAXINC LIST.
C
         IF (IWA2(JCOL) .EQ. 0) THEN
            IWA1(MAXINC) = IWA3(JCOL)
         ELSE
            IWA3(IWA2(JCOL)) = IWA3(JCOL)
            END IF
         IF (IWA3(JCOL) .GT. 0) IWA2(IWA3(JCOL)) = IWA2(JCOL)
C
C        FIND ALL COLUMNS ADJACENT TO COLUMN JCOL.
C
         IWA4(JCOL) = N
C
C        DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND
C        TO NON-ZEROES IN THE MATRIX.
C
         DO 90 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1
            IR = INDROW(JP)
C
C           FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC)
C           WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX.
C
            DO 80 IP = IPNTR(IR), IPNTR(IR+1)-1
               IC = INDCOL(IP)
C
C              ARRAY IWA4 MARKS COLUMNS WHICH ARE ADJACENT TO
C              COLUMN JCOL.
C
               IF (IWA4(IC) .LT. NUMORD) THEN
                  IWA4(IC) = NUMORD
C
C                 UPDATE THE POINTERS TO THE CURRENT INCIDENCE LISTS.
C
                  NUMINC = LIST(IC)
                  LIST(IC) = LIST(IC) + 1
                  MAXINC = MAX(MAXINC,LIST(IC))
C
C                 DELETE COLUMN IC FROM THE NUMINC LIST.
C
                  IF (IWA2(IC) .EQ. 0) THEN
                     IWA1(NUMINC) = IWA3(IC)
                  ELSE
                     IWA3(IWA2(IC)) = IWA3(IC)
                     END IF
                  IF (IWA3(IC) .GT. 0) IWA2(IWA3(IC)) = IWA2(IC)
C
C                 ADD COLUMN IC TO THE NUMINC+1 LIST.
C
                  IWA2(IC) = 0
                  IWA3(IC) = IWA1(NUMINC+1)
                  IF (IWA1(NUMINC+1) .GT. 0) IWA2(IWA1(NUMINC+1)) = IC
                  IWA1(NUMINC+1) = IC
                  END IF
   80          CONTINUE
   90       CONTINUE
C
C        END OF ITERATION LOOP.
C
         GO TO 30
  100 CONTINUE
C
C     INVERT THE ARRAY LIST.
C
      DO 110 JCOL = 1, N
         IWA2(LIST(JCOL)) = JCOL
  110    CONTINUE
      DO 120 JP = 1, N
         LIST(JP) = IWA2(JP)
  120    CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE IDO.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE NUMSRT(N,NMAX,NUM,MODE,INDEX,LAST,NEXT)
      INTEGER N,NMAX,MODE
      INTEGER NUM(N),INDEX(N),LAST(0:NMAX),NEXT(N)
C     **********.
C
C     SUBROUTINE NUMSRT
C
C     GIVEN A SEQUENCE OF INTEGERS, THIS SUBROUTINE GROUPS
C     TOGETHER THOSE INDICES WITH THE SAME SEQUENCE VALUE
C     AND, OPTIONALLY, SORTS THE SEQUENCE INTO EITHER
C     ASCENDING OR DESCENDING ORDER.
C
C     THE SEQUENCE OF INTEGERS IS DEFINED BY THE ARRAY NUM,
C     AND IT IS ASSUMED THAT THE INTEGERS ARE EACH FROM THE SET
C     0,1,...,NMAX. ON OUTPUT THE INDICES K SUCH THAT NUM(K) = L
C     FOR ANY L = 0,1,...,NMAX CAN BE OBTAINED FROM THE ARRAYS
C     LAST AND NEXT AS FOLLOWS.
C
C           K = LAST(L)
C           WHILE (K .NE. 0) K = NEXT(K)
C
C     OPTIONALLY, THE SUBROUTINE PRODUCES AN ARRAY INDEX SO THAT
C     THE SEQUENCE NUM(INDEX(I)), I = 1,2,...,N IS SORTED.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE NUMSRT(N,NMAX,NUM,MODE,INDEX,LAST,NEXT)
C
C     WHERE
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE.
C
C       NMAX IS A POSITIVE INTEGER INPUT VARIABLE.
C
C       NUM IS AN INPUT ARRAY OF LENGTH N WHICH CONTAINS THE
C         SEQUENCE OF INTEGERS TO BE GROUPED AND SORTED. IT
C         IS ASSUMED THAT THE INTEGERS ARE EACH FROM THE SET
C         0,1,...,NMAX.
C
C       MODE IS AN INTEGER INPUT VARIABLE. THE SEQUENCE NUM IS
C         SORTED IN ASCENDING ORDER IF MODE IS POSITIVE AND IN
C         DESCENDING ORDER IF MODE IS NEGATIVE. IF MODE IS 0,
C         NO SORTING IS DONE.
C
C       INDEX IS AN INTEGER OUTPUT ARRAY OF LENGTH N SET SO
C         THAT THE SEQUENCE
C
C               NUM(INDEX(I)), I = 1,2,...,N
C
C         IS SORTED ACCORDING TO THE SETTING OF MODE. IF MODE
C         IS 0, INDEX IS NOT REFERENCED.
C
C       LAST IS AN INTEGER OUTPUT ARRAY OF LENGTH NMAX + 1. THE
C         INDEX OF NUM FOR THE LAST OCCURRENCE OF L IS LAST(L)
C         FOR ANY L = 0,1,...,NMAX UNLESS LAST(L) = 0. IN
C         THIS CASE L DOES NOT APPEAR IN NUM.
C
C       NEXT IS AN INTEGER OUTPUT ARRAY OF LENGTH N. IF
C         NUM(K) = L, THEN THE INDEX OF NUM FOR THE PREVIOUS
C         OCCURRENCE OF L IS NEXT(K) FOR ANY L = 0,1,...,NMAX
C         UNLESS NEXT(K) = 0. IN THIS CASE THERE IS NO PREVIOUS
C         OCCURRENCE OF L IN NUM.
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER I,J,JINC,JL,JU,K,L
C
C     DETERMINE THE ARRAYS NEXT AND LAST.
C
      DO 10 I = 0, NMAX
         LAST(I) = 0
   10    CONTINUE
      DO 20 K = 1, N
         L = NUM(K)
         NEXT(K) = LAST(L)
         LAST(L) = K
   20    CONTINUE
      IF (MODE .EQ. 0) RETURN
C
C     STORE THE POINTERS TO THE SORTED ARRAY IN INDEX.
C
      I = 1
      IF (MODE .GT. 0) THEN
         JL = 0
         JU = NMAX
         JINC = 1
      ELSE
         JL = NMAX
         JU = 0
         JINC = -1
         END IF
      DO 50 J = JL, JU, JINC
         K = LAST(J)
   30    CONTINUE
            IF (K .EQ. 0) GO TO 40
            INDEX(I) = K
            I = I + 1
            K = NEXT(K)
            GO TO 30
   40    CONTINUE
   50    CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE NUMSRT.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,LIST,NGRP,MAXGRP,
     *               IWA)
      INTEGER N,MAXGRP
      INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(*),LIST(N),
     *        NGRP(N),IWA(N)
C     **********
C
C     SUBROUTINE SEQ
C
C     GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, THIS
C     SUBROUTINE DETERMINES A CONSISTENT PARTITION OF THE
C     COLUMNS OF A BY A SEQUENTIAL ALGORITHM.
C
C     A CONSISTENT PARTITION IS DEFINED IN TERMS OF THE LOOPLESS
C     GRAPH G WITH VERTICES A(J), J = 1,2,...,N WHERE A(J) IS THE
C     J-TH COLUMN OF A AND WITH EDGE (A(I),A(J)) IF AND ONLY IF
C     COLUMNS I AND J HAVE A NON-ZERO IN THE SAME ROW POSITION.
C
C     A PARTITION OF THE COLUMNS OF A INTO GROUPS IS CONSISTENT
C     IF THE COLUMNS IN ANY GROUP ARE NOT ADJACENT IN THE GRAPH G.
C     IN GRAPH-THEORY TERMINOLOGY, A CONSISTENT PARTITION OF THE
C     COLUMNS OF A CORRESPONDS TO A COLORING OF THE GRAPH G.
C
C     THE SUBROUTINE EXAMINES THE COLUMNS IN THE ORDER SPECIFIED
C     BY THE ARRAY LIST, AND ASSIGNS THE CURRENT COLUMN TO THE
C     GROUP WITH THE SMALLEST POSSIBLE NUMBER.
C
C     NOTE THAT THE VALUE OF M IS NOT NEEDED BY SEQ AND IS
C     THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE SEQ(N,INDROW,JPNTR,INDCOL,IPNTR,LIST,NGRP,MAXGRP,
C                      IWA)
C
C     WHERE
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW
C         INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW.
C         THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE
C         COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL.
C         THE COLUMN INDICES FOR ROW I ARE
C
C               INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1.
C
C         NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       LIST IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE ORDER TO BE USED BY THE SEQUENTIAL ALGORITHM.
C         THE J-TH COLUMN IN THIS ORDER IS LIST(J).
C
C       NGRP IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE PARTITION OF THE COLUMNS OF A. COLUMN JCOL BELONGS
C         TO GROUP NGRP(JCOL).
C
C       MAXGRP IS AN INTEGER OUTPUT VARIABLE WHICH SPECIFIES THE
C         NUMBER OF GROUPS IN THE PARTITION OF THE COLUMNS OF A.
C
C       IWA IS AN INTEGER WORK ARRAY OF LENGTH N.
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER IC,IP,IR,J,JCOL,JP
C
C     INITIALIZATION BLOCK.
C
      MAXGRP = 0
      DO 10 JP = 1, N
         NGRP(JP) = N
         IWA(JP) = 0
   10    CONTINUE
C
C     BEGINNING OF ITERATION LOOP.
C
      DO 60 J = 1, N
         JCOL = LIST(J)
C
C        FIND ALL COLUMNS ADJACENT TO COLUMN JCOL.
C
C        DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND
C        TO NON-ZEROES IN THE MATRIX.
C
         DO 30 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1
            IR = INDROW(JP)
C
C           FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC)
C           WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX.
C
            DO 20 IP = IPNTR(IR), IPNTR(IR+1)-1
               IC = INDCOL(IP)
C
C              ARRAY IWA MARKS THE GROUP NUMBERS OF THE
C              COLUMNS WHICH ARE ADJACENT TO COLUMN JCOL.
C
               IWA(NGRP(IC)) = J
   20          CONTINUE
   30       CONTINUE
C
C        ASSIGN THE SMALLEST UN-MARKED GROUP NUMBER TO JCOL.
C
         DO 40 JP = 1, MAXGRP
            IF (IWA(JP) .NE. J) GO TO 50
   40       CONTINUE
         MAXGRP = MAXGRP + 1
   50    CONTINUE
         NGRP(JCOL) = JP
   60    CONTINUE
C
C        END OF ITERATION LOOP.
C
      RETURN
C
C     LAST CARD OF SUBROUTINE SEQ.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE SETR(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA)
      INTEGER M,N
      INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(M+1),IWA(M)
C     **********
C
C     SUBROUTINE SETR
C
C     GIVEN A COLUMN-ORIENTED DEFINITION OF THE SPARSITY PATTERN
C     OF AN M BY N MATRIX A, THIS SUBROUTINE DETERMINES A
C     ROW-ORIENTED DEFINITION OF THE SPARSITY PATTERN OF A.
C
C     ON INPUT THE COLUMN-ORIENTED DEFINITION IS SPECIFIED BY
C     THE ARRAYS INDROW AND JPNTR. ON OUTPUT THE ROW-ORIENTED
C     DEFINITION IS SPECIFIED BY THE ARRAYS INDCOL AND IPNTR.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE SETR(M,N,INDROW,JPNTR,INDCOL,IPNTR,IWA)
C
C     WHERE
C
C       M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF ROWS OF A.
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW
C         INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW.
C         THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       INDCOL IS AN INTEGER OUTPUT ARRAY WHICH CONTAINS THE
C         COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       IPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH M + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL.
C         THE COLUMN INDICES FOR ROW I ARE
C
C               INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1.
C
C         NOTE THAT IPNTR(1) IS SET TO 1 AND THAT IPNTR(M+1)-1 IS
C         THEN THE NUMBER OF NON-ZERO ELEMENTS OF THE MATRIX A.
C
C       IWA IS AN INTEGER WORK ARRAY OF LENGTH M.
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER IR,JCOL,JP
C
C     STORE IN ARRAY IWA THE COUNTS OF NON-ZEROES IN THE ROWS.
C
      DO 10 IR = 1, M
         IWA(IR) = 0
   10    CONTINUE
      DO 20 JP = 1, JPNTR(N+1)-1
         IWA(INDROW(JP)) = IWA(INDROW(JP)) + 1
   20    CONTINUE
C
C     SET POINTERS TO THE START OF THE ROWS IN INDCOL.
C
      IPNTR(1) = 1
      DO 30 IR = 1, M
         IPNTR(IR+1) = IPNTR(IR) + IWA(IR)
         IWA(IR) = IPNTR(IR)
   30    CONTINUE
C
C     FILL INDCOL.
C
      DO 50 JCOL = 1, N
         DO 40 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1
            IR = INDROW(JP)
            INDCOL(IWA(IR)) = JCOL
            IWA(IR) = IWA(IR) + 1
   40       CONTINUE
   50    CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE SETR.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE SLO(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST,
     *               MAXCLQ,IWA1,IWA2,IWA3,IWA4)
      INTEGER N,MAXCLQ
      INTEGER INDROW(*),JPNTR(N+1),INDCOL(*),IPNTR(*),NDEG(N),
     *        LIST(N),IWA1(0:N-1),IWA2(N),IWA3(N),IWA4(N)
C     **********
C
C     SUBROUTINE SLO
C
C     GIVEN THE SPARSITY PATTERN OF AN M BY N MATRIX A, THIS
C     SUBROUTINE DETERMINES THE SMALLEST-LAST ORDERING OF THE
C     COLUMNS OF A.
C
C     THE SMALLEST-LAST ORDERING IS DEFINED FOR THE LOOPLESS
C     GRAPH G WITH VERTICES A(J), J = 1,2,...,N WHERE A(J) IS THE
C     J-TH COLUMN OF A AND WITH EDGE (A(I),A(J)) IF AND ONLY IF
C     COLUMNS I AND J HAVE A NON-ZERO IN THE SAME ROW POSITION.
C
C     THE SMALLEST-LAST ORDERING IS DETERMINED RECURSIVELY BY
C     LETTING LIST(K), K = N,...,1 BE A COLUMN WITH LEAST DEGREE
C     IN THE SUBGRAPH SPANNED BY THE UN-ORDERED COLUMNS.
C
C     NOTE THAT THE VALUE OF M IS NOT NEEDED BY SLO AND IS
C     THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE SLO(N,INDROW,JPNTR,INDCOL,IPNTR,NDEG,LIST,
C                      MAXCLQ,IWA1,IWA2,IWA3,IWA4)
C
C     WHERE
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       INDROW IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW
C         INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       JPNTR IS AN INTEGER INPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN INDROW.
C         THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(N+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       INDCOL IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE
C         COLUMN INDICES FOR THE NON-ZEROES IN THE MATRIX A.
C
C       IPNTR IS AN INTEGER INPUT ARRAY OF LENGTH M + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN INDCOL.
C         THE COLUMN INDICES FOR ROW I ARE
C
C               INDCOL(K), K = IPNTR(I),...,IPNTR(I+1)-1.
C
C         NOTE THAT IPNTR(M+1)-1 IS THEN THE NUMBER OF NON-ZERO
C         ELEMENTS OF THE MATRIX A.
C
C       NDEG IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE DEGREE SEQUENCE. THE DEGREE OF THE J-TH COLUMN
C         OF A IS NDEG(J).
C
C       LIST IS AN INTEGER OUTPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE SMALLEST-LAST ORDERING OF THE COLUMNS OF A. THE J-TH
C         COLUMN IN THIS ORDER IS LIST(J).
C
C       MAXCLQ IS AN INTEGER OUTPUT VARIABLE SET TO THE SIZE
C         OF THE LARGEST CLIQUE FOUND DURING THE ORDERING.
C
C       IWA1,IWA2,IWA3, AND IWA4 ARE INTEGER WORK ARRAYS OF LENGTH N.
C
C     SUBPROGRAMS CALLED
C
C       FORTRAN-SUPPLIED ... MIN
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER IC,IP,IR,JCOL,JP,MINDEG,NUMDEG,NUMORD
C
C     INITIALIZATION BLOCK.
C
      MINDEG = N
      DO 10 JP = 1, N
         IWA1(JP-1) = 0
         IWA4(JP) = N
         LIST(JP) = NDEG(JP)
         MINDEG = MIN(MINDEG,NDEG(JP))
   10    CONTINUE
C
C     CREATE A DOUBLY-LINKED LIST TO ACCESS THE DEGREES OF THE
C     COLUMNS. THE POINTERS FOR THE LINKED LIST ARE AS FOLLOWS.
C
C     EACH UN-ORDERED COLUMN IC IS IN A LIST (THE DEGREE LIST)
C     OF COLUMNS WITH THE SAME DEGREE.
C
C     IWA1(NUMDEG) IS THE FIRST COLUMN IN THE NUMDEG LIST
C     UNLESS IWA1(NUMDEG) = 0. IN THIS CASE THERE ARE
C     NO COLUMNS IN THE NUMDEG LIST.
C
C     IWA2(IC) IS THE COLUMN BEFORE IC IN THE DEGREE LIST
C     UNLESS IWA2(IC) = 0. IN THIS CASE IC IS THE FIRST
C     COLUMN IN THIS DEGREE LIST.
C
C     IWA3(IC) IS THE COLUMN AFTER IC IN THE DEGREE LIST
C     UNLESS IWA3(IC) = 0. IN THIS CASE IC IS THE LAST
C     COLUMN IN THIS DEGREE LIST.
C
C     IF IC IS AN UN-ORDERED COLUMN, THEN LIST(IC) IS THE
C     DEGREE OF IC IN THE GRAPH INDUCED BY THE UN-ORDERED
C     COLUMNS. IF JCOL IS AN ORDERED COLUMN, THEN LIST(JCOL)
C     IS THE SMALLEST-LAST ORDER OF COLUMN JCOL.
C
      DO 20 JP = 1, N
         NUMDEG = NDEG(JP)
         IWA2(JP) = 0
         IWA3(JP) = IWA1(NUMDEG)
         IF (IWA1(NUMDEG) .GT. 0) IWA2(IWA1(NUMDEG)) = JP
         IWA1(NUMDEG) = JP
   20    CONTINUE
      MAXCLQ = 0
      NUMORD = N
C
C     BEGINNING OF ITERATION LOOP.
C
   30 CONTINUE
C
C        MARK THE SIZE OF THE LARGEST CLIQUE
C        FOUND DURING THE ORDERING.
C
         IF (MINDEG + 1 .EQ. NUMORD .AND. MAXCLQ .EQ. 0)
     *       MAXCLQ = NUMORD
C
C        CHOOSE A COLUMN JCOL OF MINIMAL DEGREE MINDEG.
C
   40    CONTINUE
            JCOL = IWA1(MINDEG)
            IF (JCOL .GT. 0) GO TO 50
            MINDEG = MINDEG + 1
            GO TO 40
   50    CONTINUE
         LIST(JCOL) = NUMORD
         NUMORD = NUMORD - 1
C
C        TERMINATION TEST.
C
         IF (NUMORD .EQ. 0) GO TO 80
C
C        DELETE COLUMN JCOL FROM THE MINDEG LIST.
C
         IWA1(MINDEG) = IWA3(JCOL)
         IF (IWA3(JCOL) .GT. 0) IWA2(IWA3(JCOL)) = 0
C
C        FIND ALL COLUMNS ADJACENT TO COLUMN JCOL.
C
         IWA4(JCOL) = 0
C
C        DETERMINE ALL POSITIONS (IR,JCOL) WHICH CORRESPOND
C        TO NON-ZEROES IN THE MATRIX.
C
         DO 70 JP = JPNTR(JCOL), JPNTR(JCOL+1)-1
            IR = INDROW(JP)
C
C           FOR EACH ROW IR, DETERMINE ALL POSITIONS (IR,IC)
C           WHICH CORRESPOND TO NON-ZEROES IN THE MATRIX.
C
            DO 60 IP = IPNTR(IR), IPNTR(IR+1)-1
               IC = INDCOL(IP)
C
C              ARRAY IWA4 MARKS COLUMNS WHICH ARE ADJACENT TO
C              COLUMN JCOL.
C
               IF (IWA4(IC) .GT. NUMORD) THEN
                  IWA4(IC) = NUMORD
C
C                 UPDATE THE POINTERS TO THE CURRENT DEGREE LISTS.
C
                  NUMDEG = LIST(IC)
                  LIST(IC) = LIST(IC) - 1
                  MINDEG = MIN(MINDEG,LIST(IC))
C
C                 DELETE COLUMN IC FROM THE NUMDEG LIST.
C
                  IF (IWA2(IC) .EQ. 0) THEN
                     IWA1(NUMDEG) = IWA3(IC)
                  ELSE
                     IWA3(IWA2(IC)) = IWA3(IC)
                     END IF
                  IF (IWA3(IC) .GT. 0) IWA2(IWA3(IC)) = IWA2(IC)
C
C                 ADD COLUMN IC TO THE NUMDEG-1 LIST.
C
                  IWA2(IC) = 0
                  IWA3(IC) = IWA1(NUMDEG-1)
                  IF (IWA1(NUMDEG-1) .GT. 0) IWA2(IWA1(NUMDEG-1)) = IC
                  IWA1(NUMDEG-1) = IC
                  END IF
   60          CONTINUE
   70       CONTINUE
C
C        END OF ITERATION LOOP.
C
         GO TO 30
   80 CONTINUE
C
C     INVERT THE ARRAY LIST.
C
      DO 90 JCOL = 1, N
         IWA2(LIST(JCOL)) = JCOL
   90    CONTINUE
      DO 100 JP = 1, N
         LIST(JP) = IWA2(JP)
  100    CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE SLO.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE SRTDAT(N,NNZ,INDROW,INDCOL,JPNTR,IWA)
      INTEGER N,NNZ
      INTEGER INDROW(NNZ),INDCOL(NNZ),JPNTR(N+1),IWA(N)
C     **********
C
C     SUBROUTINE SRTDAT
C
C     GIVEN THE NON-ZERO ELEMENTS OF AN M BY N MATRIX A IN
C     ARBITRARY ORDER AS SPECIFIED BY THEIR ROW AND COLUMN
C     INDICES, THIS SUBROUTINE PERMUTES THESE ELEMENTS SO
C     THAT THEIR COLUMN INDICES ARE IN NON-DECREASING ORDER.
C
C     ON INPUT IT IS ASSUMED THAT THE ELEMENTS ARE SPECIFIED IN
C
C           INDROW(K),INDCOL(K), K = 1,...,NNZ.
C
C     ON OUTPUT THE ELEMENTS ARE PERMUTED SO THAT INDCOL IS
C     IN NON-DECREASING ORDER. IN ADDITION, THE ARRAY JPNTR
C     IS SET SO THAT THE ROW INDICES FOR COLUMN J ARE
C
C           INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C     NOTE THAT THE VALUE OF M IS NOT NEEDED BY SRTDAT AND IS
C     THEREFORE NOT PRESENT IN THE SUBROUTINE STATEMENT.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE SRTDAT(N,NNZ,INDROW,INDCOL,JPNTR,IWA)
C
C     WHERE
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF A.
C
C       NNZ IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF NON-ZERO ELEMENTS OF A.
C
C       INDROW IS AN INTEGER ARRAY OF LENGTH NNZ. ON INPUT INDROW
C         MUST CONTAIN THE ROW INDICES OF THE NON-ZERO ELEMENTS OF A.
C         ON OUTPUT INDROW IS PERMUTED SO THAT THE CORRESPONDING
C         COLUMN INDICES OF INDCOL ARE IN NON-DECREASING ORDER.
C
C       INDCOL IS AN INTEGER ARRAY OF LENGTH NNZ. ON INPUT INDCOL
C         MUST CONTAIN THE COLUMN INDICES OF THE NON-ZERO ELEMENTS
C         OF A. ON OUTPUT INDCOL IS PERMUTED SO THAT THESE INDICES
C         ARE IN NON-DECREASING ORDER.
C
C       JPNTR IS AN INTEGER OUTPUT ARRAY OF LENGTH N + 1 WHICH
C         SPECIFIES THE LOCATIONS OF THE ROW INDICES IN THE OUTPUT
C         INDROW. THE ROW INDICES FOR COLUMN J ARE
C
C               INDROW(K), K = JPNTR(J),...,JPNTR(J+1)-1.
C
C         NOTE THAT JPNTR(1) IS SET TO 1 AND THAT JPNTR(N+1)-1
C         IS THEN NNZ.
C
C       IWA IS AN INTEGER WORK ARRAY OF LENGTH N.
C
C     SUBPROGRAMS CALLED
C
C       FORTRAN-SUPPLIED ... MAX
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER I,J,K,L
C
C     STORE IN ARRAY IWA THE COUNTS OF NON-ZEROES IN THE COLUMNS.
C
      DO 10 J = 1, N
         IWA(J) = 0
   10    CONTINUE
      DO 20 K = 1, NNZ
         IWA(INDCOL(K)) = IWA(INDCOL(K)) + 1
   20    CONTINUE
C
C     SET POINTERS TO THE START OF THE COLUMNS IN INDROW.
C
      JPNTR(1) = 1
      DO 30 J = 1, N
         JPNTR(J+1) = JPNTR(J) + IWA(J)
         IWA(J) = JPNTR(J)
   30    CONTINUE
      K = 1
C
C     BEGIN IN-PLACE SORT.
C
   40 CONTINUE
         J = INDCOL(K)
         IF (K .GE. JPNTR(J)) THEN
C
C           CURRENT ELEMENT IS IN POSITION. NOW EXAMINE THE
C           NEXT ELEMENT OR THE FIRST UN-SORTED ELEMENT IN
C           THE J-TH GROUP.
C
            K = MAX(K+1,IWA(J))
         ELSE
C
C           CURRENT ELEMENT IS NOT IN POSITION. PLACE ELEMENT
C           IN POSITION AND MAKE THE DISPLACED ELEMENT THE
C           CURRENT ELEMENT.
C
            L = IWA(J)
            IWA(J) = IWA(J) + 1
            I = INDROW(K)
            INDROW(K) = INDROW(L)
            INDCOL(K) = INDCOL(L)
            INDROW(L) = I
            INDCOL(L) = J
            END IF
         IF (K .LE. NNZ) GO TO 40
      RETURN
C
C     LAST CARD OF SUBROUTINE SRTDAT.
C
      END
* ---------------------------------------------------------------------------- *
      SUBROUTINE FDJS(M,N,COL,IND,NPNTR,NGRP,NUMGRP,d,FJACD,FJAC)
      INTEGER M,N,NUMGRP
      INTEGER IND(*),NPNTR(*),NGRP(N)
      REAL D(N),FJACD(M),FJAC(*)
      LOGICAL COL
C     **********
C
C     SUBROUTINE FDJS
C
C     GIVEN A CONSISTENT PARTITION OF THE COLUMNS OF AN M BY N
C     JACOBIAN MATRIX INTO GROUPS, THIS SUBROUTINE COMPUTES
C     APPROXIMATIONS TO THOSE COLUMNS IN A GIVEN GROUP.  THE
C     APPROXIMATIONS ARE STORED INTO EITHER A COLUMN-ORIENTED
C     OR A ROW-ORIENTED PATTERN.
C
C     A PARTITION IS CONSISTENT IF THE COLUMNS IN ANY GROUP
C     DO NOT HAVE A NON-ZERO IN THE SAME ROW POSITION.
C
C     APPROXIMATIONS TO THE COLUMNS OF THE JACOBIAN MATRIX IN A
C     GIVEN GROUP CAN BE OBTAINED BY SPECIFYING A DIFFERENCE
C     PARAMETER ARRAY D WITH D(JCOL) NON-ZERO IF AND ONLY IF
C     JCOL IS A COLUMN IN THE GROUP, AND AN APPROXIMATION TO
C     JAC*D WHERE JAC DENOTES THE JACOBIAN MATRIX OF A MAPPING F.
C
C     D CAN BE DEFINED WITH THE FOLLOWING SEGMENT OF CODE.
C
C           DO 10 JCOL = 1, N
C              D(JCOL) = 0.0
C              IF (NGRP(JCOL) .EQ. NUMGRP) D(JCOL) = ETA(JCOL)
C        10    CONTINUE
C
C     IN THE ABOVE CODE NUMGRP IS THE GIVEN GROUP NUMBER,
C     NGRP(JCOL) IS THE GROUP NUMBER OF COLUMN JCOL, AND
C     ETA(JCOL) IS THE DIFFERENCE PARAMETER USED TO
C     APPROXIMATE COLUMN JCOL OF THE JACOBIAN MATRIX.
C     SUITABLE VALUES FOR THE ARRAY ETA MUST BE PROVIDED.
C
C     AS MENTIONED ABOVE, AN APPROXIMATION TO JAC*D MUST
C     ALSO BE PROVIDED. FOR EXAMPLE, THE APPROXIMATION
C
C           F(X+D) - F(X)
C
C     CORRESPONDS TO THE FORWARD DIFFERENCE FORMULA AT X.
C
C     THE SUBROUTINE STATEMENT IS
C
C       SUBROUTINE FDJS(M,N,COL,IND,NPNTR,NGRP,NUMGRP,D,FJACD,FJAC)
C
C     WHERE
C
C       M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF ROWS OF THE JACOBIAN MATRIX.
C
C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C         OF COLUMNS OF THE JACOBIAN MATRIX.
C
C       COL IS A LOGICAL INPUT VARIABLE. IF COL IS SET TRUE, THEN THE
C         JACOBIAN APPROXIMATIONS ARE STORED INTO A COLUMN-ORIENTED
C         PATTERN. IF COL IS SET FALSE, THEN THE JACOBIAN
C         APPROXIMATIONS ARE STORED INTO A ROW-ORIENTED PATTERN.
C
C       IND IS AN INTEGER INPUT ARRAY WHICH CONTAINS THE ROW
C         INDICES FOR THE NON-ZEROES IN THE JACOBIAN MATRIX
C         IF COL IS TRUE, AND CONTAINS THE COLUMN INDICES FOR
C         THE NON-ZEROES IN THE JACOBIAN MATRIX IF COL IS FALSE.
C
C       NPNTR IS AN INTEGER INPUT ARRAY WHICH SPECIFIES THE
C         LOCATIONS OF THE ROW INDICES IN IND IF COL IS TRUE, AND
C         SPECIFIES THE LOCATIONS OF THE COLUMN INDICES IN IND IF
C         COL IS FALSE. IF COL IS TRUE, THE INDICES FOR COLUMN J ARE
C
C               IND(K), K = NPNTR(J),...,NPNTR(J+1)-1.
C
C         IF COL IS FALSE, THE INDICES FOR ROW I ARE
C
C               IND(K), K = NPNTR(I),...,NPNTR(I+1)-1.
C
C         NOTE THAT NPNTR(N+1)-1 IF COL IS TRUE, OR NPNTR(M+1)-1
C         IF COL IS FALSE, IS THEN THE NUMBER OF NON-ZERO ELEMENTS
C         OF THE JACOBIAN MATRIX.
C
C       NGRP IS AN INTEGER INPUT ARRAY OF LENGTH N WHICH SPECIFIES
C         THE PARTITION OF THE COLUMNS OF THE JACOBIAN MATRIX.
C         COLUMN JCOL BELONGS TO GROUP NGRP(JCOL).
C
C       NUMGRP IS A POSITIVE INTEGER INPUT VARIABLE SET TO A GROUP
C         NUMBER IN THE PARTITION. THE COLUMNS OF THE JACOBIAN
C         MATRIX IN THIS GROUP ARE TO BE ESTIMATED ON THIS CALL.
C
C       D IS AN INPUT ARRAY OF LENGTH N WHICH CONTAINS THE
C         DIFFERENCE PARAMETER VECTOR FOR THE ESTIMATE OF
C         THE JACOBIAN MATRIX COLUMNS IN GROUP NUMGRP.
C
C       FJACD IS AN INPUT ARRAY OF LENGTH M WHICH CONTAINS
C         AN APPROXIMATION TO THE DIFFERENCE VECTOR JAC*D,
C         WHERE JAC DENOTES THE JACOBIAN MATRIX.
C
C       FJAC IS AN OUTPUT ARRAY OF LENGTH NNZ, WHERE NNZ IS THE
C         NUMBER OF ITS NON-ZERO ELEMENTS. AT EACH CALL OF FDJS,
C         FJAC IS UPDATED TO INCLUDE THE NON-ZERO ELEMENTS OF THE
C         JACOBIAN MATRIX FOR THOSE COLUMNS IN GROUP NUMGRP. FJAC
C         SHOULD NOT BE ALTERED BETWEEN SUCCESSIVE CALLS TO FDJS.
C
C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JULY 1983.
C     THOMAS F. COLEMAN, BURTON S. GARBOW, JORGE J. MORE'
C
C     **********
      INTEGER IP,IROW,JCOL,JP
C
C     COMPUTE ESTIMATES OF JACOBIAN MATRIX COLUMNS IN GROUP
C     NUMGRP. THE ARRAY FJACD MUST CONTAIN AN APPROXIMATION
C     TO JAC*D, WHERE JAC DENOTES THE JACOBIAN MATRIX AND D
C     IS A DIFFERENCE PARAMETER VECTOR WITH D(JCOL) NON-ZERO
C     IF AND ONLY IF JCOL IS A COLUMN IN GROUP NUMGRP.
c  M,N,COL,IND,NPNTR,NGRP,NUMGRP,D,FJACD,FJAC
C
      IF (COL) THEN
C
C        COLUMN ORIENTATION.
C
         DO 20 JCOL = 1, N
            IF (NGRP(JCOL) .EQ. NUMGRP) THEN
               DO 10 JP = NPNTR(JCOL), NPNTR(JCOL+1)-1
                  IROW = IND(JP)
                  FJAC(JP) = FJACD(IROW)/D(JCOL)
   10             CONTINUE
               END IF
   20       CONTINUE
      ELSE
C
C        ROW ORIENTATION.
C
         DO 50 IROW = 1, M
            DO 30 IP = NPNTR(IROW), NPNTR(IROW+1)-1
               JCOL = IND(IP)
               IF (NGRP(JCOL) .EQ. NUMGRP) THEN
c	          write (81,'(3i5, es16.8)') irow, ip, jcol, d(jcol)
                  FJAC(IP) = FJACD(IROW)/D(JCOL)
                  GO TO 40
                  END IF
   30          CONTINUE
   40       CONTINUE
   50       CONTINUE
         END IF
      RETURN
C
C     LAST CARD OF SUBROUTINE FDJS.
C
      END
* ============================================================================ *

max_stars_repo_name: Sbte/i-emic

max_issues_repo_name: Sbte/i-emic

max_forks_repo_name: nlesc-smcm/i-emic


text:
! -*- mode: F90; mode: font-lock; column-number-mode: true; vc-back-end: CVS -*-
! ------------------------------------------------------------------------------
! $Id$
! ------------------------------------------------------------------------------
! Module memory_module
! ------------------------------------------------------------------------------

!!****h* Conquest/memory_module *
!!  NAME
!!   memory_module
!!  PURPOSE
!!   Store details of memory allocated by areas of code
!!  AUTHOR
!!   D. R. Bowler
!!  CREATION DATE
!!   2006/09/27
!!  MODIFICATION HISTORY
!!   2008/02/06 08:24 dave
!!    Changed for output to file not stdout
!!   2011/12/01 L.Tong
!!    Added type_cplx, and no_bytes associated to it
!!   2014/01/17 lat
!!    Added optional what and lun arguments in:
!!      reg_alloc_mem
!!      reg_dealloc_mem
!!      write_mem_use
!!  SOURCE
!!
module memory_module

  use datatypes
  use global_module, only: io_lun

  implicit none
  save

  integer, allocatable, dimension(:) :: max_alloc_area, tot_alloc_area
  integer :: tot_alloc, max_alloc

  integer, parameter, dimension(3) :: no_bytes = (/4,8,16/)
  integer, parameter :: type_int  = 1
  integer, parameter :: type_dbl  = 2
  integer, parameter :: type_cplx = 3

  interface reg_dealloc_mem
     module procedure reg_dealloc_mem_32
     module procedure reg_dealloc_mem_64
  end interface

  ! RCS tag for object file identification
  character(len=80), private :: &
       RCSid = "$Id$"
!!***

contains

  ! -----------------------------------------------------------
  ! Subroutine reg_alloc_mem
  ! -----------------------------------------------------------

  !!****f* memory_module/reg_alloc_mem *
  !!
  !!  NAME
  !!   reg_alloc_mem
  !!  USAGE
  !!
  !!  PURPOSE
  !!   Registers the allocation of memory in a particular code area
  !!  INPUTS
  !!
  !!  USES
  !!
  !!  AUTHOR
  !!   D. R. Bowler
  !!  CREATION DATE
  !!   2006/09/27
  !!  MODIFICATION HISTORY
  !!   2014/01/17 lat
  !!    Added optional printing options
  !!  SOURCE
  !!
  subroutine reg_alloc_mem(area, amount, type, what, lun)

    use units,         only: m_units, mem_units, mem_conv
    use global_module, only: iprint_gen
    use GenComms,      only: myid

    implicit none

    ! Passed variables
    integer, intent(in) :: area
    integer, intent(in) :: amount
    integer, intent(in) :: type
    character(*), optional, intent(in) :: what
    integer,      optional, intent(in) :: lun

    if (iprint_gen>4 .and. present(what) .and. present(lun)) then
       write(lun,fmt='(10x,"Allocating ",a12," in area ",i3,f10.3," ",a2)') &
            adjustl(what), area, real(amount*no_bytes(type),double)*mem_conv, &
            mem_units(m_units)
    end if

    if (iprint_gen>4.AND.myid==0) then
       if (present(what)) then
          write(io_lun,fmt='(10x,"Allocating ",a12," in area ",i3,f10.3," ",a2)') &
               adjustl(what), area, real(amount*no_bytes(type),double)*mem_conv, &
               mem_units(m_units)
       else
          write(io_lun,fmt='(10x,"Allocating in area ",i3,f10.3," ",a2)') &
               area, real(amount*no_bytes(type),double)*mem_conv,mem_units(m_units)
       end if
    end if
    tot_alloc_area(area) = tot_alloc_area(area) + amount*no_bytes(type)
    max_alloc_area(area) = max(max_alloc_area(area),tot_alloc_area(area))
    tot_alloc = tot_alloc + amount*no_bytes(type)
    max_alloc = max(max_alloc, tot_alloc)

  end subroutine reg_alloc_mem
  !!***

  ! -----------------------------------------------------------
  ! Subroutine reg_dealloc_mem
  ! -----------------------------------------------------------

  !!****f* memory_module/reg_dealloc_mem *
  !!
  !!  NAME
  !!   reg_dealloc_mem
  !!  USAGE
  !!
  !!  PURPOSE
  !!   Registers the deallocation of memory in a particular code area
  !!  INPUTS
  !!
  !!  USES
  !!
  !!  AUTHOR
  !!   D. R. Bowler
  !!  CREATION DATE
  !!   2006/09/27
  !!  MODIFICATION HISTORY
  !!   2014/01/17 lat
  !!    Added optional printing options
  !!  SOURCE
  !!
  subroutine reg_dealloc_mem_32(area, amount, type, what, lun)

    use units, only: m_units, mem_units, mem_conv
    use global_module, only: iprint_gen

    implicit none

    ! Passed variables
    integer :: area
    integer :: amount
    integer :: type

    character(*), optional, intent(in) :: what
    integer,      optional, intent(in) :: lun

    if (iprint_gen>4 .and. present(what).and.present(lun)) then
       write(lun,fmt='(10x,"Delocating ",a12," in area ",i3,f10.3," ",a2)') &
            adjustl(what), area, real(amount*no_bytes(type),double)*mem_conv, &
            mem_units(m_units)
    end if

    tot_alloc_area(area) = tot_alloc_area(area) - amount*no_bytes(type)
    tot_alloc = tot_alloc - amount*no_bytes(type)

  end subroutine reg_dealloc_mem_32

  subroutine reg_dealloc_mem_64(area, amount, type, what, lun)

    use units, only: m_units, mem_units, mem_conv
    use global_module, only: iprint_gen

    implicit none

    ! Passed variables
    integer   :: area
    integer(8):: amount
    integer   :: type

    character(*), optional, intent(in) :: what
    integer,      optional, intent(in) :: lun

    if (iprint_gen>4 .and. present(what).and.present(lun)) then
       write(lun,fmt='(10x,"Delocating ",a12," in area ",i3,f10.3," ",a2)') &
            adjustl(what), area, real(amount*no_bytes(type),double)*mem_conv, &
            mem_units(m_units)
    end if

    tot_alloc_area(area) = tot_alloc_area(area) - amount*no_bytes(type)
    tot_alloc = tot_alloc - amount*no_bytes(type)

  end subroutine reg_dealloc_mem_64
  !!***

  ! -----------------------------------------------------------
  ! Subroutine init_reg_mem
  ! -----------------------------------------------------------

  !!****f* memory_module/init_reg_mem *
  !!
  !!  NAME
  !!   init_reg_mem
  !!  USAGE
  !!
  !!  PURPOSE
  !!   Initialises variables used in tracking memory allocation
  !!  INPUTS
  !!
  !!  USES
  !!
  !!  AUTHOR
  !!   D. R. Bowler
  !!  CREATION DATE
  !!   2006/09/27
  !!  MODIFICATION HISTORY
  !!   2014/01/17 lat
  !!    Added optional printing options
  !!  SOURCE
  !!
  subroutine init_reg_mem

    use global_module, only: n_areas
    use GenComms,      only: cq_abort

    implicit none

    ! Local variables
    integer :: stat

    allocate(max_alloc_area(n_areas), tot_alloc_area(n_areas), STAT=stat)
    if (stat /= 0) call cq_abort("Error allocating max_alloc_area: ", n_areas)
    max_alloc_area = 0
    tot_alloc_area = 0
    max_alloc = 0
    tot_alloc = 0

    return
  end subroutine init_reg_mem
  !!***
  subroutine write_mem_use(lun, area)

    use datatypes
    use global_module, only: n_areas
    use units,         only: m_units, mem_units, mem_conv
    use GenComms,      only: inode, ionode

    implicit none

    integer :: i
    integer, optional, intent(in) :: lun
    integer, optional, intent(in) :: area

    if (present(lun).and.present(area)) then
       write(lun,'(2x,"Max mem use for area ",i4," is ",f10.3," ",a2)') &
            area,real(max_alloc_area(area))*mem_conv,mem_units(m_units)
       write(io_lun,'(2x,"Max mem use for area ",i4," is ",f10.3," ",a2)') &
            area,real(max_alloc_area(area))*mem_conv,mem_units(m_units)
    else
       if(inode==ionode) then
          do i=1,n_areas
             write(io_lun,'(2x,"Max mem use for area ",i4," is ",f10.3," ",a2)') &
                  i,real(max_alloc_area(i))*mem_conv,mem_units(m_units)
          end do
          write(io_lun,'(2x,"Max total mem use is ",f10.3," ",a2)') &
               real(max_alloc)*mem_conv,mem_units(m_units)
       end if
    end if

    return
  end subroutine write_mem_use

end module memory_module

max_stars_repo_name: negf/CONQUEST-release

max_issues_repo_name: negf/CONQUEST-release

max_forks_repo_name: negf/CONQUEST-release


text:
cdis    Forecast Systems Laboratory
cdis    NOAA/OAR/ERL/FSL
cdis    325 Broadway
cdis    Boulder, CO     80303
cdis
cdis    Forecast Research Division
cdis    Local Analysis and Prediction Branch
cdis    LAPS
cdis
cdis    This software and its documentation are in the public domain and
cdis    are furnished "as is."  The United States government, its
cdis    instrumentalities, officers, employees, and agents make no
cdis    warranty, express or implied, as to the usefulness of the software
cdis    and documentation for any purpose.  They assume no responsibility
cdis    (1) for the use of the software and documentation; or (2) to provide
cdis     technical support to users.
cdis
cdis    Permission to use, copy, modify, and distribute this software is
cdis    hereby granted, provided that the entire disclaimer notice appears
cdis    in all copies.  All modifications to this software must be clearly
cdis    documented, and are solely the responsibility of the agent making
cdis    the modifications.  If significant modifications or enhancements
cdis    are made to this software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis
cdis
cdis
cdis
cdis
cdis
      FUNCTION ATAN3D(X,Y)

!       Steve Albers            1990

      real atan3d

      ATAN3D=MOD((ATAN2_DEG(X,Y)+360E0),360E0)

      RETURN

      END

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  Copyright by The HDF Group.                                              *
!  Copyright by the Board of Trustees of the University of Illinois.        *
!  All rights reserved.                                                     *
!                                                                           *
!  This file is part of HDF5.  The full HDF5 copyright notice, including    *
!  terms governing use, modification, and redistribution, is contained in   *
!   the COPYING file, which can be found at the root of the source code       *
!   distribution tree, or in https://www.hdfgroup.org/licenses.               *
!   If you do not have access to either file, you may request a copy from     *
!   help@hdfgroup.org.                                                        *
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
! This example shows how to write and read a hyperslab.
! It is used in the HDF5 Tutorial.
!

PROGRAM H5_SUBSET

  USE HDF5 ! This module contains all necessary modules

  IMPLICIT NONE

  CHARACTER(LEN=9), PARAMETER :: filename = "subset.h5"  ! File name
  CHARACTER(LEN=8), PARAMETER :: dsetname = "IntArray"   ! Dataset name

  INTEGER(HID_T) :: file_id       ! File identifier
  INTEGER(HID_T) :: dset_id       ! Dataset identifier
  INTEGER(HID_T) :: dataspace     ! Dataspace identifier
  INTEGER(HID_T) :: memspace      ! memspace identifier

  !
  ! To change the subset size, modify size of dimsm, sdata, dim0_sub,
  ! dim1_sub, and count
  !
  INTEGER(HSIZE_T), DIMENSION(1:2) :: dimsm = (/4,3/) ! Dataset dimensions
  INTEGER, DIMENSION(1:4,1:3) :: sdata                  ! Subset buffer
  INTEGER :: dim0_sub = 4
  INTEGER :: dim1_sub = 3
  INTEGER(HSIZE_T), DIMENSION(1:2) :: count = (/4,3/)  ! Size of hyperslab
  INTEGER(HSIZE_T), DIMENSION(1:2) :: offset = (/2,1/) ! Hyperslab offset
  INTEGER(HSIZE_T), DIMENSION(1:2) :: stride = (/1,1/) ! Hyperslab stride
  INTEGER(HSIZE_T), DIMENSION(1:2) :: block = (/1,1/)  ! Hyperslab block size

  INTEGER(HSIZE_T), DIMENSION(1:2) :: dimsf = (/10,8/) ! Dataset dimensions

  INTEGER, DIMENSION(1:10,1:8) :: data     ! Data to write
  INTEGER, DIMENSION(1:10,1:8) :: rdata    ! Data to read

  INTEGER :: rank = 2      ! Dataset rank ( in file )
  INTEGER :: dim0 = 10     ! Dataset size in file
  INTEGER :: dim1 = 8

  INTEGER :: i, j

  INTEGER :: error         ! Error flag
  INTEGER(HSIZE_T), DIMENSION(2) :: data_dims

  !
  ! Write data to the HDF5 file.
  !

  !
  ! Data initialization.
  !
  DO i = 1, dim0
     DO j = 1, dim1
        IF (i .LE. (dim0 / 2)) THEN
           data(i,j) = 1
        ELSE
           data(i,j) = 2
        END IF
     END DO
  END DO

  !
  ! Initialize FORTRAN interface.
  !
  CALL h5open_f(error)

  !
  ! Create a new file using default properties.
  !
  CALL h5fcreate_f(filename, H5F_ACC_TRUNC_F, file_id, error)

  !
  ! Create the data space for the  dataset.
  !
  CALL h5screate_simple_f(rank, dimsf, dataspace, error)

  !
  ! Create the dataset with default properties.
  !
  CALL h5dcreate_f(file_id, dsetname, H5T_NATIVE_INTEGER, dataspace, &
       dset_id, error)

  !
  ! Write the dataset.
  !
  data_dims(1) = dim0
  data_dims(2) = dim1
  CALL h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, data, data_dims, error)

  !
  ! Data Written to File
  !
  WRITE(*,'(/,A)') "Original Data Written to File:"
  DO i = 1, dim0
     WRITE(*,'(100(1X,I0,1X))') DATA(i,1:dim1)
  END DO

  !
  !
  ! Close the dataspace, dataset, and file.
  !
  CALL h5sclose_f(dataspace, error)
  CALL h5dclose_f(dset_id, error)
  CALL h5fclose_f(file_id, error)

  !
  ! Initialize subset data array.
  !
  sdata(1:dim0_sub,1:dim1_sub) = 5

  !
  ! Open the file.
  !
  CALL h5fopen_f(filename, H5F_ACC_RDWR_F, file_id, error)

  !
  ! Open the  dataset.
  !
  CALL h5dopen_f(file_id, dsetname, dset_id, error)

  !
  ! Get dataset's dataspace identifier and select subset.
  !
  CALL h5dget_space_f(dset_id, dataspace, error)
  CALL h5sselect_hyperslab_f(dataspace, H5S_SELECT_SET_F, &
       offset, count, error, stride, BLOCK)
  !
  ! Create memory dataspace.
  !
  CALL h5screate_simple_f(rank, dimsm, memspace, error)

  WRITE(*,'(/,A)') "Write subset to file specifying:"
  WRITE(*,'(A,/)') "   offset=2x1 stride=1x1 count=4x3 block=1x1"

  !
  ! Write subset to dataset
  !
  data_dims(1:2) = (/dim0_sub, dim1_sub/)
  CALL h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, sdata, data_dims, error, &
       memspace, dataspace)

  data_dims(1:2) = (/dim0, dim1/)
  CALL h5dread_f(dset_id, H5T_NATIVE_INTEGER, rdata, data_dims, error)

  !
  ! Read entire dataset back
  !
  WRITE(*,'(A)') "Data in File after Subset Written:"
  DO i = 1, dim0
    WRITE(*,'(100(1X,I0,1X))') rdata(i,1:dim1)
  END DO
  PRINT *, " "

  !
  ! Close everything opened.
  !
  CALL h5sclose_f(dataspace, error)
  CALL h5sclose_f(memspace, error)
  CALL h5dclose_f(dset_id, error)
  CALL h5fclose_f(file_id, error)

  !
  ! Close FORTRAN interface.
  !
  CALL h5close_f(error)

END PROGRAM H5_SUBSET


        module  OneTstep
        use global_variables
        use domain
        use InAndOut
        use netcdf
        use pulse
        use OneZstep
        use fields_var
        use particles_var
        use TestPart_var
        use output_var
        implicit none
        private
        save
        public::Get_Chi,AdvancePulse,one_step_defile,one_step_output,getlaserfield

        contains

            subroutine Get_Chi(Acom2d,outflag,time)
            implicit none
            real(RP):: ztime,delta,Chi(1:Nr+1),time
            complex(RP) :: Acom2d(1:Nr,1:Nz+1),Ac(1:Nr)
            integer :: outflag
            integer :: nzc,nzc0,f_cross

            call resetplasma(ztime,Chi,nzc,nzc0,f_cross)
            Chi2d(1:Nr,1)=Chi(1:Nr)
            Ac(:)=Acom2d(:,1);

            if(outflag==0) then
            Err2dold=Err2d; Ezz2dold=Ezz2d;Bthh2dold=Bthh2d
            end if

            do while (nzc<Nz)
                  call one_step_z(ztime,Ac,Chi,nzc,f_cross,time)!Advance Trajectory
                  if (nzc>nzc0) then
                   nzc0=nzc
                   Ac(:)=Acom2d(:,nzc+1)
                   Chi2d(1:Nr,nzc+1)=Chi(1:Nr)
                       if(outflag==1)  then
                       r2d(:,nzc+1)=r;V2d(:,nzc+1)=V;force2d(:,nzc+1)=force;gama2d(:,nzc+1)=gama
                       psi2d(:,nzc+1)=Psi(1:Nr);Er2d(:,nzc+1)=Er(1:Nr);Et2d(:,nzc+1)=Et(1:Nr);Bth2d(:,nzc+1)=Bth(1:Nr)
                       den2d(:,nzc+1)=den(1:Nr);den_e2d(:,nzc+1)=den_e(1:Nr);jr2d(:,nzc+1)=jr(1:Nr);jz2d(:,nzc+1)=jz(1:Nr)
                       else
                       Err2d(1:Nr,nzc+1)=Er(1:Nr);Ezz2d(1:Nr,nzc+1)=Et(1:Nr)
                       Bthh2d(1:Nr,nzc+1)=Bth(1:Nr)
                       end if
                  else
                   delta=(ztime-nzc*dz0)/dz0
                   Ac(:)=Acom2d(:,nzc+1)*delta+Acom2d(:,nzc)*(1-delta);
                  endif
            enddo
            end subroutine Get_Chi

            subroutine getlaserfield(Ac2d_0,Ac2d_hdt,Ac2d_dt,dt,dz,k0)
            use TestPart_var
            implicit none
            integer::i
            real(RP) :: dt,dz,k0
            complex(RP) :: Ac2d_0(1:Nr,1:Nz+1),Ac2d_hdt(1:Nr,1:Nz+1),Ac2d_dt(1:Nr,1:Nz+1)
            EL2dold=EL2d;BL2dold=BL2d;ELZ2dold=ELZ2d;

            Ac2d_hdt=(Ac2d_0+Ac2d_dt)*0.5d0

            BL2d(1:Nr,1:Nz)=ci*k0*(Ac2d_hdt(:,1:Nz)+Ac2d_hdt(:,2:Nz+1))*0.5-(Ac2d_hdt(:,2:Nz+1)-Ac2d_hdt(:,1:Nz))/dz
            EL2d(1:Nr,1:Nz)=BL2d(1:Nr,1:Nz)-0.5*(Ac2d_dt(:,2:Nz+1)+Ac2d_dt(:,1:Nz)-Ac2d_0(:,2:Nz+1)-Ac2d_0(:,1:Nz))/dt

            ELZ2d(2:Nr,1:Nz+1)=(Ac2d_hdt(2:Nr,:)-Ac2d_hdt(1:Nr-1,:))/hr

            end subroutine getlaserfield

       		subroutine AdvancePulse(Acom2d,dt)
            use fields_var
        	implicit none
        	real(RP):: dd0,dt
        	complex(RP) :: dd1(1:Nr),dd2(1:Nr),dd3(1:Nr)
    		complex(RP) :: Aold(1:Nr),Anew(1:Nr),Acom2d(1:Nr,1:Nz+1)
    		integer :: k
            dd0=1.5d0;dd1=0d0;dd2=0d0
            !Chi2d=1d0
                do k=1,Nz
                 Aold(:)=Acom2d(:,k)
                 call One_Step_Pulse(dd0,dd1,dd2,Chi2d(:,k),Aold,Anew,dz0,dt)
                 dd2(:)=dd1(:)
                 dd1(:)=Anew(:)-Aold(:)
                 Acom2d(:,k)=Anew(:)
                enddo

        	end subroutine AdvancePulse

            subroutine one_step_defile(ntc)
            implicit none
            integer :: ntc
            character(12):: filename

            write(filename,'(I5)'), ntc
            call check(nf90_create("output/"//TRIM(ADJUSTL(filename))//"Fields.nc", NF90_CLOBBER, ncidf) )
           ! call check(nf90_create("output/"//TRIM(ADJUSTL(filename))//"Particle.nc", NF90_CLOBBER, ncidp) )

            call check( nf90_def_dim(ncidf, "ztime", Nz, ztime_dim) )
            call check( nf90_def_dim(ncidf, "Rgird", Nr, r_dim) )
            !call check( nf90_def_dim(ncidp, "ztime", Nz, ztime_dim) )
            !call check( nf90_def_dim(ncidp, "Part",  Np, p_dim) )
            dims=(/r_dim,ztime_dim/)
            call check( nf90_def_var(ncidf, "Den_e", NF90_REAL, dims, var_de) )
            !call check( nf90_def_var(ncidf, "Jjj_r", NF90_REAL, dims, var_jr) )
            !call check( nf90_def_var(ncidf, "Den_r", NF90_REAL, dims, var_dr) )
            call check( nf90_def_var(ncidf, "Bth_a", NF90_REAL, dims, var_bt) )
            call check( nf90_def_var(ncidf, "Psi_i", NF90_REAL, dims, var_ps) )
            !call check( nf90_def_var(ncidf, "Jjj_z", NF90_REAL, dims, var_jz) )
            call check( nf90_def_var(ncidf, "Eee_z", NF90_REAL, dims, var_ez) )
            call check( nf90_def_var(ncidf, "Eee_r", NF90_REAL, dims, var_er) )
            call check( nf90_def_var(ncidf, "Aaa_R", NF90_REAL, dims, var_AR) )
            call check( nf90_def_var(ncidf, "Aaa_I", NF90_REAL, dims, var_AI) )
            call check( nf90_def_var(ncidf, "Byl_R", NF90_REAL, dims, var_BR) )
            call check( nf90_def_var(ncidf, "Byl_I", NF90_REAL, dims, var_BI) )
            call check( nf90_def_var(ncidf, "Exl_R", NF90_REAL, dims, var_ELR) )
            call check( nf90_def_var(ncidf, "Exl_I", NF90_REAL, dims, var_ELI) )
            call check( nf90_def_var(ncidf, "Ezl_R", NF90_REAL, dims, var_EZLR) )
            call check( nf90_def_var(ncidf, "Ezl_I", NF90_REAL, dims, var_EZLI) )

            call check(nf90_put_att(ncidf,var_de,"Grid","Half") )
            !call check(nf90_put_att(ncidf,var_jr,"Grid","Int") )
            !call check(nf90_put_att(ncidf,var_dr,"Grid","Half") )
            call check(nf90_put_att(ncidf,var_bt,"Grid","Int") )
            call check(nf90_put_att(ncidf,var_ps,"Grid","Half") )
            !call check(nf90_put_att(ncidf,var_jz,"Grid","Half") )
            call check(nf90_put_att(ncidf,var_ez,"Grid","Half") )
            call check(nf90_put_att(ncidf,var_er,"Grid","Int") )
            call check(nf90_put_att(ncidf,var_AR,"Grid","Half") )
            call check(nf90_put_att(ncidf,var_AI,"Grid","Half") )
            call check(nf90_put_att(ncidf,var_BR,"Grid","half-z,half-r") )
            call check(nf90_put_att(ncidf,var_BI,"Grid","half-z,half-r") )
            call check(nf90_put_att(ncidf,var_ELR,"Grid","half-z,half-r") )
            call check(nf90_put_att(ncidf,var_ELI,"Grid","half-z,half-r") )
            call check(nf90_put_att(ncidf,var_EZLR,"Grid","full-z,full-r") )
            call check(nf90_put_att(ncidf,var_EZLI,"Grid","full-z,full-r") )

            !dims=(/p_dim,ztime_dim/)
            !call check( nf90_def_var(ncidp, "Vvv_r", NF90_REAL, dims, var_vr) )
            !call check( nf90_def_var(ncidp, "Rrr_r", NF90_REAL, dims, var_rr) )
            !call check( nf90_def_var(ncidp, "Fff_r", NF90_REAL, dims, var_fr) )
           ! call check( nf90_def_var(ncidp, "Gam_a", NF90_REAL, dims, var_ga) )

            call check( nf90_enddef(ncidf) )
            !call check( nf90_enddef(ncidp) )

            end subroutine one_step_defile

        subroutine one_step_output
        implicit none
        call check(nf90_put_var(ncidf,var_de,den_e2d(1:Nr,1:Nz)))
        !call check(nf90_put_var(ncidf,var_jr,jr2d(1:Nr,1:Nz)))
        !call check(nf90_put_var(ncidf,var_dr,den2d(1:Nr,1:Nz)))
        call check(nf90_put_var(ncidf,var_bt,Bth2d(1:Nr,1:Nz)))
        call check(nf90_put_var(ncidf,var_ps,Psi2d(1:Nr,1:Nz)))
        !call check(nf90_put_var(ncidf,var_jz,jz2d(1:Nr,1:Nz)))
        call check(nf90_put_var(ncidf,var_ez,Et2d(1:Nr,1:Nz)))
        call check(nf90_put_var(ncidf,var_er,Er2d(1:Nr,1:Nz)))
        call check(nf90_put_var(ncidf,var_AR,real(Ac2d(1:Nr,1:Nz))*2d0))
        call check(nf90_put_var(ncidf,var_AI,AImag(Ac2d(1:Nr,1:Nz))*2d0))
        call check(nf90_put_var(ncidf,var_BR,real(BL2d(1:Nr,1:Nz))))
        call check(nf90_put_var(ncidf,var_BI,AImag(BL2d(1:Nr,1:Nz))))
        call check(nf90_put_var(ncidf,var_ELR,real(EL2d(1:Nr,1:Nz))))
        call check(nf90_put_var(ncidf,var_ELI,AImag(EL2d(1:Nr,1:Nz))))
        call check(nf90_put_var(ncidf,var_EZLR,real(ELZ2d(1:Nr,1:Nz))))
        call check(nf90_put_var(ncidf,var_EZLI,AImag(ELZ2d(1:Nr,1:Nz))))
        !call check(nf90_put_var(ncidp,var_vr,V2d(1:Np,1:Nz)))
        !call check(nf90_put_var(ncidp,var_rr,r2d(1:Np,1:Nz)))
        !call check(nf90_put_var(ncidp,var_fr,force2d(1:Np,1:Nz)))
        !call check(nf90_put_var(ncidp,var_ga,gama2d(1:Np,1:Nz)))
        call check( nf90_close(ncidf))
        !call check( nf90_close(ncidp))
        end subroutine one_step_output

        end   module  OneTstep


! { dg-do run }
! PR 50192 - on little-endian systems, this used to fail.
program main
  character(kind=4,len=2) :: c1, c2
  c1 = 4_' '
  c2 = 4_' '
  c1(1:1) = transfer(257, mold=c1(1:1))
  c2(1:1) = transfer(64, mold=c2(1:1))
  if (c1 < c2) STOP 1
end program main