
c From CMUVM.BITNET!32ZEJ7N Fri Mar 03 08:25:31 EDT 1989
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION X(16),PWR(101,15),FTBLE(101),T(202),
     *   WTBLE(101),UPTBL(101),INDUP(101),ALTBL(101),
     *   INDLW(101),INDF(101),ERROR(102),ERTST(9),ALC(101,16)
C
C THE APPROXIMATION PROBLEMS RUN WITH THIS DEMONSTRATION DRIVER
C ARE AS FOLLOWS.  NOTE THAT P(1) MEANS THE FIRST NUMERATOR
C COEFFICIENT, WHILE P(1.0) (RESP. P(Z)) MEANS THE VALUE OF THE
C NUMERATOR AT 1.0 (RESP. AT Z), AND SIMILARLY FOR Q.
C
C IN JOB 1, THE GRID IS A 101 POINT EQUALLY SPACED SUBDIVISION OF
C   THE CLOSED INTERVAL (0.0,2.0).  F(Z) IS 1.0 FOR 0.0 .LE. Z .LE. 1.0,
C   UNDEFINED FOR 1.0 .LT. Z .LT. 1.2, AND 0.0 FOR 1.2 .LE. Z .LE. 2.0.
C   (P/Q)(Z) IS (P(1)+P(2)*Z)/(Q(1)+Q(2)*Z+Q(3)*Z**2).  WE REQUIRE
C   (P/Q)(Z) .GE. 0.9999 FOR 0.0 .LE. Z .LE. 0.3, (P/Q)(Z) .GE. 0.0 FOR
C   1.0 .LT. Z .LE. 2.0, AND Q(Z) .GE. 1000.0*B**(-T) FOR ALL Z,
C   WHERE T IS THE NUMBER OF BASE B DIGITS IN THE MANTISSA OF
C   FLOATING POINT NUMBERS.  (THIS LAST REQUIREMENT IS USED IN THE
C   FIRST FOUR JOBS TO AVOID THE POSSIBILITY THAT A BAD DENOMINATOR
C   VALUE AT A POINT WHERE F IS UNDEFINED, SO THE DIFFERENTIAL
C   CORRECTION INEQUALITIES DO NOT PROVIDE ANY CONTROL ON THE
C   DENOMINATOR, MIGHT PREMATURELY TERMINATE THE PROGRAM.)
C
C IN JOB 2, THE GRID IS A 21 POINT EQUALLY SPACED SUBDIVISION OF
C   THE CLOSED INTERVAL (0.0,PI).  F(Z) IS 1.0 FOR 0.0 .LE. Z .LE. PI/2,
C   UNDEFINED FOR PI/2 .LT. Z .LT. 3*PI/5, AND 0.0 FOR 3*PI/5 .LE. Z
C   .LE. PI.  (P/Q)(Z) IS (P(1)+P(2)*COSZ+...+P(7)*COS6Z)/(Q(1)+Q(2)*Z+
C   ...+Q(8)*COS7Z).  THE WEIGHT FUNCTION (FOR WEIGHTED ERROR CURVE
C   APPROXIMATION) IS 0.5 FOR 0.0 .LE. Z .LE. PI/2, UNDEFINED FOR
C   PI/2 .LT. Z .LT. 3*PI/5, AND 1.0 FOR 3*PI/5 .LE. Z .LE. PI.
C   WE REQUIRE (P/Q)(Z) .GE. 0.0 FOR PI/2 .LT. Z .LE. PI AND
C   Q(Z) .GE. 1000.0*B**(-T) FOR ALL Z.
C
C JOB 3 IS THE SAME AS JOB 2 EXCEPT THE GRID IS A 41 POINT EQUALLY
C   SPACED SUBDIVISION OF THE CLOSED INTERVAL (0.0,PI), AND WE USE
C   THE P/Q PRODUCED IN JOB 2 FOR INITIALIZATION.
C
C IN JOB 4 THE GRID IS THE SET OF ORDERED PAIRS (Z1,Z2) WHERE Z1
C   AND Z2 ARE NONNEGATIVE INTEGERS WITH Z1**2+Z2**2 .LE. 16.0.
C   F(Z1,Z2) IS SQRT(16.0-Z1**2-Z2**2).  (P/Q)(Z1,Z2) IS (P(1)+
C   P(2)*Z1**2+P(3)*Z2**2)/(Q(1)+Q(2)*(Z1*Z2)**2).  WE REQUIRE
C   Q(Z1,Z2) .GE. 1000.0*B**(-T) FOR ALL (Z1,Z2).
C
C IN JOB 5 THE GRID IS A 21 POINT EQUALLY SPACED SUBDIVISION OF THE
C   CLOSED INTERVAL (0.0,2.0).  F(Z) IS 1.0 FOR Z=0.0 AND 0.0 FOR
C   0.0 .LT. Z .LE. 2.0.  (P/Q)(Z) IS P(1)/(Q(1)+Q(2)*Z).  THIS IS AN
C   EXAMPLE WHERE NO BEST APPROXIMATION EXISTS.
C
C JOB 6 IS THE SAME AS JOB 5 EXCEPT WE REQUIRE Q(Z) .GE. 10.0**(-7) ON
C   THE GRID.  HERE A BEST APPROXIMATION DOES EXIST.
C
C IN JOB 7 THE GRID IS A 101 POINT EQUALLY SPACED SUBDIVISION OF THE
C   CLOSED INTERVAL (0.0,2.0).  F(Z) IS 3.0/(1.0+Z) + G(Z), WHERE
C   G HAS VALUES 0.2, -0.2, 0.2, AND -0.2 AT Z = 0.0, 0.2, 1.0,
C   AND 2.0 RESPECTIVELY, AND G IS LINEAR BETWEEN THESE POINTS.
C   (P/Q)(Z) IS (P(1)+P(2)*Z)/(Q(1)+Q(2)*Z+Q(3)*Z**2).  THIS IS AN
C   EXAMPLE WHERE THE BEST APPROXIMATION IS DEGENERATE.
C
C IN JOB 8 THE GRID IS A 101 POINT EQUALLY SPACED SUBDIVISION OF
C   THE CLOSED INTERVAL (0.0,2.0), F(Z) IS EXP(Z) FOR ALL Z, AND
C   (P/Q)(Z) IS (P(1)+P(2)*Z+P(3)*Z**2)/Q(1).  WE REQUIRE Q(1)=1.0,
C   (P/Q)(1.0) = F(1.0), AND 2.0*P(3) .GE. 3.0.  THIS EXAMPLE
C   ILLUSTRATES THE USE OF EQUALITY CONSTRAINTS (INCLUDING AN
C   INTERPOLATORY CONSTRAINT) AND AN ADDITIONAL INEQUALITY
C   CONSTRAINT FROM OPTION 4.  NOTE THAT IT IS NOT NECESSARY TO
C   EXPLICITLY IMPOSE THE CONSTRAINT Q(1)=1.0 FOR (GENERALIZED)
C   POLYNOMIAL APPROXIMATION AS LONG AS NO CONSTRAINTS FROM
C   OPTION 4 HAVE NONZERO RIGHT SIDES, BUT IT MAKES LESS WORK FOR
C   THE LINEAR PROGRAMMING SUBROUTINE IF THE CONSTRAINT IS IMPOSED.
C
C JOB 9 IS THE SAME AS JOB 1 EXCEPT THAT THE CONSTRAINTS Q(Z) .GE.
C   1000.0*B**(-T) FOR ALL Z ARE REPLACED BY THE CONSTRAINTS Q(Z) .GE.
C   0.5 FOR ALL Z, AND THE CONSTRAINTS (P/Q)(Z) .GE. 0.9999 FOR
C   0.0 .LE. Z .LE. 0.3 ARE REPLACED BY THE CONSTRAINT P(1) =
C   (Q(1)+Q(2))/2.0 + 1.0.
C
C THE DATA CARDS NEEDED FOR THESE JOBS ARE AS FOLLOWS (WITH THE
C CS IN THE FIRST COLUMN REPLACED BY BLANKS).
C   1    2    3    1  101
C     1011
C                0.0                 2.0
C   1    7    8    1   21
C    21011
C                0.03.141592653589793238
C   1    7    8    1   41
C   121011
C                0.03.141592653589793238
C   2    3    2    5    5
C     1001
C                0.0                 0.0                 4.0
C   1    1    2    1   21
C        0
C                0.0                 2.0
C   1    1    2    1   21
C        1
C                0.0                 2.0
C   1    2    3    1  101
C        0
C                0.0                 2.0
C   1    3    1    1  101
C        2
C   2    1
C                0.0                 2.0
C   1    2    3    1  101
C     1013
C   1    0
C                0.0                 2.0
C   0    0    0    0    0
C
C SET MACHINE DEPENDENT PARAMETERS FOR THE DEMONSTRATION DRIVER.
C SET INPUT AND OUTPUT UNIT NUMBERS.
      NREAD=I1MACH(1)
      NWRIT=I1MACH(2)
C SET PRECISION DEPENDENT CONSTANTS.
      ONE=1.0D0
      ZERO=ONE-ONE
      TWO=ONE+ONE
      THREE=ONE+TWO
      FOUR=TWO+TWO
      TEN=FOUR+FOUR+TWO
      SIXTN=TEN+THREE+THREE
C SET SPCMN=B**(-T), WHERE B IS THE BASE AND T IS THE NUMBER
C OF BASE B DIGITS IN THE MANTISSA.  SPCMN IS THE MINIMUM
C RELATIVE SPACING ABS((X1-X2)/X2) BETWEEN TWO SUCCESSIVE
C FLOATING POINT NUMBERS, SO IT IS THE SPACING BETWEEN TWO
C SUCCESSIVE FLOATING POINT NUMBERS IN THE CLOSED INTERVAL
C (0.1,1.0).  WE ALSO HAVE SPCMN=10.0**(-T*(LOG10)(B))=
C 10.0**(-TNMAN), WHERE TNMAN IS THE BASE 10 EQUIVALENT OF
C THE LENGTH OF THE MANTISSA.
      SPCMN=D1MACH(3)
C SET THE ERROR NORM TEST TOLERANCE ERTOL=MAX(10.0**(-11),
C 1000.0*SPCMN).
      ERTOL=SPCMN*TEN**3
      IF(ERTOL-TEN**(-11))10,20,20
   10 ERTOL=TEN**(-11)
C END OF SETTING MACHINE AND PRECISION DEPENDENT CONSTANTS FOR THE
C DEMONSTRATION DRIVER.
C
   20 JOBNM=1
C SET THE ERROR NORMS COMPUTED WITH THE IBM 3090 AT CENTRAL MICHIGAN
C UNIVERSITY FOR COMPARISON PURPOSES.  THIS MACHINE HAS 14 BASE 16
C DIGITS IN DOUBLE PRECISION.
      ERTST(1)=0.369559497889575D0
      ERTST(2)=0.8792063120*TEN**(-5)
      ERTST(3)=0.54401686050*TEN**(-4)
      ERTST(4)=0.447875655532295D0
      ERTST(5)=0.22*TEN**(-13)
      ERTST(6)=0.999998000*TEN**(-6)
      ERTST(7)=0.200000000000000D0
      ERTST(8)=0.148781030293754D0
      ERTST(9)=0.380060239141926D0
C
C READ NUMDM, NUMP=NUMBER OF COEFFICIENTS IN THE NUMERATOR,
C NUMQ=NUMBER OF COEFFICIENTS IN THE DENOMINATOR, NRWGR=
C NUMBER OF ROWS IN THE GRID, NCLGR=NUMBER OF
C COLUMNS IN THE GRID.
   30 READ(NREAD, 40 )NUMDM,NUMP,NUMQ,NRWGR,NCLGR
   40 FORMAT(5I5)
C
      NMPP1=NUMP+1
      NMPPQ=NUMP+NUMQ
      IF(NUMDM) 60 , 50 , 60
