      SUBROUTINE SSPLX1(TA,ITA,JTA,B0,C,M,N0,IND,IBASIS,R,Z,
     1    ITER,MXITER,EPS0,RERRMN,RERRMX,RERR,NUMLE,NUMGE,
     2    BI,XB,Y,BASIS,INDEX)
C     -----------------
C     NSTEP = 1   ELIMINATE THE NEGATIVE VARIABLES
C     NSTEP = 2   PHASE 1 OF THE SIMPLEX ALGORITHM
C     NSTEP = 3   PHASE 2 OF THE SIMPLEX ALGORITHM
C     -----------------
C     MXITER = THE MAXIMUM NUMBER OF ITERATIONS PERMITTED
C     ITER = THE NUMBER OF THE CURRENT ITERATION
C     ICOUNT = THE NUMBER OF ITERATIONS SINCE THE LAST INVERSION
C     -----------------
C     NUMLE = THE NUMBER OF .LE. CONSTRAINTS
C     NUMGE = THE NUMBER OF .GE. CONSTRAINTS
C     -----------------
C     THE ROUTINE ASSUMES THAT THE .LE. CONSTRAINTS PRECEDE THE .GE.
C     CONSTRAINTS AND THAT THE .EQ. CONSTRAINTS COME LAST. THERE ARE
C     M CONSTRAINTS. X(N0+I) IS THE SLACK, SURPLUS, OR ARTIFICIAL
C     VARIABLE FOR THE I-TH CONSTRAINT (I=1,...,M).
C     -----------------
C     N0 = THE NUMBER OF ORGINAL VARIABLES
C     NS = THE NUMBER OF ORGINAL AND SLACK VARIABLES
C     N  = THE NUMBER OF ORGINAL, SLACK, AND SURPLUS VARIABLES
C     NUM = THE TOTAL NUMBER OF VARIABLES
C     -----------------
C     RERRMN = THE SMALLEST RELATIVE ERROR TOLERANCE THAT IS USED
C     RERRMX = THE MINIMUM RELATIVE ERROR REQUIRING REINVERSION
C     RERR   = THE CURRENT RELATIVE ERROR TOLERANCE BEING EMPLOYED
C     -----------------
C     LET A DENOTE AN MXN0 MATRIX. THE TRANSPOSE OF A IS STORED IN
C     SPARSE FORM IN TA,ITA,JTA. ASSUME THAT
C         B0 = (B0(1),...,B0(M))
C         C  = (C(1),...,C(N0))
C         Z  = C(1)*X(1)+...+C(N0)*X(N0)
C     THE PROBLEM IS TO MAXIMIZE Z SUBJECT TO
C         AX(LE,EQ,GE)B0
C         X.GE.0
C     -----------------
C     ON INPUT IND CAN HAVE THE VALUES
C         IND = 0   NO BEGINNING BASIS IS PROVIDED BY THE USER
C         IND = 1   THE ARRAY IBASIS HAS BEEN SET BY THE USER
C     ON OUTPUT IND IS ASSIGNED ONE OF THE VALUES
C         IND = 0   Z WAS SUCCESSFULLY MAXIMIZED
C         IND = 1   THE PROBLEM HAS NO FEASIBLE SOLUTION
C         IND = 2   MXITER ITERATIONS WERE PERFORMED
C         IND = 3   SUFFICIENT ACCURACY CANNOT BE MAINTAINED
C         IND = 4   THE PROBLEM HAS AN UNBOUNDED SOLUTION
C         IND = 5   THERE IS AN INPUT ERROR
C         IND = 6   Z WAS POSSIBLY MAXIMIZED
C     -----------------
C     BASIS IS AN INTEGER ARRAY OF DIMENSION N0+M. FOR J.LE.N
C         BASIS(J) = 1  IF X(J) IS A BASIC VARIABLE
C         BASIS(J) = 0  IF X(J) IS NOT A BASIC VARIABLE
C     IF THE BASIC VARIABLES ARE X(I1),...,X(IM) THEN
C         IBASIS = (I1,...,IM)
C     ALSO XB(1),...,XB(M) ARE THE CORRESPONDING VALUES OF THE
C     BASIC VARIABLES.
C     -----------------
C     BI IS AN MXM ARRAY CONTAINING THE INVERSE OF THE BASIS MATRIX.
C     -----------------
C     R IS AN ARRAY OF DIMENSION N. ON OUTPUT R CONTAINS THE CURRENT
C     VALUE OF X. DURING COMPUTATION R NORMALLY CONTAINS THE REDUCED
C     COSTS USED FOR THE SELECTION OF THE VARIABLE TO BE MADE BASIC.
C     -----------------
      DIMENSION TA(1),ITA(1),JTA(1)
      DIMENSION B0(M),C(N0)
      DIMENSION BI(M,M),XB(M),Y(M),R(1)
      INTEGER   IBASIS(M),BASIS(1)
      INTEGER   BFLAG,INDEX(M)
      DOUBLE PRECISION DZERO,DONE,DSUM,DSUMP,DSUMN,DT
      DATA ZERO/0.0/,DZERO/0.D0/,ONE/1.0/,DONE/1.D0/
