MODULE constants_NSWC
! Contains the NSWC functions SPMPAR, DPMPAR, EPSLN, DEPSLN, EXPARG & DXPARG
!-----------------------------------------------------------------------
!     WRITTEN using F90 intrinsics by
!        Alan Miller
!        CSIRO Mathematical & Information Sciences
!        CLAYTON, VICTORIA, AUSTRALIA 3169
!     Latest revision - 1 February 1997
!-----------------------------------------------------------------------

IMPLICIT NONE
INTEGER, PARAMETER     :: dp = SELECTED_REAL_KIND(15, 60)

CONTAINS

FUNCTION spmpar (i) RESULT(fn_val)
!-----------------------------------------------------------------------

!     SPMPAR PROVIDES THE SINGLE PRECISION MACHINE CONSTANTS FOR
!     THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT
!     I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE
!     SINGLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND
!     ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN

!        SPMPAR(1) = B**(1 - M), THE MACHINE PRECISION,

!        SPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE,

!        SPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE.
!-----------------------------------------------------------------------

IMPLICIT NONE
INTEGER, INTENT(IN) :: i
REAL                :: fn_val

! Local variable
REAL                :: one = 1.0

SELECT CASE (i)
  CASE (1)
    fn_val = EPSILON(one)
  CASE (2)
    fn_val = TINY(one)
  CASE (3)
    fn_val = HUGE(one)
END SELECT

RETURN
END FUNCTION spmpar



FUNCTION dpmpar (i) RESULT(fn_val)
!-----------------------------------------------------------------------

!     DPMPAR PROVIDES THE DOUBLE PRECISION MACHINE CONSTANTS FOR
!     THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT
!     I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE
!     DOUBLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND
!     ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN

!        DPMPAR(1) = B**(1 - M), THE MACHINE PRECISION,

!        DPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE,

!        DPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE.
!-----------------------------------------------------------------------

IMPLICIT NONE
INTEGER, INTENT(IN) :: i
REAL (dp)           :: fn_val

! Local variable
REAL (dp)    :: one = 1._dp

SELECT CASE (i)
  CASE (1)
    fn_val = EPSILON(one)
  CASE (2)
    fn_val = TINY(one)
  CASE (3)
    fn_val = HUGE(one)
END SELECT

RETURN
END FUNCTION dpmpar


FUNCTION epsln () RESULT(fn_val)
!--------------------------------------------------------------------
!     THE EVALUATION OF LN(EPS) WHERE EPS IS THE SMALLEST NUMBER
!     SUCH THAT 1.0 + EPS .GT. 1.0 .  L IS A DUMMY ARGUMENT.
!--------------------------------------------------------------------
IMPLICIT NONE
REAL                :: fn_val

! Local variable
REAL                :: one = 1.0

fn_val = LOG( EPSILON(one) )
RETURN
END FUNCTION epsln


FUNCTION exparg (l) RESULT(fn_val)
!--------------------------------------------------------------------
!     IF L = 0 THEN  EXPARG(L) = THE LARGEST POSITIVE W FOR WHICH
!     EXP(W) CAN BE COMPUTED.
!
!     IF L IS NONZERO THEN  EXPARG(L) = THE LARGEST NEGATIVE W FOR
!     WHICH THE COMPUTED VALUE OF EXP(W) IS NONZERO.
!
!     NOTE... ONLY AN APPROXIMATE VALUE FOR EXPARG(L) IS NEEDED.
!--------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: l
REAL                :: fn_val

! Local variable
REAL                :: one = 1.0

IF (l == 0) THEN
  fn_val = LOG( HUGE(one) )
ELSE
  fn_val = LOG( TINY(one) )
END IF
RETURN
END FUNCTION exparg


FUNCTION depsln () RESULT(fn_val)
!--------------------------------------------------------------------
!     THE EVALUATION OF LN(EPS) WHERE EPS IS THE SMALLEST NUMBER
!     SUCH THAT 1.D0 + EPS .GT. 1.D0 .  L IS A DUMMY ARGUMENT.
!--------------------------------------------------------------------
IMPLICIT NONE
REAL (dp)           :: fn_val

! Local variable
REAL (dp)    :: one = 1._dp

fn_val = LOG( EPSILON(one) )
RETURN
END FUNCTION depsln


FUNCTION dxparg (l) RESULT(fn_val)
!--------------------------------------------------------------------
!     IF L = 0 THEN  DXPARG(L) = THE LARGEST POSITIVE W FOR WHICH
!     DEXP(W) CAN BE COMPUTED.
!
!     IF L IS NONZERO THEN  DXPARG(L) = THE LARGEST NEGATIVE W FOR
!     WHICH THE COMPUTED VALUE OF DEXP(W) IS NONZERO.
!
!     NOTE... ONLY AN APPROXIMATE VALUE FOR DXPARG(L) IS NEEDED.
!--------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: l
REAL (dp)           :: fn_val

! Local variable
REAL (dp)    :: one = 1._dp

IF (l == 0) THEN
  fn_val = LOG( HUGE(one) )
ELSE
  fn_val = LOG( TINY(one) )
END IF
RETURN
END FUNCTION dxparg

END MODULE constants_NSWC



SUBROUTINE smplx (a, b0, c, ka, m, n0, ind, ibasis, x, z, iter, mxiter,   &
                  numle, numge, bi, rerr)