C A DATA CARD WITH 0 FOR NUMDM IS USED AS A SIGNAL TO STOP.
   50 STOP
   60 WRITE(NWRIT, 70 )JOBNM
   70 FORMAT(///30H ***********THIS IS JOB NUMBER,I4)
      NUMGR=NRWGR*NCLGR
C
C READ THE OPTION SWITCH IOPT.  INSTRUCTIONS FOR SETTING IT
C ARE IN THE INITIAL COMMENTS OF CDFCOR.
      READ(NREAD, 80)IOPT
   80 FORMAT(I10)
C
      IONE=IOPT-(IOPT/10)*10
C
C IF THE ONES DIGIT OF IOPT IS .GT. 1, READ THE NUMBER MEQ OF
C EQUALITY CONSTRAINTS AND THE NUMBER MINEQ OF EXTRA INEQUALITY
C CONSTRAINTS (NOT INCLUDING THE DENOMINATOR OR RESTRICTED
C RANGE CONSTRAINTS SET ELSEWHERE).
      IF(IONE-1)87,87,83
   83 READ(NREAD,85)MEQ,MINEQ
   85 FORMAT(2I5)
C
C READ THE ENDPOINTS INTO ASWP1, ASWP2 IN THE ONE DIMENSIONAL CASE.
C IN THE TWO DIMENSIONAL CASE ASWP1, ASWP2 AND ANEP1, ANEP2 ARE THE
C COORDINATES OF THE LOWER LEFT AND UPPER RIGHT CORNERS OF
C THE GRID RESPECTIVELY.
   87 IF(NUMDM-1)90,90,110
   90 READ(NREAD,100)ASWP1,ASWP2
  100 FORMAT(2F20.10)
      GO TO 130
  110 READ(NREAD, 120)ASWP1,ASWP2,ANEP1,ANEP2
  120 FORMAT(4F20.10)
C PUT THE GRID POINTS INTO T.
  130 CALL STCOMP(NUMDM,NRWGR,NCLGR,ASWP1,ASWP2,ANEP1,
     *ANEP2,T)
C
C SET UP PWR AND FTBLE.
      I=0
      IDBM=-1
      IDB=0
  140 I=I+1
      IF(I-NUMGR) 150, 150, 540
  150 IDBM=IDBM+2
      IDB=IDB+2
C
C COMPUTE THE BASIS FUNCTIONS AND F AND, IF NECESSARY,
C WTBLE, INDUP, UPTBL, INDLW, ALTBL, AND INDF.
C IN THE ONE DIMENSIONAL CASE THE COORDINATE WILL BE T(I),
C WHILE IN THE TWO DIMENSIONAL CASE THE COORDINATES WILL BE
C (T(IDBM),T(IDB)).
      GO TO ( 160, 230, 230, 300, 160, 160, 160, 160, 160),JOBNM
C
C INPUT OF BASIS FUNCTIONS FOR JOBS 1, 5, 6, 7, 8, AND 9.
  160 PWR(I,1)=ONE
      IF(NUMP-1) 190, 190, 170
  170 DO 180  J=2,NUMP
        PWR(I,J)=(T(I))**(J-1)
  180   CONTINUE
  190 PWR(I,NMPP1)=ONE
      IF(NUMQ-1) 220, 220, 200
  200 DO 210  J=2,NUMQ
        JJ=NUMP+J
        PWR(I,JJ)=(T(I))**(J-1)
  210   CONTINUE
  220 GO TO ( 310, 310, 310, 310, 460, 460, 490, 533, 310),JOBNM
C
C INPUT OF BASIS FUNCTIONS FOR JOBS 2 AND 3.
  230 PWR(I,1)=ONE
      IF(NUMP-1) 260, 260, 240
  240 DO 250  J=2,NUMP
        PWR(I,J)=COS((J-1)*T(I))
  250   CONTINUE
  260 PWR(I,NMPP1)=ONE
      IF(NUMQ-1) 290, 290, 270
  270 DO 280  J=2,NUMQ
        JJ=NUMP+J
        PWR(I,JJ)=COS((J-1)*T(I))
  280   CONTINUE
  290 CONTINUE
      GO TO 380
C
C INPUT OF BASIS FUNCTIONS FOR JOB 4.
  300 PWR(I,1)=ONE
      PWR(I,2)=T(IDBM)**2
      PWR(I,3)=T(IDB)**2
      PWR(I,4)=ONE
      PWR(I,5)=PWR(I,2)*PWR(I,3)
      GO TO 430
C
C INPUT OF F, INDF, DENOMINATOR LOWER BOUND, AND LOWER CONSTRAINTS
C FOR JOBS 1 AND 9, AND ADDITIONAL EQUALITY CONSTRAINT FOR JOB 9.
  310 NUM1=3*(NUMGR-1)/20+1
      NUM2=(NUMGR-1)/2+1
      NUM3=3*(NUMGR-1)/5+1
C IF I=1 SET THE DENOMINATOR LOWER BOUND FOR JOBS 1 AND 9 AND THE
C ADDITIONAL EQUALITY CONSTRAINT FOR JOB 9.  NOTE THAT THESE THINGS
C NEED BE SET ONLY ONCE IN ANY ONE JOB.
      IF(I-1)311,311,317
  311 IF(JOBNM-1)312,312,315
  312 EPSQ=SPCMN*TEN**3
      GO TO 317
  315 EPSQ=ONE/TWO
      ALC(1,1)=ONE
      ALC(1,2)=ZERO
      ALC(1,3)=-ONE/TWO
      ALC(1,4)=-ONE/TWO
      ALC(1,5)=ZERO
      ALC(1,6)=ONE
C SET RESTRICTED RANGE CONSTRAINTS FOR JOBS 1 AND 9.
  317 IF(I-NUM1) 318, 318, 330
  318 IF(JOBNM-1)320,320,340
  320 INDF(I)=1
      FTBLE(I)=ONE
      INDLW(I)=1
      ALTBL(I)=ONE-TEN**(-4)
      GO TO 140
  330 IF(I-NUM2) 340, 340, 350
  340 INDF(I)=1
      FTBLE(I)=ONE
      INDLW(I)=0
      GO TO 140
  350 IF(I-NUM3) 360, 370, 370
  360 INDF(I)=0
      INDLW(I)=1
      ALTBL(I)=ZERO
      GO TO 140
  370 INDF(I)=1
      FTBLE(I)=ZERO
      INDLW(I)=1
      ALTBL(I)=ZERO
      GO TO 140
C
C INPUT OF F, INDF, DENOMINATOR LOWER BOUND, WEIGHT FUNCTION, AND
C LOWER CONSTRAINTS FOR JOBS 2 AND 3.
  380 NUM1=(NUMGR-1)/2+1
      NUM2=3*(NUMGR-1)/5+1
      EPSQ=SPCMN*TEN**3
      IF(I-NUM1) 390, 390, 400
  390 INDF(I)=1
      FTBLE(I)=ONE
      WTBLE(I)=ONE/TWO
      INDLW(I)=0
      GO TO 140
  400 IF(I-NUM2) 410, 420, 420
  410 INDF(I)=0
      INDLW(I)=1
      ALTBL(I)=ZERO
      GO TO 140
  420 INDF(I)=1
      FTBLE(I)=ZERO
      WTBLE(I)=ONE
      INDLW(I)=1
      ALTBL(I)=ZERO
      GO TO 140
C
C INPUT OF F, INDF, AND DENOMINATOR LOWER BOUND FOR JOB 4.
  430 EPSQ=SPCMN*TEN**3
C WE USE 16.0001 INSTEAD OF 16.0 IN THE NEXT TEST TO AVOID TROUBLE
C DUE TO ROUND OFF ERROR.
      IF(T(IDBM)**2+T(IDB)**2-(SIXTN+TEN**(-4))) 440, 440, 450
  440 INDF(I)=1
      IF(SIXTN-T(IDBM)**2-T(IDB)**2)447,447,443
  443 FTBLE(I)=SQRT(SIXTN-T(IDBM)**2-T(IDB)**2)
      GO TO 140
  447 FTBLE(I)=ZERO
      GO TO 140
  450 INDF(I)=0
      GO TO 140
C
C INPUT OF F FOR JOBS 5 AND 6.  WE ALSO INPUT EPSQ HERE, ALTHOUGH
C IT IS NOT ACTUALLY USED IN JOB 5.
  460 IF(I-1)470,470,480
  470 FTBLE(I)=ONE
      EPSQ=TEN**(-7)
      GO TO 140
  480 FTBLE(I)=ZERO
      GO TO 140
C
C INPUT OF F FOR JOB 7.
  490 NUM1=(NUMGR-1)/10+1
      NUM2=(NUMGR-1)/2+1
      IF(I-NUM1) 500, 500, 510
  500 FTBLE(I)=THREE/(ONE+T(I))-TWO*T(I)+TWO/TEN
      GO TO 140
  510 IF(I-NUM2) 520, 520, 530
  520 FTBLE(I)=THREE/(ONE+T(I))+T(I)/TWO-THREE/TEN
      GO TO 140
  530 FTBLE(I)=THREE/(ONE+T(I))-FOUR*T(I)/TEN+(THREE+THREE)/TEN
      GO TO 140
C
C INPUT OF F AND EQUALITY CONSTRAINTS FOR JOB 8.  WE NEED SET
C ALC ONLY ONCE.
  533 FTBLE(I)=EXP(T(I))
      IF(I-1)535,535,140
  535 DO 537 J=1,NUMP
        ALC(1,J)=ZERO
        ALC(2,J)=ONE
        ALC(3,J)=ZERO
  537   CONTINUE
      ALC(1,NMPP1)=ONE
      ALC(1,NMPP1+1)=ONE
      ALC(2,NMPP1)=ZERO
      ALC(2,NMPP1+1)=EXP(ONE)
      ALC(3,NUMP)=-TWO
      ALC(3,NMPP1)=ZERO
      ALC(3,NMPP1+1)=-THREE
C END COMPUTING OF BASIS FUNCTIONS AND F.
      GO TO 140
C
  540 CALL CDFCOR(NUMP,NUMQ,NUMGR,FTBLE,PWR,IOPT,EPSQ,INDLW,ALTBL,
     *INDUP,UPTBL,INDF,WTBLE,MEQ,MINEQ,ALC,MEQA,X,ERROR,QMIN,IQMIN,
     *ITER,INDLP,IFLAG)
C CHECK THE OUTPUT.  FIRST MAKE SURE IFLAG IS NONPOSITIVE.
      IF(IFLAG)570,570,550
  550 WRITE(NWRIT,560)
  560 FORMAT(/42H THE PROGRAM FAILED IN THE INITIALIZATION,,
     *28H SO WE GIVE A FULL PRINTOUT.)
      GO TO 620
C CHECK TO SEE IF THE COMPUTED ERROR NORM IS CORRECT WITHIN ERTOL.
  570 IF(ABS(ERROR(NUMGR+1)-ERTST(JOBNM))-ERTOL)580,580,600
  580 WRITE(NWRIT,590)
  590 FORMAT(/33H THE PROGRAM HAS WORKED NORMALLY.)
C
C**********MODIFICATION  IF A FULL PRINTOUT IS DESIRED EVEN IF THE
C PROGRAM SUCCEEDS, THE NEXT STATEMENT SHOULD BE GO TO 620
C OTHERWISE THE NEXT STATEMENT SHOULD BE GO TO 790
      GO TO 620
C**********END MODIFICATION
C
  600 WRITE(NWRIT,610)ERTOL
  610 FORMAT(/36H THE ERROR NORM WAS OFF BY MORE THAN,E12.2,
     *30H   SO WE GIVE A FULL PRINTOUT.)
C
  620 WRITE(NWRIT, 630)NUMDM
  630 FORMAT(/33H WE ARE DEALING WITH FUNCTIONS OF,I4,
     *12H   VARIABLES)
C
      WRITE(NWRIT, 640)NUMP,NUMQ,NUMGR,NRWGR,NCLGR
  640 FORMAT(/6H P HAS,I4,24H   COEFFICIENTS    Q HAS,I4,
     *15H   COEFFICIENTS//13H THE GRID HAS,I5,
     *18H   POINTS ARRANGED,I5,5H   BY,I5)
C
      WRITE(NWRIT, 650)IOPT
  650 FORMAT(/8H IOPT IS,I10)
C
C IF THE ONES DIGIT OF IOPT EXCEEDS 1, WRITE THE NUMBERS OF
C EQUALITY AND EXTRA INEQUALITY CONSTRAINTS, AND THE NUMBER OF
C REDUNDANT EQUALITY CONSTRAINTS.
      IF(IONE-1)657,657,653
  653 MEQDF=MEQ-MEQA
      WRITE(NWRIT,655)MEQ,MEQDF,MINEQ
  655 FORMAT(/10H THERE ARE,I5,31H  EQUALITY CONSTRAINTS OF WHICH,
     *I5,15H  ARE REDUNDANT//10H THERE ARE,I5,
     *30H  EXTRA INEQUALITY CONSTRAINTS)
C
  657 IF(NUMDM-2) 660, 680, 660
  660 WRITE(NWRIT, 670)ASWP1,ASWP2
  670 FORMAT(/22H THE LEFT END POINT IS,E25.15,
     *//23H THE RIGHT END POINT IS,E25.15)
      GO TO 700
  680 WRITE(NWRIT, 690)ASWP1,ASWP2,ANEP1,ANEP2
  690 FORMAT(/40H THE COORDINATES OF THE SW AND NE POINTS,
     *4H ARE/2E25.15,6H   AND/2E25.15)
C
C WRITE IFLAG AND INDLP.
  700 WRITE(NWRIT, 710)IFLAG,INDLP
  710 FORMAT(/9H IFLAG IS,I5,13H     INDLP IS,I4)
C
C WRITE THE DENOMINATOR MINIMUM AND THE INDEX OF THE POINT
C WHERE IT OCCURRED.
      WRITE(NWRIT,720)QMIN,IQMIN
  720 FORMAT(/40H THE MINIMUM VALUE OF THE DENOMINATOR IS,
     *E25.15/23H ACHIEVED AT GRID POINT,I5)
C
C WRITE ITER AND DELTK.
      DELTK=ERROR(NUMGR+1)
      WRITE(NWRIT, 730)ITER,DELTK
  730 FORMAT(/21H THE ERROR NORM AFTER,I5,13H   ITERATIONS,
     *3H IS,E25.15)
C
C WRITE THE COEFFICIENTS OF THE NUMERATOR AND DENOMINATOR.
      WRITE(NWRIT, 740)(X(I),I=1,NUMP)
  740 FORMAT(/26H THE COEFFICIENTS OF P ARE/(/3E25.15))
C
      WRITE(NWRIT, 750)(X(I),I=NMPP1,NMPPQ)
  750 FORMAT(/26H THE COEFFICIENTS OF Q ARE/(/3E25.15))
C
C WRITE THE ERRORS.
      WRITE(NWRIT, 760)
  760 FORMAT(/15H THE ERRORS ARE)
      DO 780  I=1,NRWGR
        IND1=1+(I-1)*NCLGR
        IND2=I*NCLGR
        WRITE(NWRIT, 770)(ERROR(IT),IT=IND1,IND2)
  770   FORMAT(/4E20.10)
  780   CONTINUE
C
      JOBNM=JOBNM+1
      GO TO 30
      END
      SUBROUTINE CDFCOR(NUMP,NUMQ,NUMGR,FTBLE,PWR,IOPT,EPSQ,INDLW,
     *ALTBL,INDUP,UPTBL,INDF,WTBLE,MEQ,MINEQ,ALC,MEQA,X,ERROR,QMIN,
     *IQMIN,ITER,INDLP,IFLAG)
C***BEGIN PROLOGUE  CDFCOR
C***DATE WRITTEN   720120   (YYMMDD)
C***REVISION DATE  880226   (YYMMDD)
C***CATEGORY NO.  E2C1
C***KEYWORDS  APPROXIMATION,CHEBYSHEV,RATIONAL,UNIFORM,BEST,
C             GENERALIZED,DIFFERENTIAL CORRECTION,MULTIVARIATE,
C             CONSTRAINTS,RESTRICTED RANGE,WEIGHTED
C***AUTHOR  KAUFMAN, EDWIN H. JR.,
C             DEPARTMENT OF MATHEMATICS
C             CENTRAL MICHIGAN UNIVERSITY
C             MOUNT PLEASANT, MICHIGAN 48859
C           TAYLOR, GERALD D.,
C             DEPARTMENT OF MATHEMATICS
C             COLORADO STATE UNIVERSITY
C             FORT COLLINS, COLORADO 80523
C***PURPOSE  THIS SUBROUTINE COMPUTES A BEST UNIFORM GENERALIZED
C            RATIONAL APPROXIMATION P/Q TO A FUNCTION F DEFINED
C            ON A FINITE SET OF GRID POINTS, WITH THE POSSIBILITY
C            OF INCLUDING A WEIGHT FUNCTION AND CONSTRAINTS ON
C            P/Q AND Q, AND ADDITIONAL ARBITRARY LINEAR CONSTRAINTS
C            ON THE COEFFICIENTS OF P AND Q.
C***DESCRIPTION
C
C THE VARIABLES IN THE CALLING SEQUENCE ARE AS FOLLOWS.  NONE OF
C THE STRICTLY INPUT VARIABLES IS CHANGED BY THIS PROGRAM.
C
C  NUMP   (INPUT) THIS IS THE NUMBER OF BASIS FUNCTIONS FOR THE
C         NUMERATOR.
C
C  NUMQ   (INPUT) THIS IS THE NUMBER OF BASIS FUNCTIONS FOR THE
C         DENOMINATOR.
C
C  NUMGR  (INPUT) THIS IS THE NUMBER OF GRID POINTS.
C
C  FTBLE  (INPUT) FTBLE(I) IS THE VALUE OF THE FUNCTION F AT GRID
C         POINT I (I=1,...,NUMGR).
C
C  PWR    (INPUT) PWR(I,J) IS THE VALUE OF THE JTH BASIS FUNCTION
C         (WHERE THE NUMERATOR BASIS FUNCTIONS PRECEDE THE
C         DENOMINATOR BASIS FUNCTIONS) AT GRID POINT I (J=1,...,
C         NUMP+NUMQ, I=1,...,NUMGR).
C
C  IOPT   (INPUT)  THIS IS THE OPTION SWITCH.  SETTING IOPT=0 TURNS
C         OFF ALL EXTRA OPTIONS.  THE PROPER SETTING OF IOPT TO
C         ACTIVATE ONE OR MORE OF THE EXTRA OPTIONS IS DISCUSSED IN
C         THE LONG DESCRIPTION BELOW.
C
C  EPSQ, INDLW, ALTBL, INDUP, UPTBL, INDF, WTBLE, MEQ, MINEQ, ALC,
C         MEQA  (OPTIONAL INPUT OR OUTPUT)
C         THESE ELEVEN VARIABLES ARE USED FOR INPUT OR OUTPUT FOR THE
C         EXTRA OPTIONS DESCRIBED IN THE LONG DESCRIPTION BELOW.  NONE
C         OF THEM NEED BE ASSIGNED VALUES IF IOPT=0.
C
C  X      (OUTPUT AND OPTIONAL INPUT) THE PROGRAM WILL PLACE THE
C         COMPUTED COEFFICIENT OF THE JTH BASIS FUNCTION IN X(J)
C         (J=1,...,NUMP+NUMQ).  X IS ALSO USED IF THE USER WISHES
C         TO SUPPLY AN INITIAL APPROXIMATION, IN WHICH CASE IOPT
C         MUST BE SET APPROPRIATELY (SEE LONG DESCRIPTION BELOW).
C
C  ERROR  (OUTPUT) THE PROGRAM WILL PLACE THE ERROR AT GRID POINT
C         I, WHICH WILL BE (F-P/Q)(GRID POINT I) IF NEITHER WEIGHT
C         FUNCTION OPTION IS USED (SEE LONG DESCRIPTION BELOW FOR
C         THE WEIGHT FUNCTION OPTIONS), IN ERROR(I) (I=1,...,NUMGR).
C         THE UNIFORM ERROR NORM WILL BE PLACED IN ERROR(NUMGR+1).
C
C  QMIN   (OUTPUT) THE PROGRAM WILL PLACE THE MINIMUM VALUE OF THE
C         DENOMINATOR ON THE GRID IN QMIN.
C
C  IQMIN  (OUTPUT) THE PROGRAM WILL PLACE THE INDEX OF THE GRID POINT
C         AT WHICH THE MINIMUM VALUE OF THE DENOMINATOR IS ACHIEVED
C         IN IQMIN.
C
C  ITER   (OUTPUT) THE PROGRAM WILL PLACE THE NUMBER OF ITERATIONS
C         (NOT COUNTING INITIALIZATION) REQUIRED IN ITER.
C
C  INDLP  (OUTPUT) THE PROGRAM WILL PLACE INFORMATION ABOUT THE
C         PERFORMANCE OF THE LINEAR PROGRAMMING SUBROUTINE SLNPRO
C         IN THE TWO DIGIT VARIABLE INDLP (SEE COMMENTS AT THE END OF
C         THIS SUBROUTINE AND THE BEGINNING OF SUBROUTINE SLNPRO).
C         THE USER NORMALLY NEED NOT CONSIDER INDLP UNLESS IFLAG (SEE
C         BELOW) OR SOMETHING ELSE INDICATES SOMETHING IS WRONG.  IN
C         THAT CASE A CAREFUL STUDY OF THE INDICATED COMMENTS MAY
C         PROVIDE A CLUE AS TO WHAT HAPPENED.  FOR EXAMPLE, INDLP=44
C         USUALLY INDICATES INCORRECT INPUT TO THE PROGRAM.
C
C  IFLAG  (OUTPUT) IFLAG IS THE ERROR FLAG FOR THE PROGRAM.  IFLAG=0
C         MEANS THE PROGRAM APPEARED TO WORK NORMALLY, WITH NONE OF
C         THE CONDITIONS DESCRIBED BELOW FOR IFLAG .NE. O OCCURRING.
C         IFLAG POSITIVE IS A WARNING OF FAILURE IN THE INITIALIZATION.
C         IFLAG NEGATIVE INDICATES CAUTION SHOULD BE USED, ALTHOUGH
C         THE RETURNED APPROXIMATION MAY WELL BE ACCEPTABLE.  FOR
C         DETAILS CONCERNING SPECIFIC NONZERO VALUES OF IFLAG, SEE
C         THE LONG DESCRIPTION BELOW.
C
C THE DIMENSIONS OF THE ARRAYS IN THE CALLING SEQUENCE FOR
C CDFCOR ARE
C
C FTBLE(101), PWR(101,15), INDLW(101), ALTBL(101), INDUP(101),
C UPTBL(101), INDF(101), WTBLE(101), X(16), ERROR(102), ALC(101,16).
C
C IN ADDITION, IF SUBROUTINE STCOMP IS TO BE USED,
C DIMENSION T(202) MUST BE ADDED.
C
C THESE DIMENSIONS ARE SET SO THAT THE PROGRAM CAN BE APPLIED TO
C ANY PROBLEM WITH UP TO 15 BASIS FUNCTIONS (SO NUMP+NUMQ .LE. 15),
C UP TO 101 GRID POINTS (SO NUMGR .LE. 101), AND UP TO 101
C ADDITIONAL LINEAR CONSTRAINTS FROM OPTION 4 BELOW (SO
C MEQ+MINEQ .LE. 101), EVEN IF ALL THE EXTRA OPTIONS DESCRIBED
C IN THE LONG DESCRIPTION BELOW ARE ACTIVATED.  SEE THE LONG
C DESCRIPTION FOR THE DIMENSION PATTERN REQUIRED IN THE DRIVER
C PROGRAM AND SUBROUTINES IF MORE BASIS FUNCTIONS, GRID POINTS,
C OR ADDITIONAL LINEAR CONSTRAINTS ARE DESIRED, AND FOR THE
C DIMENSIONS OF THE INTERNAL WORK ARRAYS.
C
C***LONG DESCRIPTION
C
C THE USER CAN, BY SETTING IOPT AND THE CORRESPONDING OPTIONAL
C VARIABLES, USE ANY COMBINATION OF THE SEVEN EXTRA OPTIONS
C DESCRIBED BELOW.  AN OPTIONAL VARIABLE NEED BE ASSIGNED A VALUE
C BY THE USER ONLY IF ITS OPTION IS TO BE USED.  SETTING IOPT=0
C TURNS OFF ALL THE EXTRA OPTIONS.
C
C  OPTION 1  CONSTRAINTS (DENOMINATOR)
C
C            TO ACTIVATE THIS OPTION, SET THE ONES DIGIT OT IOPT
C            TO 1 (IF OPTION 4 IS NOT USED) OR TO 3 (IF OPTION 4 IS
C            USED) AND SET EPSQ AS FOLLOWS.
C
C  EPSQ   (INPUT) THIS IS THE CONSTANT LOWER BOUND ON THE VALUES
C         OF THE DENOMINATOR. THUS THE PROGRAM WILL SET UP
C         EXTRA CONSTRAINTS IN SUBROUTINE SLNPRO TO FORCE
C         Q .GE. EPSQ EVERYWHERE ON THE GRID.  THIS OPTION IS
C         USEFUL BECAUSE THE COMPUTED APPROXIMATION P/Q MAY BE OF
C         MORE VALUE IF ITS DENOMINATOR IS KEPT FARTHER AWAY FROM
C         ZERO, AND ALSO BECAUSE IN DIFFICULT CASES THIS OPTION
C         MAY ALLEVIATE NUMERICAL DIFFICULTIES IN THE PROGRAM.
C
C  OPTION 2  CONSTRAINTS (LOWER RESTRICTED RANGE)
C
C            TO ACTIVATE THIS OPTION, SET THE TENS DIGIT OF IOPT
C            TO 1 AND SET INDLW AND ALTBL AS FOLLOWS.
C
C  INDLW  (INPUT) INDLW(I)=1 IF THERE IS A LOWER RESTRICTED
C         RANGE CONSTRAINT AT GRID POINT I, AND INDLW(I)=0
C         OTHERWISE (I=1,...,NUMGR).
C
C  ALTBL  (INPUT) ALTBL(I) CONTAINS THE LOWER BOUND AT GRID POINT I
C         (THUS WE ARE FORCING (P/Q)(GRID POINT I) .GE. ALTBL(I))
C         IF THERE IS ONE, AND ALTBL(I) NEED NOT BE ASSIGNED A
C         VALUE OTHERWISE (I=1,...,NUMGR).
C
C  OPTION 3  CONSTRAINTS (UPPER RESTRICTED RANGE)
C
C            TO ACTIVATE THIS OPTION, SET THE HUNDREDS DIGIT OF
C            IOPT TO 1 AND SET INDUP AND UPTBL AS FOLLOWS.
C
C  INDUP  (INPUT) INDUP(I)=1 IF THERE IS AN UPPER RESTRICTED
C         RANGE CONSTRAINT AT GRID POINT I, AND INDUP(I)=0
C         OTHERWISE (I=1,...,NUMGR).
C
C  UPTBL  (INPUT) UPTBL(I) CONTAINS THE UPPER BOUND AT GRID POINT I
C         (THUS WE ARE FORCING (P/Q)(GRID POINT I) .LE. UPTBL(I))
C         IF THERE IS ONE, AND UPTBL(I) NEED NOT BE ASSIGNED A VALUE
C         OTHERWISE (I=1,...,NUMGR).
C
C  OPTION 4  CONSTRAINTS (ADDITIONAL LINEAR CONSTRAINTS)
C
C            TO ACTIVATE THIS OPTION, SET THE ONES DIGIT OF
C            IOPT TO 2 (IF OPTION 1 IS NOT USED) OR TO 3
C            (IF OPTION 1 IS USED) AND SET MEQ, MINEQ, AND
C            ALC AS FOLLOWS.
C
C  MEQ    (INPUT)  THIS IS THE NUMBER OF ADDITIONAL EQUALITY
C         CONSTRAINTS.
C
C  MINEQ  (INPUT)  THIS IS THE NUMBER OF ADDITIOINAL INEQUALITY
C         CONSTRAINTS.
C
C  ALC    (INPUT)  FOR I=1,...,MEQ, THE ITH ROW OF ALC SPECIFIES
C         THE EQUALITY CONSTRAINT
C           ALC(I,1)*X(1)+...+ALC(I,NUMP+NUMQ)*X(NUMP+NUMQ) =
C           ALC(I,NUMP+NUMQ+1)
C         AND FOR I=MEQ+1,...,MEQ+MINEQ, THE ITH ROW OF ALC
C         SPECIFIES THE INEQUALITY CONSTRAINT
C           ALC(I,1)*X(1)+...+ALC(I,NUMP+NUMQ)*X(NUMP+NUMQ) .LE.
C           ALC(I,NUMP+NUMQ+1).
C
C  OPTION 5  IGNORE F AT SPECIFIED POINTS
C
C            TO ACTIVATE THIS OPTION, SET THE THOUSANDS DIGIT OF
C            IOPT TO 1 AND SET INDF AS FOLLOWS.
C
C  INDF   (INPUT) INDF(I)=1 IF F IS TO BE APPROXIMATED AT GRID
C         POINT I, AND INDF(I)=0 IF F IS TO BE IGNORED AT GRID
C         POINT I (ALTHOUGH OPTIONS 1, 2, AND 3 CAN BE USED AT
C         SUCH POINTS) (I=1,...,NUMGR).  NEITHER FTBLE(I) NOR
C         WTBLE(I) NEED BE ASSIGNED A VALUE IF INDF(I)=0.
C
C            PROGRAMMING HINT... IF OPTION 5 IS USED WITH A NONTRIVIAL
C            DENOMINATOR, IT IS A GOOD IDEA TO ALSO USE OPTION 1 OR 4
C            TO IMPOSE A POSITIVE LOWER BOUND ON THE DENOMINATOR AT
C            LEAST AT THE POINTS WHERE F IS IGNORED, SINCE THE
C            DIFFERENTIAL CORRECTION INQUALITIES ARE NOT PRESENT TO
C            CONTROL THE DENOMINATOR AT SUCH POINTS, AND A DENOMINATOR
C            VALUE WHICH IS SMALL IN ABSOLUTE VALUE OR NEGATIVE COULD
C            TERMINATE THE PROGRAM BEFORE A BEST APPROXIMATION IS FOUND.
C            A REASONABLE LOWER BOUND IS 1000.0*B**(-T) (WHERE T IS THE
C            NUMBER OF BASE B DIGITS IN THE MANTISSA).  NOTE THAT AT
C            POINTS WHERE F IS NOT IGNORED, IF
C            (ABS(F*Q-WT*P)-DELTK*Q)/QK .LE. W .LT. 0.0, THEN
C            Q .GE. (ABS(F*Q-WT*P)+ABS(W)*QK)/DELTK .GT. 0.0
C
C  OPTION 6  WEIGHTED APPROXIMATION
C
C            TO ACTIVATE STRAIGHT WEIGHTED APPROXIMATION, THAT IS,
C            TO MINIMIZE MAX(ABS(F-WT*P/Q)), SET THE TEN THOUSANDS
C            DIGIT OF IOPT TO 1 AND SET WTBLE AS BELOW.
C
C            TO ACTIVATE WEIGHTED ERROR CURVE APPROXIMATION, THAT IS,
C            TO MINIMIZE MAX(ABS(WT*(F-P/Q))), SET THE TEN THOUSANDS
C            DIGIT OF IOPT TO 2 AND SET WTBLE AS BELOW.
C
C  WTBLE  (INPUT) WTBLE(I) IS THE VALUE OF THE WEIGHT FUNCTION WT
C         (WHICH NEED NOT BE POSITIVE) AT GRID POINT I (I=1,...,NUMGR).
C
C  OPTION 7  USER SUPPLIED INITIAL APPROXIMATION
C
C            TO ACTIVATE THIS OPTION, SET THE HUNDRED THOUSANDS DIGIT
C            OF IOPT TO 1 AND PLACE THE NUMP+NUMQ INITIAL
C            COEFFICIENTS IN X (WHERE X IS DEFINED IN THE
C            DESCRIPTION SECTION ABOVE).  SOMETIMES ROUND OFF ERROR
C            CAN BE ALLEVIATED, AND COMPUTER TIME CAN BE SAVED, BY
C            FIRST RUNNING A PROBLEM WITH A COARSE GRID AND THEN
C            USING THE RESULT TO INITIALIZE THE PROBLEM ON THE
C            DESIRED GRID.
C
C THE MEANINGS OF SPECIFIC NONZERO VALUES OF IFLAG ARE AS FOLLOWS.
C
C   IFLAG POSITIVE (INITIALIZATION FAILURE) RARELY OCCURS IF THE
C     DATA FOR THE PROGRAM IS ENTERED CORRECTLY.  THERE ARE THREE
C     POSSIBILITIES.
C
C       IFLAG=7777 MEANS OPTION 4 ABOVE WAS USED, AND THE
C         LINEAR INEQUALITY CONSTRAINTS APPEARED TO BE
C         INCONSISTENT.
C
C       IFLAG=3333 MEANS THE APPROXIMATION RETURNED TO THE USER IS
C         THE INITIAL APPROXIMATION P/Q=0.0/PWR(.,NUMP+1) (THAT IS,
C         ALL COEFFICIENTS ARE 0.0 EXCEPT THE FIRST DENOMINATOR
C         COEFFICIENT, WHICH IS 1.0).  THIS APPROXIMATION IS LIKELY
C         TO BE OF LITTLE VALUE.
C
C       IFLAG=6666 MEANS EVEN THE APPROXIMATION P/Q=0.0/PWR(.,NUMP+1)
C         WAS REJECTED BY SUBROUTINE SPQERD BECAUSE ITS DENOMINATOR
C         WAS VERY SMALL IN ABSOLUTE VALUE OR NEGATIVE AT SOME GRID
C         POINT, OR ELSE THE APPROXIMATION VIOLATED SOME CONSTRAINT
C         BY MORE THAN EPRES = 10000.0*B**(-T).  NO APPROXIMATION
C         WAS RETURNED, ALL THE COEFFICIENTS OF P AND Q WERE SET
C         TO 0.0, AND ITER WAS SET TO -1.
C
C   IN CASE IFLAG IS NEGATIVE (INDICATING CAUTION), IT MAY BE
C     ADVISABLE TO EXAMINE THE ERRORS AT THE GRID POINTS, AS
C     NORMALLY THE NUMBER OF POINTS WHERE THE ERROR NORM IS
C     ACHIEVED OR A USER SUPPLIED CONSTRAINT IS SATISFIED WITH
C     EQUALITY, PLUS THE NUMBER OF DENOMINATOR COEFFICIENTS
C     WITH ABSOLUTE VALUE 1.0, WILL BE AT LEAST NUMP+NUMQ+1.
C     THE SPECIFIC POSSIBILITIES FOR IFLAG NEGATIVE, MORE THAN
C     ONE OF WHICH CAN OCCUR, ARE AS FOLLOWS.
C
C       IF THE ONES DIGIT OF IFLAG IS 1, THE INITIAL APPROXIMATION
C         WAS NOT IMPROVED, AND WAS RETURNED TO THE USER.  THIS
C         FREQUENTLY HAPPENS NATURALLY (AND WITH NO ILL
C         CONSEQUENCES) WHEN ONE IS DOING GENERALIZED POLYNOMIAL
C         APPROXIMATION, THAT IS, WHEN ONE HAS SET NUMQ=1 AND
C         PWR(I,NUMP+1)=1.0 FOR I=1,...,NUMGR.
C
C       IF THE ONES DIGIT OF IFLAG IS 2, THE PROGRAM WAS TERMINATED
C         BECAUSE ITLIM ITERATIONS WERE PERFORMED, WHERE CURRENTLY
C         ITLIM=100.  ONE MAY WISH TO RESTART THE PROGRAM,
C         INITIALIZING WITH THE CURRENT APPROXIMATION, TO SEE IF THIS
C         APPROXIMATION CAN BE IMPROVED.
C
C       IF THE TENS DIGIT OF IFLAG IS 1, THEN DURING THE COMPUTATION
C         OF THE RETURNED APPROXIMATION IT WAS NECESSARY TO ADJUST
C         THE TOLERANCES IN SUBROUTINE SLNPRO TO AVOID FAILURE.
C         THUS THE DANGER THAT THE RETURNED APPROXIMATION HAS BEEN
C         SERIOUSLY AFFECTED BY ROUND OFF ERROR IS GREATER.
C
C       IF THE HUNDREDS DIGIT OF IFLAG IS 1, A USER SUPPLIED
C         CONSTRAINT WAS VIOLATED BY MORE THAN EPRES, WHERE
C         EPRES=10000.0*B**(-T).  (B IS THE BASE FOR FLOATING POINT
C         NUMBERS FOR THE MACHINE AND T IS THE NUMBER OF BASE B
C         DIGITS IN THE MANTISSA.)
C
C IF IT IS DESIRED TO CHANGE FROM DOUBLE PRECISION TO SINGLE PRECISION
C (FOR EXAMPLE), ONLY THREE THINGS NEED TO BE DONE.
C
C    (1)  REPLACE ALL THE STATEMENTS  ONE=1.0D0  BY  ONE=1.0
C
C    (2)  REPLACE ALL THE STATEMENTS  IMPLICIT REAL*8 (A-H,O-Z)  BY
C         C       IMPLICIT REAL*8 (A-H,O-Z)
C
C    (3  REPLACE ALL OCCURRENCES OF  D1MACH  BY  R1MACH
C        (OR MODIFY THIS PROGRAM OR THE SUBPROGRAM D1MACH SUPPLIED BY
C        NETLIB SO THAT WHEREVER  THE STATEMENT  SPCMN=D1MACH(3)
C        APPEARS IN THE CURRENT VERSION, SPCMN WILL BE COMPUTED AS
C        B**(-T), WHERE B IS THE BASE FOR FLOATING POINT NUMBERS AND
C        T IS THE NUMBER OF BASE B DIGITS IN THE MANTISSA)
C        (NOTE THAT THE ONLY VALUES WHICH THIS VERSION OF THIS PACKAGE
C        USES FROM THE NETLIB SUPPLIED SUBPROGRAMS I1MACH, R1MACH, AND
C        D1MACH ARE D1MACH(3), I1MACH(1), AND I1MACH(2), WITH THE LATTER
C        TWO VALUES USED ONLY IN THE SAMPLE DRIVER PROGRAM.)
C
C THE PATTERN FOR DIMENSION DECLARATIONS WHICH WILL SET ASIDE
C SUFFICIENT SPACE TO SOLVE ANY PROBLEM WITH UP TO MFN BASIS
C FUNCTIONS (THUS NUMP+NUMQ .LE. MFN), MGR GRID POINTS (THUS
C NUMGR .LE. MGR), AND MALC ADDITIONAL LINEAR CONSTRAINTS IN
C OPTION 4 ABOVE (THUS MEQ+MINEQ .LE. MALC) IS
C
C FTBLE(MGR), PWR(MGR,MFN), INDLW(MGR), ALTBL(MGR), INDUP(MGR),
C UPTBL(MGR), INDF(MGR), WTBLE(MGR), X(MFN+1), ERROR(MGR+1),
C T(2*MGR), ALC(MALC,MFN+1), XTEMP(MFN+1), FTBWT(MGR), IYRCT(MAXCN),
C IYSAV(MAXCN), V(MAXCN+1,MFN+2), IYCCT(MFN+1), IXRCT(MAXCN), Y(MAXCN),
C EQ(MFN,MFN+1), NSCOL(MFN), IRESL(2*MFN),
C
C WHERE MAXCN=5*MGR+2*MFN+MALC-2.  NOTE THAT THE LAST 11 ARRAYS LISTED
C ARE INTERNAL WORK ARRAYS AND NEED NOT BE DIMENSIONED IN THE
C USERS DRIVER PROGRAM.  WE NOTE THAT FOR EACH OF OPTIONS 2 OR 3
C WHICH IS NOT USED, MAXCN COULD BE REDUCED BY MGR, AND MAXCN COULD BE
C REDUCED BY MALC IF OPTION 4 IS NOT USED.
C
C THE METHOD USED BY THIS PROGRAM IS THE DIFFERENTIAL CORRECTION
C ALGORITHM, WHICH REQUIRES AT EACH ITERATION THE SOLUTION OF THE
C LINEAR PROGRAMMING PROBLEM
C     MINIMIZE MAX(ABS(F*Q-WT*P)-DELTK*Q)/QK,
C     SUBJECT TO ABS(Q(J)) .LE. 1.0 FOR J=1,...,NUMQ,
C       (AND THE OTHER INEQUALITY CONSTRAINTS)
C WHERE QK AND DELTK ARE THE DENOMINATOR AND THE UNIFORM ERROR
C NORM FOR THE PREVIOUS APPROXIMATION, WT IS THE WEIGHT FUNCTION
C (NOTE THAT WT WILL BE 1.0 FOR UNWEIGHTED APPROXIMATION, AND F
C WILL HAVE BEEN MULTIPLIED BY WT FOR WEIGHTED ERROR CURVE
C APPROXIMATION), AND Q(J) IS THE JTH COEFFICIENT OF THE NEW Q.
C HERE MAX MEANS THE MAXIMUM COMPUTED OVER THE GRID.
C THE EQUALITY CONSTRAINTS (IF ANY, SEE OPTION 4 ABOVE) ARE
C HANDLED BY SOLVING THEM FOR MEQA OF THE COEFFICIENTS IN
C TERMS OF THE REMAINING COEFFICIENTS (USING A TOTAL
C PIVOTING STRATEGY) AND IN EFFECT SUBSTITUTING TO REDUCE
C THE NUMBER OF VARIABLES IN THE REST OF THE PROBLEM BY MEQA.
C IN THEORY, THE ERROR NORMS OF THE APPROXIMATIONS PRODUCED BY
C THE ALGORITHM WILL CONVERGE TO THE INFIMUM OF ALL POSSIBLE
C ERROR NORMS.
C
C IN ORDER TO SPEED UP CONVERGENCE FAR FROM THE SOLUTION, FOR
C ITERATIONS BEYOND THE INITIALIZATION AND FIRST MAIN ITERATION
C THE DELTK IN THE ABOVE EXPRESSION IS REPLACED BY 0.5*DELTK
C UNTIL THE PROGRAM FAILS TO PRODUCE A BETTER APPROXIMATION,
C AFTER WHICH THE PROGRAM REVERTS TO ORDINARY DIFFERENTIAL
C CORRECTION FROM THEN ON.  IN THEORY THIS PROCESS, KNOWN AS
C PARAMETRIC BISECTION, WILL RESULT IN MORE THAN HALVING THE
C ERROR NORM AT EACH ITERATION IN WHICH THE ERROR NORM IS
C MORE THAN TWICE THE OPTIMAL ERROR NORM.  THIS IDEA IS DUE
C TO ALI SAZEGARI.
C
C IN THE EVENT OF FAILURE OF THE LINEAR PROGRAMMING SUBROUTINE
C SLNPRO TO PRODUCE AN APPROXIMATION P/Q BETTER THAN THE PREVIOUS
C ONE (IN THE ERROR NORM SENSE), THE PROGRAM TRIES TWICE MORE, ONCE
C WITH LOOSENED AND ONCE WITH TIGHTENED TOLERANCES IN SLNPRO.  AN
C APPROXIMATION PRODUCED BY SLNPRO IS TESTED IN SUBROUTINE SPQERD TO
C SEE IF ITS DENOMINATOR IS NEGATIVE OR VERY SMALL IN ABSOLUTE VALUE
C (THAT IS, LESS THAN 100.0*B**-(T)) AT ANY GRID POINT.  IF SO, THE
C APPROXIMATION IS REJECTED (EXCEPT IN THE INITIALIZATION, WHERE
C THE PROGRAM RESETS THE DENOMINATOR AT THE BAD POINTS AND
C ATTEMPTS TO CONTINUE, AND IF A BETTER APPROXIMATION IS NOT
C PRODUCED, THEN THIS APPROXIMATION IS REJECTED).
C TO SAVE TIME, THE TWO EXTRA LINEAR PROGRAMMING ATTEMPTS ARE
C NOT DONE IF THE FAILURE OCCURS DURING THE PARAMETRIC BISECTION
C PHASE SINCE IN THAT PHASE THE PARAMETRIC BISECTION ATTEMPT WAS
C THE PROBABLE CAUSE OF THE FAILURE, AND THE PROGRAM WILL HAVE
C ANOTHER CHANCE USING ORDINARY DIFFERENTIAL CORRECTION.
C
C THE INITIAL APPROXIMATION FOR THE ALGORITHM IS OBTAINED BY ONE
C OF THREE METHODS, WITH THE USER BEING ABLE TO SELECT THE FIRST
C OR (BY DEFAULT) THE SECOND.  IF A METHOD FAILS, THE PROGRAM
C WILL TRY THE NEXT METHOD IN THE SEQUENCE.  FAILURE OF A METHOD
C MEANS EITHER NO RATIONAL APPROXIMATION P/Q COULD BE PRODUCED,
C OR ONE WAS PRODUCED WHICH VIOLATED SOME CONSTRAINT BY MORE THAN
C EPRES, OR WHICH HAD A DENOMINATOR WHICH WAS VERY SMALL IN
C ABSOLUTE VALUE OR NEGATIVE AT SOME POINT (AND THE MAIN
C ALGORITHM WAS UNABLE TO CORRECT THIS PROBLEM).
C THE THREE METHODS ARE
C
C  METHOD 1  THE USER SUPPLIES AN INITIAL APPROXIMATION (SEE
C            OPTION 7 ABOVE).
C
C  METHOD 2  ONE ITERATION OF LOEBS ALGORITHM IS PERFORMED, THAT
C            IS, SUBROUTINE SLNPRO IS USED TO SOLVE THE PROBLEM
C                MINIMIZE MAX(ABS(F*Q-WT*P))
C                SUBJECT TO Q .GE. MAX(0.00001,10000.0*B**(-T))
C                EVERYWHERE ON THE GRID, ABS(Q(J)) .LE. 1.0
C                FOR J=1,...,NUMQ (AND THE OTHER
C                INEQUALITY CONSTRAINTS).
C            IF WE ARE PLACING A LOWER BOUND ON THE DENOMINATOR
C            IN THE MAIN ITERATIONS, THE LOWER BOUND HERE IS
C            INCREASED (IF NECESSARY) TO MATCH THAT BOUND.
C            IF SLNPRO FAILS TO PRODUCE AN APPROXIMATION P/Q WE
C            TRY WITH BOTH LOOSENED AND TIGHTENED TOLERANCES.
C
C  METHOD 3  THE INITIAL APPROXIMATION IS TAKEN TO BE P/Q=
C            0.0/PWR(.,NUMP+1).
C
C THE ALGORITHM IS TERMINATED WHEN A BETTER APPROXIMATION CANNOT
C BE PRODUCED OR ITLIM ITERATIONS HAVE BEEN PERFORMED, AND THE
C BEST APPROXIMATON FOUND SO FAR IS THEN RETURNED TO THE USER.
C IF IT WAS NECESSARY TO ADJUST THE TOLERANCES IN SLNPRO TO
C OBTAIN AN APPROXIMAITON P/Q, THEN VIOLATION OF A USER SUPPLIED
C CONSTRAINT BY MORE THAN 10000.0*B**(-T) WILL CAUSE THAT
C APPROXIMATION TO BE REJECTED SINCE THE ADJUSTED TOLERANCES
C MAY HAVE ALLOWED REDUCTION OF THE UNIFORM ERROR NORM AT THE
C EXPENSE OF A SERIOUS VIOLATION OF THE USER SUPPLIED CONSTRAINTS.
C
C***REFERENCES  CHENEY, E. W. AND LOEB, H. L., TWO NEW ALGORITHMS
C                 FOR RATIONAL APPROXIMATION, NUMER. MATH. 4,
C                 PP. 124-127, 1962.
C               BARRODALE, I., POWELL, M. J. D., AND ROBERTS,
C                 F. D. K., THE DIFFERENTIAL CORRECTION ALGORITHM
C                 FOR RATIONAL L-INFINITY APPROXIMATION, SIAM J.
C                 NUMER. ANAL. 9, PP. 493-504, 1972.
C               KAUFMAN, EDWIN H. JR., LEEMING, D. J., AND TAYLOR,
C                 G. D., UNIFORM RATIONAL APPROXIMATION BY
C                 DIFFERENTIAL CORRECTION AND REMES-DIFFERENTIAL
C                 CORRECTION, INT. J. FOR NUMER. METH. IN ENGRG.
C                 17, PP. 1273-1278, 1981.
C***REFER TO  (NONE)
C***ROUTINES CALLED  SLNPRO,SPQERD,SSETV,EQELM,SCRCH,EXPND,
C                    ICNCK,NRMLZ
C***END PROLOGUE  CDFCOR
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION FTBLE(101),PWR(101,15),V(635,17),WTBLE(101),
     *   UPTBL(101),ALTBL(101),INDUP(101),INDLW(101),
     *   ERROR(102),IYRCT(634),IYSAV(634),XTEMP(16),X(16),
     *   INDF(101),FTBWT(101),ALC(101,16),EQ(15,16),NSCOL(15),
     *   IRESL(30)
C
C SET MACHINE DEPENDENT PARAMETERS FOR SUBROUTINE CDFCOR.
C SET SPCMN=B**(-T), WHERE B IS THE BASE AND T IS THE NUMBER
C OF BASE B DIGITS IN THE MANTISSA.  SPCMN IS THE MINIMUM
C RELATIVE SPACING ABS((X1-X2)/X2) BETWEEN TWO SUCCESSIVE
C FLOATING POINT NUMBERS, SO IT IS THE SPACING BETWEEN TWO
C SUCCESSIVE FLOATING POINT NUMBERS IN THE CLOSED INTERVAL
C (0.1,1.0).  WE ALSO HAVE SPCMN=10.0**(-T*(LOG10)(B))=
C 10.0**(-TNMAN), WHERE TNMAN IS THE BASE 10 EQUIVALENT OF
C THE LENGTH OF THE MANTISSA.
C
C***FIRST EXECUTABLE STATEMENT  CDFCOR
      SPCMN=D1MACH(3)
C SET PRECISION DEPENDENT CONSTANTS FOR CDFCOR.
      ONE=1.0D0
      ZERO=ONE-ONE
      TWO=ONE+ONE
      FOUR=TWO+TWO
      TEN=FOUR+FOUR+TWO
      TENTH=ONE/TEN
C SET EPRES=10000.0*SPCMN.  THUS EPRES=10.0**(-(TNMAN-4)).
C AN APPROXIMATION WILL BE REJECTED IF A USER SUPPLIED
C CONSTRAINT IS VIOLATED BY MORE THAN EPRES AND IT WAS NECESSARY
C TO ADJUST THE TOLERANCES IN SUBROUTINE SLNPRO WHILE COMPUTING
C THE APPROXIMATION.
      EPRES=SPCMN*TEN**4
C SET TOLEL=100.0*SPCMN.  IN SUBROUTINE EQELM, NUMBERS WITH
C ABSOLUTE VALUE .LE. TOLEL WILL BE TREATED AS ZEROES.
C THIS WILL ALSO BE THE CASE IN SUBROUTINE NRMLZ.
      TOLEL=TEN*TEN*SPCMN
C END OF SETTING MACHINE AND PRECISION DEPENDENT PARAMETERS FOR CDFCOR.
C
      ITLIM=100
      ITER=0
      NMPP1=NUMP+1
      NMPPQ=NUMP+NUMQ
      N=NMPPQ+1
      NP1=N+1
      NLP=N
      NLP1=NP1
      IADJS=0
      INDIC=0
      IRET1=0
      ILBFL=0
      ILST1=0
      QMIN=ZERO
      IQMIN=0
      INITL=0
      IFLAG=0
      IPBIT=1
      IONE=IOPT-(IOPT/10)*10
C IF THE ONES DIGIT OF IOPT IS .LE. 1, SET MEQ AND MINEQ TO ZERO.
      IF(IONE-1)100,100,200
  100 MEQ=0
      MINEQ=0
  200 MEQA=0
C
C IF F IS NOT TO BE IGNORED AT ANY OF THE GRID POINTS,
C SET INDF IDENTICALLY EQUAL TO 1.
      IF(IOPT-(IOPT/10000)*10000-1000)300,500,500
  300 DO 400 I=1,NUMGR
        INDF(I)=1
  400   CONTINUE
C
C IF WE ARE DEALING WITH UNWEIGHTED APPROXIMATION SET THE
C WEIGHT FUNCTION IDENTICALLY EQUAL TO 1.
  500 IF(IOPT-(IOPT/100000)*100000-10000)600,800,800
  600 DO 700 I=1,NUMGR
        WTBLE(I)=ONE
  700   CONTINUE
C
C COPY FTBLE INTO FTBWT FOR USE BY THE REST OF THIS PROGRAM
C UNLESS WEIGHTED ERROR CURVE APPROXIMATION IS DESIRED.  IN
C THAT CASE SET FTBWT(I)=WTBLE(I)*FTBLE(I) AT THOSE POINTS
C WHERE F IS TO BE APPROXIMATED, SO WE WILL HAVE
C FTBWT(I)-WTBLE(I)*P/Q=WTBLE(I)*(FTBLE(I)-P/Q).  THE VECTOR
C FTBWT IS INTRODUCED SO THAT FTBLE NEED NOT BE CHANGED BY
C THE PROGRAM.
  800 IF(IOPT-(IOPT/100000)*100000-20000)900,1200,1200
  900 DO 1100 I=1,NUMGR
        IF(INDF(I))1100,1100,1000
 1000   FTBWT(I)=FTBLE(I)
 1100   CONTINUE
      GO TO 1500
 1200 DO 1400 I=1,NUMGR
        IF(INDF(I))1400,1400,1300
 1300   FTBWT(I)=WTBLE(I)*FTBLE(I)
 1400   CONTINUE
C
C IF THERE ARE EQUALITY CONSTRAINTS, ATTEMPT TO SOLVE FOR MEQA
C OF THE COEFFICIENTS IN TERMS OF THE REMAINING NUMP+NUMQ-MEQA
C COEFFICIENTS, WHERE MEQA IS THE NUMBER OF INDEPENDENT
C EQUALITY CONSTRAINTS.
 1500 IF(MEQ)2000,2000,1600
C HERE THERE ARE MEQ EQUALITY CONSTRAINTS, AND WE COPY THE FIRST
C MEQ ROWS OF ALC INTO EQ.
 1600 DO 1800 I=1,MEQ
        DO 1700 J=1,N
          EQ(I,J)=ALC(I,J)
 1700     CONTINUE
 1800   CONTINUE
      CALL EQELM(NMPPQ,MEQ,TOLEL,EQ,NSCOL,IRESL,MEQA,IFLAG)
C HERE IFLAG WILL BE EITHER 0 OR 7777.  IN THE LATTER CASE THE
C EQUALITY CONSTRAINTS WERE INCONSISTENT AND WE RETURN.
      IF(IFLAG)1900,1900,9600
C RESET N AND NP1 FOR USE IN SLNPRO.
 1900 NLP=N-MEQA
      NLP1=NP1-MEQA
C
C IF THE 100000 DIGIT OF IOPT IS 1, THEN THE
C INITIALIZATION BY LOEBS ALGORITHM WILL BE SKIPPED, AND
C THE COEFFICIENTS PLACED IN X BY THE USER WILL BE TAKEN FOR
C THE INITIAL APPROXIMATION.
 2000 IF(IOPT/100000)2200,2200,2100
C SET INITL=-1 TO INDICATE WE ARE USING A USER SUPPLIED INITIAL
C APPROXIMATION.
 2100 INITL=-1
      GO TO 3400
C
C HERE WE ATTEMPT TO COMPUTE AN INITIAL APPROXIMATION USING
C LOEBS ALGORITHM.
 2200 INITL=0
C WE TRY TO FIND P, Q SUCH THAT Q IS POSITIVE (AND NOT TOO
C SMALL) AT EVERY GRID POINT, P AND Q SATISFY THE CONSTRAINTS
C (INCLUDING THE CONSTRAINTS MAX(ABS(Q(J))) .LE. 1.0), AND
C MAX(ABS(F*Q-WT*P)) IS MINIMIZED AT THOSE GRID POINTS WHERE
C F IS TO BE APPROXIMATED.
 2300 CALL SSETV(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,INDLW,ALTBL,
     *INDUP,UPTBL,INDF,WTBLE,MEQ,MINEQ,ALC,0,V,M)
C SET IYRCT(1)=-1 AS A SIGNAL TO DO THE INITIAL EXCHANGES
C IN SLNPRO STRICTLY ACCORDING TO A PIVOTING STRATEGY.
      IYRCT(1)=-1
C REDUCE V IF MEQA .GT. 0
      IF(MEQA)2500,2500,2400
 2400 CALL SCRCH(M,NMPPQ,MEQA,NSCOL,IRESL,EQ,V)
C IF SLNPRO FAILED DURING THE INITIALIZATION WITH NORMAL TOLERANCES,
C AS INDICATED BY IADJS=1, WE SET M=-M AS A SIGNAL TO LOOSEN THE
C TOLERANCES IN SLNPRO (I. E. INCREASE REA AND REA1) AND TRY AGAIN.
C SLNPRO WILL RESET M TO ITS ORIGINAL VALUE.  ANY LOSS OF ACCURACY
C DUE TO THIS LOOSENING WILL PROBABLY BE CORRECTED IN LATER ITERATIONS.
 2500 M=M*(1-4*IADJS+2*IADJS**2)
C IF SLNPRO FAILED DURING THE INITIALIZATION WITH BOTH NORMAL
C TOLERANCES AND LOOSENED TOLERANCES, AS INDICATED BY IADJS=2,
C WE SET NLP=-NLP AS A SIGNAL TO TIGHTEN THE TOLERANCES IN SLNPRO
C (I. E. DECREASE REA AND REA1) AND TRY ONCE AGAIN.  SLNPRO WILL
C RESET NLP TO ITS ORIGINAL VALUE.
      NLP=NLP*(1+IADJS-IADJS**2)
      CALL SLNPRO(V,M,NLP,IYRCT,X,INDIC)
C SAVE INDIC IN CASE THE APPROXIMATION JUST COMPUTED IS THE ONE
C RETURNED.
      IRET1=INDIC
C IF THIS WAS THE FIRST LOEB ATTEMPT WE SAVE INDIC IN ILBFL FOR
C OUTPUT PURPOSES IN CASE BOTH LOEB AND THE FIRST MAIN
C ITERATION FAIL.
      IF(IADJS)2600,2600,2700
 2600 ILBFL=INDIC
C IF INDIC IS POSITIVE SOMETHING WENT WRONG IN SLNPRO.
 2700 IF(INDIC)3200,3200,2800
 2800 IF(IADJS-1)2900,2900,3000
C BEFORE CONCEDING DEFEAT WE WILL TRY AGAIN WITH CHANGED
C TOLERANCES.
 2900 IADJS=IADJS+1
      GO TO 2300
C
C HERE THE LOEB INITIALIZATION HAS FAILED, AND WE TAKE OUR INITIAL
C APPROXIMATION TO BE P/Q=0.0/PWR(.,NUMP+1).  WE ALSO SET INITL=1
C AS A WARNING THAT THIS HAS BEEN DONE.
 3000 INITL=1
      DO 3100 J=1,NMPPQ
        X(J)=ZERO
 3100   CONTINUE
      X(NMPP1)=ONE
      GO TO 3500
C
C HERE INDIC WAS NONPOSITIVE, SO SLNPRO PRODUCED AN APPROXIMATION.
C IF MEQA .GT. 0, EXPAND X TO INCLUDE THE MEQA VARIABLES
C SOLVED FOR IN SUBROUTINE EQELM.
 3200 IF(MEQA)3400,3400,3300
 3300 CALL EXPND(NMPPQ,MEQA,NSCOL,IRESL,EQ,X)
C NORMALIZE P AND Q IF POSSIBLE.
 3400 CALL NRMLZ(NUMP,NUMQ,MEQ,MINEQ,NSCOL,IRESL,ALC,TOLEL,X)
C
C COMPUTE THE VALUES OF P AND Q AT THE GRID POINTS, THE
C ERRORS AT THE GRID POINTS, AND THE UNIFORM ERROR NORM.
 3500 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X,
     *V,QMIN,IQMIN)
      DELTK=V(2*NUMGR+1,NP1)
      DLOLD=DELTK
C RESET IADJS.
      IADJS=0
C PUT THE ERRORS AND ERROR NORM INTO ERROR.
      DO 3600 I=1,NUMGR
        ERROR(I)=V(I,N)
 3600   CONTINUE
      ERROR(NUMGR+1)=DELTK
C
C IF THE DENOMINATOR WAS VERY SMALL IN ABSOLUTE VALUE OR
C NEGATIVE AT SOME POINT (WHICH IS INDICATED BY
C V(2*NUMGR+1,NP1) BEING NEGATIVE), WE RESET
C V(2*NUMGR+1,NP1) TO THE ERROR NORM COMPUTED WITH THE BAD
C POINTS IGNORED FOR USE IN SSETV2, AND ATTEMPT TO CONTINUE.
C NOTE THAT DELTK, DLOLD, AND ERROR(NUMGR+1) WILL REMAIN
C NEGATIVE AS A WARNING UNLESS AN ACCEPTABLE APPROXIMATION WITH
C SIGNIFICANTLY POSITIVE DENOMINATOR IS FOUND BELOW.
C THIS SITUATION SOMETIMES OCCURS WHEN AN APPROXIMATION
C COMPUTED ON A COARSE GRID IS USED AS A STARTING
C APPROXIMATION ON A FINE GRID.
      IF(DELTK)3700,3800,3800
 3700 V(2*NUMGR+1,NP1)=-V(2*NUMGR+1,NP1)-TWO
C SET ICON0=1 TO INDICATE THIS IS NOT AN ACCEPTABLE
C APPROXIMATION.
      ICON0=1
      GO TO 3900
C
C HERE DELTK .GE. 0.0.  MAKE ICON0=0 IF THE CONSTRAINTS ARE
C SATISFIED WITHIN EPRES AND ICON0=1 OTHERWISE.
 3800 ICON0=ICNCK(IOPT,NUMP,NUMQ,NUMGR,INDLW,ALTBL,INDUP,
     *UPTBL,EPSQ,X,V,MEQ,MINEQ,ALC,EPRES)
C
C END OF INITIALIZATION PHASE.
C
C
C THE NEXT GROUP OF STATEMENTS (THROUGH STATEMENT 5400) ARE
C FOR THE FIRST MAIN DIFFERENTIAL CORRECTION ITERATION.
C IF THIS ITERATION FAILS TO PRODUCE AN ACCEPTABLE
C APPROXIMATION AND THE INITIAL APPROXIMATION IS ALSO
C UNACCEPTABLE WE QUIT, OTHERWISE WE TRY THE NEXT
C INITIALIZATION STAGE.
C
C SET UP THE LINEAR PROGRAMMING PROBLEM.
 3900 CALL SSETV(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,INDLW,
     *ALTBL,INDUP,UPTBL,INDF,WTBLE,MEQ,MINEQ,ALC,1,V,M)
C SET IYRCT(1)=-1 AS A SIGNAL TO DO THE INITIAL EXCHANGES IN
C SLNPRO STRICTLY ACCORDING TO A PIVOTING STRATEGY RATHER
C THAN ATTEMPTING TO USE INFORMATION ABOUT A PREVIOUS VERTEX.
      IYRCT(1)=-1
C REDUCE V IF MEQA .GT. 0
      IF(MEQA)4100,4100,4000
 4000 CALL SCRCH(M,NMPPQ,MEQA,NSCOL,IRESL,EQ,V)
C SET M=-M AS A SIGNAL TO LOOSEN THE TOLERANCES IN SLNPRO IF
C IADJS=2 AND LEAVE M UNCHANGED IF IADJS=0 OR 1.  SLNPRO WILL
C RESET M TO ITS ORIGINAL VALUE IF IT IS CHANGED HERE.
 4100 M=M*(1+IADJS-IADJS**2)
C SET NLP=-NLP AS A SIGNAL TO TIGHTEN THE TOLERANCES IN SLNPRO IF
C IADJS=1 AND LEAVE NLP UNCHANGED IF IADJS=0 OR 2.  SLNPRO WILL
C RESET NLP TO ITS ORIGINAL VALUE IF IT IS CHANGED HERE.
      NLP=NLP*(1-4*IADJS+2*IADJS**2)
      CALL SLNPRO(V,M,NLP,IYRCT,XTEMP,INDIC)
C IF THIS WAS THE FIRST SLNPRO CALL IN THIS ITERATION, SO THE
C TOLERANCES IN SLNPRO WERE NOT ADJUSTED BY CDFCOR, FOR
C POSSIBLE OUTPUT PURPOSES WE SAVE INDIC IN ILST1.
      IF(IADJS)4200,4200,4300
 4200 ILST1=INDIC
C IF INDIC IS POSITIVE SOMETHING WENT WRONG IN SLNPRO SO WE
C TRY AGAIN WITH BOTH TIGHTER AND LOOSER TOLERANCES IF THIS
C HAS NOT ALREADY BEEN TRIED.
 4300 IF(INDIC)4900,4900,4400
C
C HERE WE DO NOT HAVE AN ACCEPTABLE APPROXIMATION.
 4400 IF(IADJS-1)4500,4500,4600
 4500 IADJS=IADJS+1
C CALL SPQERD WITH THE COEFFICIENTS FROM THE INITIALIZATION
C (WHICH ARE IN X) TO RECOMPUTE THE VALUES OF P AND Q AND
C PUT THEM IN V FOR USE BY SSETV.
      CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X,V,
     *QMIN,IQMIN)
      GO TO 3900
C
C HERE WE HAVE TRIED AND FAILED WITH NORMAL, TIGHTENED, AND
C LOOSENED TOLEARANCES SO THE FIRST MAIN DIFFERENTIAL
C CORRECTION ITERATION HAS FAILED TO PRODUCE AN ACCEPTABLE
C APPROXIMATION.  IF THE APPROXIMATION FROM THE INITIALIZATION
C PHASE WAS ACCEPTABLE, WE SET IFLAG AND INDLP, AND RETURN.
 4600 IF(ICON0)8500,8500,4700
C HERE THE INITIAL APPROXIMATION WAS UNACCEPTABLE.  IF IT WAS
C USER SUPPLIED WE RESET IADJS AND TRY LOEB, IF IT CAME
C FROM LOEB WE TRY P/Q = 0.0/1.0, AND IF IT WAS P/Q = 0.0/1.0
C WE SET IFLAG AND INDLP, AND RETURN.
 4700 IF(INITL)4800,3000,8500
 4800 IADJS=0
      GO TO 2200
C END OF BLOCK WITHIN FIRST MAIN DIFFERENTIAL CORRECTION
C ITERATION TO DEAL WITH AN UNACCEPTABLE APPROXIMATION.
C
C HERE INDIC WAS NONPOSITIVE, SO SLNPRO PLACED AN APPROXIMATION
C IN XTEMP, WHICH WE NOW CHECK FOR ACCEPTABLENESS.  FIRST
C EXPAND XTEMP IF MEQA .GT. O
 4900 IF(MEQA)5100,5100,5000
 5000 CALL EXPND(NMPPQ,MEQA,NSCOL,IRESL,EQ,XTEMP)
C CALL SPQERD TO COMPUTE P, Q, THE ERRORS, AND THE
C UNIFORM ERROR NORM (AND TO PLACE THESE VALUES IN V FOR
C USE BY OTHER ROUTINES).
 5100 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,
     *XTEMP,V,QMIN,IQMIN)
      DELTK=V(2*NUMGR+1,NP1)