C     -----------------
C     ****** XMAX IS A MACHINE DEPENDENT CONSTANT. XMAX IS THE
C            LARGEST POSITIVE FLOATING POINT NUMBER.
C
                 XMAX = SPMPAR(3)
C     -----------------
      ITER=0
      ICOUNT=0
      MCHECK=MIN0(5,1+M/15)
      Z=ZERO
C
C                CHECK FOR INPUT ERRORS
C
      MS=NUMLE+NUMGE
      IF (M.LT.2.OR.N0.LT.2.OR.MS.GT.M) GO TO 12
      DO 10 I=1,M
      IF (B0(I)) 12,10,10
   10 XB(I)=ZERO
      RTOL=XMAX
      DO 11 I=1,N0
      IF (C(I).NE.ZERO) RTOL=AMIN1(ABS(C(I)),RTOL)
   11 CONTINUE
      RERR0=0.1*RERRMX
      RTOL=RERR0*RTOL
      GO TO 20
   12 IND=5
      RETURN
C
C     FORMATION OF THE IBASIS AND BASIS ARRAYS. (IF IND=1
C     THEN THE IBASIS ARRAY IS DEFINED BY THE USER.)
C
   20 NS=N0+NUMLE
      N=NS+NUMGE
      IF (IND.EQ.0) GO TO 30
      NUM=N
      DO 21 I=1,M
      IF (IBASIS(I).GT.N) NUM=NUM+1
   21 CONTINUE
      GO TO 32
   22 IF (IND.EQ.0) GO TO 590
      IND=0
C
   30 NUM=N0+M
      DO 31 I=1,M
   31 IBASIS(I)=N0+I
   32 BFLAG=0
      DO 33 I=1,N
   33 BASIS(I)=0
      DO 34 I=1,M
      KI=IBASIS(I)
   34 BASIS(KI)=1
      IF (IND.EQ.1) GO TO 100
C
C          CALCULATION OF XB AND BI WHEN IND=0
C
      RERR=RERRMN
      DO 41 J=1,M
      XB(J)=B0(J)
      DO 40 I=1,M
   40 BI(I,J)=ZERO
   41 BI(J,J)=ONE
      IF (NUMGE.EQ.0) GO TO 630
      JMIN=NUMLE+1
      DO 42 J=JMIN,MS
      XB(J)=-XB(J)
   42 BI(J,J)=-1.0
      GO TO 601
C
C                  REORDER THE BASIS
C
  100 IBEG=1
      IEND=M
      DO 102 I=1,M
      IF (IBASIS(I).LE.N0) GO TO 101
      INDEX(IBEG)=IBASIS(I)
      IBEG=IBEG+1
      GO TO 102
  101 INDEX(IEND)=IBASIS(I)
      IEND=IEND-1
  102 CONTINUE
      IF (IEND.EQ.M) GO TO 22
      DO 103 I=1,M
  103 IBASIS(I)=INDEX(I)