!-----------------------------------------------------------------------
!     SIMPLEX PROCEDURE FOR SOLVING LINEAR PROGRAMMING PROBLEMS
!-----------------------------------------------------------------------
! Finds non-negative x's to maximize:
!     c1.x1 + ... + cn.xn
!
! Subject to the constraints:
!     a11.x1 + ... + a1n.xn {<=,=,>=} b1
!              ...
!     am1.x1 + ... + amn.xn {<=,=,>=} bm
!
! The constraints must be ordered so that the <= constraints, if any, come
! first, then any >= constraints, with the equality constraints last.
!
! Arguments:
! a(ka,n0)   INPUT   coefficients in the constraints.
! b0(m)      INPUT   right-hand side of the constraints.
! c(n0)      INPUT   vector of `costs' in the objective function.
! ka         INPUT   first dimension of array a.
! m          INPUT   dimension of array b0.
! n0         INPUT   2nd dimension of array a.
! ind        IN/OUT  If IND = 0 on input, the routine selects its own initial
!                    basis, else if IND = 1 then the indices of the initial
!                    basis should be in IBASIS.
!                    On output:
!                    IND = 0 the problem was solved
!                    IND = 1 the problem has no solution
!                    IND = 2 MXITER iterations were performed; more needed
!                    IND = 3 sufficient accuracy could not be maintained to
!                            solve the problem
!                    IND = 4 the problem has an unbounded solution
!                    IND = 5 input error detected
!                    IND = 6 the solution may have been obtained
! ibasis(m)  IN/OUT  indices of the variables in the basis.
! x()        OUTPUT  Dimension must be >= n + numle + numge.
!                    If IND = 0 or 6, it contains the values of the original,
!                    slack and surplus variables.
! z          OUTPUT  If IND = 0 or 6, contains the value of the objective.
! iter       OUTPUT  number of iterations used.
! mxiter     INPUT   maximum number of iterations.
! numle      INPUT   number of <= constraints.
! numge      INPUT   number of >= constraints.
! bi(m,m)    OUTPUT  the inverse of the basis matrix.
! rerr       OUTPUT  the estimated relative error achieved.
! N.B. The last two arguments in the NSWC routine have been eliminated.
!      These were work spaces.   They are now declared as automatic arrays.
!-----------------------------------------------------------------------
!     WRITTEN BY ALFRED H. MORRIS JR.
!        NAVAL SURFACE WEAPONS CENTER
!        DAHLGREN, VIRGINIA
!------------------------
!     INITIAL VERSION  DEC  1977
!     LAST UPDATE      OCT  1990
!-----------------------------------------------------------------------
!     Converted using F90 intrinsics by
!        Alan Miller
!        CSIRO Mathematical & Information Sciences
!        CLAYTON, VICTORIA, AUSTRALIA 3169
!     Latest revision - 5 February 1997
!-----------------------------------------------------------------------

USE constants_NSWC
IMPLICIT NONE

INTEGER, INTENT(IN)                    :: ka, m, n0, mxiter, numle, numge
INTEGER, INTENT(IN OUT)                :: ind
INTEGER, DIMENSION(:), INTENT(IN OUT)  :: ibasis
INTEGER, INTENT(OUT)                   :: iter
REAL (dp), DIMENSION(:,:), INTENT(IN)  :: a
REAL (dp), DIMENSION(:), INTENT(IN)    :: b0, c
REAL (dp), INTENT(OUT)                 :: z, rerr
REAL (dp), DIMENSION(:), INTENT(OUT)   :: x
REAL (dp), DIMENSION(:,:), INTENT(OUT) :: bi

INTERFACE
  SUBROUTINE smplx1 (a, b0, c, ka, m, n0, ind, ibasis, r, z, iter, mxiter,  &
                     eps0, rerrmn, rerrmx, rerr, numle, numge, bi)
    USE constants_NSWC
    IMPLICIT NONE
    INTEGER, INTENT(IN)                    :: ka, m, n0, mxiter, numle, numge
    INTEGER, INTENT(OUT)                   :: ind, iter
    REAL (dp), DIMENSION(:,:), INTENT(IN)  :: a
    REAL (dp), DIMENSION(:), INTENT(IN)    :: b0, c
    REAL (dp), DIMENSION(:,:), INTENT(OUT) :: bi
    REAL (dp), DIMENSION(:), INTENT(OUT)   :: r
    REAL (dp), INTENT(OUT)                 :: z, rerr
    REAL (dp), INTENT(IN)                  :: eps0, rerrmn, rerrmx
    INTEGER, DIMENSION(:), INTENT(IN OUT)  :: ibasis
  END SUBROUTINE smplx1
END INTERFACE

!------------------------
!     DIMENSION X(N0+NUMLE+NUMGE)
!------------------------

!     ********** EPS0 IS A MACHINE DEPENDENT PARAMETER. ASSIGN EPS0
!                THE VALUE U WHERE U IS THE SMALLEST POSITIVE FLOATING
!                POINT NUMBER SUCH THAT 1.0 + U .GT. 1.0.

!     Local variables
REAL (dp) :: eps0, rerrmn, rerrmx

eps0 = dpmpar(1)

!------------------------
rerrmn = 10.0_dp*eps0
rerrmx = 1.d-4
IF (eps0 < 1.d-13) rerrmx = 1.d-5

CALL smplx1(a, b0, c, ka, m, n0, ind, ibasis, x, z, iter, mxiter, eps0,   &
        rerrmn, rerrmx, rerr, numle, numge, bi)
RETURN
END SUBROUTINE smplx



SUBROUTINE smplx1 (a, b0, c, ka, m, n0, ind, ibasis, r, z, iter, mxiter,  &
                   eps0, rerrmn, rerrmx, rerr, numle, numge, bi)
!----------------------
!     NSTEP = 1   ELIMINATE THE NEGATIVE VARIABLES
!     NSTEP = 2   PHASE 1 OF THE SIMPLEX ALGORITHM
!     NSTEP = 3   PHASE 2 OF THE SIMPLEX ALGORITHM
!----------------------
!     MXITER = THE MAXIMUM NUMBER OF ITERATIONS PERMITTED
!     ITER = THE NUMBER OF THE CURRENT ITERATION
!     ICOUNT = THE NUMBER OF ITERATIONS SINCE THE LAST INVERSION
!----------------------
!     NUMLE = THE NUMBER OF .LE. CONSTRAINTS
!     NUMGE = THE NUMBER OF .GE. CONSTRAINTS
!----------------------
!     THE ROUTINE ASSUMES THAT THE .LE. CONSTRAINTS PRECEDE THE .GE.
!     CONSTRAINTS AND THAT THE .EQ. CONSTRAINTS COME LAST. THERE ARE
!     M CONSTRAINTS. X(N0+I) IS THE SLACK, SURPLUS, OR ARTIFICIAL
!     VARIABLE FOR THE I-TH CONSTRAINT (I=1, ..., M).
!----------------------
!     N0 = THE NUMBER OF ORGINAL VARIABLES
!     NS = THE NUMBER OF ORGINAL AND SLACK VARIABLES
!     N  = THE NUMBER OF ORGINAL, SLACK, AND SURPLUS VARIABLES
!     NUM = THE TOTAL NUMBER OF VARIABLES
!----------------------
!     RERRMN = THE SMALLEST RELATIVE ERROR TOLERANCE USED
!     RERRMX = THE LARGEST RELATIVE ERROR TOLERACE USED
!     RERR   = THE ESTIMATED CURRENT RELATIVE ERROR
!----------------------
!     ASSUME THAT
!         B0 = (B0(1), ..., B0(M))
!         C  = (C(1), ..., C(N0))
!         Z  = C(1)*X(1)+...+C(N0)*X(N0)
!     THE PROBLEM IS TO MAXIMIZE Z SUBJECT TO
!         AX(LE,EQ,GE)B0
!         X.GE.0
!----------------------
!     ON INPUT IND CAN HAVE THE VALUES
!         IND = 0   NO BEGINNING BASIS IS PROVIDED BY THE USER
!         IND = 1   THE ARRAY IBASIS HAS BEEN SET BY THE USER
!     ON OUTPUT IND IS ASSIGNED ONE OF THE VALUES
!         IND = 0   Z WAS SUCCESSFULLY MAXIMIZED
!         IND = 1   THE PROBLEM HAS NO FEASIBLE SOLUTION
!         IND = 2   MXITER ITERATIONS WERE PERFORMED
!         IND = 3   SUFFICIENT ACCURACY CANNOT BE MAINTAINED
!         IND = 4   THE PROBLEM HAS AN UNBOUNDED SOLUTION
!         IND = 5   THERE IS AN INPUT ERROR
!         IND = 6   Z WAS POSSIBLY MAXIMIZED
!----------------------
!     BASIS IS AN INTEGER ARRAY OF DIMENSION N0+M. FOR J.LE.N
!         BASIS(J) = 1  IF X(J) IS A BASIC VARIABLE
!         BASIS(J) = 0  IF X(J) IS NOT A BASIC VARIABLE
!     IF THE BASIC VARIABLES ARE X(I1), ..., X(IM) THEN
!         IBASIS = (I1, ..., IM)
!     ALSO XB(1), ..., XB(M) ARE THE CORRESPONDING VALUES OF THE
!     BASIC VARIABLES.
!----------------------
!     BI IS AN MXM ARRAY CONTAINING THE INVERSE OF THE BASIS MATRIX.
!----------------------
!     R IS AN ARRAY OF DIMENSION N. ON OUTPUT R CONTAINS THE CURRENT
!     VALUE OF X. DURING COMPUTATION R NORMALLY CONTAINS THE REDUCED
!     COSTS USED FOR THE SELECTION OF THE VARIABLE TO BE MADE BASIC.
!----------------------
USE constants_NSWC
IMPLICIT NONE
INTEGER, INTENT(IN)                    :: ka, m, n0, mxiter, numle, numge
INTEGER, INTENT(OUT)                   :: ind, iter
REAL (dp), DIMENSION(:,:), INTENT(IN)  :: a
REAL (dp), DIMENSION(:), INTENT(IN)    :: b0, c
REAL (dp), DIMENSION(:,:), INTENT(OUT) :: bi
REAL (dp), DIMENSION(:), INTENT(OUT)   :: r
REAL (dp), INTENT(OUT)                 :: z, rerr
REAL (dp), INTENT(IN)                  :: eps0, rerrmn, rerrmx
INTEGER, DIMENSION(:), INTENT(IN OUT)  :: ibasis
!----------------------

INTERFACE
  SUBROUTINE crout1(a, ka, n, iend, indx, temp, ierr)
    USE constants_NSWC
    IMPLICIT NONE
    INTEGER, INTENT(IN)                     :: ka, n, iend
    INTEGER, INTENT(OUT)                    :: ierr
    REAL (dp), DIMENSION(:), INTENT(IN OUT) :: a, temp
    INTEGER, DIMENSION(:), INTENT(IN OUT)   :: indx
  END SUBROUTINE crout1
END INTERFACE

!     Local variables
INTEGER   :: i, ibeg, icount, iend, ierr, ii, il, imin, iout, ip, j, jj,  &
             jmin, jp, k, ki, kj, l, ll, lrow, m0, mcheck, ms, n, npos,   &
             nrow, ns, nstep, num, bflag
REAL (dp) :: xb(m), y(m)
INTEGER   :: basis(m+n0), indx(m)
REAL (dp) :: amax, binorm, bmax, bmin, bnorm, cmin, const, eps, epsi, ratio, &
             rerr1, rmin, rtol, s, sgn, t, tol, total, w, xmax, zero = 0._dp, &
             dsum, dsump, dsumn, dt

!     ****** XMAX IS A MACHINE DEPENDENT CONSTANT. XMAX IS THE
!            LARGEST POSITIVE FLOATING POINT NUMBER.

xmax = dpmpar(3)

!----------------------
iter = 0
icount = 0
mcheck = MIN(5, 1 + m/15)
z = zero

!                CHECK FOR INPUT ERRORS

ms = numle + numge
IF (m < 2 .OR. n0 < 2 .OR. ms > m .OR. ka < m)  &
GO TO 12
DO i = 1, m
  IF (b0(i) < zero) GO TO 12
  xb(i) = zero
END DO
rtol = xmax
DO i = 1, n0
  IF (c(i) /= zero) rtol = MIN(ABS(c(i)), rtol)
END DO
rtol = rerrmx*rtol
GO TO 20

12 ind = 5
RETURN

!     FORMATION OF THE IBASIS AND BASIS ARRAYS.
!     (IF IND = 1  THEN THE IBASIS ARRAY IS DEFINED BY THE USER.)

20 ns = n0 + numle
n = ns + numge
IF (ind == 0) GO TO 30
num = n
DO i = 1, m
  IF (ibasis(i) > n) num = num + 1
END DO
GO TO 32
22 IF (ind == 0) GO TO 590
ind = 0

30 num = n0 + m
DO i = 1, m
  ibasis(i) = n0 + i
END DO
32 bflag = 0
basis(1:n) = 0
DO i = 1, m
  ki = ibasis(i)
  basis(ki) = 1
END DO
IF (ind == 1) GO TO 100

!          CALCULATION OF XB AND BI WHEN IND = 0

rerr = rerrmn
DO j = 1, m
  xb(j) = b0(j)
  bi(1:m,j) = zero
  bi(j,j) = 1.0_dp
END DO
IF (numge == 0) GO TO 630
jmin = numle + 1
DO j = jmin, ms
  xb(j) = -xb(j)
  bi(j,j) = -1.0_dp
END DO
GO TO 601

!                  REORDER THE BASIS

100 ibeg = 1
iend = m
DO i = 1, m
  IF (ibasis(i) > n0) THEN
    indx(ibeg) = ibasis(i)
    ibeg = ibeg + 1
  ELSE
    indx(iend) = ibasis(i)
    iend = iend - 1
  END IF
END DO
IF (iend == m) GO TO 22
ibasis(1:m) = indx(1:m)

!            REINVERSION OF THE BASIS MATRIX

DO j = 1, m
  kj = ibasis(j)
  IF (kj <= n0) GO TO 110
  IF (kj <= ns) GO TO 120
  IF (kj <= n) GO TO 130
  GO TO 120
  
  110 bi(1:m,j) = a(1:m,kj)
  CYCLE
  
  120 l = kj - n0
  bi(1:m,j) = zero
  bi(l,j) = 1.0_dp
  CYCLE
  
  130 l = kj - n0
  bi(1:m,j) = zero
  bi(l,j) = -1.0_dp
END DO

icount = 0
CALL crout1 (bi(1:,1), m, m, iend, indx, y, ierr)
IF (ierr /= 0) GO TO 580

!         CHECK THE ACCURACY OF BI AND RESET RERR

bnorm = zero
DO j = 1, m
  kj = ibasis(j)
  IF (kj <= n0) GO TO 140
  total = 1.0_dp
  GO TO 142
  140 total = SUM( ABS(a(1:m,kj)) )
  142 bnorm = MAX(bnorm, total)
END DO

binorm = zero
DO j = 1, m
  total = SUM( ABS(bi(1:m,j)) )
  binorm = MAX(binorm, total)
END DO
rerr = MAX(rerrmn, eps0*bnorm*binorm)
IF (rerr > 1.d-2) GO TO 580
bflag = 0

!                 RECALCULATION OF XB

DO i = 1, m
  dsump = zero
  dsumn = zero
  DO l = 1, m
    dt = bi(i,l)*b0(l)
    IF (dt > zero) THEN
      dsump = dsump + dt
    ELSE
      dsumn = dsumn + dt
    END IF
  END DO
  xb(i) = dsump + dsumn
  s = dsump
  t = dsumn
  tol = rerrmx*MAX(s, -t)
  IF (ABS(xb(i)) < tol) xb(i) = zero
END DO
GO TO 601

!     FIND THE NEXT VECTOR A(--, JP) TO BE INSERTED INTO THE BASIS

200 jp = 0
rmin = zero
IF (nstep == 3) rmin = -rtol
DO j = 1, n0
  IF (basis(j) /= 0) CYCLE
  IF (r(j) >= rmin) CYCLE
  jp = j
  rmin = r(j)
END DO
IF (n0 == n) GO TO 203
jmin = n0 + 1
rmin = rmin*1.1_dp
DO j = jmin, n
  IF (basis(j) /= 0) CYCLE
  IF (r(j) >= rmin) CYCLE
  jp = j
  rmin = r(j)
END DO
203 IF (jp /= 0) GO TO 300
IF (nstep < 2) THEN
   GO TO 800
ELSE IF (nstep == 2) THEN
  GO TO 230
ELSE
  GO TO 250
END IF

!     INSERT THE VALUES OF THE ORGINAL, SLACK, AND SURPLUS
!             VARIABLES INTO R, THEN TERMINATE.

220 r(1:n) = zero
DO i = 1, m
  ki = ibasis(i)
  IF (ki <= n) r(ki) = xb(i)
END DO
RETURN

!             COMPLETION OF THE NSTEP = 2 CASE

230 DO i = 1, m
  IF (ibasis(i) <= n) CYCLE
  IF (xb(i) > zero) GO TO 800
END DO
GO TO 680

240 IF (icount >= 5) GO TO 100
ind = 1
GO TO 220

!             COMPLETION OF THE NSTEP = 3 CASE

250 IF (rerr > 1.d-2) GO TO 251
ind = 0
GO TO 800
251 IF (icount >= 5) GO TO 100
ind = 6
GO TO 800

!     IF MXITER ITERATIONS HAVE NOT BEEN PERFORMED THEN BEGIN THE
!     NEXT ITERATION. COMPUTE THE JP-TH COLUMN OF BI*A AND STORE IT IN Y.

300 IF (iter < mxiter) GO TO 301
ind = 2
GO TO 220
301 iter = iter + 1
icount = icount + 1
IF (jp > ns) GO TO 330
IF (jp > n0) GO TO 320

nrow = 0
amax = zero
DO i = 1, m
  IF (a(i,jp) == zero) CYCLE
  nrow = nrow + 1
  indx(nrow) = i
  amax = MAX(ABS(a(i,jp)), amax)
END DO
IF (nrow /= 0) GO TO 310
ind = 4
GO TO 220

310 rerr1 = rerrmx*amax
DO i = 1, m
  dsum = zero
  DO ll = 1, nrow
    l = indx(ll)
    dsum = dsum + bi(i,l)*a(l,jp)
  END DO
  y(i) = dsum
  IF (ABS(y(i)) >= 5.d-3) CYCLE
  bmax = zero
  DO l = 1, m
    bmax = MAX(ABS(bi(i,l)), bmax)
  END DO
  tol = rerr1*bmax
  IF (ABS(y(i)) < tol) y(i) = zero
END DO
GO TO 350

320 l = jp - n0
y(1:m) = bi(1:m,l)
GO TO 350

330 l = jp - n0
y(1:m) = -bi(1:m,l)

350 DO i = 1, m
  IF (y(i) /= zero) GO TO 360
END DO
r(jp) = zero
iter = iter - 1
icount = icount - 1
GO TO 200

360 IF (nstep == 2) GO TO 430
IF (nstep > 2) GO TO 440

!     FINDING THE VARIABLE XB(IP) TO BE MADE NONBASIC FOR THE NSTEP = 1 CASE

npos = 0
ip = 0
eps = zero
epsi = xmax
DO i = 1, m
  IF (xb(i) < zero .OR. y(i) <= zero) CYCLE
  ratio = xb(i)/y(i)
  IF (ratio < epsi) THEN
    epsi = ratio
    npos = 1
    indx(1) = i
    CYCLE
  ELSE IF (ratio > epsi) THEN
    CYCLE
  END IF
  npos = npos + 1
  indx(npos) = i
END DO
IF (npos == 0) GO TO 420
IF (epsi == zero) GO TO 460

DO i = 1, m
  IF (xb(i) >= zero .OR. y(i) >= zero) CYCLE
  ratio = xb(i)/y(i)
  IF (ratio > epsi) CYCLE
  IF (ratio < eps) CYCLE
  eps = ratio
  ip = i
END DO
IF (ip /= 0) GO TO 500
GO TO 460

420 DO i = 1, m
  IF (xb(i) >= zero .OR. y(i) >= zero) CYCLE
  ratio = xb(i)/y(i)
  IF (ratio < eps) CYCLE
  eps = ratio
  ip = i
END DO
GO TO 500

!     FINDING THE VARIABLE XB(IP) TO BE MADE NONBASIC FOR THE NSTEP = 2 CASE

430 npos = 0
epsi = xmax
DO i = 1, m
  IF (y(i) <= zero) CYCLE
  ratio = xb(i)/y(i)
  IF (ratio < epsi) THEN
    epsi = ratio
    npos = 1
    indx(1) = i
  ELSE IF (ratio > epsi) THEN
    CYCLE
  END IF
  npos = npos + 1
  indx(npos) = i
END DO
GO TO 450

!     FINDING THE VARIABLE XB(IP) TO BE MADE NONBASIC FOR THE NSTEP = 3 CASE

440 npos = 0
epsi = xmax
DO i = 1, m
  IF (y(i) < zero) THEN
    IF (ibasis(i) <= n) CYCLE
    ip = i
    GO TO 500
  ELSE IF (y(i) > zero) THEN
    ratio = xb(i)/y(i)
    IF (ratio < epsi) THEN
      epsi = ratio
      npos = 1
      indx(1) = i
    ELSE IF (ratio > epsi) THEN
      CYCLE
    END IF
    npos = npos + 1
    indx(npos) = i
  END IF
END DO

450 IF (npos /= 0) GO TO 460
IF (icount >= 5) GO TO 100
ind = 4
GO TO 220

!              TIE BREAKING PROCEDURE

460 ip = indx(1)
IF (npos == 1) GO TO 500
ip = 0
bmin = xmax
cmin = xmax
DO ii = 1, npos
  i = indx(ii)
  l = ibasis(i)
  IF (l > n0) GO TO 461
  IF (c(l) <= zero) cmin = MIN(zero, cmin)
  IF (c(l) > cmin) CYCLE
  imin = i
  cmin = c(l)
  CYCLE
  461 IF (l <= n) GO TO 462
  ip = i
  GO TO 500
  462 lrow = l - n0
  s = b0(lrow)
  IF (lrow <= numle) THEN
    IF (s > bmin) CYCLE
    ip = i
    bmin = s
  ELSE
    s = -s
    bmin = MIN(zero, bmin)
    IF (s > bmin) CYCLE
    ip = i
    bmin = s
  END IF
END DO
IF (cmin <= zero .OR. ip == 0) ip = imin

!               TRANSFORMATION OF XB

500 IF (xb(ip) == zero) GO TO 510
const = xb(ip)/y(ip)
DO i = 1, m
  s = xb(i)
  xb(i) = xb(i) - const*y(i)
  IF (xb(i) >= zero) CYCLE
  IF (s >= zero .OR. xb(i) >= rerrmx*s) xb(i) = zero
END DO
xb(ip) = const

!               TRANSFORMATION OF BI

510 DO j = 1, m
  IF (bi(ip,j) == zero) CYCLE
  const = bi(ip,j)/y(ip)
  bi(1:m,j) = bi(1:m,j) - const*y(1:m)
  bi(ip,j) = const
END DO

!             UPDATING IBASIS AND BASIS

iout = ibasis(ip)
ibasis(ip) = jp
basis(iout) = 0
basis(jp) = 1
IF (iout > n) num = num - 1

!        CHECK THE ACCURACY OF BI AND RESET RERR

IF (rerr > 1.d-2) GO TO 530
k = 0
DO j = 1, m
  kj = ibasis(j)
  IF (kj > n0) CYCLE
  total = DOT_PRODUCT( bi(j,1:m), a(1:m, kj) )
  rerr = MAX(rerr, ABS(1.0_dp - total))
  k = k + 1
  IF (k >= mcheck) GO TO 522
END DO
522 IF (rerr <= 1.d-2) GO TO 600

!        THE ACCURACY CRITERIA ARE NOT SATISFIED

530 IF (icount < 5) GO TO 600
bflag = 1
GO TO 100

580 IF (iter == 0) GO TO 12
IF (bflag == 0) GO TO 590
bflag = 0
DO ip = 1, m
  IF (jp == ibasis(ip)) GO TO 582
END DO
582 ibasis(ip) = iout
basis(jp) = 0
basis(iout) = 1
IF (iout > n) num = num + 1
GO TO 100

590 ind = 3
GO TO 220

!       SET UP THE R ARRAY FOR THE NSTEP = 1 CASE

600 IF (nstep == 2) GO TO 630
IF (nstep > 2) GO TO 700
601 DO j = 1, m
  IF (xb(j) < zero) GO TO 610
END DO
GO TO 630

610 nstep = 1
m0 = 0
DO l = 1, m
  IF (xb(l) >= zero) CYCLE
  m0 = m0 + 1
  indx(m0) = l
END DO

DO j = 1, m
  dsump = zero
  dsumn = zero
  DO ll = 1, m0
    l = indx(ll)
    IF (bi(l,j) < zero) THEN
      dsumn = dsumn + bi(l,j)
    ELSE IF (bi(l,j) > zero) THEN
      dsump = dsump + bi(l,j)
    END IF
  END DO
  y(j) = dsump + dsumn
  s = dsump
  t = dsumn
  tol = rerrmx*MAX(s, -t)
  IF (ABS(y(j)) < tol) y(j) = zero
END DO
GO TO 650

!       SET UP THE R ARRAY FOR THE NSTEP = 2 CASE

630 IF (n == num) GO TO 680
nstep = 2
m0 = 0
DO l = 1, m
  IF (ibasis(l) <= n) CYCLE
  m0 = m0 + 1
  indx(m0) = l
END DO

DO j = 1, m
  dsump = zero
  dsumn = zero
  DO ll = 1, m0
    l = indx(ll)
    IF (bi(l,j) < zero) THEN
      dsumn = dsumn + bi(l,j)
    ELSE IF (bi(l,j) > zero) THEN
      dsump = dsump + bi(l,j)
    END IF
  END DO
  y(j) = -(dsump + dsumn)
  s = dsump
  t = dsumn
  tol = rerrmx*MAX(s, -t)
  IF (ABS(y(j)) < tol) y(j) = zero
END DO

650 DO j = 1, n0
  IF (basis(j) == 0) THEN
    r(j) = DOT_PRODUCT( y(1:m), a(1:m,j) )
  ELSE
    r(j) = zero
  END IF
END DO

660 IF (n0 == ns) GO TO 670
jmin = n0 + 1
DO j = jmin, ns
  r(j) = zero
  IF (basis(j) /= 0) CYCLE
  jj = j - n0
  r(j) = y(jj)
END DO

670 IF (ns == n) GO TO 200
jmin = ns + 1
DO j = jmin, n
  r(j) = zero
  IF (basis(j) /= 0) CYCLE
  jj = j - n0
  r(j) = -y(jj)
END DO
GO TO 200

!      SET UP A NEW R ARRAY FOR THE NSTEP = 3 CASE

680 nstep = 3
DO j = 1, m
  dsum = zero
  DO l = 1, m
    il = ibasis(l)
    IF (il <= n0) dsum = dsum + c(il)*bi(l,j)
  END DO
  y(j) = dsum
END DO

DO j = 1, n0
  r(j) = zero
  IF (basis(j) /= 0) CYCLE
  dsum = -c(j) + DOT_PRODUCT( y(1:m), a(1:m,j) )
  r(j) = dsum
  IF (r(j) >= zero) CYCLE
  tol = rerrmx*ABS(c(j))
  IF (ABS(r(j)) < tol) r(j) = zero
END DO
GO TO 660

!       UPDATE THE R ARRAY FOR THE NSTEP = 3 CASE

700 const = r(jp)
DO j = 1, n0
  IF (basis(j) /= 0) THEN
    r(j) = zero
  ELSE
    total = DOT_PRODUCT( bi(ip,1:m), a(1:m,j) )
    r(j) = r(j) - const*total
    IF (r(j) >= zero) CYCLE
    tol = rerrmx*ABS(c(j))
    IF (ABS(r(j)) < tol) r(j) = zero
  END IF
END DO

IF (n0 == ns) GO TO 720
jmin = n0 + 1
DO j = jmin, ns
  IF (basis(j) /= 0) THEN
    r(j) = zero
  ELSE
    jj = j - n0
    r(j) = r(j) - const*bi(ip,jj)
  END IF
END DO

720 IF (ns == n) GO TO 200
jmin = ns + 1
DO j = jmin, n
  IF (basis(j) /= 0) THEN
    r(j) = zero
  ELSE
    jj = j - n0
    r(j) = r(j) + const*bi(ip,jj)
  END IF
END DO
GO TO 200
!-----------------------------------------------------------------------
!               REFINE XB AND STORE THE RESULT IN Y
!-----------------------------------------------------------------------
800 y(1:m) = zero

m0 = 0
DO j = 1, m
  kj = ibasis(j)
  IF (kj <= n0) GO TO 810
  IF (kj <= ns) GO TO 820
  IF (kj <= n) GO TO 830
  GO TO 820
  
  810 m0 = m0 + 1
  indx(m0) = j
  CYCLE
  
  820 l = kj - n0
  y(l) = xb(j)
  CYCLE
  
  830 l = kj - n0
  y(l) = -xb(j)
END DO

IF (m0 == 0) THEN
  r(1:m) = b0(1:m) - y(1:m)
ELSE
  DO i = 1, m
    dsum = y(i)
    DO jj = 1, m0
      j = indx(jj)
      kj = ibasis(j)
      dsum = dsum + a(i,kj)*xb(j)
    END DO
    r(i) = b0(i) - dsum
  END DO
END IF

rerr1 = MIN(rerrmx, rerr)
DO i = 1, m
  y(i) = zero
  IF (xb(i) < zero) THEN
    sgn = -1.0_dp
    dsump = zero
    dsumn = xb(i)
  ELSE IF (xb(i) > zero) THEN
    sgn = 1.0_dp
    dsump = xb(i)
    dsumn = zero
  ELSE
    CYCLE
  END IF
  DO l = 1, m
    dt = bi(i,l)*r(l)
    IF (dt > zero) THEN
      dsump = dsump + dt
    ELSE
      dsumn = dsumn + dt
    END IF
  END DO
  w = dsump + dsumn
  IF (w == zero) CYCLE
  IF (sgn /= SIGN(1.0_dp, w)) CYCLE
  s = dsump
  t = dsumn
  tol = rerr1*MAX(s, -t)
  IF (ABS(w) > tol) y(i) = w
END DO
IF (nstep == 2) GO TO 870
IF (nstep > 2) GO TO 880

!         CHECK THE REFINEMENT (NSTEP = 1)

DO i = 1, m
  IF (y(i) >= zero) GO TO 861
  IF (y(i) < -rerrmx) GO TO 240
  y(i) = zero
  861 xb(i) = y(i)
END DO
GO TO 630

!         CHECK THE REFINEMENT (NSTEP = 2)

870 DO i = 1, m
  IF (ibasis(i) <= n) GO TO 871
  IF (y(i) > rerrmx) GO TO 240
  y(i) = zero
  871 xb(i) = y(i)
END DO
GO TO 680

!              COMPUTE Z  (NSTEP = 3)

880 dsum = zero
DO i = 1, m
  ki = ibasis(i)
  IF (ki > n0) GO TO 881
  dsum = dsum + c(ki)*y(i)
  881 xb(i) = y(i)
END DO
z = dsum
GO TO 220
END SUBROUTINE smplx1



SUBROUTINE crout1(a, ka, n, iend, indx, temp, ierr)
!     ******************************************************************
!     CROUT PROCEDURE FOR INVERTING MATRICES
!     ******************************************************************
!     A IS A MATRIX OF ORDER N WHERE N IS GREATER THAN OR EQUAL TO 1.
!     THE INVERSE OF A IS COMPUTED AND STORED IN A.

!     KA = LENGTH OF THE COLUMNS OF THE ARRAY A

!     IEND MAY BE 0, 1, ..., N-1.  IT IS ASSUMED THAT EACH OF THE FIRST
!     IEND COLUMNS OF THE MATRIX A CONTAINS ONLY ONE NONZERO ELEMENT
!     AND THAT THE NONZERO ELEMENT IS 1 OR -1.

!     indx IS AN ARRAY OF DIMENSION N-1 OR LARGER THAT IS USED BY THE
!     ROUTINE FOR KEEPING TRACK OF THE ROW INTERCHANGES THAT ARE MADE.

!     TEMP IS A TEMPORARY STORAGE ARRAY.

!     IERR REPORTS THE STATUS OF THE RESULTS. IF A IS NONSINGULAR THEN
!     THE INVERSE OF A IS COMPUTED AND IERR=0. OTHERWISE IF A IS FOUND
!     TO BE SINGULAR THEN IERR=1 AND THE ROUTINE TERMINATES.
!     --------------------
USE constants_NSWC
IMPLICIT NONE
INTEGER, INTENT(IN)                     :: ka, n, iend
INTEGER, INTENT(OUT)                    :: ierr
REAL (dp), DIMENSION(:), INTENT(IN OUT) :: a, temp
INTEGER, DIMENSION(:), INTENT(IN OUT)   :: indx

!     Local variables
INTEGER   :: i, ibeg, ij, ik, il, j, jcol, jj, k, kcol, kcount, kj,  &
             kj0, kk, kl, km1, kp1, l, lj, lj0, lk, lmin, maxdim, mcol,  &
             ncol, nk, nm1, nmj, nmk, nn
REAL (dp) :: zero = 0._dp
REAL (dp) :: dsum, c, pmin, s

maxdim = ka*n
mcol = iend*ka
IF (iend == 0) GO TO 100

!           PROCESS THE FIRST IEND COLUMNS OF A

kcol = 0
DO k = 1, iend
  kk = kcol + k
  nk = kcol + n
  DO lk = kk, nk
    IF (a(lk) < zero) GO TO 20
    IF (a(lk) > zero) GO TO 30
  END DO
  GO TO 300
  
  20 l = lk - kcol
  lj0 = mcol + l
  DO lj = lj0, maxdim, ka
    a(lj) = -a(lj)
  END DO
  
  30 l = lk - kcol
  indx(k) = l
  IF (k == l) GO TO 32
  lj = lk
  DO kj = kk, maxdim, ka
    c = a(kj)
    a(kj) = a(lj)
    a(lj) = c
    lj = lj + ka
  END DO
  32 kcol = kcol + ka
END DO

!           PROCESS THE REMAINING COLUMNS OF A

100 nm1 = n - 1
ierr = 0
pmin = zero
ibeg = iend + 1
IF (ibeg == n) GO TO 190

k = ibeg
km1 = iend
kp1 = k + 1
kcol = mcol
kk = kcol + k
DO kcount = ibeg, nm1
  
!     SEARCH FOR THE K-TH PIVOT ELEMENT (K=IBEG, ..., N-1)
  
  l = k
  s = ABS(a(kk))
  DO i = kp1, n
    ik = kcol + i
    c = ABS(a(ik))
    IF (s >= c) CYCLE
    l = i
    s = c
  END DO
  
  IF (k > ibeg .AND. s >= pmin) GO TO 120
  pmin = s
  IF (s == zero) GO TO 300
  
!              INTERCHANGING ROWS K AND L
  
  120 indx(k) = l
  IF (k == l) GO TO 130
  kj0 = mcol + k
  lj  = mcol + l
  DO kj = kj0, maxdim, ka
    c = a(kj)
    a(kj) = a(lj)
    a(lj) = c
    lj = lj + ka
  END DO
  
!       COMPUTE THE K-TH ROW OF U (K=IBEG, ..., N-1)
  
  130 c = a(kk)
  IF (k > ibeg) GO TO 140
  kj0 = kk + ka
  DO kj = kj0, maxdim, ka
    a(kj) = a(kj)/c
  END DO
  GO TO 160
  
  140 kl = mcol + k
  DO l = ibeg, km1
    temp(l) = a(kl)
    kl = kl + ka
  END DO
  
  kj0 = kk + ka
  DO kj = kj0, maxdim, ka
    jcol = kj - k
    dsum = -a(kj)
    DO l = ibeg, km1
      lj = jcol + l
      dsum = dsum + temp(l)*a(lj)
    END DO
    a(kj) = -dsum/c
  END DO
  
!      COMPUTE THE K-TH COLUMN OF L (K=IBEG+1, ..., N)
  
  160 km1 = k
  k = kp1
  kp1 = k + 1
  kcol = kcol + ka
  kk = kcol + k
  DO l = ibeg, km1
    lk = kcol + l
    temp(l) = a(lk)
  END DO
  
  DO i = k, n
    il = mcol + i
    dsum = zero
    DO l = ibeg, km1
      dsum = dsum + a(il)*temp(l)
      il = il + ka
    END DO
    a(il) = a(il) - dsum
  END DO
END DO

!           CHECK THE N-TH PIVOT ELEMENT

190 ncol = maxdim - ka
nn = ncol + n
c = ABS(a(nn))
IF (c > pmin) GO TO 200
IF (c == zero) GO TO 300

!          REPLACE L WITH THE INVERSE OF L

200 IF (ibeg == n) GO TO 213
jj = mcol + ibeg
i = ka + 1
DO j = ibeg, nm1
  a(jj) = 1.0_dp/a(jj)
  temp(j) = a(jj)
  kj = jj
  DO km1 = j, nm1
    k = km1 + 1
    kj = kj + 1
    dsum = zero
    kl = kj
    DO l = j, km1
      dsum = dsum + a(kl)*temp(l)
      kl = kl + ka
    END DO
    a(kj) = -dsum/a(kl)
    temp(k) = a(kj)
  END DO
  jj = jj + i
END DO
213 a(nn) = 1.0_dp/a(nn)
IF (n == 1) RETURN

!       SOLVE UX = Y WHERE Y IS THE INVERSE OF L

DO nmk = 1, nm1
  k = n - nmk
  lmin = MAX(ibeg, k+1)
  kl = (lmin-1)*ka + k
  DO l = lmin, n
    temp(l) = a(kl)
    a(kl) = zero
    kl = kl + ka
  END DO
  
  kj0 = mcol + k
  DO kj = kj0, maxdim, ka
    dsum = -a(kj)
    lj = (kj - k) + lmin
    DO l = lmin, n
      dsum = dsum + temp(l)*a(lj)
      lj = lj + 1
    END DO
    a(kj) = -dsum
  END DO
END DO

!                 COLUMN INTERCHANGES

jcol = ncol - ka
DO nmj = 1, nm1
  j = n - nmj
  k = indx(j)
  IF (j == k) GO TO 251
  ij = jcol
  ik = (k-1)*ka
  DO i = 1, n
    ij = ij + 1
    ik = ik + 1
    c = a(ij)
    a(ij) = a(ik)
    a(ik) = c
  END DO
  251 jcol = jcol - ka
END DO
RETURN

!                    ERROR RETURN

300 ierr = 1
RETURN
END SUBROUTINE crout1