C IF DELTK IS NEGATIVE (INDICATING A VERY SMALL OR NEGATIVE
C DENOMINATOR VALUE AT SOME POINT) WE REJECT THIS
C APPROXIMATION AND TRY AGAIN.
      IF(DELTK)4400,5200,5200
C HERE DELTK IS NONNEGATIVE, AND WE NORMALIZE XTEMP AND CHECK
C THE CONSTRAINTS.
 5200 CALL NRMLZ(NUMP,NUMQ,MEQ,MINEQ,NSCOL,IRESL,ALC,TOLEL,XTEMP)
      ICON1=ICNCK(IOPT,NUMP,NUMQ,NUMGR,INDLW,ALTBL,INDUP,
     *UPTBL,EPSQ,XTEMP,V,MEQ,MINEQ,ALC,EPRES)
C IF THE APPROXIMATION IS UNACCEPTABLE WE TRY AGAIN.
      IF(ICON1)5300,5300,4400
C HERE THE APPROXIMATION COMPUTED BY THE FIRST MAIN CDFCOR
C ITERATION SATISFIES THE CONSTRAINTS WITHIN EPRES, AND WE WILL
C ACCEPT IT IF THE APPROXIMATION FROM THE INITIALIZATION WAS
C UNACCEPTABLE OR DOES NOT HAVE A SMALLER ERROR NORM.
C OTHERWISE WE WILL REJECT IT AND TRY AGAIN SINCE SOMETHING
C IS FISHY.
 5300 IF(ICON0)5400,5400,5500
 5400 IF(DELTK-DLOLD)5500,5500,4400