C
C            REINVERSION OF THE BASIS MATRIX
C
      DO 132 J=1,M
      KJ=IBASIS(J)
      IF (KJ.LE.N0) GO TO 110
      IF (KJ.LE.NS) GO TO 120
      IF (KJ.LE.N) GO TO 130
      GO TO 120
C
  110 DO 111 I=1,M
  111 BI(I,J)=ZERO
      LMIN=ITA(KJ)
      LMAX=ITA(KJ+1)-1
      IF (LMIN.GT.LMAX) GO TO 132
      DO 112 LL=LMIN,LMAX
      L=JTA(LL)
  112 BI(L,J)=TA(LL)
      GO TO 132
C
  120 L=KJ-N0
      DO 121 I=1,M
  121 BI(I,J)=ZERO
      BI(L,J)=ONE
      GO TO 132
C
  130 L=KJ-N0
      DO 131 I=1,M
  131 BI(I,J)=ZERO
      BI(L,J)=-1.0
  132 CONTINUE
C
      ICOUNT=0
      CALL CROUT1(BI,M,M,IEND,INDEX,Y,JCOL,IERR)
      IF (IERR.NE.0) GO TO 580
C
C         CHECK THE ACCURACY OF BI AND RESET RERR
C
      BNORM=ZERO
      DO 142 J=1,M
      KJ=IBASIS(J)
      IF (KJ.LE.N0) GO TO 140
      SUM=ONE
      GO TO 142
  140 SUM=ZERO
      LMIN=ITA(KJ)
      LMAX=ITA(KJ+1)-1
      DO 141 LL=LMIN,LMAX
  141 SUM=SUM+ABS(TA(LL))
  142 BNORM=AMAX1(BNORM,SUM)
C
      BINORM=ZERO
      DO 151 J=1,M
      SUM=ZERO
      DO 150 I=1,M
  150 SUM=SUM+ABS(BI(I,J))
  151 BINORM=AMAX1(BINORM,SUM)
      RERR=AMAX1(RERRMN,EPS0*BNORM*BINORM)
      IF (RERR.GT.1.E-2) GO TO 580
      BFLAG=0
C
C                 RECALCULATION OF XB
C
  180 RERR1=AMIN1(RERR,RERRMX)
      DO 184 I=1,M
      DSUMP=DZERO
      DSUMN=DZERO
      DO 183 L=1,M
      DT=BI(I,L)*B0(L)
      IF (DT) 181,183,182
  181 DSUMN=DSUMN+DT
      GO TO 183
  182 DSUMP=DSUMP+DT
  183 CONTINUE
      XB(I)=DSUMP+DSUMN
      S=DSUMP
      T=DSUMN
      TOL=RERR1*AMAX1(S,-T)
      IF (ABS(XB(I)).LE.TOL) XB(I)=ZERO
  184 CONTINUE
      GO TO 601
C
C     FIND THE NEXT VECTOR A(--,JP) TO BE INSERTED INTO
C                       THE BASIS
C
  200 JP=0
      RMIN=ZERO
      IF (NSTEP.EQ.3) RMIN=-RTOL
      DO 201 J=1,N0
      IF (BASIS(J).EQ.1) GO TO 201
      IF (R(J).GE.RMIN) GO TO 201
      JP=J
      RMIN=R(J)
  201 CONTINUE
      IF (N0.EQ.N) GO TO 203
      JMIN=N0+1
      RMIN=RMIN*1.1
      DO 202 J=JMIN,N
      IF (BASIS(J).EQ.1) GO TO 202
      IF (R(J).GE.RMIN) GO TO 202
      JP=J
      RMIN=R(J)
  202 CONTINUE
  203 IF (JP.NE.0) GO TO 300
      IF (NSTEP-2) 800,230,250
C
C     INSERT THE VALUES OF THE ORGINAL, SLACK, AND SURPLUS
C             VARIABLES INTO R. THEN TERMINATE.
C
  220 DO 221 J=1,N
  221 R(J)=ZERO
      DO 222 I=1,M
      KI=IBASIS(I)
      IF (KI.LE.N) R(KI)=XB(I)
  222 CONTINUE
      RETURN