C
C END OF GROUP OF STATEMENTS FOR FIRST MAIN DIFFERENTIAL
C CORRECTION ITERATION.  FROM NOW ON ONLY STATEMENTS BELOW
C THIS POINT IN CDFCOR WILL BE EXECUTED.
C
C
C HERE THE PRESENT APPROXIMATION WILL BE KEPT, AND WE PREPARE
C TO COMPUTE ANOTHER ONE.
C COPY XTEMP INTO X.
 5500 DO 5600 I=1,NMPPQ
        X(I)=XTEMP(I)
 5600   CONTINUE
C PUT THE ERRORS AND ERROR NORM INTO ERROR.
      DO 5700 I=1,NUMGR
        ERROR(I)=V(I,N)
 5700   CONTINUE
      ERROR(NUMGR+1)=DELTK
      ITER=ITER+1
      IF(ITER-ITLIM)5800,6000,6000
 5800 DLOLD=DELTK
C SAVE INDIC IN CASE THE APPROXIMATION JUST COMPUTED IS THE ONE
C RETURNED.
      IRET1=INDIC
C RESET IADJS.
      IADJS=0
C IF IPBIT=1 WE ARE STILL IN THE PARAMETRIC BISECTION PHASE,
C AND WE CUT V(2*NUMGR+1,N+1) IN HALF IN ORDER TO TRY FOR AT
C LEAST A HALVING OF THE ERROR NORM IN THE NEXT STEP.
      IF(IPBIT)6300,6300,5900
 5900 V(2*NUMGR+1,NP1)=V(2*NUMGR+1,NP1)/TWO
      GO TO 6300
C
C HERE ITLIM ITERATIONS HAVE BEEN PERFORMED AND WE WILL RETURN
C AFTER SETTING IFLAG AND INDLP (SEE BELOW).
 6000 INDLP=5-INDIC
      IF(INDIC)6100,6200,6200
 6100 IFLAG=-12
      GO TO 9400
 6200 IFLAG=-2
      GO TO 9400
C
C SAVE IYRCT IN IYSAV SO WE CAN LATER USE INFORMATION FROM THIS
C SUCCESSFUL SLNPRO ITERATION TO RESTART SLNPRO.
 6300 DO 6400 I=1,M
        IYSAV(I)=IYRCT(I)
 6400   CONTINUE
C CALL SSETV TO SET UP FOR THE NEXT ITERATION.
 6500 CALL SSETV(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,INDLW,
     *ALTBL,INDUP,UPTBL,INDF,WTBLE,MEQ,MINEQ,ALC,1,V,M)
C REDUCE V IF MEQA .GT. 0
      IF(MEQA)6700,6700,6600
 6600 CALL SCRCH(M,NMPPQ,MEQA,NSCOL,IRESL,EQ,V)
C IF IADJS=1, SET NLP=-NLP (AND LEAVE M ALONE) AS A SIGNAL TO
C TIGHTEN THE TOLERANCES IN SLNPRO, AND IF IADJS=2 SET M=-M
C (AND LEAVE NLP ALONE) AS A SIGNAL TO LOOSEN THE TOLERANCES
C IN SLNPRO.  IF IADJS=0, LEAVE THE TOLERANCES ALONE.  M AND
C NLP WILL BE RESET TO THEIR ORIGINAL VALUES BY SLNPRO IF THEY
C ARE CHANGED HERE.
 6700 M=M*(1+IADJS-IADJS**2)
      NLP=NLP*(1-4*IADJS+2*IADJS**2)
      CALL SLNPRO(V,M,NLP,IYRCT,XTEMP,INDIC)
C IF THIS WAS THE FIRST SLNPRO CALL IN THIS ITERATION, SO THE
C TOLERANCES IN SLNPRO WERE NOT ADJUSTED BY CDFCOR, FOR
C POSSIBLE OUTPUT PURPOSES WE SAVE INDIC IN ILST1.
      IF(IADJS)6800,6800,6900
 6800 ILST1=INDIC
C IF INDIC IS POSITIVE SOMETHING WENT WRONG IN SLNPRO SO WE GO
C INTO THE BLOCK OF STATEMENTS DESIGNED TO HANDLE FAILURE TO
C PRODUCE AN ACCEPTABLE APPROXIMATION.
 6900 IF(INDIC)7600,7600,7000
C
C HERE WE DID NOT GET AN ACCEPTABLE APPROXIMATION.  IF WE ARE
C STILL DOING PARAMETRIC BISECTION, TO SAVE TIME WE DO NOT TRY
C SLNPRO AGAIN WITH ADJUSTED TOLERANCES SINCE THE PARAMETRIC
C BISECTION PROBABLY CAUSED THE FAILURE, BUT RATHER WE SWITCH
C PERMANENTLY INTO ORDINARY DIFFERENTIAL CORRECTION.
 7000 IF(IPBIT)7200,7200,7100
 7100 IPBIT=0
      GO TO 7400
C IF WE WERE NOT DOING PARAMETRIC BISECTION WE TRY AGAIN WITH
C TIGHTER AND LOOSER TOLERANCES, IF THIS HAS NOT ALREADY BEEN
C TRIED.  IF IT HAS BEEN TRIED, WE SET IFLAG AND INDLP AND RETURN.
 7200 IF(IADJS-1)7300,7300,8500
 7300 IADJS=IADJS+1
C CALL SPQERD WITH THE COEFFICIENTS FROM THE PREVIOUS ITERATION
C (WHICH ARE IN X) TO RECOMPUTE THE VALUES OF P AND Q AND PUT
C THEM IN V FOR USE BY SSETV.
 7400 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X,V,
     *QMIN,IQMIN)
C RESTORE IYRCT TO ITS VALUE AT THE PREVIOUS ITERATION.
      DO 7500 I=1,M
        IYRCT(I)=IYSAV(I)
 7500   CONTINUE
      GO TO 6500
C END OF BLOCK OF STATEMENTS WITHIN THE SECOND OR LATER MAIN
C DIFFERENTIAL CORRECTION ITERATION TO HANDLE FAILURE TO
C PRODUCE AN ACCEPTABLE APPROXIMATION.
C
C HERE INDIC WAS NONPOSITIVE, SO SLNPRO PRODUCED AN
C APPROXIMATION IN XTEMP WHICH WE NOW CHECK FOR
C ACCEPTABILITY.  FIRST EXPAND XTEMP IF MEQA .GT. 0
 7600 IF(MEQA)7800,7800,7700
 7700 CALL EXPND(NMPPQ,MEQA,NSCOL,IRESL,EQ,XTEMP)
C CALL SPQERD TO COMPUTE P, Q, THE ERRORS, AND THE UNIFORM
C ERROR NORM (AND PLACE THEM IN V FOR USE BY OTHER ROUTINES).
 7800 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,
     *XTEMP,V,QMIN,IQMIN)
      DELTK=V(2*NUMGR+1,NP1)
C IF DELTK IS NEGATIVE (INDICATING A VERY SMALL OR NEGATIVE
C DENOMINATOR VALUE AT SOME POINT) WE REJECT THIS
C APPROXIMATION AND TRY AGAIN.
      IF(DELTK)7000,7900,7900
C IF INDIC IS NEGATIVE HERE IT WAS NECESSARY
C TO ADJUST TOLERANCES IN SLNPRO, SO THE APPROXIMATION
C COMPUTED IS MORE LIKELY TO HAVE BEEN AFFECTED
C BY ROUNDOFF ERROR.  THUS IN THIS CASE WE CHECK THE
C CONSTRAINTS, AND REJECT THE APPROXIMATION IF ANY OF THEM ARE
C VIOLATED BY MORE THAN EPRES.  NOTE THAT INDIC CANNOT BE
C POSITIVE HERE SINCE SLNPRO PRODUCED AN APPROXIMATION.
 7900 IF(INDIC)8000,8100,8100
 8000 ICON=ICNCK(IOPT,NUMP,NUMQ,NUMGR,INDLW,ALTBL,INDUP,
     *UPTBL,EPSQ,XTEMP,V,MEQ,MINEQ,ALC,EPRES)
      IF(ICON)8100,8100,7000
C REJECT THIS APPROXIMATION AND TRY AGAIN IF THE ERROR NORM HAS
C NOT DECREASED.
 8100 IF(DELTK-DLOLD)8200,7000,7000
C RARELY AN APPROXIMATION WITH ALL COEFFICIENTS VERY SMALL
C WILL BE PRODUCED.  TO BE SAFE, IF ALL DENOMINATOR
C COEFFICIENTS ARE .LE. 0.1 IN ABSOLUTE VALUE WE ATTEMPT TO
C NORMALIZE THE APPROXIMATION AND CHECK THE CONSTRAINTS.  IF
C THEY ARE SATISFIED WITHIN EPRES WE ACCEPT THIS
C APPROXIMATION, AND OTHERWISE WE TRY AGAIN.
 8200 DO 8400 J=NMPP1,NMPPQ
        IF(XTEMP(J)-TENTH)8300,8300,5500
 8300   IF(XTEMP(J)+TENTH)5500,8400,8400
 8400   CONTINUE
      CALL NRMLZ(NUMP,NUMQ,MEQ,MINEQ,NSCOL,IRESL,ALC,TOLEL,
     *XTEMP)
      ICON=ICNCK(IOPT,NUMP,NUMQ,NUMGR,INDLW,ALTBL,INDUP,
     *UPTBL,EPSQ,XTEMP,V,MEQ,MINEQ,ALC,EPRES)
      IF(ICON)5500,5500,7000
C END OF COMPUTATION OF APPROXIMATIONS.
C
C
C HERE WE ARE FINISHED COMPUTING APPROXIMATIONS, SO WE SET INDLP,
C IFLAG, QMIN, AND IQMIN, AND WE RETURN.
C WE CARRY INFORMATION BACK TO THE CALLING PROGRAM ABOUT TWO
C SLNPRO CALLS BY SETTING THE TWO DIGIT VARIABLE INDLP WITH
C TENS DIGIT 5-ILAST AND ONES DIGIT 5-IRET, WHERE IN MOST
C CASES ILAST AND IRET GIVE INFORMATION ON THE LAST (REJECTED)
C SLNPRO ATTEMPT AND THE RETURNED APPROXIMATION RESPECTIVELY.
C TO BE EXACT, IF THE PROGRAM WAS TERMINATED BECAUSE ITLIM
C ITERATIONS WERE PERFORMED SUCCESSFULLY, ILAST WILL BE 5 (AND
C SO 5-ILAST WILL BE 0).  OTHERWISE, ILAST WILL BE THE VALUE OF
C INDIC PRODUCED BY SLNPRO AT THE LAST (FAILED) MAIN ITERATION,
C WITH THE TOLERANCES IN SLNPRO NOT HAVING BEEN ADJUSTED BY
C CDFCOR YET.
C IF THE PROGRAM FAILED IN THE INITIALIZATION (SO IFLAG WAS
C POSITIVE), IRET WILL BE THE VALUE OF INDIC PRODUCED BY
C SLNPRO DURING THE (FAILED) LOEB INITIALIZATION, WITH THE
C TOLERANCES IN SLNPRO NOT HAVING BEEN ADJUSTED BY CDFCOR
C YET.  IF A USER SUPPLIED INITIAL APPROXIMATION WAS RETURNED
C TO THE USER, IRET WILL BE O.  OTHERWISE, IRET WILL BE THE
C VALUE OF INDIC PRODUCED BY SLNPRO DURING THE COMPUTATION
C OF THE APPROXIMATION WHICH WAS ACTUALLY RETURNED TO THE USER.
C THUS IF THE APPROXIMATION RETURNED WAS COMPUTED BY SLNPRO
C (WHICH WILL BE THE CASE IF EITHER THE ONES DIGIT OF IFLAG IS
C 0 OR 2, OR IT IS 1 AND THE RETURNED APPROXIMATION IS NOT A USER
C SUPPLIED INITIAL APPROXIMATION) AND THE ONES DIGIT OF INDLP IS 5,
C THEN SLNPRO APPEARED TO WORK NORMALLY DURING THE COMPUTATION OF
C THE RETURNED APPROXIMATION.
 8500 INDLP=10*(5-ILST1)+5-IRET1
C
C WE NOW SET UP IFLAG FOR THE INFORMATION OF THE USER.
      IF(ITER)8600,8600,9200
C
C HERE ITER=0, SO THE INITIAL APPROXIMATION WAS NOT IMPROVED.
 8600 IF(INITL)9100,9100,8700
C HERE THE UNIMPROVED INITIAL APPROXIMATION WAS P/Q=
C 0.0/PWR(.,NUMP+1).  IF IT WAS FOUND BY SPQERD TO HAVE A
C DENOMINATOR WHICH IS VERY SMALL IN ABSOLUTE VALUE OR
C NEGATIVE AT SOME POINT, OR IF IT VIOLATED A CONSTRAINT (EITHER OF
C WHICH WOULD BE INDICATED BY ICON0=1), THEN SET IFLAG=6666,
C ITER=-1, AND ZERO OUT THE ERRORS AND COEFFICIENTS AND RETURN.
C OTHERWISE SET IFLAG=3333 AND RETURN.
C FIRST RESET INDLP SO THAT 5-(ONES DIGIT OF INDIC) IS THE VALUE
C OF INDIC THE FIRST TIME SLNPRO WAS CALLED IN THE (FAILED) LOEB
C INITIALIZATION.
 8700 INDLP=10*(5-ILST1)+5-ILBFL
      IF(ICON0)9000,9000,8800
 8800 IFLAG=6666
      ITER=-1
C NOTE THAT HERE X(I) ALREADY EQUALS 0.0 FOR I=1,...,NMPPQ,
C I.NE.NUMP+1.
      X(NMPP1)=ZERO
      DO 8900 I=1,NUMGR
        ERROR(I)=ZERO
 8900   CONTINUE
      RETURN
 9000 IFLAG=3333
      RETURN
C
C HERE THE UNIMPROVED INITIAL APPROXIMATION WAS USER SUPPLIED OR
C COMPUTED BY THE LOEB INITIALIZATION.  SET THE ONES DIGIT OF
C IFLAG TO INDICATE THE INITIAL APPROXIMATION WAS NOT IMPROVED.
 9100 IFLAG=-1
C
C IF IRET1 .LT. 0 SET THE TENS DIGIT OF IFLAG TO INDICATE
C THAT THE TOLERANCES IN SLNPRO HAD TO BE ADJUSTED DURING THE
C COMPUTATION OF THE APPROXIMATION RETURNED.
 9200 IF(IRET1)9300,9400,9400
 9300 IFLAG=IFLAG-10
C NOW WE CHECK THE USER SUPPLIED CONSTRAINTS (IF ANY).
C CALL SPQERD TO RECOMPUTE THE VALUES OF P AND Q.  QMIN AND
C IQMIN ARE ALSO RECOMPUTED BY SPQERD.
 9400 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X,
     *V,QMIN,IQMIN)
      ICON=ICNCK(IOPT,NUMP,NUMQ,NUMGR,INDLW,ALTBL,INDUP,
     *UPTBL,EPSQ,X,V,MEQ,MINEQ,ALC,EPRES)
      IF(ICON)9600,9600,9500
 9500 IFLAG=IFLAG-100
 9600 RETURN
      END
      SUBROUTINE STCOMP(NUMDM,NRWGR,NCLGR,ASWP1,ASWP2,
     *ANEP1,ANEP2,T)
C***BEGIN PROLOGUE  STCOMP
C***DATE WRITTEN   720120   (YYMMDD)
C***REVISION DATE  880226   (YYMMDD)
C***CATEGORY NO.  M2,Q
C***KEYWORDS  GRID CONSTRUCTION
C***AUTHOR  KAUFMAN, EDWIN H. JR.,
C             DEPARTMENT OF MATHEMATICS
C             CENTRAL MICHIGAN UNIVERSITY
C             MOUNT PLEASANT, MICHIGAN 48859
C           TAYLOR, GERALD D.,
C             DEPARTMENT OF MATHEMATICS
C             COLORADO STATE UNIVERSITY
C             FORT COLLINS, COLORADO 80523
C***PURPOSE  THIS SUBROUTINE, WHICH IS CALLED ONLY BY THE
C            USER, CAN BE USED TO COMPUTE THE COORDINATES
C            OF THE GRID POINTS IF THEY ARE AN EQUALLY
C            SPACED (IN EACH DIRECTION) SUBSET OF AN INTERVAL
C            OR RECTANGLE.
C***DESCRIPTION
C
C  NUMDM  (INPUT) THIS IS THE NUMBER OF DIMENSIONS (1 OR 2).
C
C  NRWGR  (INPUT) THIS NEED NOT BE ASSIGNED A VALUE IN THE
C         1 DIMENSIONAL CASE.  IN THE 2 DIMENSIONAL CASE IT
C         IS THE NUMBER OF ROWS OF THE GRID.
C
C  NCLGR  (INPUT) THIS IS THE NUMBER OF GRID POINTS IN THE 1
C         DIMENSIONAL CASE.  IN THE 2 DIMENSIONAL CASE IT
C         IS THE NUMBER OF COLUMNS OF THE GRID.
C
C  ASWP1  (INPUT) THIS IS THE LEFT END POINT IN THE 1 DIMENSIONAL
C         CASE.  IN THE 2 DIMENSIONAL CASE IT IS THE FIRST
C         COORDINATE OF THE LOWER LEFT CORNER POINT.
C
C  ASWP2  (INPUT) THIS IS THE RIGHT END POINT IN THE 1 DIMENSIONAL
C         CASE.  IN THE 2 DIMENSIONAL CASE IT IS THE SECOND
C         COORDINATE OF THE LOWER LEFT CORNER POINT.
C
C  ANEP1  (INPUT) THIS NEED NOT BE ASSIGNED A VALUE IN THE 1
C         DIMENSIONAL CASE.  IN THE 2 DIMENSIONAL CASE IT IS
C         THE FIRST COORDINATE OF THE UPPER RIGHT CORNER POINT.
C
C  ANEP2  (INPUT) THIS NEED NOT BE ASSIGNED A VALUE IN THE 1
C         DIMENSIONAL CASE.  IN THE 2 DIMENSIONAL CASE IT IS
C         THE SECOND COORDINATE OF THE UPPER RIGHT CORNER POINT.
C
C  T      (OUTPUT) IN THE 1 DIMENSIONAL CASE THE COORDINATE
C         OF THE ITH GRID POINT WILL BE PLACED IN T(I). IN
C         THE 2 DIMENSIONAL CASE THE COORDINATES OF THE ITH
C         GRID POINT (READING ROW BY ROW, TOP TO BOTTOM) WILL
C         BE PLACED IN (T(2*I-1),T(2*I)).
C***REFERENCES  (NONE)
C***ROUTINES CALLED  (NONE)
C***END PROLOGUE  STCOMP
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION T(202)
C
C SET PRECISION DEPENDENT CONSTANTS FOR STCOMP.
C***FIRST EXECUTABLE STATEMENT  STCOMP
      ONE=1.0D0
      ZERO=ONE-ONE
C END SETTING PRECISION DEPENDENT CONSTANTS FOR STCOMP.
C
      IF(NUMDM-2)1,3,1
C
C THIS IS THE 1 DIMENSIONAL CASE.
    1 HSPAC=(ASWP2-ASWP1)/FLOAT(NCLGR-1)
      DO 2 I=1,NCLGR
        T(I)=ASWP1+HSPAC*FLOAT(I-1)
    2   CONTINUE
      RETURN
C
C THIS IS THE 2 DIMENSIONAL CASE.
    3 NCGDB=2*NCLGR
      IF(NCLGR-1)4,5,4
    4 HSPAC=(ANEP1-ASWP1)/FLOAT(NCLGR-1)
      GO TO 6
    5 HSPAC=ZERO
    6 IF(NRWGR-1)7,8,7
    7 VSPAC=(ANEP2-ASWP2)/FLOAT(NRWGR-1)
      GO TO 9
    8 VSPAC=ZERO
    9 DO 11 INDEX=1,NRWGR
        YCOR=ANEP2-VSPAC*FLOAT(INDEX-1)
        IND1=2+NCGDB*(INDEX-1)
        IND2=IND1+NCGDB-2
C FILL IN THE Y COORDINATES IN ROW INDEX.
        DO 10 J=IND1,IND2,2
          T(J)=YCOR
   10     CONTINUE
   11   CONTINUE
      DO 13 INDEX=1,NCLGR
        XCOR=ASWP1+HSPAC*FLOAT(INDEX-1)
        IND1=2*INDEX-1
        IND2=IND1+NCGDB*(NRWGR-1)
C FILL IN THE X COORDINATES IN COLUMN INDEX.
        DO 12 J=IND1,IND2,NCGDB
          T(J)=XCOR
   12     CONTINUE
   13   CONTINUE
      RETURN
      END
      SUBROUTINE SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X,
     *V,QMIN,IQMIN)
C***BEGIN PROLOGUE  SPQERD
C***REFER TO CDFCOR
C***ROUTINES CALLED  (NONE)
C***PURPOSE  THIS SUBROUTINE COMPUTES THE VALUES OF P, Q, AND
C            THE ERROR F-WT*P/Q, AS WELL AS THE UNIFORM ERROR NORM
C            AND THE MINIMUM OF THE DENOMINATOR, FOR THE
C            DIFFERENTIAL CORRECTION ALGORITHM.
C***END PROLOGUE  SPQERD
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION FTBWT(101),PWR(101,15),V(635,17),WTBLE(101),
     *   X(16),INDF(101)
C
C FOR EACH GRID POINT I THIS SUBROUTINE COMPUTES P,Q, AND
C THE ERROR F-WT*P/Q.  THE UNIFORM ERROR NORM
C MAX ABS(F-WT*P/Q) IS ALSO COMPUTED.  RECALL THAT IF WE ARE
C DOING WEIGHTED ERROR CURVE APPROXIMATION, F WILL ACTUALLY
C BE WT*F.
C
C SET MACHINE DEPENDENT PARAMETERS FOR SUBROUTINE SPQERD.
C SET SPCMN=B**(-T), WHERE B IS THE BASE AND T IS THE NUMBER
C OF BASE B DIGITS IN THE MANTISSA.  SPCMN IS THE MINIMUM
C RELATIVE SPACING ABS((X1-X2)/X2) BETWEEN TWO SUCCESSIVE
C FLOATING POINT NUMBERS, SO IT IS THE SPACING BETWEEN TWO
C SUCCESSIVE FLOATING POINT NUMBERS IN THE CLOSED INTERVAL
C (0.1,1.0).  WE ALSO HAVE SPCMN=10.0**(-T*(LOG10)(B))=
C 10.0**(-TNMAN), WHERE TNMAN IS THE BASE 10 EQUIVALENT OF
C THE LENGTH OF THE MANTISSA.
C
C***FIRST EXECUTABLE STATEMENT  SPQERD
      SPCMN=D1MACH(3)
C SET PRECISION DEPENDENT CONSTANTS FOR SPQERD.
      ONE=1.0D0
      ZERO=ONE-ONE
      TWO=ONE+ONE
      FOUR=TWO+TWO
      TEN=FOUR+FOUR+TWO
C SET EPS2=100.0*SPCMN.  THUS EPS2=10.0**(-(TNMAN-2)).
C THE DENOMINATOR WILL NOT BE ALLOWED TO BE SMALLER THAN
C EPS2 AT ANY GRID POINT.  THIS WILL AVOID A DIVIDE FAULT
C AND ALSO REJECT SOME BAD APPROXIMATIONS.
C
      EPS2=TEN*TEN*SPCMN
C END OF SETTING MACHINE AND PRECISION DEPENDENT PARAMETERS FOR SPQERD.
C
      NMPP1=NUMP+1
      NMPPQ=NUMP+NUMQ
      N=NMPPQ+1
      NP1=N+1
      NN=2*NUMGR+1
      I=0
      IDBM=-1
      IDB=0
      V(NN,NP1)=ZERO
      ABAD=ONE
    1 I=I+1
      IF(I-NUMGR)3,3,1001
C
C REPLACE V(NN,NP1) BY -V(NN,NP1)-2.0 IF THE DENOMINATOR AT
C SOME POINT WAS VERY SMALL IN ABSOLUTE VALUE OR NEGATIVE.
 1001 V(NN,NP1)=ABAD*V(NN,NP1)+ABAD-ABS(ABAD)
      RETURN
C
    3 IDBM=IDBM+2
      IDB=IDB+2
C
C COMPUTE THE VALUE OF P AT GRID POINT I AND PUT IT IN
C V(2*I-1,N+1).
      V(IDBM,NP1)=ZERO
      DO 4 J=1,NUMP
        V(IDBM,NP1)=V(IDBM,NP1)+X(J)*PWR(I,J)
    4   CONTINUE
C
C COMPUTE THE VALUE OF Q AT GRID POINT I AND PUT IT IN
C V(2*I,N+1).
      V(IDB,NP1)=ZERO
      DO 5 J=NMPP1,NMPPQ
        V(IDB,NP1)=V(IDB,NP1)+X(J)*PWR(I,J)
    5   CONTINUE
C KEEP TRACK OF THE MINIMUM OF THE DENOMINATOR IN QMIN AND THE
C INDEX OF THE GRID POINT WHERE IT OCCURS IN IQMIN.
      IF(I-1)10010,10010,10000
10000 IF(V(IDB,NP1)-QMIN)10010,1003,1003
10010 QMIN=V(IDB,NP1)
      IQMIN=I
 1003 IF(V(IDB,NP1)-EPS2)1005,10020,10020
C
C NOW COMPUTE THE ERROR AT GRID POINT I AND PUT IT IN
C V(I,N).  KEEP TRACK OF THE ERROR NORM IN V(2*NUMGR+1,N+1).
C IF INDF(I)=0, THEN F IS TO BE IGNORED AT GRID POINT I,
C SO SET V(I,N)=0.0 AND LOOK AT THE NEXT GRID POINT.
10020 IF(INDF(I))10030,10030,1004
10030 V(I,N)=ZERO
      GO TO 1
 1004 V(I,N)=FTBWT(I)-WTBLE(I)*V(IDBM,NP1)/V(IDB,NP1)
      ABERR=ABS(V(I,N))
      IF(ABERR-V(NN,NP1))1,1,6
    6 V(NN,NP1)=ABERR
      GO TO 1
C
C HERE V(2*I,N+1) IS VERY SMALL IN ABSOLUTE VALUE OR NEGATIVE.
C WE SET THE ERROR EQUAL TO ZERO AND REPLACE THE FINAL ERROR
C NORM BY ITS NEGATIVE, MINUS 2.0, AS A WARNING.  THIS
C SITUATION SOMETIMES OCCURS IF THE PREVIOUS APPROXIMATION
C WAS NEAR BEST, OR IF AN APPROXIMATION COMPUTED ON A COARSE
C GRID WAS USED TO INITIALIZE THE MAIN ALGORITHM ON A FINE
C GRID.
 1005 ABAD=-ONE
      V(I,N)=ZERO
C TO INCREASE THE PROBABILITY THAT THE DENOMINATOR OF THE
C NEXT COMPUTED APPROXIMATION WILL BE .GE. EPS2 IF THE MAIN
C PROGRAM CONTINUES ON, WE RESET THE DENOMINATOR OF THE
C PRESENT APPROXIMATION AT THE BAD POINTS.
      V(IDB,NP1)=TEN*EPS2
      GO TO 1
      END
      SUBROUTINE SLNPRO(V,M,N,IYRCT,X,INDIC)
C***BEGIN PROLOGUE  SLNPRO
C***REFER TO  CDFCOR
C***ROUTINES CALLED  SJELIM
C***PURPOSE  THIS SUBROUTINE SOLVES THE LINEAR PROGRAMMING PROBLEM
C              MAXIMIZE Z = -V(M+1,1)*X(1)-...-V(M+1,N)*X(N)
C              WHERE X(1),...,X(N) ARE FREE, SUBJECT TO
C              V(I,1)*X(1)+...+V(I,N)*X(N) .LE. V(I,N+1), FOR I=1,..,N.
C            (INFORMATION CONCERNING TOLERANCES AND BASIC VARIABLES
C            IS ALSO TRANSMITTED USING M, N, AND IYRCT.)
C***REFERENCES  AVDEYEVA, L. I. AND ZUKHOVITSKIY, S. I.,
C                 LINEAR AND CONVEX PROGRAMMING,
C                 SAUNDERS, PHILADELPHIA, 1966.
C***END PROLOGUE  SLNPRO
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION V(635,17),IYRCT(634),X(16),Y(634),
     *   IXRCT(634),IYCCT(16)
C
C GIVEN INTEGERS M AND N (WITH M.GE.N) AND A MATRIX V,
C THIS SUBROUTINE SOLVES THE LINEAR PROGRAMMING PROBLEM
C    MAXIMIZE Z=-V(M+1,1)X(1)-...-V(M+1,N)X(N)+V(M+1,N+1)
C SUBJECT TO THE CONSTRAINTS
C    V(I,1)X(1)+...+V(I,N)X(N).LE.V(I,N+1), I=1,...,M
C USING ESSENTIALLY THE METHOD IN THE BOOK BY AVDEYEVA AND
C ZUKHOVITSKIY.  Y(I)=-V(I,1)X(1)-...-V(I,N)X(N)+V(I,N+1),
C I=1,...,M ARE SLACK VARIABLES.  THE METHOD HAS 4 PHASES.
C
C FIRST, XS ARE EXCHANGED FOR YS TO GET A PROBLEM
C INVOLVING ONLY NONNEGATIVE VARIABLES, THE YS BEING
C SELECTED IN THE ORDER DETERMINED BY IYRCT AND A PIVOTING
C STRATEGY.  AT THE BEGINNING OF THIS ROUTINE WE MUST HAVE
C IYRCT(1) NONPOSITIVE, OR IYRCT MUST CONTAIN SOME
C PERMUTATION OF THE INTEGERS 1,...,M (SEE BELOW).
C SECOND, THE SLACK CONSTANTS OF THE DUAL PROBLEM ARE MADE
C (SIGNIFICANTLY) NONNEGATIVE.
C THIRD, THE COST COEFFICIENTS OF THE DUAL PROBLEM ARE MADE
C (SIGNIFICANTLY) NONNEGATIVE.
C FINALLY, THE SOLUTION VECTOR IS COMPUTED.
C
C THE VARIABLE INDIC WILL BE GIVEN VALUE
C 0, IF A SOLUTION WAS FOUND NORMALLY
C 1, IF THERE WAS TROUBLE IN PHASE 1
C 2, IF THERE WAS TROUBLE IN PHASE 2 (EITHER ROUND OFF
C   ERROR, OR THE ORIGINAL PROBLEM WAS INCONSISTENT OR
C   UNBOUNDED)
C 3, IF THERE WAS TROUBLE IN PHASE 3 (EITHER ROUND OFF
C   ERROR, OR THE ORIGINAL CONSTRAINTS WERE INCONSISTENT)
C 4, IF 300 MODIFIED JORDAN ELIMINATIONS WERE USED IN
C   PHASES 2 AND 3
C -1, IF A SOLUTION WAS FOUND BUT IN ORDER TO OVERCOME
C   TROUBLE IN PHASE 2 OR 3 IT WAS NECESSARY TO TEMPORARILY
C   RELAX THE RESTRICTION ON DIVISION BY NUMBERS WITH SMALL
C   ABSOLUTE VALUE.  THUS THERE IS AN INCREASED CHANCE OF
C   SERIOUS ROUNDOFF ERROR IN THE RESULTS.
C -2, IF A SOLUTION WAS FOUND NORMALLY, EXCEPT THAT
C   THE PARAMETERS REA AND REA1 WERE INCREASED BY A SIGNAL
C   FROM THE CALLING PROGRAM (NAMELY, M=-M).  THE INCREASED
C   TOLERANCES MAY HAVE ALLOWED MORE ERROR THAN USUAL.
C -3, IF IN ORDER TO COMPLETE PHASE 1 IT WAS NECESSARY TO
C   TEMPORARILY RELAX THE RESTRICTION ON DIVISION BY NUMBERS
C   WITH SMALL ABSOLUTE VALUE.  THUS THERE IS AN INCREASED
C   CHANCE OF SERIOUS ROUNDOFF ERROR IN THE RESULTS.
C -4, IF A SOLUTION WAS FOUND NORMALLY, EXCEPT THAT REA AND REA1
C   WERE DECREASED BY A SIGNAL FROM THE CALLING PROGRAM (NAMELY
C   N=-N) IN ORDER TO TRY FOR MORE ACCURACY.  THIS INCREASES THE
C   CHANCES OF SERIOUS ROUNDOFF ERROR IN THE RESULTS.
C NOTE THAT INDIC=-3 WILL OVERWRITE (AND THUS CONCEAL) INDIC=-2
C   OR INDIC=-4, AND INDIC=-1 WILL OVERWRITE INDIC=-2, -3, OR -4
C
C SET MACHINE DEPENDENT PARAMETERS FOR SUBROUTINE SLNPRO.
C SET SPCMN=B**(-T), WHERE B IS THE BASE AND T IS THE NUMBER
C OF BASE B DIGITS IN THE MANTISSA.  SPCMN IS THE MINIMUM
C RELATIVE SPACING ABS((X1-X2)/X2) BETWEEN TWO SUCCESSIVE
C FLOATING POINT NUMBERS, SO IT IS THE SPACING BETWEEN TWO
C SUCCESSIVE FLOATING POINT NUMBERS IN THE CLOSED INTERVAL
C (0.1,1.0).  WE ALSO HAVE SPCMN=10.0**(-T*(LOG10)(B))=
C 10.0**(-TNMAN), WHERE TNMAN IS THE BASE 10 EQUIVALENT OF
C THE LENGTH OF THE MANTISSA.
C
C***FIRST EXECUTABLE STATEMENT  SLNPRO
      SPCMN=D1MACH(3)