C
C              COMPLETION OF THE NSTEP=2 CASE
C
  230 DO 231 I=1,M
      IF (IBASIS(I).LE.N) GO TO 231
      IF (XB(I)) 231,231,800
  231 CONTINUE
      GO TO 680
C
  240 IF (ICOUNT.GE.5) GO TO 100
      IND=1
      GO TO 220
C
C              COMPLETION OF THE NSTEP=3 CASE
C
  250 IF (RERR.GT.RERRMX) GO TO 251
      IND=0
      GO TO 800
  251 IF (ICOUNT.GE.5) GO TO 100
      IND=6
      GO TO 800
C
C     IF MXITER ITERATIONS HAVE NOT BEEN PERFORMED THEN
C     BEGIN THE NEXT ITERATION. COMPUTE THE JP-TH COLUMN
C               OF BI*A AND STORE IT IN Y.
C
  300 IF (ITER.LT.MXITER) GO TO 301
      IND=2
      GO TO 220
  301 ITER=ITER+1
      ICOUNT=ICOUNT+1
      IF (JP.GT.NS) GO TO 330
      IF (JP.GT.N0) GO TO 320
C
      LMIN=ITA(JP)
      LMAX=ITA(JP+1)-1
      IF (LMIN.LE.LMAX) GO TO 305
      IND=4
      GO TO 220
  305 AMAX=ZERO
      DO 306 LL=LMIN,LMAX
  306 AMAX=AMAX1(ABS(TA(LL)),AMAX)
C
  310 RERR1=RERR0*AMAX
      DO 313 I=1,M
      DSUM=DZERO
      DO 311 LL=LMIN,LMAX
      L=JTA(LL)
  311 DSUM=DSUM+DBLE(BI(I,L)*TA(LL))
      Y(I)=DSUM
      IF (ABS(Y(I)).GE.1.E-3) GO TO 313
      BMAX=ZERO
      DO 312 L=1,M
  312 BMAX=AMAX1(ABS(BI(I,L)),BMAX)
      TOL=RERR1*BMAX
      IF (ABS(Y(I)).LE.TOL) Y(I)=ZERO
  313 CONTINUE
      GO TO 350
C
  320 L=JP-N0
      DO 321 I=1,M
  321 Y(I)=BI(I,L)
      GO TO 340
C
  330 L=JP-N0
      DO 331 I=1,M
  331 Y(I)=-BI(I,L)
C
  340 DO 342 I=1,M
      IF (ABS(Y(I)).GE.1.E-3) GO TO 342
      BMAX=ZERO
      DO 341 L=1,M
  341 BMAX=AMAX1(ABS(BI(I,L)),BMAX)
      TOL=RERR0*BMAX
      IF (ABS(Y(I)).LE.TOL) Y(I)=ZERO
  342 CONTINUE
C
  350 DO 352 I=1,M
      IF (Y(I)) 351,352,351
  351 IF (NSTEP-2) 400,430,440
  352 CONTINUE
      R(JP)=ZERO
      ITER=ITER-1
      ICOUNT=ICOUNT-1
      GO TO 200
C
C     FINDING THE VARIABLE XB(IP) TO BE MADE NONBASIC
C               FOR THE NSTEP=1 CASE
C
  400 NPOS=0
      IP=0
      EPS=ZERO
      EPSI=XMAX
      DO 403 I=1,M
      IF (XB(I).LT.ZERO.OR.Y(I).LE.ZERO) GO TO 403
      RATIO=XB(I)/Y(I)
      IF (RATIO-EPSI) 401,402,403
  401 EPSI=RATIO
      NPOS=1
      INDEX(1)=I
      GO TO 403
  402 NPOS=NPOS+1
      INDEX(NPOS)=I
  403 CONTINUE
      IF (NPOS.EQ.0) GO TO 420
      IF (EPSI.EQ.ZERO) GO TO 460