C SET PRECISION DEPENDENT CONSTANTS FOR SLNPRO.
      ONE=1.0D0
      ZERO=ONE-ONE
      TWO=ONE+ONE
      FOUR=TWO+TWO
      TEN=FOUR+TWO+TWO
C SET REA (ROUND OFF ERROR ADJUSTMENT) =
C MAX(10.0**(-8),100.0*SPCMN).  THUS REA=10.0**(-EXREA),
C WHERE EXREA=MIN(8,TNMAN-2).
C DIVISION BY NUMBERS .LE. REA IN ABSOLUTE VALUE WILL NOT BE
C PERMITTED.
      REA=TEN*TEN*SPCMN
      IF(REA-TEN**(-8))10000,10010,10010
10000 REA=TEN**(-8)
C SET REA1=10.0*SPCMN.  THUS REA1=10.0**(-(TNMAN-1)).
C NUMBERS IN ROW M+1 OR COLUMN N+1 WHICH ARE .LE. REA1 IN
C ABSOLUTE VALUE WILL BE TREATED AS ZEROES.  SLNPRO ASSUMES
C THAT 0.0.LT.REA1.LE.REA.
10010 REA1=TEN*SPCMN
C END OF INITIAL SETTING OF MACHINE DEPENDENT PARAMETERS FOR
C SLNPRO.  THESE PARAMETERS MAY BE ADJUSTED BY A COMMAND FROM
C CDFCOR.
C
      INDIC=0
C M BEING NEGATIVE IS A SIGNAL TO INCREASE REA AND REA1,
C THUS TREATING MORE NUMBERS WITH SMALL ABSOLUTE VALUES AS
C ZEROES.  THIS MAY GIVE THIS ROUTINE A BETTER CHANCE TO
C SUCCEED, BUT MAY ALSO CAUSE LARGER ERRORS.
      IF(M)1001,10001,10001
C RESET M.
 1001 M=-M
      REA=SQRT(REA)
      REA1=SQRT(REA1)
      INDIC=-2
C N BEING NEGATIVE IS A SIGNAL TO DECREASE REA AND REA1 TO TRY
C FOR MORE ACCURACY.  AMONG OTHER THINGS, THIS MAKES IT MORE
C LIKELY THAT THE PREVIOUS VERTEX WILL BE RETAINED IN PHASE 1
C BELOW, BUT IT ALSO COULD INCREASE ROUND OFF ERROR.
10001 IF(N)10002,1002,1002
C RESET N.
10002 N=-N
      REA=REA1
      REA1=REA1/(TEN*TEN)
      INDIC=-4
C PRESERVE REA IN CASE IT MUST BE TEMPORARILY RELAXED.
C IRLAX=0 INDICATES REA IS NOT RELAXED AT THIS STAGE.
 1002 REAKP=REA
      IRLAX=0
C IN COLUMN N+1, NUMBERS .LE. REA2 IN ABSOLUTE VALUE WILL BE
C TREATED AS ZEROES.
      REA2=REA1
      NP1=N+1
      MP1=M+1
      KTJOR=0
      IBACK=0
C SET V(MP1,NP1)=0.0 SO THE DESCRIPTIONS IN AND FOLLOWING THE
C PROLOGUE WILL AGREE.
      V(MP1,NP1)=ZERO
C THE ONLY REASON FOR THE FOLLOWING THREE STATEMENTS IS TO
C AVOID THE ERROR MESSAGE ON SOME MACHINES THAT THESE
C VARIABLES HAVE NOT BEEN ASSIGNED A VALUE.  THEY WILL BE
C REASSIGNED A VALUE BEFORE THE PROGRAM REACHES A SPOT WHERE
C THEY WILL ACTUALLY BE USED.
      DIST=ONE
      AMPRV=ONE
      AMPR2=ONE
C SET IXRCT.  IXRCT(I)=0 MEANS SOME Y IS IN ROW I, WHILE
C IXRCT(I)=K.NE.0 MEANS X(K) IS IN ROW I.
      DO 1 I=1,M
        IXRCT(I)=0
    1   CONTINUE
C
C EXCHANGE THE XS AT THE TOP OF THE TABLE FOR YS.
C IF IYRCT(1) IS NONPOSITIVE, WE SET IYRCT AND CHOOSE THE
C LARGEST POSSIBLE RESOLVENTS FOR THE EXCHANGES.  IF
C IYRCT(1) IS POSITIVE, IYRCT WILL HAVE BEEN PREVIOUSLY SET
C AND WE TRY TO EXCHANGE IN ROWS IYRCT(1),...,IYRCT(N),
C STILL EMPLOYING A PIVOTING STRATEGY, BUT IF WE CANNOT, WE
C EXCHANGE IN ROWS IYRCT(N+1),...,IYRCT(M).
      IF(IYRCT(1))1003,1003,1005
 1003 I10=1
      I20=M
C IF WE HAVE NO INFORMATION FROM A PREVIOUS VERTEX, WE GIVE
C UP A LITTLE ACCURACY IN COLUMN N+1 TO HAVE A BETTER CHANCE
C OF SUCCESS.
      REA2=REA
C THIS ROUTINE HAS A BACKTRACKING OPTION WHICH SOMETIMES
C INCREASES ACCURACY BUT SOMETIMES LEADS TO FAILURE DUE TO
C CYCLING.  IT IS SUGGESTED THAT THIS OPTION BE EMPLOYED IF
C INFORMATION ABOUT A STARTING VERTEX IS AVAILABLE, AND
C OTHERWISE BE DISABLED BY SETTING IBACK=1.
      IBACK=1
      DO 1004 I=1,M
        IYRCT(I)=I
 1004   CONTINUE
      GO TO 1006
 1005 I10=1
      I20=N
 1006 J=0
C SET THE LOWER BOUND ON THE ABSOLUTE VALUE OF A RESOLVENT IN
C PHASE 1.  ALSO SET IFAIL=0 TO INDICATE THE RESOLVENT SEARCH
C IN THIS COLUMN HAS NOT FAILED.
      REA3=REA
      IFAIL=0
    2 J=J+1
      IF(J-N)1007,1007,9
C SET I1, I2 ACCORDING TO THE STRATEGY WE ARE USING.
 1007 I1=I10
      I2=I20
      AMAX=ZERO
C SEARCH FOR A RESOLVENT IN ROWS IYRCT(I1),...,IYRCT(I2).
10003 DO 1012 I=I1,I2
        IYTMP=IYRCT(I)
        IF(IXRCT(IYTMP))1012,1009,1012
 1009   ABSV=ABS(V(IYTMP,J))
        IF(ABSV-AMAX)1012,1012,1011
 1011   IYRI=IYTMP
        AMAX=ABSV
 1012   CONTINUE
C CHECK TO SEE IF THE PROSPECTIVE RESOLVENT IS LARGE ENOUGH
C IN ABSOLUTE VALUE.
      IF(AMAX-REA3)1013,1013,7
C EXCHANGE X(J) FOR Y(IYRI).
    7 CALL SJELIM(MP1,1,NP1,IYRI,J,V)
      IXRCT(IYRI)=J
      IYCCT(J)=IYRI
C IYCCT(J)=IYRI MEANS Y(IYRI) IS IN COLUMN J.
C RESET REA3 AND IFAIL SINCE WE SUCCESSFULLY FOUND A RESOLVENT IN
C THIS COLUMN, AND THE RESOLVENT SEARCH IN THE NEXT COLUMN HAS
C NOT FAILED.
      REA3=REA
      IFAIL=0
      GO TO 2
C WE HAVE NOT FOUND A SUITABLE RESOLVENT IN ROWS IYRCT(I1),
C ...IYRCT(I2).  IF I2.LT.M WE SEARCH THE REST OF COLUMN J.
 1013 IF(I2-M)1014,10004,10004
 1014 I1=I2+1
      I2=M
      GO TO 10003
C HERE WE FAILED TO FIND A RESOLVENT IN COLUMN J WITH ABSOLUTE
C VALUE .GT. REA3.  IF IFAIL=0 WE SET INDIC=-3 AND TRY AGAIN
C WITH REA3 REDUCED.  IF THIS HAS ALREADY BEEN TRIED WE SET
C INDIC=1 AND RETURN.
10004 IF(IFAIL)10005,10005,8
10005 IFAIL=1
      INDIC=-3
      REA3=REA1
      GO TO 1007
    8 INDIC=1
      RETURN
C
C REARRANGE THE ROWS OF V SO THAT X(1),...,X(N) COME FIRST
C IN THAT ORDER.  REDEFINE IYRCT SO THAT AFTER THE
C REARRANGEMENT IS DONE, IYRCT(I)=K WILL MEAN Y(K) IS IN
C ROW I (FOR I GREATER THAN N).
    9 DO 10 I=1,M
        IYRCT(I)=I
   10   CONTINUE
      IROW=0
   11 IROW=IROW+1
      IF(IROW-M)12,12,20
   12 IF(IXRCT(IROW))13,11,13
   13 IF(IXRCT(IROW)-IROW)14,11,14
C NOW X(L) IS IN ROW IROW, BUT WE WANT IT IN ROW L.
   14 L=IXRCT(IROW)
      LL=IXRCT(L)
      IF(LL)15,16,15
C X(L) IS IN ROW IROW, WHILE X(LL) IS IN ROW L.
   15 IXRCT(IROW)=LL
      IXRCT(L)=L
      GO TO 17
C X(L) IS IN ROW IROW, WHILE Y(IYRCT(L)) IS IN ROW L.
   16 IXRCT(IROW)=0
      IYRCT(IROW)=IYRCT(L)
      IXRCT(L)=L
C NOW EXCHANGE THE CONTENTS OF ROWS IROW AND L.
   17 DO 18 J=1,NP1
        TEMP=V(IROW,J)
        V(IROW,J)=V(L,J)
        V(L,J)=TEMP
   18   CONTINUE
      IF(IXRCT(IROW))19,11,19
   19 IF(IXRCT(IROW)-IROW)14,11,14
C NOW IXRCT IS NO LONGER NEEDED, SO STORE THE PRESENT IYCCT
C IN IT.
   20 DO 21 I=1,N
        IXRCT(I)=IYCCT(I)
   21   CONTINUE
C END OF PHASE 1.
C
C THE FIRST N ROWS OF V GIVE THE XS IN TERMS OF CERTAIN
C YS.  THESE ROWS WILL NOT BE CHANGED BY LATER OPERATIONS.
C
C WE NOW ATTACK THE MAXIMIZATION PROBLEM BY LOOKING AT THE
C DUAL PROBLEM OF MINIMIZING A FORM GIVEN BY THE
C COEFFICIENTS IN V(N+1,N+1) THROUGH V(M,N+1) WITH SLACK
C TERMS IN THE BOTTOM ROW OF V.
C SEARCH ROW M+1 FOR A SIGNIFICANTLY NEGATIVE ELEMENT.  IF
C THERE ARE NONE, PROCEED TO THE ACTUAL MINIMIZATION
C PROBLEM.  STICK WITH COLUMN JJ UNTIL V(M+1,JJ).GE.-REA1.
      JJ=0
   22 JJ=JJ+1
      IF(JJ-N)1015,1015,1035
 1015 IF(V(MP1,JJ)+REA1)24,22,22
C
C WE HAVE V(M+1,JJ) SIGNIFICANTLY NEGATIVE.  SEARCH COLUMN
C JJ FOR A POSITIVE ELEMENT, TREATING A VERY SMALL V(I,J)
C AS A ZERO.  IF THERE ARE NO POSITIVE ELEMENTS THE DUAL
C CONSTRAINTS WERE INCONSISTENT, SO THE ORIGINAL PROBLEM WAS
C INCONSISTENT OR UNBOUNDED.
   24 I=N
      INAMP=0
   25 I=I+1
      IF(I-M)1016,1016,1020
 1016 IF(V(I,JJ)-REA)25,25,1017
C
C NOW V(I,JJ).GT.REA.  WE SEARCH ROW I FOR INDICES K SUCH
C THAT V(M+1,K).GE.0.0.OR.K.LT.JJ, AND V(I,K).LT.-REA, AND
C FIND THE MAXIMUM RATIO (I.E. THE RATIO WITH SMALLEST
C ABSOLUTE VALUE, IF V(M+1,K).GE.0.0) V(M+1,K)/V(I,K).  IF
C THERE IS NO SUCH K WE LOOK AT POSITIVE V(I,K) BELOW.
 1017 INDST=0
      DO 32 J=1,N
        IF(V(MP1,J))1018,28,28
 1018   IF(J-JJ)28,32,32
   28   IF(V(I,J)+REA)29,32,32
   29   DIST1=V(MP1,J)/V(I,J)
        IF(INDST)31,31,30
   30   IF(DIST1-DIST)32,32,31
   31   DIST=DIST1
        INDST=1
        K=J
   32   CONTINUE
      IF(INDST)35,35,33
C
C WE NOW COMPUTE V(I,JJ)*DIST AND GO ON TO LOOK AT OTHER
C ROWS TO MINIMIZE THIS QUANTITY (I.E. TO MAXIMIZE ITS
C ABSOLUTE VALUE, IF V(M+1,K).GE.0.0).  THIS IS THE NEGATIVE
C OF THE CHANGE IN V(M+1,JJ).
   33 BMPRV=V(I,JJ)*DIST
      IF(INAMP)34,34,1019
 1019 IF(BMPRV-AMPRV)34,25,25
   34 AMPRV=BMPRV
      INAMP=1
      KPMP1=I
      KPMP2=K
C (KPMP1,KPMP2) GIVES THE POSITION OF THE BEST PROSPECTIVE
C RESOLVENT FOUND SO FAR.
      GO TO 25
C
C IF THERE WAS NO INDEX K SUCH THAT V(M+1,K).GE.0.0.OR.K.LT.
C JJ, AND V(I,K).LT.-REA, WE LOOK FOR THE SMALLEST (I.E.
C LARGEST IN ABSOLUTE VALUE) RATIO V(M+1,K)/V(I,K) FOR
C V(I,K).GT.REA AND V(M+1,K).LT.0.0, AND PERFORM AN
C ELIMINATION WITH RESOLVENT V(I,K).  THERE IS AT LEAST ONE
C SUCH K, NAMELY JJ.
C THIS WILL FINISH PHASE 2 UNLESS BACKTRACKING IS NECESSARY.
   35 DIST=ONE
      DO 39 J=1,N
        IF(V(MP1,J))36,39,39
   36   IF(V(I,J)-REA)39,39,37
   37   DIST1=V(MP1,J)/V(I,J)
        IF(DIST1-DIST)38,39,39
   38   DIST=DIST1
        K=J
   39   CONTINUE
      GO TO 49
 1020 IF(INAMP)1021,1021,1023
C AT THIS POINT INAMP IS POSITIVE IFF THERE WAS AT LEAST ONE
C ELEMENT .GT. REA IN COLUMN JJ.  IF THERE WERE NONE, WE
C TEMPORARILY RELAX REA AND TRY AGAIN.
 1021 IF(IRLAX)1022,1022,41
 1022 IRLAX=1
      INDIC=-1
      REA=REA1
      GO TO 24
   41 INDIC=2
      RETURN
C CHECK TO SEE IF V(MP1,KPMP2) IS VERY SMALL IN ABSOLUTE
C VALUE OR NEGATIVE.  THIS INDICATES DEGENERACY.
 1023 IF(V(MP1,KPMP2)-REA)1024,1024,43
C DO AN ELIMINATION WITH RESOLVENT V(KPMP1,KPMP2).
   43 I=KPMP1
      K=KPMP2
      GO TO 49
C
C WE ARE NOW STUCK IN DEGENERATE COLUMN KPMP2.  WE SEARCH
C EACH DEGENERATE COLUMN IN WHICH WE ARE STUCK FOR A
C RESOLVENT WHICH WILL KEEP US FROM GETTING STUCK IN THIS
C COLUMN NEXT TIME, AND TO REDUCE THE ROUND-OFF ERROR WE
C TAKE THE SMALLEST OF THESE (I.E. LARGEST IN ABSOLUTE
C VALUE) AS OUR ACTUAL RESOLVENT.
 1024 AMIN=ONE
      DO 1034 J=1,N
C COLUMN J MAY BE DEGENERATE IF 0.0.LE.V(M+1,J).LE.REA,
C OR V(M+1,J).LT.0.0.AND.J.LT.JJ.
        IF(V(MP1,J))1025,1026,1026
 1025   IF(J-JJ)1027,1034,1034
 1026   IF(V(MP1,J)-REA)1027,1027,1034
C WE WILL BE STUCK IN COLUMN J IFF THERE IS AN INDEX ID FOR
C WHICH V(ID,JJ).GT.REA AND V(ID,J).LT.-REA.  IF THIS IS THE
C CASE, CHOOSING SUCH AN ID SO THAT V(ID,JJ)/V(ID,J) IS
C MINIMIZED (I.E. MAXIMIZED IN ABSOLUTE VALUE) AND TAKING
C V(ID,J) AS THE RESOLVENT WILL INSURE THAT WE DONT GET
C STUCK IN COLUMN J NEXT TIME.
 1027   DIST=ONE
        DO 1031 I=NP1,M
          IF(V(I,JJ)-REA)1031,1031,1028
 1028     IF(V(I,J)+REA)1029,1031,1031
 1029     DIST1=V(I,JJ)/V(I,J)
          IF(DIST1-DIST)1030,1031,1031
 1030     DIST=DIST1
          ID=I
 1031     CONTINUE
        IF(DIST-ONE/TWO)1032,1034,1034
C WE HAVE NOW DETERMINED THAT WE ARE STUCK IN COLUMN J.
C IF V(ID,J).LT.AMIN THEN V(ID,J) IS THE BEST RESOLVENT
C FOUND SO FAR.
 1032   IF(V(ID,J)-AMIN)1033,1034,1034
 1033   AMIN=V(ID,J)
        KPMP1=ID
        KPMP2=J
 1034   CONTINUE
C THE BEST RESOLVENT IS V(KPMP1,KPMP2), SO WE DO AN
C ELIMINATION.
      GO TO 43
   49 KTJOR=KTJOR+1
      IF(KTJOR-300)50,50,73
   50 CALL SJELIM(MP1,NP1,NP1,I,K,V)
      ITEMP=IYRCT(I)
      IYRCT(I)=IYCCT(K)
      IYCCT(K)=ITEMP
C RESET REA AND IRLAX.
      REA=REAKP
      IRLAX=0
C IF NOW V(M+1,JJ) HAS BEEN MADE NOT SIGNIFICANTLY NEGATIVE,
C WE GO TO THE NEXT COLUMN.  OTHERWISE WE TRY AGAIN IN
C COLUMN JJ.
      IF(V(MP1,JJ)+REA1)24,22,22
C IN THE UNLIKELY EVENT THAT SOME V(M+1,J) IS STILL VERY
C SIGNIFICANTLY NEGATIVE WE BACKTRACK TO COLUMN J.  THIS
C COULD NOT HAPPEN IF THERE WERE NO ROUNDOFF ERROR AND WE
C COULD ALLOW DIVISION BY NUMBERS WITH VERY SMALL ABSOLUTE
C VALUE.  OMIT BACKTRACKING IF IBACK=1.
 1035 IF(IBACK)1036,1036,51
 1036 DO 1038 J=1,N
        IF(V(MP1,J)+REA)1037,1037,1038
 1037   JJ=J
        GO TO 24
 1038   CONTINUE
C END OF PHASE 2.
C
   51 I=N
      KKK=0
C
C SEARCH FOR A SIGNIFICANTLY NEGATIVE ELEMENT BETWEEN
C V(N+1,N+1) AND V(N+1,M).  IF THERE ARE NONE WE HAVE THE
C MINIMAL POINT OF THE DUAL PROBLEM (AND THUS THE MAXIMAL
C POINT OF THE DIRECT PROBLEM) ALREADY.
   52 I=I+1
      IF(I-M)1039,1039,1043
 1039 IF(V(I,NP1)+REA2)1040,52,52
C
C SEARCH FOR A NEGATIVE ELEMENT IN ROW I, TREATING A NUMBER
C WHICH IS VERY SMALL IN ABSOLUTE VALUE AS A ZERO.  IF THERE
C ARE NO NEGATIVE ELEMENTS THE DUAL PROBLEM WAS UNBOUNDED
C BELOW, SO THE ORIGINAL CONSTRAINTS WERE INCONSISTENT.
C FIND THE INDEX K OF THE LARGEST (I.E. SMALLEST IN ABSOLUTE
C VALUE, IF V(M+1,K).GE.0.0) RATIO V(M+1,K)/V(I,K) WITH
C V(I,K).LT.-REA.
 1040 INDST=0
      DO 58 J=1,N
        IF(V(I,J)+REA)55,58,58
   55   DIST1=V(MP1,J)/V(I,J)
        IF(INDST)57,57,56
   56   IF(DIST1-DIST)58,58,57
   57   K=J
        INDST=1
        DIST=DIST1
   58   CONTINUE
      IF(INDST)1041,1041,60
C RELAX REA AND LOOK FOR NEGATIVE ELEMENTS WITH SMALLER
C ABSOLUTE VALUE.
 1041 IF(IRLAX)1042,1042,59
 1042 IRLAX=1
      INDIC=-1
      REA=REA1
      GO TO 1040
   59 INDIC=3
      RETURN
C
C COMPUTE THE IMPROVEMENT DIST*V(I,N+1) IN THE VALUE OF THE
C FORM USING V(I,K) AS THE RESOLVENT.  SET KKK=1 TO INDICATE
C A SIGNIFICANTLY NEGATIVE V(I,N+1) WAS FOUND, AND LOOK AT
C THE OTHER ROWS TO FIND THE RESOLVENT GIVING THE LARGEST
C IMPROVEMENT.
   60 BMPR2=DIST*V(I,NP1)