C
      DO 410 I=1,M
      IF (XB(I).GE.ZERO.OR.Y(I).GE.ZERO) GO TO 410
      RATIO=XB(I)/Y(I)
      IF (RATIO.GT.EPSI) GO TO 410
      IF (RATIO.LT.EPS) GO TO 410
      EPS=RATIO
      IP=I
  410 CONTINUE
      IF (IP.NE.0) GO TO 500
      GO TO 460
C
  420 DO 421 I=1,M
      IF (XB(I).GE.ZERO.OR.Y(I).GE.ZERO) GO TO 421
      RATIO=XB(I)/Y(I)
      IF (RATIO.LT.EPS) GO TO 421
      EPS=RATIO
      IP=I
  421 CONTINUE
      GO TO 500
C
C     FINDING THE VARIABLE XB(IP) TO BE MADE NONBASIC
C               FOR THE NSTEP=2 CASE
C
  430 NPOS=0
      EPSI=XMAX
      DO 433 I=1,M
      IF (Y(I).LE.ZERO) GO TO 433
      RATIO=XB(I)/Y(I)
      IF (RATIO-EPSI) 431,432,433
  431 EPSI=RATIO
      NPOS=1
      INDEX(1)=I
      GO TO 433
  432 NPOS=NPOS+1
      INDEX(NPOS)=I
  433 CONTINUE
      GO TO 450
C
C     FINDING THE VARIABLE XB(IP) TO BE MADE NONBASIC
C               FOR THE NSTEP=3 CASE
C
  440 NPOS=0
      EPSI=XMAX
      DO 445 I=1,M
      IF (Y(I)) 441,445,442
  441 IF (IBASIS(I).LE.N) GO TO 445
      IP=I
      GO TO 500
  442 RATIO=XB(I)/Y(I)
      IF (RATIO-EPSI) 443,444,445
  443 EPSI=RATIO
      NPOS=1
      INDEX(1)=I
      GO TO 445
  444 NPOS=NPOS+1
      INDEX(NPOS)=I
  445 CONTINUE
C
  450 IF (NPOS.NE.0) GO TO 460
      IF (ICOUNT.GE.5) GO TO 100
      IND=4
      GO TO 220
C
C              TIE BREAKING PROCEDURE
C
  460 IP=INDEX(1)
      IF (NPOS.EQ.1) GO TO 500
      IP = 0
      BMIN=XMAX
      CMIN=XMAX
      DO 464 II=1,NPOS
      I=INDEX(II)
      L=IBASIS(I)
      IF (L.GT.N0) GO TO 461
      IF (C(L).LE.ZERO) CMIN=AMIN1(ZERO,CMIN)
      IF (C(L).GT.CMIN) GO TO 464
      IMIN=I
      CMIN=C(L)
      GO TO 464
  461 IF (L.LE.N) GO TO 462
      IP=I
      GO TO 500
  462 LROW=L-N0
      S=B0(LROW)
      IF (LROW.GT.NUMLE) GO TO 463
      IF (S.GT.BMIN) GO TO 464
      IP=I
      BMIN=S
      GO TO 464
  463 S=-S
      BMIN=AMIN1(ZERO,BMIN)
      IF (S.GT.BMIN) GO TO 464
      IP=I
      BMIN=S
  464 CONTINUE
      IF (CMIN.LE.ZERO.OR.IP.EQ.0) IP=IMIN
C
C               TRANSFORMATION OF XB
C
  500 IF (XB(IP).EQ.ZERO) GO TO 510
      RERR1=AMIN1(10.0*RERR,RERR0)
      CONST=XB(IP)/Y(IP)
      DO 501 I=1,M
      S=XB(I)
      T=CONST*Y(I)
      XB(I)=S-T
      IF (S.GE.ZERO.AND.XB(I).LT.ZERO) XB(I)=ZERO
      TOL=RERR1*AMAX1(ABS(S),ABS(T))
      IF (ABS(XB(I)).LE.TOL) XB(I)=ZERO
  501 CONTINUE
      XB(IP)=CONST
C
C               TRANSFORMATION OF BI
C
  510 DO 512 J=1,M
      IF (BI(IP,J).EQ.ZERO) GO TO 512
      CONST=BI(IP,J)/Y(IP)
      DO 511 I=1,M
  511 BI(I,J)=DBLE(BI(I,J))-DBLE(CONST*Y(I))
      BI(IP,J)=CONST
  512 CONTINUE
C
C             UPDATING IBASIS AND BASIS
C
      IOUT=IBASIS(IP)
      IBASIS(IP)=JP
      BASIS(IOUT)=0
      BASIS(JP)=1
      IF (IOUT.GT.N) NUM=NUM-1
C
C        CHECK THE ACCURACY OF BI AND RESET RERR
C
      IF (RERR.GT.RERRMX) GO TO 530
      K=0
      DO 521 J=1,M
      KJ=IBASIS(J)
      IF (KJ.GT.N0) GO TO 521
      SUM=ZERO
      LMIN=ITA(KJ)
      LMAX=ITA(KJ+1)-1
      DO 520 LL=LMIN,LMAX
      L=JTA(LL)
  520 SUM=SUM+BI(J,L)*TA(LL)
      RERR=AMAX1(RERR,ABS(ONE-SUM))
      K=K+1
      IF (K.GE.MCHECK) GO TO 522
  521 CONTINUE
  522 IF (RERR.LE.RERRMX) GO TO 600
C
C        THE ACCURACY CRITERIA ARE NOT SATISFIED
C
  530 IF (ICOUNT.LT.5) GO TO 600
      BFLAG=1
      GO TO 100
C
  580 IF (ITER.EQ.0) GO TO 12
      IF (BFLAG.EQ.0) GO TO 590
      BFLAG=0
      DO 581 IP=1,M
      IF (JP.EQ.IBASIS(IP)) GO TO 582
  581 CONTINUE
  582 IBASIS(IP)=IOUT
      BASIS(JP)=0
      BASIS(IOUT)=1
      IF (IOUT.GT.N) NUM=NUM+1
      GO TO 100
C
  590 IND=3
      GO TO 220
C
C       SET UP THE R ARRAY FOR THE NSTEP=1 CASE
C
  600 IF (NSTEP-2) 601,630,700
  601 DO 602 J=1,M
      IF (XB(J)) 610,602,602
  602 CONTINUE
      GO TO 630
C
  610 NSTEP=1
      M0=0
      DO 612 L=1,M
      IF (XB(L)) 611,612,612
  611 M0=M0+1
      INDEX(M0)=L
  612 CONTINUE
C
      RERR1=AMIN1(RERR,RERRMX)
      DO 623 J=1,M
      DSUMP=DZERO
      DSUMN=DZERO
      DO 622 LL=1,M0
      L=INDEX(LL)
      IF (BI(L,J)) 620,622,621
  620 DSUMN=DSUMN+DBLE(BI(L,J))
      GO TO 622
  621 DSUMP=DSUMP+DBLE(BI(L,J))
  622 CONTINUE
      Y(J)=DSUMP+DSUMN
      S=DSUMP
      T=DSUMN
      TOL=RERR1*AMAX1(S,-T)
      IF (ABS(Y(J)).LE.TOL) Y(J)=ZERO
  623 CONTINUE
      GO TO 650
C
C       SET UP THE R ARRAY FOR THE NSTEP=2 CASE
C
  630 IF (N.EQ.NUM) GO TO 680
      NSTEP=2
      M0=0
      DO 631 L=1,M
      IF (IBASIS(L).LE.N) GO TO 631
      M0=M0+1
      INDEX(M0)=L
  631 CONTINUE
C
      RERR1=AMIN1(RERR,RERRMX)
      DO 643 J=1,M
      DSUMP=DZERO
      DSUMN=DZERO
      DO 642 LL=1,M0
      L=INDEX(LL)
      IF (BI(L,J)) 640,642,641
  640 DSUMN=DSUMN+DBLE(BI(L,J))
      GO TO 642
  641 DSUMP=DSUMP+DBLE(BI(L,J))
  642 CONTINUE
      Y(J)=-(DSUMP+DSUMN)
      S=DSUMP
      T=DSUMN
      TOL=RERR1*AMAX1(S,-T)
      IF (ABS(Y(J)).LE.TOL) Y(J)=ZERO
  643 CONTINUE