C RESET IRLAX SO THAT THE NEXT ROW WHICH NEEDS RELAXING DOES
C NOT TERMINATE THE ROUTINE.  REA WILL REMAIN RELAXED UNTIL
C AFTER THE NEXT ELIMINATION.
      IRLAX=0
      IF(KKK)62,62,61
   61 IF(BMPR2-AMPR2)52,52,62
   62 KKK=1
      KEEP=I
      KEEP1=K
      AMPR2=BMPR2
      GO TO 52
C KKK=0 HERE IFF NONE OF THE COST COEFFICIENTS ARE
C SIGNIFICANTLY NEGATIVE.
 1043 IF(KKK)1048,1044,1048
C CHECK TO SEE IF ANY OF THE NUMBERS IN THE BOTTOM ROW HAVE
C BECOME VERY SIGNIFICANTLY NEGATIVE.  IF SO, WE MUST
C BACKTRACK TO PHASE 2 (SEE COMMENT ABOVE STATEMENT 1035).
C OMIT BACKTRACKING IF IBACK=1.
 1044 IF(IBACK)1045,1045,74
 1045 DO 1047 J=1,N
        IF(V(MP1,J)+REA)1046,1046,1047
 1046   JJ=J
        GO TO 24
 1047   CONTINUE
      GO TO 74
C CHECK TO SEE IF V(MP1,KEEP1) IS VERY SMALL IN ABSOLUTE
C VALUE OR NEGATIVE.  THIS INDICATES DEGENERACY.
 1048 IF(V(MP1,KEEP1)-REA)1049,1049,65
C DO AN ELIMINATION WITH RESOLVENT V(KEEP,KEEP1).
   65 I=KEEP
      K=KEEP1
      GO TO 71
C
C WE ARE NOW STUCK IN DEGENERATE COLUMN KEEP1.  WE SEARCH
C EACH DEGENERATE COLUMN IN WHICH WE ARE STUCK FOR A
C RESOLVENT WHICH WILL KEEP US FROM GETTING STUCK IN THIS
C COLUMN NEXT TIME.  IF WE ARE NOT USING THE OPTION
C DESCRIBED IN THE COMMENTS PRECEDING STATEMENT 1055, WE
C TAKE THE SMALLEST OF THESE (I.E. THE LARGEST IN ABSOLUTE
C VALUE) AS OUR ACTUAL RESOLVENT IN ORDER TO REDUCE THE
C GROWTH OF ROUND-OFF ERROR.
 1049 AMIN=ONE
      MXRKN=NP1
      DO 1072 J=1,N
C COLUMN J MAY BE DEGENERATE IF V(M+1,J).LE.REA.
        IF(V(MP1,J)-REA)1050,1050,1072
C WE WILL BE STUCK IN COLUMN J IFF THERE IS AN INDEX ID FOR
C WHICH V(ID,N+1).LT.-REA2 AND V(ID,J).LT.-REA.  IF THIS
C IS THE CASE, CHOOSING SUCH AN ID SO THAT V(ID,N+1)/V(ID,J)
C IS MAXIMIZED AND TAKING V(ID,J) AS THE RESOLVENT WILL
C INSURE THAT WE DONT GET STUCK IN COLUMN J NEXT TIME.
 1050   DIST=-ONE
        DO 1054 I=NP1,M
          IF(V(I,NP1)+REA2)1051,1054,1054
 1051     IF(V(I,J)+REA)1052,1054,1054
 1052     DIST1=V(I,NP1)/V(I,J)
          IF(DIST1-DIST)1054,1054,1053
 1053     DIST=DIST1
          ID=I
 1054     CONTINUE
        IF(DIST+ONE/TWO)1072,1072,1055
C WE HAVE NOW DETERMINED THAT WE ARE STUCK IN COLUMN J.
C THE FOLLOWING STATEMENTS ATTEMPT TO BREAK DEGENERACY
C FASTER BY LOOKING ONE ITERATION INTO THE FUTURE, THAT IS,
C BY CHOOSING FROM THE PROSPECTIVE RESOLVENTS FOUND ABOVE
C THAT ONE WHICH MINIMIZES THE MINIMUM NUMBER OF STICKING
C PLACES IN ANY ROW AT THE NEXT STAGE.
C BECAUSE OF COMPUTER TIME AND THE POSSIBLE LOSS OF ACCURACY
C DUE TO LESSENED PIVOTING (EVEN THOUGH TIES ARE ALWAYS
C BROKEN IN FAVOR OF THE RESOLVENT WITH GREATEST ABSOLUTE
C VALUE), IT IS SUGGESTED THAT THIS OPTION BE OMITTED IF
C INFORMATION WAS AVAILABLE FROM A PREVIOUS VERTEX.  THIS
C WILL BE THE CASE IFF THE BACKTRACKING OPTION WAS USED,
C THAT IS, IFF IBACK=0.
 1055   IF(IBACK)1070,1070,1056
C COMPUTE WHAT THE NEW BOTTOM ROW WOULD BE (EXCEPT FOR
C POSITION J) IF V(ID,J) WERE USED AS THE RESOLVENT, AND
C PUT THE RESULTS INTO Y.
 1056   ROWQ=V(MP1,J)/V(ID,J)
        DO 1058 L=1,N
          IF(L-J)1057,1058,1057
 1057     Y(L)=V(MP1,L)-V(ID,L)*ROWQ
 1058     CONTINUE
        LRKNT=-1
C WE LOOK FOR A ROW WHICH WILL HAVE A SIGNIFICANTLY NEGATIVE
C LAST ELEMENT BUT A MINIMUM NUMBER OF PLACES WHERE WE WILL
C BE STUCK IN DEGENERATE COLUMNS.  LRKNT=-1 MEANS WE HAVE
C NOT YET FOUND A ROW WHICH WILL HAVE A SIGNIFICANTLY
C NEGATIVE LAST ELEMENT.
        DO 1068 II=NP1,M
          IF(II-ID)1059,1068,1059
 1059     ROWQ=V(II,J)/V(ID,J)
          RTCOL=V(II,NP1)-V(ID,NP1)*ROWQ
          IF(RTCOL+REA2)1060,1068,1068
C IF WE HAVE ALREADY LOCATED A RESOLVENT WHICH WILL FINISH
C THE ROUTINE, BUT THE PRESENT PROSPECTIVE RESOLVENT WOULD
C GIVE A ROW WITH A SIGNIFICANTLY NEGATIVE LAST ELEMENT, WE
C LOOK AT THE NEXT PROSPECTIVE RESOLVENT FOR PIVOTING
C PURPOSES.
 1060     IF(MXRKN+1)1061,1072,1061
 1061     LRKNT=0
C NOW COUNT THE NUMBER (LRKNT) OF STICKING PLACES IN ROW II
C AT THE NEXT ITERATION.
          DO 1065 JJ=1,N
            IF(JJ-J)1062,1065,1062
 1062       IF(Y(JJ)-REA)1063,1063,1065
 1063       IF(V(II,JJ)-V(ID,JJ)*ROWQ+REA)1064,1065,1065
 1064       LRKNT=LRKNT+1
            IF(LRKNT-MXRKN)1065,1065,1068
 1065       CONTINUE
          IF(LRKNT-MXRKN)1067,1066,1068
 1066     IF(V(ID,J)-AMIN)1067,1068,1068
 1067     MXRKN=LRKNT
          AMIN=V(ID,J)
          KEEP=ID
          KEEP1=J
 1068     CONTINUE
C LRKNT=-1 HERE WOULD MEAN THIS RESOLVENT WOULD FINISH THE
C ROUTINE.  IF LRKNT.GE.0 THEN MXRKN.GE.0 ALSO, SO WE WILL
C NOT HAVE EARLIER FOUND A RESOLVENT WHICH WILL FINISH THE
C ROUTINE.
        IF(LRKNT+1)1072,1069,1072
 1069   IF(MXRKN+1)1071,1070,1071
 1070   IF(V(ID,J)-AMIN)1071,1072,1072
 1071   MXRKN=-1
        AMIN=V(ID,J)
        KEEP=ID
        KEEP1=J
 1072   CONTINUE
C THE BEST RESOLVENT IS V(KEEP,KEEP1), SO WE DO AN
C ELIMINATION.
      GO TO 65
   71 KTJOR=KTJOR+1
      IF(KTJOR-300)72,72,73
   72 CALL SJELIM(MP1,NP1,NP1,I,K,V)
      ITEMP=IYRCT(I)
      IYRCT(I)=IYCCT(K)
      IYCCT(K)=ITEMP
C RESET REA AND IRLAX.
      REA=REAKP
      IRLAX=0
      GO TO 51
   73 INDIC=4
      RETURN
C END OF PHASE 3.  WE NOW HAVE THE VERTEX WE ARE SEEKING.
C
C READ OFF THE Y VALUES FOR THIS VERTEX.
   74 DO 75 J=1,N
        IYCJ=IYCCT(J)
        Y(IYCJ)=ZERO
   75   CONTINUE
      DO 76 I=NP1,M
        IYRI=IYRCT(I)
        Y(IYRI)=V(I,NP1)
   76   CONTINUE
C COMPUTE THE XS FROM THE YS.  RECALL THAT IXRCT CONTAINS
C THE FORMER IYCCT.
      DO 78 I=1,N
        X(I)=V(I,NP1)
        DO 77 J=1,N
          IXRJ=IXRCT(J)
          X(I)=X(I)-V(I,J)*Y(IXRJ)
   77     CONTINUE
   78   CONTINUE
C
C NOW PUT THE VALUES IN IYCCT INTO THE FIRST N POSITIONS OF
C IYRCT IN DECREASING ORDER SO THAT WHEN SLNPRO IS CALLED
C AGAIN THE YS AT THE PRESENT MINIMIZING VERTEX WILL BE
C SCANNED FIRST, BEGINNING WITH THOSE CORRESPONDING TO
C -1.0.LE.Q(J).LE.1.0.  THUS THESE WILL PROBABLY APPEAR IN
C THE INITIAL VERTEX AFTER EXCHANGE OF XS AND YS.
C TO ACCOMPLISH THIS, MAKE IXRCT(I)=-1 IF IYCCT(J)=I FOR
C SOME J, THEN SCAN IXRCT BACKWARDS.
      DO 79 J=1,N
        IYCJ=IYCCT(J)
        IXRCT(IYCJ)=-1
   79   CONTINUE
      K=1
      I=MP1
   80 I=I-1
      IF(I)83,83,81
   81 IF(IXRCT(I)+1)80,82,80
   82 IYRCT(K)=I
      K=K+1
      GO TO 80
C NOW FILL IN THE REST OF IYRCT BY SCANNING IXRCT AGAIN.
   83 I=MP1
   84 I=I-1
      IF(I)87,87,85
   85 IF(IXRCT(I))84,86,86
   86 IYRCT(K)=I
      K=K+1
      GO TO 84
   87 RETURN
      END
      SUBROUTINE SJELIM(L,LL,K,IR,IS,V)
C***BEGIN PROLOGUE  SJELIM
C***REFER TO  SLNPRO
C***ROUTINES CALLED  (NONE)
C***PURPOSE  THIS SUBROUTINE PERFORMS A MODIFIED JORDAN
C            ELIMINATION ON THE L-LL+1 BY K MATRIX
C            CONSISTING OF ROWS LL THROUGH L OF V AND
C            COLUMNS 1 THROUGH K OF V.  THE RESOLVENT
C            IS V(IR,IS).
C***END PROLOGUE  SJELIM
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION V(635,17)
C
C SET PRECISION DEPENDENT CONSTANTS FOR SJELIM.
C***FIRST EXECUTABLE STATEMENT  SJELIM
      ONE=1.0D0
C EMD OF SETTING PRECISION DEPENDENT CONSTANTS FOR SJELIM.
C
C DIVIDE THE ENTRIES IN THE RESOLVENT ROW (EXCEPT FOR THE
C RESOLVENT) BY THE RESOLVENT.
      RESOL=V(IR,IS)
      DO 2 J=1,K
        IF(J-IS)1001,2,1001
 1001   V(IR,J)=V(IR,J)/RESOL
    2   CONTINUE
C SWEEP OUT IN ALL BUT ROW IR AND COLUMN IS.
      DO 6 I=LL,L
        IF(I-IR)1002,6,1002
 1002   FACT=-V(I,IS)
        DO 5 J=1,K
          IF(J-IS)1003,5,1003
 1003     V(I,J)=V(I,J)+V(IR,J)*FACT
    5     CONTINUE
    6   CONTINUE
C DIVIDE THE ENTRIES IN THE RESOLVENT COLUMN (EXCEPT FOR THE
C RESOLVENT) BY THE NEGATIVE OF THE RESOLVENT.
      DO 8 I=LL,L
        IF(I-IR)1004,8,1004
 1004   V(I,IS)=-V(I,IS)/RESOL
    8   CONTINUE
C REPLACE THE RESOLVENT BY ITS RECIPROCAL.
      V(IR,IS)=ONE/RESOL
      RETURN
      END
      SUBROUTINE SSETV(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,
     *INDLW,ALTBL,INDUP,UPTBL,INDF,WTBLE,MEQ,MINEQ,ALC,MODE,V,M)
C***BEGIN PROLOGUE  SSETV
C***REFER TO  CDFCOR
C***ROUTINES CALLED  (NONE)
C***PURPOSE  THIS SUBROUTINE SETS UP THE LINEAR PROGRAMMING PROBLEM
C            FOR THE INITIAL (IF MODE=0) OR MAIN (IF MODE=1)
C            ITERATIONS OF THE DIFFERENTIAL CORRECTION ALGORITHM.
C***END PROLOGUE  SSETV2
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION FTBWT(101),PWR(101,15),V(635,17),WTBLE(101),
     *   UPTBL(101),ALTBL(101),INDUP(101),INDLW(101),
     *   INDF(101),ALC(101,16)
C
C IN THE INITIALIZATION WE MINIMIZE W WITH THE RESTRICTIONS
C   ABS(F*Q-WT*P) .LE. W AT ALL GRID POINTS WHERE F IS TO BE
C   APPROXIMATED,
C   Q .GE. EPS AT ALL GRID POINTS,
C   ABS(Q(J)) .LE. 1.0 FOR ALL J,
C   PLUS THE OTHER CONSTRAINTS (IF ANY).
C IN THE MAIN ITERTATIONS WE MINIMIZE W WITH THE RESTRICTIONS
C   (ABS(F*Q-WT*P)-DELTK*Q)/QK .LE. W AT ALL GRID POINTS WHERE
C   F IS TO BE APPROXIMATED,
C   ABS(Q(J)) .LE. 1.0 FOR ALL J,
C   PLUS THE OTHER CONSTRAINTS (IF ANY).
C EACH ABSOLUTE VALUE CONSTRAINT IS BROKEN INTO TWO CONSTRAINTS
C TO ACCOUNT FOR THE ABSOLUTE VALUE.  IN THE MAIN ITERATIONS
C THE SUBROUTINE USES THE FACT THAT THE VALUE OF QK AT GRID
C POINT I IS STORED IN V(2*I,N+1) AND DELTK IS IN
C V(2*NUMGR+1,N+1) FROM A PREVIOUS CALL TO SPQERD.
C
C***FIRST EXECUTABLE STATEMENT  SSETV2
      IONE=IOPT-(IOPT/10)*10
C
C SET MACHINE AND PRECISION DEPENDENT CONSTANTS FOR
C SUBROUTINE SSETV.
      ONE=1.0D0
      ZERO=ONE-ONE
      TWO=ONE+ONE
      FOUR=TWO+TWO
      TEN=FOUR+FOUR+TWO
C IF WE ARE IN THE MAIN PART OF DIFFERENTIAL CORRECTION AND
C IONE IS ODD WE SET EPS=EPSQ.  IN THIS CASE WE WILL REQUIRE
C THAT THE DENOMINATOR BE .GE. EPS AT EVERY GRID POINT.
      IF(MODE)300,300,100
  100 IF(IONE-(IONE/2)*2)800,800,200
  200 EPS=EPSQ
      GO TO 800
C HERE WE ARE IN THE INITIALIZATION AND WE AT LEAST TEMPORARILY
C COMPUTE EPS=MAX(0.0001,1000.0*SPCMN)
  300 SPCMN=D1MACH(3)
      EPS=SPCMN*TEN**3
      IF(EPS-TEN**(-4))400,500,500
  400 EPS=TEN**(-4)
C IF IONE IS ODD HERE WE INCREASE EPS TO EPSQ IF NECESSARY.
  500 IF(IONE-(IONE/2)*2)800,800,600
  600 IF(EPS-EPSQ)700,800,800
  700 EPS=EPSQ
C END OF SETTING MACHINE AND PRECISION DEPENDENT CONSTANTS
C FOR SSETV.
C
  800 NMPP1=NUMP+1
      NMPPQ=NUMP+NUMQ
      N=NMPPQ+1
      NP1=N+1
      DELTK=V(2*NUMGR+1,NP1)
C M WILL KEEP TRACK OF WHICH ROW OF V IS BEING SET.  AT THE
C CONCLUSION OF SSETV M WILL BE THE TOTAL NUMBER OF
C CONSTRAINTS.
      M=0
C
C SET THE UPPER CONSTRAINTS, IF THERE ARE ANY, EXCEPT FOR
C THE LAST COLUMN.  THEY WILL BE OF THE FORM P-Q*UP.LE.0.0
      IF(IOPT-(IOPT/1000)*1000-100)1400,900,900
  900 DO 1300 I=1,NUMGR
C INDUP(I)=1 IF THERE IS AN UPPER CONSTRAINT AT GRID POINT
C I, AND INDUP(I)=0 OTHERWISE.
        IF(INDUP(I))1300,1300,1000
 1000   M=M+1
        DO 1100 J=1,NUMP
          V(M,J)=PWR(I,J)
 1100     CONTINUE
        DO 1200 J=NMPP1,NMPPQ
          V(M,J)=-UPTBL(I)*PWR(I,J)
 1200     CONTINUE
        V(M,N)=ZERO
 1300   CONTINUE
C
C SET THE LOWER CONSTRAINTS, IF THERE ARE ANY, EXCEPT FOR
C THE LAST COLUMN.  THEY WILL BE OF THE FORM -P+Q*LW.LE.0.0
 1400 IF(IOPT-(IOPT/100)*100-10)2000,1500,1500
 1500 DO 1900 I=1,NUMGR
C INDLW(I)=1 IF THERE IS A LOWER CONSTRAINT AT GRID POINT I,
C AND INDLW(I)=0 OTHERWISE.
        IF(INDLW(I))1900,1900,1600
 1600   M=M+1
        DO 1700 J=1,NUMP
          V(M,J)=-PWR(I,J)
 1700     CONTINUE
        DO 1800 J=NMPP1,NMPPQ
          V(M,J)=ALTBL(I)*PWR(I,J)
 1800     CONTINUE
        V(M,N)=ZERO
 1900   CONTINUE
C
 2000 IF(MODE)2100,2100,2700
C
C HERE WE ARE IN THE INITIALIZATION PHASE, AND FOR EACH GRID
C POINT AT WHICH F IS TO BE APPROXIMATED WE SET THE CONSTRAINTS
C -WT*P+F*Q-W .LE. 0.0 AND WT*P-F*Q-W .LE. 0.0
 2100 DO 2500 I=1,NUMGR
C INDF(I)=1 IF F IS TO BE APPROXIMATED AT GRID POINT I, AND
C INDF(I)=0 IF F IS TO BE IGNORED AT GRID POINT I.
        IF(INDF(I))2500,2500,2200
 2200   M=M+2
        MM1=M-1
        DO 2300 J=1,NUMP
          V(MM1,J)=-WTBLE(I)*PWR(I,J)
          V(M,J)=-V(MM1,J)
 2300     CONTINUE
        DO 2400 J=NMPP1,NMPPQ
          V(MM1,J)=FTBWT(I)*PWR(I,J)
          V(M,J)=-V(MM1,J)
 2400     CONTINUE
        V(MM1,N)=-ONE
        V(M,N)=-ONE
 2500   CONTINUE
C SINCE WE ARE IN THE INITIALIZATION, COLUMN N+1 IS NOT NEEDED
C FOR STORAGE OF QK, SO ZERO IT OUT FOR THOSE ROWS ALREADY SET.
      DO 2600 I=1,M
        V(I,NP1)=ZERO
 2600   CONTINUE
      GO TO 3300
C
C HERE WE ARE IN THE MAIN DIFFERENTIAL CORRECTION PHASE, AND
C FOR EACH GRID POINT AT WHICH F IS TO BE APPROXIMATED SET
C THE CONSTRAINTS -WT*P+(-DELTK+F)*Q-QK*W.LE.0.0
C AND WT*P+(-DELTK-F)*Q-QK*W.LE.0.0
C (EXCEPT FOR THE LAST COLUMN).
 2700 DO 3100 I=1,NUMGR
C INDF(I)=1 IF F IS TO BE APPROXIMATED AT GRID POINT I,
C AND INDF(I)=0 IF F IS TO BE IGNORED AT GRID POINT I.
        IF(INDF(I))3100,3100,2800
 2800   M=M+2
        MM1=M-1
        DO 2900 J=1,NUMP
          V(MM1,J)=-WTBLE(I)*PWR(I,J)
          V(M,J)=-V(MM1,J)
 2900     CONTINUE
        DEL1=-DELTK+FTBWT(I)
        DEL2=-DELTK-FTBWT(I)
        DO 3000 J=NMPP1,NMPPQ
          V(MM1,J)=DEL1*PWR(I,J)
          V(M,J)=DEL2*PWR(I,J)
 3000     CONTINUE
C QK AT GRID POINT I IS STORED IN V(2*I,N+1) FROM AN
C EARLIER CALL TO SPQERD.
        IDB=2*I
        V(MM1,N)=-V(IDB,NP1)
        V(M,N)=-V(IDB,NP1)
 3100   CONTINUE
C NOW THE LAST COLUMN IS NO LONGER NEEDED FOR THE STORAGE OF
C QK, SO ZERO IT OUT FOR THOSE ROWS ALREADY SET.
      DO 3200 I=1,M
        V(I,NP1)=ZERO
 3200   CONTINUE
C IF THE ONES DIGIT OF IOPT IS ODD HERE WE FORCE Q .GE. EPS.
      IF(IONE-(IONE/2)*2)3700,3700,3300
C
C HERE WE FORCE Q .GE. EPS.
 3300 DO 3600 I=1,NUMGR
        M=M+1
        DO 3400 J=1,N
          V(M,J)=ZERO
 3400     CONTINUE
        DO 3500 J=NMPP1,NMPPQ
          V(M,J)=-PWR(I,J)
 3500     CONTINUE
        V(M,NP1)=-EPS
 3600   CONTINUE