C
  650 DO 652 J=1,N0
      SUM=ZERO
      IF (BASIS(J).NE.0) GO TO 652
      LMIN=ITA(J)
      LMAX=ITA(J+1)-1
      IF (LMIN.GT.LMAX) GO TO 652
      DO 651 LL=LMIN,LMAX
      L=JTA(LL)
  651 SUM=SUM+Y(L)*TA(LL)
  652 R(J)=SUM
C
  660 IF (N0.EQ.NS) GO TO 670
      JMIN=N0+1
      DO 661 J=JMIN,NS
      R(J)=ZERO
      IF (BASIS(J).NE.0) GO TO 661
      JJ=J-N0
      R(J)=Y(JJ)
  661 CONTINUE
C
  670 IF (NS.EQ.N) GO TO 200
      JMIN=NS+1
      DO 671 J=JMIN,N
      R(J)=ZERO
      IF (BASIS(J).NE.0) GO TO 671
      JJ=J-N0
      R(J)=-Y(JJ)
  671 CONTINUE
      GO TO 200
C
C      SET UP A NEW R ARRAY FOR THE NSTEP=3 CASE
C
  680 NSTEP=3
      RERR1=AMIN1(RERR,RERRMX)
      DO 684 J=1,M
      DSUMP=DZERO
      DSUMN=DZERO
      DO 683 L=1,M
      IL=IBASIS(L)
      IF (IL.GT.N0) GO TO 683
      DT=C(IL)*BI(L,J)
      IF (DT) 681,683,682
  681 DSUMN=DSUMN+DT
      GO TO 683
  682 DSUMP=DSUMP+DT
  683 CONTINUE
      Y(J)=DSUMP+DSUMN
      S=DSUMP
      T=DSUMN
      TOL=RERR1*AMAX1(S,-T)
      IF (ABS(Y(J)).LE.TOL) Y(J)=ZERO
  684 CONTINUE
C
      DO 691 J=1,N0
      R(J)=ZERO
      IF (BASIS(J).NE.0) GO TO 691
      DSUM=-C(J)
      LMIN=ITA(J)
      LMAX=ITA(J+1)-1
      R(J)=-C(J)
      IF (LMIN.GT.LMAX) GO TO 691
      DO 690 LL=LMIN,LMAX
      L=JTA(LL)
  690 DSUM=DSUM+DBLE(Y(L)*TA(LL))
      R(J)=DSUM
      IF (R(J).GE.ZERO) GO TO 691
      TOL=RERR0*ABS(C(J))
      IF (ABS(R(J)).LE.TOL) R(J)=ZERO
  691 CONTINUE
      GO TO 660
C
C       UPDATE THE R ARRAY FOR THE NSTEP=3 CASE
C
  700 CONST=R(JP)
      DO 703 J=1,N0
      IF (BASIS(J).EQ.0) GO TO 701
      R(J)=ZERO
      GO TO 703
  701 SUM=ZERO
      LMIN=ITA(J)
      LMAX=ITA(J+1)-1
      IF (LMIN.GT.LMAX) GO TO 703
      DO 702 LL=LMIN,LMAX
      L=JTA(LL)
  702 SUM=SUM+BI(IP,L)*TA(LL)
      R(J)=R(J)-CONST*SUM
      IF (R(J).GE.ZERO) GO TO 703
      TOL=RERR0*ABS(C(J))
      IF (ABS(R(J)).LE.TOL) R(J)=ZERO
  703 CONTINUE
C
  710 IF (N0.EQ.NS) GO TO 720
      JMIN=N0+1
      DO 712 J=JMIN,NS
      IF (BASIS(J).EQ.0) GO TO 711
      R(J)=ZERO
      GO TO 712
  711 JJ=J-N0
      R(J)=R(J)-CONST*BI(IP,JJ)
  712 CONTINUE