C
C NOW INCLUDE THE ADDITIONAL INEQUALITY CONSTRAINTS FROM ALC,
C IF THERE ARE ANY.
 3700 IF(MINEQ)4100,4100,3800
 3800 IALC=MEQ
      DO 4000 I=1,MINEQ
        M=M+1
        IALC=IALC+1
        DO 3900 J=1,NMPPQ
          V(M,J)=ALC(IALC,J)
 3900     CONTINUE
        V(M,N)=ZERO
        V(M,NP1)=ALC(IALC,N)
 4000   CONTINUE
C
C NOW SET THE CONSTRAINTS OF THE FORM -Q(J).LE.1.0
C AND Q(J).LE.1.0
 4100 DO 4300 I=1,NUMQ
        M=M+2
        MM1=M-1
        DO 4200 J=1,N
          V(MM1,J)=ZERO
          V(M,J)=ZERO
 4200     CONTINUE
        NNN=NUMP+I
        V(MM1,NNN)=-ONE
        V(M,NNN)=ONE
        V(MM1,NP1)=ONE
        V(M,NP1)=ONE
 4300   CONTINUE
C
C NOW SET THE BOTTOM ROW.  TO MINIMIZE W WE MAXIMIZE -W.
      MP1=M+1
      DO 4400 J=1,NP1
        V(MP1,J)=ZERO
 4400   CONTINUE
      V(MP1,N)=ONE
      RETURN
      END
      SUBROUTINE EQELM(NVAR,MEQ,TOLEL,EQ,NSCOL,IRESL,MEQA,IFLAG)
C***BEGIN PROLOGUE  EQELM
C***REFER TO  CDFCOR
C***ROUTINES CALLED  (NONE)
C***PURPOSE  THIS SUBROUTINE ATTEMPTS TO SOLVE FOR SOME OF THE
C            VARIABLES IN THE MEQ EQUALITY CONSTRAINTS IN TERMS
C            OF THE OTHER VARIABLES.
C***END PROLOGUE  EQELM
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION EQ(15,16),NSCOL(15),IRESL(30)
C
C GIVEN THE EQUATIONS
C   EQ(I,1)*X(1)+...+EQ(I,NVAR)*X(NVAR) = EQ(I,NVAR+1)
C FOR I=1,...,MEQ, WE CHOOSE RESOLVENTS BY TOTAL PIVOTING IN THE
C LEFT MEQ BY NVAR PORTION OF EQ, STORE THE LTH RESOLVENT
C POSITION IN (IRESL(2*L-1),IRESL(2*L)), AND APPLY ROW
C OPERATIONS TO MAKE 1.0 IN THE RESOLVENT POSITION AND 0.0
C ELSEWHERE IN THE RESOLVENT COLUMN.  ONLY RESOLVENTS WITH ABSOLUTE
C VALUE .GT. TOLEL ARE ACCEPTED.  MEQA (WHICH STANDS FOR MEQ ACTUAL)
C WILL BE THE NUMBER OF ACCEPTABLE RESOLVENTS FOUND.  IF
C MEQA .LT. MEQ, AND IF SOME ROW NOT CONTAINING A RESOLVENT
C HAS LAST ELEMENT .GT. TOLEL IN ABSOLUTE VALUE, WE SET
C IFLAG=7777 AS AN INCONSISTENCY WARNING AND RETURN.  IF
C MEQA .LT. MEQ BUT THERE IS NO INCONSISTENCY THEN THE MEQ-MEQA
C ROWS OF EQ NOT CONTAINING RERSOLVENTS WILL BE IGNORED (THIS IS
C THE REDUNDANT CONSTRAINTS SITUATION).
C NSCOL (WHICH STANDS FOR NON SOLVED FOR COLUMN) IS DEFINED BY
C NSCOL(K) = THE SUBSCRIPT OF THE KTH UNSOLVED FOR VARIABLE
C (THAT IS, VARIABLE WITH NO RESOLVENT IN ITS COLUMN) FOR K=1,...,
C NVAR-MEQA.  THUS FOR K=1,...,MEQA WE WILL HAVE
C   X(IRESL(2*K)) = -EQ(IRESL(2*K-1),NSCOL(1))*X(NSCOL(1))-...-
C   EQ(IRESL(2*K-1),NSCOL(NVAR-MEQA))*X(NSCOL(NVAR-MEQA))+
C   EQ(IRESL(2*K-1),NVAR+1).
C
C SET PRECISION DEPENDENT CONSTANTS FOR EQELM.
C***FIRST EXECUTABLE STATEMENT  EQELM
      ONE=1.0D0
      ZERO=ONE-ONE
C END OF SETTING PRECISION DEPENDENT CONSTANTS FOR EQELM.
C
      N=NVAR+1
      NRESL=0
      MEQA=MEQ
C IF MEQ IS NOT POSITIVE THERE IS NOTHING TO DO HERE AND WE RETURN.
      IF(MEQ)30,30,70
   30 RETURN
C ZERO OUT NSCOL.  LATER EACH ELEMENT OF NSCOL CORRESPONDING TO A
C COLUMN CONTAINING A RESOLVENT WILL BE SET TO 1, AND NSCOL WILL
C BE READJUSTED TO ITS FINAL STATE AT THE END OF THIS SUBROUTINE
C IF THE EQUATIONS IN EQ ARE CONSISTENT.
   70 DO 100 J=1,NVAR
        NSCOL(J)=0
  100   CONTINUE
C
C WE ATTEMPT TO FIND MEQ RESOLVENTS USING TOTAL PIVOTING.
      DO 1400 L=1,MEQ
        AMAX=ZERO
        DO 800 I=1,MEQ
          DO 700 J=1,NVAR
            IF(NRESL)500,500,200
C WHILE SEARCHING, IGNORE ROWS AND COLUMNS CONTAINING AN
C EARLIER RESOLVENT.
  200       DO 400 K=1,NRESL
              IF(I-IRESL(2*K-1))300,800,300
  300         IF(J-IRESL(2*K))400,700,400
  400         CONTINUE
C HERE LOCATION (I,J) IS NOT IN A ROW OR COLUMN WITH AN
C EARLIER RESOLVENT.
  500       IF(ABS(EQ(I,J))-AMAX)700,700,600
  600       AMAX=ABS(EQ(I,J))
            IMAX=I
            JMAX=J
  700       CONTINUE
  800     CONTINUE
        IF(AMAX-TOLEL)1500,1500,900
C
C HERE WE HAVE FOUND AN ACCEPTABLE RESOLVENT.
  900   NRESL=NRESL+1
        IRESL(2*NRESL-1)=IMAX
        IRESL(2*NRESL)=JMAX
        NSCOL(JMAX)=1
C NOW DO ROW OPERATIONS TO MAKE 1.0 IN THE (IMAX,JMAX) POSITION
C AND 0.0 ELSEWHERE IN COLUMN JMAX (THUS SOLVING FOR X(JMAX)
C IN TERMS OF THE OTHER VARIABLES).
C FIRST DIVIDE ROW IMAX BY THE RESOLVENT.
        RESOL=EQ(IMAX,JMAX)
        DO 1000 J=1,N
          EQ(IMAX,J)=EQ(IMAX,J)/RESOL
 1000     CONTINUE
        DO 1300 I=1,MEQ
          IF(I-IMAX)1100,1300,1100
 1100     FACT=EQ(I,JMAX)
          DO 1200 J=1,N
            EQ(I,J)=EQ(I,J)-FACT*EQ(IMAX,J)
 1200       CONTINUE
 1300     CONTINUE
 1400   CONTINUE
      GO TO 2100
C
C HERE WE WERE UNABLE TO FIND AN ACCEPTABLE RESOLVENT.
C RESET MEQA TO THE NUMBER OF RESOLVENTS FOUND AND CHECK TO
C SEE IF ANY OF THE RIGHT SIDES OF ROWS NOT CONTAINING A
C RESOLVENT ARE .GT. TOLEL IN ABSOLUTE VALUE.  IF SO, THE
C EQUATIONS IN EQ WERE APPARENTLY INCONSISTENT, AND IF NOT,
C THEY WERE APPARENTLY REDUNDANT.
 1500 MEQA=NRESL
      DO 1900 I=1,MEQ
        IF(NRESL)1800,1800,1600
 1600   DO 1700 K=1,NRESL
          IF(I-IRESL(2*K-1))1700,1900,1700
 1700     CONTINUE
C HERE ROW I DOES NOT CONTAIN A RESOLVENT.
 1800   IF(ABS(EQ(I,N))-TOLEL)1900,1900,2000
 1900   CONTINUE
      GO TO 2100
C
C HERE THE CONSTRAINTS WERE INCONSISTENT.
 2000 IFLAG=7777
      RETURN
C
C HERE THE CONSTRAINTS WERE CONSISTENT, AND WE SET NSCOL AND
C RETURN.
 2100 IARG=0
      DO 2300 J=1,NVAR
        IF(NSCOL(J))2200,2200,2300
 2200   IARG=IARG+1
        NSCOL(IARG)=J
 2300   CONTINUE
      RETURN
      END
      SUBROUTINE SCRCH(M,NVAR,MEQA,NSCOL,IRESL,EQ,V)
C***BEGIN PROLOGUE  SCRCH
C***REFER TO  CDFCOR
C***ROUTINES CALLED  (NONE)
C***PURPOSE  THIS SUBROUTINE USES THE NONREDUNDANT EQUALITY
C            CONSTRAINTS TO ELIMINATE MEQA VARIABLES, THUS
C            SCRUNCHING V FROM MP1 BY NP1 TO MP1 BY NP1-MEQA.
C***END PROLOGUE  SCRCH
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION NSCOL(15),IRESL(30),EQ(15,16),V(635,17)
C
C IF MEQA=0 THERE IS NO SCRUNCHING TO BE DONE, AND WE RETURN.
C***FIRST EXECUTABLE STATEMENT  SCRCH
      IF(MEQA)100,100,200
  100 RETURN
C
  200 N=NVAR+1
      NP1=N+1
      MP1=M+1
      NS=NVAR-MEQA
C
      DO 800 L=1,MEQA
C THE RESOLVENT IS AT POSITION (LROW,LCOL) OF EQ.
        LROW=IRESL(2*L-1)
        LCOL=IRESL(2*L)
C IF ALL THE VARIABLES X(1),...,X(NVAR) ARE SOLVED FOR WE GO
C TO COLUMN NP1.
        IF(NS)600,600,300
  300   DO 500 J=1,NS
          JCOL=NSCOL(J)
          FACT=EQ(LROW,JCOL)
C WE REPLACE COLUMN JCOL OF V BY ITSELF MINUS EQ(LROW,JCOL)
C TIMES COLUMN LCOL OF V.
          DO 400 I=1,MP1
            V(I,JCOL)=V(I,JCOL)-FACT*V(I,LCOL)
  400       CONTINUE
  500     CONTINUE
C
C NOW DO COLUMN NP1 OF V.
  600   FACT=EQ(LROW,N)
        DO 700 I=1,MP1
          V(I,NP1)=V(I,NP1)-FACT*V(I,LCOL)
  700     CONTINUE
  800   CONTINUE
C
C NOW COUMNS IRESL(2),...,IRESL(2*MEQA) OF V ARE NO LONGER
C NEEDED, AND WE COMPACTIFY V.
      IF(NS)1300,1300,900
  900 DO 1200 J=1,NS
        JCOL=NSCOL(J)
C COLUMN JCOL OF V IS NOW MOVED TO THE COLUMN J POSITION IF IT
C IS NOT ALREADY THERE.
        IF(JCOL-J)1000,1200,1000
 1000   DO 1100 I=1,MP1
          V(I,J)=V(I,JCOL)
 1100     CONTINUE
 1200   CONTINUE
C NOW MOVE THE LAST TWO COLUMNS OF V TO THE LEFT MEQA SPACES.
C NOTE THAT THE VARIABLE W CORRESPONDING TO COLUMN N = NVAR+1
C OF V IS AUTOMATICALLY AN UNSOLVED FOR VARIABLE SINCE IT IS
C NOT INCLUDED IN EQ.
 1300 NMEQ=N-MEQA
      NPMEQ=NMEQ+1
      DO 1400 I=1,MP1
        V(I,NMEQ)=V(I,N)
        V(I,NPMEQ)=V(I,NP1)
 1400   CONTINUE
      RETURN
      END
      SUBROUTINE EXPND(NVAR,MEQA,NSCOL,IRESL,EQ,X)
C***BEGIN PROLOGUE  EXPND
C***REFER TO  CDFCOR
C***ROUTINES CALLED  (NONE)
C***PURPOSE  THIS SUBROUTINE COMPUTES THE VALUES OF THE MEQA
C            SOLVED FOR VARIABLES AND SPLICES THEM INTO X.
C***END PROLOGUE  EXPND
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION NSCOL(15),IRESL(30),EQ(15,16),X(16)
C
C IF MEQA=0 NO EXPANDING NEEDS TO BE DONE AND WE RETURN.
      IF(MEQA)100,100,200
  100 RETURN
C
  200 N=NVAR+1
C SHIFT THE LAST ELEMENT IN X INTO ITS PROPER PLACE.  NOTE
C THAT THIS ELEMENT CORRESPONDS TO THE VARIABLE W, WHICH IS
C UNSOLVED FOR SINCE IT IS NOT INCLUDED IN EQ.
      NMEQ=N-MEQA
      X(N)=X(NMEQ)
      NS=NVAR-MEQA
C SHIFT THE NS UNSOLVED FOR VARIABLES (NOT COUNTING W) INTO THE
C PROPER PLACES IN X, IF THERE ARE ANY.
      IF(NS)500,500,300
  300 DO 400 J=1,NS
        JJ=NS-J+1
        X(NSCOL(JJ))=X(JJ)
  400   CONTINUE
C
C NOW COMPUTE THE MEQA SOLVED FOR VARIABLES.
  500 DO 800 L=1,MEQA
C THE LTH RESOLVENT POSITION IN EQ IS (LROW,LCOL).  WE
C COMPUTE X(LCOL).
        LROW=IRESL(2*L-1)
        LCOL=IRESL(2*L)
        X(LCOL)=EQ(LROW,N)
        IF(NS)800,800,600
  600   DO 700 J=1,NS
          JCOL=NSCOL(J)
          X(LCOL)=X(LCOL)-EQ(LROW,JCOL)*X(JCOL)
  700     CONTINUE
  800   CONTINUE
      RETURN
      END
      SUBROUTINE NRMLZ(NUMP,NUMQ,MEQ,MINEQ,NSCOL,IRESL,ALC,
     *TOLNR,X)
C***BEGIN PROLOGUE  NRMLZ
C***REFER TO  CDFCOR
C***ROUTINES CALLED  (NONE)
C***PURPOSE  THIS SUBROUTINE ATTEMPTS TO NORMALIZE X SO THAT
C            MAX(ABS(Q(J)))=1.0, OR SOME INEQUALITY CONSTRAINT
C            IN ALC WITH LEFT AND RIGHT SIDES .GT. TOLNR HOLDS
C            WITH EQUALITY.  THE NORMALIZATION IS SKIPPED IF IT
C            WOULD REQUIRE MULTIPLYING BY A FACTOR LESS THAN 1.0,
C            OR IF SOME EQUALITY CONSTRAINT HAS RIGHT SIDE WITH
C            ABSOLUTE VALUE .GT. TOLNR.
C***END PROLOGUE  NRMLZ
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION NSCOL(15),IRESL(30),ALC(101,16),X(16)
C
C SET PRECISION DEPENDENT CONSTANTS FOR NRMLZ.
C***FIRST EXECUTABLE STATEMENT  NRMLZ
      ONE=1.0D0
C END OF SETTING PRECISION DEPENDENT CONSTANTS FOR NRMLZ.
C
      NMPPQ=NUMP+NUMQ
      N=NMPPQ+1
      IF(MEQ)300,300,100
C IF THE RIGHT SIDE OF SOME EQUALITY CONSTRAINT IS .GT. TOLNR,
C WE SKIP THE NORMALIZATION.
  100 DO 200 I=1,MEQ
        IF(ABS(ALC(I,N))-TOLNR)200,200,1500
  200   CONTINUE
C
C COMPUTE AMAX = MAX(TOLNR,ABS(Q(1)),...,ABS(Q(NUMQ))).
  300 AMAX=TOLNR
      DO 500 J=1,NUMQ
        NPJ=NUMP+J
        IF(ABS(X(NPJ))-AMAX)500,500,400
  400   AMAX=ABS(X(NPJ))
  500   CONTINUE
      FACT=ONE/AMAX
      IF(MINEQ)1200,1200,600
C
C CHECK THE INEQUALITY CONSTRAINTS IN ALC WITH SIGNIFICANTLY
C POSITIVE LEFT AND RIGHT SIDES.  MULTIPLYING ALL COEFFICIENTS
C BY A FACTOR .GT. 1.0 WILL NOT AFFECT INEQUALITY CONSTRAINTS
C WITH LEFT AND RIGHT SIDES OF OPPOSITE SIGN OR ZERO RIGHT SIDES
C AND CAN ONLY HELP INEQUALITY CONSTRAINTS WITH NEGATIVE LEFT
C AND RIGHT SIDES.
  600 DO 1100 I=1,MINEQ
        MPI=MEQ+I
        IF(ALC(MPI,N)-TOLNR)1100,1100,700
  700   ALFT=ALC(MPI,1)*X(1)
        DO 800 J=2,NMPPQ
          ALFT=ALFT+ALC(MPI,J)*X(J)
  800     CONTINUE
        IF(ALFT-TOLNR)1100,1100,900
  900   IF(ALC(MPI,N)/ALFT-FACT)1000,1100,1100
 1000   FACT=ALC(MPI,N)/ALFT
 1100   CONTINUE
C
C IF FACT .GT. 1.0, MULTIPLY ALL COEFFICIENTS BY FACT.
 1200 IF(FACT-ONE)1500,1500,1300
 1300 DO 1400 J=1,NMPPQ
        X(J)=FACT*X(J)
 1400   CONTINUE
 1500 RETURN
      END
      FUNCTION ICNCK(IOPT,NUMP,NUMQ,NUMGR,INDLW,ALTBL,INDUP,
     *UPTBL,EPSQ,X,V,MEQ,MINEQ,ALC,EPRES)
C***BEGIN PROLOGUE  ICNCK
C***REFER TO  CDFCOR
C***ROUTINES CALLED  (NONE)
C***PURPOSE  THIS SUBPROGRAM CHECKS ALL CONSTRAINTS, RETURNING
C            WITH VALUE 0 IF ALL ARE OK AND VALUE 1 IF ANY ARE
C            NOT SATISFIED WITHIN EPRES.
C***END PROLOGUE  ICNCK
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION INDLW(101),ALTBL(101),INDUP(101),UPTBL(101),
     *X(16),V(635,17),ALC(101,16)
C
C SET PRECISION DEPENDENT CONSTANTS FOR FUNCTION ICNCK.
C***FIRST EXECUTABLE STATEMENT  ICNCK
      ONE=1.0D0
C END OF SETTING PRECISION DEPENDENT CONSTANTS FOR ICNCK.
C
      NMPPQ=NUMP+NUMQ
      N=NMPPQ+1
      NP1=N+1
      IONE=IOPT-(IOPT/10)*10
C
C IF V(2*NUMGR+1,N+1) IS NEGATIVE, THEN SPQERD DETECTED A POINT
C WHERE THE DENOMINATOR WAS VERY SMALL OR NEGATIVE, SO WE SET
C ICNCK=1 AND RETURN.
      IF(V(2*NUMGR+1,NP1))2200,100,100
C
C CHECK LOWER RESTRICTED RANGE CONSTRAINTS, IF THERE ARE ANY.
  100 IF(IOPT-(IOPT/100)*100-10)500,200,200
  200 DO 400 I=1,NUMGR
        IF(INDLW(I))400,400,300
  300   IF(V(2*I-1,NP1)/V(2*I,NP1)-ALTBL(I)+EPRES)2200,400,400
  400   CONTINUE
C
C CHECK UPPER RESTRICTED RANGE CONSTRAINTS, IF THERE ARE ANY.
  500 IF(IOPT-(IOPT/1000)*1000-100)900,600,600
  600 DO 800 I=1,NUMGR
        IF(INDUP(I))800,800,700
  700   IF(V(2*I-1,NP1)/V(2*I,NP1)-UPTBL(I)-EPRES)800,800,2200
  800   CONTINUE
C
C CHECK DENOMINATOR LOWER BOUND CONSTRAINTS, IF THERE ARE ANY.
  900 IF(IONE-(IONE/2)*2)1200,1200,1000
 1000 DO 1100 I=1,NUMGR
        IF(V(2*I,NP1)-EPSQ+EPRES)2200,1100,1100
 1100   CONTINUE
C
C CHECK ADDITIONAL EQUALITY CONSTRAINTS, IF THERE ARE ANY.
 1200 IF(MEQ)1600,1600,1300
 1300 DO 1500 L=1,MEQ
        CHK=-ALC(L,N)
        DO 1400 J=1,NMPPQ
          CHK=CHK+ALC(L,J)*X(J)
 1400     CONTINUE
        IF(ABS(CHK)-EPRES)1500,1500,2200
 1500   CONTINUE
C
C CHECK ADDITIONAL INEQUALITY CONSTRAINTS, IF THERE ARE ANY.
 1600 IF(MINEQ)2000,2000,1700
 1700 DO 1900 L=1,MINEQ
        MEQL=MEQ+L
        CHK=-ALC(MEQL,N)
        DO 1800 J=1,NMPPQ
          CHK=CHK+ALC(MEQL,J)*X(J)
 1800     CONTINUE
        IF(CHK-EPRES)1900,1900,2200
 1900   CONTINUE
C
C CHECK THE DENOMINATOR COEFFICIENT BOUNDS.
 2000 OPE=ONE+EPRES
      DO 2100 J=1,NUMQ
        NPJ=NUMP+J
        IF(ABS(X(NPJ))-OPE)2100,2100,2200
 2100   CONTINUE
C
C IF WE REACH THIS SPOT, THEN ALL CONSTRAINTS ARE OK WITHIN EPRES.
      ICNCK=0
      RETURN
C
C HERE AT LEAST ONE CONSTRAINT WAS NOT SATISFIED WITHIN EPRES.
 2200 ICNCK=1
      RETURN
      END
    1    2    3    1  101
      1011
                 0.0                 2.0
    1    7    8    1   21
     21011
                 0.03.141592653589793238
    1    7    8    1   41
    121011
                 0.03.141592653589793238
    2    3    2    5    5
      1001
                 0.0                 0.0                 4.0                 4.
0
    1    1    2    1   21
         0
                 0.0                 2.0
    1    1    2    1   21
         1
                 0.0                 2.0
    1    2    3    1  101
         0
                 0.0                 2.0
    1    3    1    1  101
         2
    2    1
                 0.0                 2.0
    1    2    3    1  101
      1013
    1    0
                 0.0                 2.0
    0    0    0    0    0