C
  720 IF (NS.EQ.N) GO TO 200
      JMIN=NS+1
      DO 722 J=JMIN,N
      IF (BASIS(J).EQ.0) GO TO 721
      R(J)=ZERO
      GO TO 722
  721 JJ=J-N0
      R(J)=R(J)+CONST*BI(IP,JJ)
  722 CONTINUE
      GO TO 200
C     ------------------------------------------------------------------
C
C        REFINE XB AND STORE THE RESULT IN Y
C
  800 DO 801 I=1,M
      R(I)=ZERO
  801 Y(I)=ZERO
C
      DO 831 J=1,M
      KJ=IBASIS(J)
      IF (KJ.LE.N0) GO TO 810
      IF (KJ.LE.NS) GO TO 820
      IF (KJ.LE.N) GO TO 830
      GO TO 820
C
  810 LMIN=ITA(KJ)
      LMAX=ITA(KJ+1)-1
      DO 811 LL=LMIN,LMAX
      L=JTA(LL)
      DT=DBLE(R(L))+DBLE(Y(L))
      DT=DT+DBLE(TA(LL)*XB(J))
      R(L)=DT
  811 Y(L)=DT-DBLE(R(L))
      GO TO 831
C
  820 L=KJ-N0
      DT=DBLE(R(L))+DBLE(Y(L))
      DT=DT+DBLE(XB(J))
      R(L)=DT
      Y(L)=DT-DBLE(R(L))
      GO TO 831
C
  830 L=KJ-N0
      DT=DBLE(R(L))+DBLE(Y(L))
      DT=DT-DBLE(XB(J))
      R(L)=DT
      Y(L)=DT-DBLE(R(L))
  831 CONTINUE
C
      DO 840 I=1,M
      DT=DBLE(R(I))+DBLE(Y(I))
  840 R(I)=DBLE(B0(I))-DT
C
  850 RERR1=AMIN1(RERR0,RERR)
      DO 858 I=1,M
      IF (XB(I)) 851,857,852
  851 DSUMP=DZERO
      DSUMN=XB(I)
      GO TO 853
  852 DSUMP=XB(I)
      DSUMN=DZERO
  853 DO 856 L=1,M
      DT=BI(I,L)*R(L)
      IF (DT) 854,856,855
  854 DSUMN=DSUMN+DT
      GO TO 856
  855 DSUMP=DSUMP+DT
  856 CONTINUE
      Y(I)=DSUMP+DSUMN
      IF (XB(I).GT.ZERO.AND.Y(I).LT.ZERO) GO TO 857
      IF (XB(I).LT.ZERO.AND.Y(I).GT.ZERO) GO TO 857
      S=DSUMP
      T=DSUMN
      TOL=RERR1*AMAX1(S,-T)
      IF (ABS(Y(I)).GT.TOL) GO TO 858
  857 Y(I)=ZERO
  858 CONTINUE
      IF (NSTEP-2) 860,870,880
C
C          CHECK THE REFINEMENT (NSTEP=1)
C
  860 DO 861 I=1,M
      IF (Y(I).GE.ZERO) GO TO 861
      IF (Y(I).LT.-RERR0) GO TO 240
      Y(I)=ZERO
  861 XB(I)=Y(I)
      GO TO 630
C
C          CHECK THE REFINEMENT (NSTEP=2)
C
  870 DO 871 I=1,M
      IF (IBASIS(I).LE.N) GO TO 871
      IF (Y(I).GT.RERR0) GO TO 240
      Y(I)=ZERO
  871 XB(I)=Y(I)
      GO TO 680
C
C               COMPUTE Z  (NSTEP=3)
C
  880 DSUM=DZERO
      DO 881 I=1,M
      KI=IBASIS(I)
      IF (KI.GT.N0) GO TO 881
      DSUM=DSUM+DBLE(C(KI)*Y(I))
  881 XB(I)=Y(I)
      Z=DSUM
      GO TO 220
      END
