Subject: DIFCOR

From: [email-masked] (G. D. TAYLOR, MATH, COLORADO STATE UNIV.)
Date: Thu, 19 Feb 87 11:45:41 MST

Jack, attached is the code SDIFCOR complete with driver and input
at end of code.  I ran it on our CYBER 830 before sending it to
be safe and it ran.  (This the most recent version of the code and
I had just received it from Kaufman via a tape.)  Let me if I can
be of any help.  Jerry
 
      PROGRAM DRIVR(INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT)
      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(7)
C
C THE APPROXIMATION PROBLEMS RUN WITH THIS DEMONSTRATION DRIVER
C ARE AS FOLLOWS.
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 AND (P/Q)(Z) .GE. 0.0 FOR
C   1.0 .LT. Z .LE. 2.0.
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 .LE. PI.
C   (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.  WE REQUIRE
C   (P/Q)(Z) .GE. 0.0 FOR PI/2 .LT. Z .LE. PI.
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).
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 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     1010
C                0.0                 2.0
C   1    7    8    1   21
C    21010
C                0.03.141592653589793238
C   1    7    8    1   41
C   121010
C                0.03.141592653589793238
C   2    3    2    5    5
C     1000
C                0.0                 0.0                 4.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   0    0    0    0    0
C SET MACHINE DEPENDENT PARAMETERS FOR THE DEMONSTRATION DRIVER.
C SET INPUT AND OUTPUT UNIT NUMBERS.
      NREAD=I1MACH(1)
      NWRIT=I1MACH(2)
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=R1MACH(3)
C SET THE ERROR NORM TEST TOLERANCE ERTOL=MAX(10.0**(-12),
C 100.0*SPCMN).
      ERTOL=100.0*SPCMN
      IF(ERTOL-1.0E-12)10,20,20
   10 ERTOL=1.0E-12
   20 JOBNM=1
C SET THE ERROR NORMS COMPUTED WITH THE CDC CYBER 172 AT
C CENTRAL MICHIGAN UNIVERSITY FOR COMPARISON PURPOSES.
      ERTST(1)=0.36955949788958
      ERTST(2)=0.879206364E-5
      ERTST(3)=0.5440168621E-4
      ERTST(4)=0.44787565553231
      ERTST(5)=0.463E-11
      ERTST(6)=0.99999800E-6
      ERTST(7)=0.20000000000001
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)
      NMPP1=NUMP+1
      NMPPQ=NUMP+NUMQ
      NUMG1=NMPPQ
      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 READ THE OPTION SWITCH IOPT.  INSTRUCTIONS FOR SETTING IT
C ARE IN THE INITIAL COMMENTS OF SDFCOR.
      READ(NREAD, 80)IOPT
   80 FORMAT(I10)
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.
      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,
     1ANEP2,T)
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 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),JOBNM
C INPUT OF BASIS FUNCTIONS FOR JOBS 1, 5, 6, AND 7.
  160 PWR(I,1)=1.0
      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)=1.0
      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),JOBNM
C INPUT OF BASIS FUNCTIONS FOR JOBS 2 AND 3.
  230 PWR(I,1)=1.0
      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)=1.0
      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 INPUT OF BASIS FUNCTIONS FOR JOB 4.
  300 PWR(I,1)=1.0
      PWR(I,2)=T(IDBM)**2
      PWR(I,3)=T(IDB)**2
      PWR(I,4)=1.0
      PWR(I,5)=PWR(I,2)*PWR(I,3)
      GO TO 430
C INPUT OF F, INDF, AND LOWER CONSTRAINTS FOR JOB 1.
  310 NUM1=3*(NUMGR-1)/20+1
      NUM2=(NUMGR-1)/2+1
      NUM3=3*(NUMGR-1)/5+1
      IF(I-NUM1) 320, 320, 330
  320 INDF(I)=1
      FTBLE(I)=1.0
      INDLW(I)=1
      ALTBL(I)=0.9999
      GO TO 140
  330 IF(I-NUM2) 340, 340, 350
  340 INDF(I)=1
      FTBLE(I)=1.0
      INDLW(I)=0
      GO TO 140
  350 IF(I-NUM3) 360, 370, 370
  360 INDF(I)=0
      INDLW(I)=1
      ALTBL(I)=0.0
      GO TO 140
  370 INDF(I)=1
      FTBLE(I)=0.0
      INDLW(I)=1
      ALTBL(I)=0.0
      GO TO 140
C INPUT OF F, INDF, WEIGHT FUNCTION, AND LOWER CONSTRAINTS
C FOR JOBS 2 AND 3.
  380 NUM1=(NUMGR-1)/2+1
      NUM2=3*(NUMGR-1)/5+1
      IF(I-NUM1) 390, 390, 400
  390 INDF(I)=1
      FTBLE(I)=1.0
      WTBLE(I)=0.5
      INDLW(I)=0
      GO TO 140
  400 IF(I-NUM2) 410, 420, 420
  410 INDF(I)=0
      INDLW(I)=1
      ALTBL(I)=0.0
      GO TO 140
  420 INDF(I)=1
      FTBLE(I)=0.0
      WTBLE(I)=1.0
      INDLW(I)=1
      ALTBL(I)=0.0
      GO TO 140
C INPUT OF F AND INDF FOR JOB 4.
C WE USE 16.0001 INSTEAD OF 16.0 IN THE NEXT TEST TO AVOID TROUBLE
C DUE TO ROUND OFF ERROR.
  430 IF(T(IDBM)**2+T(IDB)**2-16.0001) 440, 440, 450
  440 INDF(I)=1
      FTBLE(I)=SQRT(16.0-T(IDBM)**2-T(IDB)**2)
      GO TO 140
  450 INDF(I)=0
      GO TO 140
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 EPSQ=1.0E-7
      IF(I-1)470,470,480
  470 FTBLE(I)=1.0
      GO TO 140
  480 FTBLE(I)=0.0
      GO TO 140
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)=3.0/(1.0+T(I))-2.0*T(I)+0.2
      GO TO 140
  510 IF(I-NUM2) 520, 520, 530
  520 FTBLE(I)=3.0/(1.0+T(I))+0.5*T(I)-0.3
      GO TO 140
  530 FTBLE(I)=3.0/(1.0+T(I))-0.4*T(I)+0.6
C END COMPUTING OF BASIS FUNCTIONS AND F.
      GO TO 140
  540 CALL SDFCOR(NUMP,NUMQ,NUMGR,FTBLE,PWR,IOPT,EPSQ,INDLW,ALTBL,
     *INDUP,UPTBL,INDF,WTBLE,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.)
  620 WRITE(NWRIT, 630)NUMDM
  630 FORMAT(/33H WE ARE DEALING WITH FUNCTIONS OF,I4,
     112H   VARIABLES)
      WRITE(NWRIT, 640)NUMP,NUMQ,NUMGR,NRWGR,NCLGR
  640 FORMAT(/6H P HAS,I4,24H   COEFFICIENTS    Q HAS,I4,
     131H   COEFFICIENTS    THE GRID HAS,I5,
     218H   POINTS ARRANGED,I5,5H   BY,I5)
      WRITE(NWRIT, 650)IOPT
  650 FORMAT(/8H IOPT IS,I10)
      IF(NUMDM-2) 660, 680, 660
  660 WRITE(NWRIT, 670)ASWP1,ASWP2
  670 FORMAT(/22H THE LEFT END POINT IS,E24.14,
     125H   THE RIGHT END POINT IS,E24.14)
      GO TO 700
  680 WRITE(NWRIT, 690)ASWP1,ASWP2,ANEP1,ANEP2
  690 FORMAT(/40H THE COORDINATES OF THE SW AND NE POINTS,
     14H ARE,2E15.5,6H   AND,2E15.5)
C WRITE IFLAG AND INDLP.
  700 WRITE(NWRIT, 710)IFLAG,INDLP
  710 FORMAT(/9H IFLAG IS,I5,13H     INDLP IS,I4)
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,25H,  ACHIEVED AT GRID POINT,I5)
C WRITE ITER AND DELTK.
      DELTK=ERROR(NUMGR+1)
      WRITE(NWRIT, 730)ITER,DELTK
  730 FORMAT(/21H THE ERROR NORM AFTER,I5,13H   ITERATIONS,
     13H IS,E25.15)
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/(/5E24.14))
      WRITE(NWRIT, 750)(X(I),I=NMPP1,NMPPQ)
  750 FORMAT(/26H THE COEFFICIENTS OF Q ARE/(/5E24.14))
      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(/6E20.10)
  780   CONTINUE
  790 JOBNM=JOBNM+1
      GO TO 30
      END
      FUNCTION I1MACH(I)
C
C THIS SUBPROGRAM IS NOT TO BE INCLUDED IN THE SLATEC LIBRARY.
C IT APPEARS HERE FOR TESTING PURPOSES TO MIMIC THE BEHAVIOR OF
C A FUNCTION ALREADY IN THE LIBRARY FOR ASSIGNING MACHINE
C DEPENDENT PARAMETERS, IN THE CASE OF THE CDC CYBER 172 AT
C CENTRAL MICHIGAN UNIVERSITY.
C
C I1MACH(1) IS THE INPUT UNIT NUMBER, AND I1MACH(2) IS THE OUTPUT
C UNIT NUMBER.
      IF(I-1)1,1,2
    1 I1MACH=5
      RETURN
    2 I1MACH=6
      RETURN
      END
      FUNCTION R1MACH(I)
C
C THIS SUBPROGRAM IS NOT TO BE INCLUDED IN THE SLATEC LIBRARY.
C IT APPEARS HERE FOR TESTING PURPOSES TO MIMIC THE BEHAVIOR
C OF A FUNCTION ALREADY IN THE LIBRARY FOR ASSIGNING MACHINE
C DEPENDENT PARAMETERS, IN THE CASE OF THE CDC CYBER 172 AT
C CENTRAL MICHIGAN UNIVERSITY.
C
C R1MACH(3) IS B**(-T), WHERE B IS THE BASE FOR FLOATING POINT NUMBERS
C AND T IS THE NUMBER OF BASE B DIGITS IN THE MANTISSA.
      R1MACH=2.0**(-48)
      RETURN
      END
      SUBROUTINE SDFCOR(NUMP,NUMQ,NUMGR,FTBLE,PWR,IOPT,EPSQ,INDLW,
     *ALTBL,INDUP,UPTBL,INDF,WTBLE,X,ERROR,QMIN,IQMIN,ITER,INDLP,IFLAG)
C***BEGIN PROLOGUE  SDFCOR
C***DATE WRITTEN   720120   (YYMMDD)
C***REVISION DATE  860328   (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.
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  (OPTIONAL INPUT)
C         THESE SEVEN VARIABLES ARE USED FOR INPUT FOR THE EXTRA
C         OPTIONS DESCRIBED IN THE LONG DESCRIPTION BELOW.  NONE OF
C         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.
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 SDFCOR 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).
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 AND UP TO 101 GRID POINTS (SO NUMGR .LE. 101), EVEN IF ALL THE
C EXTRA OPTIONS DESCRIBED IN THE LONG DESCRIPTION BELOW ARE
C ACTIVATED.  SEE THE LONG DESCRIPTION FOR THE DIMENSION
C PATTERN REQUIRED IN THE DRIVER PROGRAM AND SUBROUTINES IF
C MORE BASIS FUNCTIONS OR GRID POINTS 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 SIX 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 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  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  OPTION 5  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 6  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 TWO
C     POSSIBILITIES.
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.  NO APPROXIMATION WAS RETURNED, ALL THE COEFFICIENTS
C         OF P AND Q WERE SET 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 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) AND MGR GRID POINTS (THUS
C NUMGR .LE. MGR) 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), XTEMP(MFN+1), FTBWT(MGR), IYRCT(MAXCN), IYSAV(MAXCN),
C V(MAXCN+1,MFN+2), IYCCT(MFN+1), IXRCT(MAXCN), Y(MAXCN),
C
C WHERE MAXCN=5*MGR+2*MFN-2.  NOTE THAT THE LAST 8 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.
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 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 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.
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 HAD A DENOMINATOR WHICH WAS VERY
C SMALL IN ABSOLUTE VALUE OR NEGATIVE AT SOME POINT (AND THE MAIN
C ALGORITHM WAS UNABLE TO CORRECT THIS PROBLEM).  THE THREE
C METHODS ARE
C
C  METHOD 1  THE USER SUPPLIES AN INITIAL APPROXIMATION (SEE
C            OPTION 6 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, AND Q(1)=1.0
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, AND
C            THE ABSOLUTE VALUES OF THE FREE DENOMINATOR
C            COEFFICIENTS ARE BOUNDED BY 1.0
C            IF SLNPRO FAILS TO PRODUCE AN APPROXIMATION P/Q WE
C            TRY WITH BOTH LOOSENED AND TIGHTENED TOLERANCES,
C            AND ALSO WITH Q(1)=-1.0
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***ROUTINES CALLED  SINTPQ,SLNPRO,SPQERD,SSETV1,SSETV2
C***END PROLOGUE  SDFCOR
      DIMENSION FTBLE(101),PWR(101,15),V(534,17),WTBLE(101),
     *   UPTBL(101),ALTBL(101),INDUP(101),INDLW(101),
     *   ERROR(102),IYRCT(533),IYSAV(533),XTEMP(16),X(16),
     *   INDF(101),FTBWT(101)
C
C THE BASIC PATTERN OF STATEMENT NUMBERS IN THIS PROGRAM IS THAT
C STATEMENTS NUMBERED 1, 2,... ARE FROM A VERSION OF THE PROGRAM
C PUBLISHED PRIOR TO THE THIRD REFERENCE CITED ABOVE, STATEMENT
C NUMBERS 1001, 1002,... ARE FROM THE VERSION OF THAT PAPER, AND
C STATEMENT NUMBERS 10000 OR HIGHER REPRESENT LATER MODIFICATIONS.
C
C SET MACHINE DEPENDENT PARAMETERS FOR SUBROUTINE SDFCOR.
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  SDFCOR
      SPCMN=R1MACH(3)
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=10000.0*SPCMN
C END OF SETTING MACHINE DEPENDENT PARAMETERS FOR SDFCOR.
C
      ITLIM=100
      ITER=0
      NMPP1=NUMP+1
      NMPPQ=NUMP+NUMQ
      N=NMPPQ+1
      NP1=N+1
      IADJS=0
      INDIC=0
      IRET1=0
      ILBFL=0
      ILST1=0
      QMIN=0.0
      IQMIN=0
      INITL=0
      ICHK=0
      IFLAG=0
      IPBIT=1
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)1001,1003,1003
 1001 DO 1002 I=1,NUMGR
        INDF(I)=1
 1002   CONTINUE
C
C IF WE ARE DEALING WITH UNWEIGHTED APPROXIMATION SET THE
C WEIGHT FUNCTION IDENTICALLY EQUAL TO 1.
 1003 IF(IOPT-(IOPT/100000)*100000-10000)1004,1006,1006
 1004 DO 1005 I=1,NUMGR
        WTBLE(I)=1.0
 1005   CONTINUE
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.
 1006 IF(IOPT-(IOPT/100000)*100000-20000)10000,1007,1007
10000 DO 10020 I=1,NUMGR
        IF(INDF(I))10020,10020,10010
10010   FTBWT(I)=FTBLE(I)
10020   CONTINUE
      GO TO 1010
 1007 DO 1009 I=1,NUMGR
        IF(INDF(I))1009,1009,1008
 1008   FTBWT(I)=WTBLE(I)*FTBLE(I)
 1009   CONTINUE
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.
 1010 IF(IOPT/100000)1012,1012,10030
C SET INITL=-1 TO INDICATE WE ARE USING A USER SUPPLIED INITIAL
C APPROXIMATION.
10030 INITL=-1
      GO TO 1024
C
C DURING THE LOEB INITIALIZATION X(N) WILL BE USED TO STORE
C THE VALUE ASSIGNED TO Q(1) (EITHER 1.0 OR -1.0).
 1012 X(N)=1.0
C WE TRY TO FIND P,Q SUCH THAT Q(1)=1.0, Q IS POSITIVE (AND
C NOT TOO SMALL) AT EVERY GRID POINT, P/Q SATISFIES THE
C RESTRICTED RANGE CONSTRAINTS, AND MAX ABS(F*Q-WT*P) IS
C MINIMIZED AT THOSE GRID POINTS WHERE F IS TO BE
C APPROXIMATED.
 1013 CALL SSETV1(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,INDLW,ALTBL,
     *INDUP,UPTBL,INDF,WTBLE,X,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 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.
      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 NMPPQ=-NMPPQ AS A SIGNAL TO TIGHTEN THE TOLERANCES IN SLNPRO
C (I. E. DECREASE REA AND REA1) AND TRY ONCE AGAIN.  SLNPRO WILL
C RESET NMPPQ TO ITS ORIGINAL VALUE.
      NMPPQ=NMPPQ*(1+IADJS-IADJS**2)
      CALL SLNPRO(V,M,NMPPQ,IYRCT,X,INDIC)
      IRET1=INDIC
C IF THIS WAS THE FIRST SLNPRO CALL IN THE LOEB INITIALIZATION,
C SO WE HAD X(N)=1.0 AND THE TOLERANCES IN SLNPRO WERE NOT
C ADJUSTED BY SDFCOR, FOR POSSIBLE OUTPUT PURPOSES WE SAVE
C INDIC IN ILBFL.
      IF(X(N))10060,10060,10040
10040 IF(IADJS)10050,10050,10060
10050 ILBFL=INDIC
C IF INDIC IS POSITIVE SOMETHING WENT WRONG IN SLNPRO.
C OTHERWISE WE PUT THE COEFFICIENTS OF Q INTO THEIR PROPER
C PLACES IN X AND NORMALIZE P AND Q.
10060 IF(INDIC)1023,1023,1016
C IF WE CANNOT FIND P/Q (WITH Q POSITIVE) SATISFYING THE
C CONSTRAINTS WITH Q(1)=1.0, WE TRY WITH Q(1)=-1.0.
 1016 IF(X(N))1018,1018,1017
 1017 X(N)=-1.0
      GO TO 1013
 1018 IF(IADJS-1)1019,1019,1020
C BEFORE CONCEDING DEFEAT WE WILL TRY AGAIN WITH CHANGED
C TOLERANCES.
 1019 IADJS=IADJS+1
      GO TO 1012
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.
 1020 INITL=1
      DO 10070 J=1,NMPPQ
        X(J)=0.0
10070   CONTINUE
      X(NMPP1)=1.0
C RESET PARAMETERS.
      IADJS=0
      INDIC=0
      IRET1=0
      QMIN=0.0
      IQMIN=0
      GO TO 1024
C
C HERE INDIC WAS NONPOSITIVE, SO SLNPRO PRODUCED AN APPROXIMATION.
C PUT THE COEFFICIENTS OF P AND Q IN THE PROPER PLACES IN X AND
C NORMALIZE Q.
 1023 CALL SINTPQ(NUMP,NUMQ,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.
 1024 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X,
     *V,QMIN,IQMIN)
      DELTK=V(2*NUMGR+1,NP1)
      DLOLD=DELTK
      IADJS=0
C PUT THE ERRORS AND ERROR NORM INTO ERROR.
      DO 6 I=1,NUMGR
        ERROR(I)=V(I,N)
    6   CONTINUE
      ERROR(NUMGR+1)=DELTK
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 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)1025,1026,1026
 1025 V(2*NUMGR+1,NP1)=-V(2*NUMGR+1,NP1)-2.0
C END OF INITIALIZATION PHASE.
C
C
C SET UP AND SOLVE THE LINEAR PROGRAMMING PROBLEM TO GET THE
C NEXT APPROXIMATION.
 1026 CALL SSETV2(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,INDLW,
     *ALTBL,INDUP,UPTBL,INDF,WTBLE,V,M)
C IF WE ARE IN THE FIRST STEP OF THE MAIN ALGORITHM, SET
C IYRCT(1)=-1 AS A SIGNAL TO DO THE INITIAL EXCHANGES IN
C SLNPRO STRICTLY ACCORDING TO A PIVOTING STRATEGY.
C HENCEFORTH WHEN SLNPRO IS CALLED THE OLD VERTEX WILL BE
C TAKEN AS THE STARTING POINT.
      IF(ITER)1027,1027,1028
 1027 IYRCT(1)=-1
C SAVE IYRCT IN IYSAV.
 1028 DO 1029 I=1,M
        IYSAV(I)=IYRCT(I)
 1029   CONTINUE
C SET M=-M AS A SIGNAL TO LOOSEN THE TOLERANCES IN SLNPRO (I.E.
C INCREASE REA AND REA1) IF IADJS=1, AND LEAVE M UNCHANGED IF
C IADJS=0 OR IADJS=2.  SLNPRO WILL AUTOMATICALLY RESET M IF IT
C IS CHANGED HERE.
      M=M*(1-4*IADJS+2*IADJS**2)
C SET N=-N AS A SIGNAL TO TIGHTEN THE TOLERANCES IN SLNPRO (I.E.
C DECREASE REA AND REA1) IF IADJS=2, AND LEAVE N UNCHANGED IF
C IADJS=0 OR IADJS=1.  SLNPRO WILL AUTOMATICALLY RESET N IF IT
C IS CHANGED HERE.
      N=N*(1+IADJS-IADJS**2)
      CALL SLNPRO(V,M,N,IYRCT,XTEMP,INDIC)
C IF THIS WAS THE FIRST SLNPRO CALL IN THIS ITERATION, SO THE
C TOLERANCES IN SLNPRO WERE NOT ADJUSTED BY SDFCOR, FOR POSSIBLE
C OUTPUT PURPOSES WE SAVE INDIC IN ILST1.
      IF(IADJS)10080,10080,10090
10080 ILST1=INDIC
C IF INDIC IS POSITIVE SOMETHING WENT WRONG IN SLNPRO, SO
C WE RETURN WITH THE X AND ERROR COMPUTED AT THE PREVIOUS
C ITERATION AFTER TRYING AGAIN WITH BOTH TIGHTER AND LOOSER
C TOLERANCES, IF THIS HAS NOT ALREADY BEEN TRIED.
10090 IF(INDIC)1030,1030,1045
C
C HERE INDIC WAS NONPOSITIVE, SO SLNPRO PRODUCED AN APPROXIMATION.
C COMPUTE THE VALUES OF P AND Q AT THE GRID POINTS, THE
C ERRORS AT THE GRID POINTS, AND THE UNIFORM ERROR NORM.
 1030 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,XTEMP,
     *V,QMIN,IQMIN)
      DELTK=V(2*NUMGR+1,NP1)
C REJECT THIS APPROXIMATION IF THE DENOMINATOR WAS VERY
C SMALL IN ABSOLUTE VALUE OR NEGATIVE AT SOME POINT.
C THIS IS INDICATED BY DELTK BEING NEGATIVE.
      IF(DELTK)1045,1031,1031
C IF INDIC IS NEGATIVE HERE IT WAS NECESSARY TO ADJUST THE
C TOLERANCES IN SLNPRO TO GET AN ACCEPTABLE APPROXIMATION, SO THE
C APPROXIMATION COMPUTED IS MORE LIKELY TO HAVE BEEN
C AFFECTED BY ROUND OFF ERROR.  THIS PRESENTS NO DANGER IF NEITHER
C RESTRICTED RANGE NOR DENOMINATOR CONSTRAINTS WERE USED, SINCE AN
C APPROXIMATION WITH ERROR NORM WORSE THAN THE PREVIOUS ONE
C WILL BE REJECTED BELOW ANYWAY, BUT IF THERE ARE RESTRICTED
C RANGE OR DENOMINATOR CONSTRAINTS THE ERROR NORM COULD HAVE
C BEEN REDUCED AT THE EXPENSE OF VIOLATING THE CONSTRAINTS.
C THUS WE REJECT THE APPROXIMATION IF ANY CONSTRAINT IS VIOLATED
C BY MORE THAN EPRES.
 1031 IF(INDIC)1032,1040,1040
C
C CHECK UPPER CONSTRAINTS, IF THERE ARE ANY.
 1032 IF(IOPT-(IOPT/1000)*1000-100)1036,1033,1033
 1033 DO 1035 I=1,NUMGR
        IF(INDUP(I))1035,1035,1034
 1034   IF(V(2*I-1,NP1)/V(2*I,NP1)-UPTBL(I)-EPRES)1035,1035,10140
 1035   CONTINUE
C CHECK LOWER CONSTRAINTS, IF THERE ARE ANY.
 1036 IF(IOPT-(IOPT/100)*100-10)10100,1037,1037
 1037 DO 1039 I=1,NUMGR
        IF(INDLW(I))1039,1039,1038
 1038   IF(V(2*I-1,NP1)/V(2*I,NP1)-ALTBL(I)+EPRES)10140,1039,1039
 1039   CONTINUE
C CHECK DENOMINATOR CONSTRAINTS, IF THERE ARE ANY.
10100 IF(IOPT-(IOPT/10)*10)10130,10130,10110
10110 DO 10120 I=1,NUMGR
        IF(V(2*I,NP1)-EPSQ+EPRES)10140,10120,10120
10120   CONTINUE
C HERE NO USER SUPPLIED CONSTRAINT WAS VIOLATED BY MORE THEN EPRES.
C IF WE ARE IN THE FINAL CHECK WE RETURN.
10130 IF(ICHK)1040,1040,15
C HERE AT LEAST ONE USER SUPPLIED CONSTRAINT WAS VIOLATED BY MORE
C THAN EPRES.  IF WE ARE IN THE FINAL CHECK WE ADJUST IFLAG AND
C RETURN.
10140 IF(ICHK)1045,1045,10150
10150 IFLAG=IFLAG-100
      RETURN
C
C IF DLOLD IS NEGATIVE, THEN WE WILL HAVE ITER=0, AND THE
C INITIAL APPROXIMATION WILL HAVE A DENOMINATOR WHICH IS
C VERY SMALL IN ABSOLUTE VALUE OR NEGATIVE AT SOME POINT.
C THUS THE NEW APPROXIMATION, WHICH HAS A SIGNIFICANTLY
C POSITIVE DENOMINATOR EVERYWHERE ON THE GRID, WILL BE
C PREFERRED.
 1040 IF(DLOLD)10,1041,1041
C IF THE ERROR NORM HAS NOT DECREASED, RETURN WITH THE X
C AND ERROR COMPUTED AT THE PREVIOUS ITERATION AFTER TRYING
C AGAIN WITH CHANGED TOLERANCES, IF BOTH TIGHTER AND LOOSER
C TOLERANCES HAVE NOT ALREADY BEEN TRIED.
 1041 IF(DELTK-DLOLD)1042,1045,1045
C IN RARE INSTANCES AN APPROXIMATION WITH ALL COEFFICIENTS
C SMALL BUT WITH A GOOD COMPUTED ERROR NORM WILL BE
C PRODUCED.  TO AVOID THIS WE REJECT AN APPROXIMATION IF
C MAX(ABS(Q(J))).LE.0.1
 1042 DO 1043 J=NMPP1,NMPPQ
        IF(ABS(XTEMP(J))-0.1)1043,1043,10
 1043   CONTINUE
      GO TO 1045
C
C HERE THE PRESENT APPROXIMATION WILL BE KEPT, AND WE PREPARE
C TO COMPUTE ANOTHER ONE.
C COPY XTEMP INTO X.
   10 DO 11 I=1,NMPPQ
        X(I)=XTEMP(I)
   11   CONTINUE
C PUT THE ERRORS AND ERROR NORM INTO ERROR.
      DO 12 I=1,NUMGR
        ERROR(I)=V(I,N)
   12   CONTINUE
      ERROR(NUMGR+1)=DELTK
      ITER=ITER+1
      IF(ITER-ITLIM)14,10160,10160
   14 DLOLD=DELTK
      IRET1=INDIC
      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)1026,1026,10155
10155 V(2*NUMGR+1,NP1)=0.5*V(2*NUMGR+1,NP1)
      GO TO 1026
C HERE ITLIM ITERATIONS HAVE BEEN PERFORMED AND WE WILL RETURN
C AFTER SETTING IFLAG AND INDLP (SEE BELOW).
10160 INDLP=5-INDIC
      IF(INDIC)10170,10180,10180
10170 IFLAG=-12
      GO TO 10300
10180 IFLAG=-2
      GO TO 10300
C
C HERE SLNPRO FAILED TO PRODUCED AN APPROXIMATION WHICH IS BETTER
C THAN THE PREVIOUS APPROXIMATION.  IF WE HAVE NOT ALREADY
C DONE SO, WE TRY AGAIN WITH TIGHTER TOLERANCES.  IF THIS
C HAS BEEN TRIED, WE TRY ONCE MORE WITH LOOSER TOLERANCES.
C THE ONLY EXCEPTION TO THIS IS THAT IF WE ARE STILL IN THE
C PARAMETRIC BISECTION PHASE AND ITER .GT. 0, TO SAVE TIME WE
C DO NOT TRY WITH CHANGED TOLERANCES, SINCE THE FAILURE WAS
C PROBABLY DUE TO THE PARAMETRIC BISECTION.  INSTEAD WE TRY
C WITH ORDINARY DIFFERENTIAL CORRECTION, AND USE ORDINARY
C DIFFERENTIAL CORRECTION FROM NOW ON.
 1045 IF(ITER*IPBIT)10185,10185,10183
10183 IPBIT=0
      GO TO 10187
10185 IF(IADJS-1)1046,10190,1046
 1046 IADJS=2-IADJS/2
C RESTORE THE PREVIOUS IYRCT FOR THE NEXT ATTEMPT.
10187 DO 1047 I=1,M
        IYRCT(I)=IYSAV(I)
 1047   CONTINUE
C CALL SPQERD WITH THE PREVIOUS APPROXIMATION (WHICH IS IN X)
C TO RECOMUTE THE PREVIOUS VALUES OF P AND Q AND PUT THEM
C IN V FOR USE BY SSETV2.
      CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X,
     *V,QMIN,IQMIN)
      GO TO 1026
C
C HERE THE APPROXIMATION HAS NOT BEEN IMPROVED EVEN BY ADJUSTING
C TOLERANCES TWICE IN SLNPRO.  IF IT IS AN INITIAL APPROXIMATION
C WHICH WAS FOUND BY SPQERD TO HAVE A DENOMINATOR WHICH IS VERY
C SMALL IN ABSOLUTE VALUE OR NEGATIVE AT SOME POINT (INDCIATED BY
C DLOLD NEGATIVE), AND IF IT IS A USER SUPPLIED APPROXIMATION
C (INDICATED BY INITL=-1), THEN WE ATTEMPT A LOEB REINITIALIZATION.
C IF IT WAS A LOEB APPROXIMATION (INDICATED BY INITL=0) WITH
C UNACCEPTABLE DENOMINATOR, THEN WE TRY P/Q=0.0/PWR(.,NUMP+1).
10190 IF(DLOLD)10200,1048,1048
10200 IF(INITL)10210,1020,1048
C SET UP TO TRY THE LOEB INITIALIZATION.
10210 INITL=0
      IADJS=0
      INDIC=0
      IRET1=0
      QMIN=0.0
      IQMIN=0
      GO TO 1012
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 SDFCOR 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 SDFCOR
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.
 1048 INDLP=10*(5-ILST1)+5-IRET1
C
C WE NOW SET UP IFLAG FOR THE INFORMATION OF THE USER.
      IF(ITER)10220,10220,10280
C HERE ITER=0, SO THE INITIAL APPROXIMATION WAS NOT IMPROVED.
10220 IF(INITL)10270,10270,10230
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, THEN SET IFLAG=6666, ITER=-1, AND
C ZERO OUT THE ERRORS (WHICH WERE INVALID ANYWAY) AND
C COEFFICIENTS AND RETURN.  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.
10230 INDLP=10*(5-ILST1)+5-ILBFL
      IF(DLOLD)10240,10260,10260
10240 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)=0.0
      DO 10250 I=1,NUMGR
        ERROR(I)=0.0
10250   CONTINUE
      RETURN
10260 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.
10270 IFLAG=-1
C IF IRET1.LT.0 SET THE TENS DIGIT OF IFLAG TO INDICATE THAT
C THE TOLERANCES IN SLNPRO HAD TO BE ADJUSTED DURING THE
C COMPUTATION OF THE APPROXIMATION RETURNED.
10280 IF(IRET1)10290,10300,10300
10290 IFLAG=IFLAG-10
C NOW WE GO BACK AND CHECK THE USER SUPPLIED CONSTRAINTS (IF ANY).
C SET ICHK=1 AS A SIGNAL AND CALL SPQERD TO RECOMPUTE THE VALUES
C OF P AND Q.  QMIN AND IQMIN ARE ALSO RECOMPUTED BY SPQERD.
10300 ICHK=1
      CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X,
     *V,QMIN,IQMIN)
      GO TO 1032
   15 RETURN
      END
      SUBROUTINE STCOMP(NUMDM,NRWGR,NCLGR,ASWP1,ASWP2,
     *ANEP1,ANEP2,T)
C***BEGIN PROLOGUE  STCOMP
C***DATE WRITTEN   720120   (YYMMDD)
C***REVISION DATE  840301   (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
      DIMENSION T(501)
C
C***FIRST EXECUTABLE STATEMENT  STCOMP
      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=0.0
    6 IF(NRWGR-1)7,8,7
    7 VSPAC=(ANEP2-ASWP2)/FLOAT(NRWGR-1)
      GO TO 9
    8 VSPAC=0.0
    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 SSETV1(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,INDLW,
     *ALTBL,INDUP,UPTBL,INDF,WTBLE,X,V,M)
C***BEGIN PROLOGUE  SSETV1
C***REFER TO  SDFCOR
C***ROUTINES CALLED  (NONE)
C***PURPOSE  THIS SUBROUTINE SETS UP THE LINEAR PROGRAMMING PROBLEM
C            FOR THE LOEB INITIALIZATION OF THE DIFFERENTIAL
C            CORRECTION ALGORITHM.
C***END PROLOGUE  SSETV1
      DIMENSION FTBWT(101),PWR(101,15),V(534,17),WTBLE(101),
     *   UPTBL(101),ALTBL(101),INDUP(101),INDLW(101),X(16),
     *   INDF(101)
C
C THIS SUBROUTINE SETS UP V FOR THE LINEAR PROGRAMMING
C PROBLEM OF MINIMIZING W WITH THE RESTRICTIONS
C ABS(F*Q-WT*P).LE.W, WITH Q(1) FIXED AS X(NMPPQ+1).
C
C SET MACHINE DEPENDENT PARAMETERS FOR SUBROUTINE SSETV1.
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  SSETV1
      SPCMN=R1MACH(3)
C SET EPS=MAX(0.0001,10000.0*SPCMN).  THUS EPS=10.0**(-EXEPS),
C WHERE EXEPS=MIN(4,TNMAN-4).
C WE WILL REQUIRE Q.GE.EPS AT EVERY GRID POINT.
      EPS=10000.0*SPCMN
      IF(EPS-0.0001)10000,10010,10010
10000 EPS=0.0001
C IF WE ARE PLACING A LOWER BOUND ON THE DENOMINATOR IN THE MAIN
C ITERATIONS (INDICATED BY THE ONES DIGIT OF IOPT BEING 1) WE
C INCREASE EPS, IF NECESSARY, TO MATCH THIS LOWER BOUND.
10010 IF(IOPT-(IOPT/10)*10)10040,10040,10020
10020 IF(EPSQ-EPS)10040,10040,10030
10030 EPS=EPSQ
C END OF SETTING MACHINE DEPENDENT PARAMETERS FOR SSETV1.
C
10040 NMPP1=NUMP+1
      N1=NUMP+NUMQ
      NMPQ1=N1-1
      N1P1=N1+1
C M WILL KEEP TRACK OF WHICH ROW OF V IS BEING SET.  AT THE
C CONCLUSION OF SSETV1 M WILL BE THE TOTAL NUMBER OF
C CONSTRAINTS.
      M=0
C
C SET THE UPPER CONSTRAINTS, IF THERE ARE ANY.  THEY WILL BE
C OF THE FORM P-Q*UP.LE.0.0, WITH THE CONSTANT Q(1)*UP TERM
C TRANSPOSED.
      IF(IOPT-(IOPT/1000)*1000-100)1008,1001,1001
 1001 DO 1007 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))1007,1007,1002
 1002   M=M+1
        DO 1003 J=1,NUMP
          V(M,J)=PWR(I,J)
 1003     CONTINUE
        IF(NUMQ-1)1006,1006,1004
 1004   DO 1005 J=NMPP1,NMPQ1
          JP1=J+1
          V(M,J)=-UPTBL(I)*PWR(I,JP1)
 1005     CONTINUE
 1006   V(M,N1)=0.0
        V(M,N1P1)=X(N1P1)*UPTBL(I)*PWR(I,NMPP1)
 1007   CONTINUE
C
C SET THE LOWER CONSTRAINTS, IF THERE ARE ANY.  THEY WILL BE
C OF THE FORM -P+Q*LW.LE.0.0, WITH THE CONSTANT Q(1)*LW TERM
C TRANSPOSED.
 1008 IF(IOPT-(IOPT/100)*100-10)1016,1009,1009
 1009 DO 1015 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))1015,1015,1010
 1010   M=M+1
        DO 1011 J=1,NUMP
          V(M,J)=-PWR(I,J)
 1011     CONTINUE
        IF(NUMQ-1)1014,1014,1012
 1012   DO 1013 J=NMPP1,NMPQ1
          JP1=J+1
          V(M,J)=ALTBL(I)*PWR(I,JP1)
 1013     CONTINUE
 1014   V(M,N1)=0.0
        V(M,N1P1)=-X(N1P1)*ALTBL(I)*PWR(I,NMPP1)
 1015   CONTINUE
C
C FOR EACH GRID POINT AT WHICH F IS TO BE APPROXIMATED SET
C THE CONSTRAINTS -WT*P+F*Q-W.LE.0.0 AND WT*P-F*Q-W.LE.0.0,
C WITH THE CONSTANT Q(1)*F TERM TRANSPOSED.
 1016 DO 1022 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))1022,1022,1017
 1017   M=M+2
        MM1=M-1
        DO 1018 J=1,NUMP
          V(MM1,J)=-WTBLE(I)*PWR(I,J)
          V(M,J)=-V(MM1,J)
 1018     CONTINUE
        IF(NUMQ-1)1021,1021,1019
 1019   DO 1020 J=NMPP1,NMPQ1
          JP1=J+1
          V(MM1,J)=FTBWT(I)*PWR(I,JP1)
          V(M,J)=-V(MM1,J)
 1020     CONTINUE
 1021   V(MM1,N1)=-1.0
        V(M,N1)=-1.0
        V(MM1,N1P1)=-X(N1P1)*FTBWT(I)*PWR(I,NMPP1)
        V(M,N1P1)=-V(MM1,N1P1)
 1022   CONTINUE
C
C NOW SET THE CONSTRAINTS OF THE FORM -Q.LE.-EPS, WITH THE
C CONSTANT Q(1) TERM TRANSPOSED.
      DO 1027 I=1,NUMGR
        M=M+1
        DO 1023 J=1,NUMP
          V(M,J)=0.0
 1023     CONTINUE
        IF(NUMQ-1)1026,1026,1024
 1024   DO 1025 J=NMPP1,NMPQ1
          JP1=J+1
          V(M,J)=-PWR(I,JP1)
 1025     CONTINUE
 1026   V(M,N1)=0.0
        V(M,N1P1)=X(N1P1)*PWR(I,NMPP1)-EPS
 1027   CONTINUE
C
C IF WE ARE PLACING A LOWER BOUND ON THE DENOMINATOR (INDICATED
C BY THE ONES COEFFICIENT OF IOPT BEING 1), THEN WE BOUND THE
C FREE DENOMINATOR COEFFICIENTS (IF THERE ARE ANY) SO THAT WHEN
C THE DENOMINATOR IS NORMALIZED IN SINTPQ, THE DENOMINATOR BOUND
C WILL NOT BE VIOLATED.
      IF(IOPT-(IOPT/10)*10)10090,10090,10050
10050 IF(NUMQ-1)10090,10090,10060
C SET THE CONSTRAINTS OF THE FORM -Q(J) .LE. 1.0 AND Q(J)
C .LE. 1.0 FOR J=2,...,NUMQ.
10060 DO 10080 I=2,NUMQ
        M=M+2
        MM1=M-1
        DO 10070 J=1,N1
          V(MM1,J)=0.0
          V(M,J)=0.0
10070     CONTINUE
        NNN1=NUMP+I-1
        V(MM1,NNN1)=-1.0
        V(M,NNN1)=1.0
        V(MM1,N1P1)=1.0
        V(M,N1P1)=1.0
10080   CONTINUE
C
C NOW SET THE BOTTOM ROW.  TO MINIMIZE W WE MAXIMIZE -W.
10090 MP1=M+1
      DO 6 J=1,N1P1
        V(MP1,J)=0.0
    6   CONTINUE
      V(MP1,N1)=1.0
      RETURN
      END
      SUBROUTINE SINTPQ(NUMP,NUMQ,X)
C***BEGIN PROLOGUE  SINTPQ
C***REFER TO  SDFCOR
C***ROUTINES CALLED  (NONE)
C***PURPOSE  THIS SUBROUTINE PUTS THE DENOMINATOR COEFFICIENTS
C            INTO THE PROPER LOCATIONS IN X AND NORMALIZES P/Q
C            DURING THE LOEB INITIALIZATION OF THE DIFFERENTIAL
C            CORRECTION ALGORITHM.
C***END PROLOGUE  SINTPQ
      DIMENSION X(16)
C
C***FIRST EXECUTABLE STATEMENT  SINTPQ
      NMPP1=NUMP+1
      NMPPQ=NUMP+NUMQ
      NMQ1=NUMQ-1
C PUT THE COEFFICIENTS OF Q IN THEIR PROPER PLACES IN X.
      IF(NUMQ-1)1002,1002,1001
 1001 DO 1 J=1,NMQ1
        K=NMPPQ+1-J
        X(K)=X(K-1)
    1   CONTINUE
 1002 X(NMPP1)=X(NMPPQ+1)
C
C WE DIVIDE BY THE LARGEST ABSOLUTE VALUE OF THE
C COEFFICIENTS OF Q TO NORMALIZE P AND Q.
      AMAX=0.0
      DO 12 J=NMPP1,NMPPQ
        IF(ABS(X(J))-AMAX)12,12,11
   11   AMAX=ABS(X(J))
   12   CONTINUE
      DO 13 J=1,NMPPQ
        X(J)=X(J)/AMAX
   13   CONTINUE
      RETURN
      END
      SUBROUTINE SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X,
     *V,QMIN,IQMIN)
C***BEGIN PROLOGUE  SPQERD
C***REFER TO SDFCOR
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
      DIMENSION FTBWT(101),PWR(101,15),V(534,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=R1MACH(3)
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=100.0*SPCMN
C END OF SETTING MACHINE 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)=0.0
      ABAD=1.0
    1 I=I+1
      IF(I-NUMGR)3,3,1001
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
    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)=0.0
      DO 4 J=1,NUMP
        V(IDBM,NP1)=V(IDBM,NP1)+X(J)*PWR(I,J)
    4   CONTINUE
C COMPUTE THE VALUE OF Q AT GRID POINT I AND PUT IT IN
C V(2*I,N+1).
      V(IDB,NP1)=0.0
      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 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)=0.0
      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=-1.0
      V(I,N)=0.0
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)=10.0*EPS2
      GO TO 1
      END
      SUBROUTINE SLNPRO(V,M,N,IYRCT,X,INDIC)
C***BEGIN PROLOGUE  SLNPRO
C***REFER TO  SDFCOR
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
      DIMENSION V(534,17),IYRCT(533),X(16),Y(533),
     *   IXRCT(533),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 DE(ENDENT 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=R1MACH(3)
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=100.0*SPCMN
      IF(REA-1.0E-8)10000,10010,10010
10000 REA=1.0E-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=10.0*SPCMN
C END OF INITIAL SETTING OF MACHINE DEPENDENT PARAMETERS FOR
C SLNPRO.  THESE PARAMETERS MAY BE ADJUSTED BY A COMMAND FROM
C SDFCOR.
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/100.0
      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)=0.0
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=1.0
      AMPRV=1.0
      AMPR2=1.0
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=0.0
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=1.0
      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=1.0
      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=1.0
        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-0.5)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=1.0
      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=-1.0
        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+0.5)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)=0.0
   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
      DIMENSION V(534,17)
C
C DIVIDE THE ENTRIES IN THE RESOLVENT ROW (EXCEPT FOR THE
C RESOLVENT) BY THE RESOLVENT.
C
C***FIRST EXECUTABLE STATEMENT  SJELIM
      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)=1.0/RESOL
      RETURN
      END
      SUBROUTINE SSETV2(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,
     *INDLW,ALTBL,INDUP,UPTBL,INDF,WTBLE,V,M)
C***BEGIN PROLOGUE  SSETV2
C***REFER TO  SDFCOR
C***ROUTINES CALLED  (NONE)
C***PURPOSE  THIS SUBROUTINE SETS UP THE LINEAR PROGRAMMING PROBLEM
C            FOR THE MAIN ITERATIONS OF THE DIFFERENTIAL CORRECTION
C            ALGORITHM.
C***END PROLOGUE  SSETV2
      DIMENSION FTBWT(101),PWR(101,15),V(534,17),WTBLE(101),
     *   UPTBL(101),ALTBL(101),INDUP(101),INDLW(101),
     *   INDF(101)
C
C THIS SUBROUTINE SETS UP V FOR THE LINEAR PROGRAMMING
C PROBLEM OF MINIMIZING W WITH THE RESTRICTIONS
C   (ABS(F*Q-WT*P)-DELTK*Q)/QK.LE.W AND
C   ABS(Q(J)).LE.1.0 FOR ALL J.
C EACH RESTRICTION IS BROKEN INTO TWO CONSTRAINTS TO ACCOUNT
C FOR THE ABSOLUTE VALUE.  THE SUBROUTINE USES THE FACT THAT
C THE VALUE OF QK AT GRID POINT I IS STORED IN V(2*I,N+1)
C AND DELTK IS IN V(2*NUMGR+1,N+1) FROM A PREVIOUS CALL TO
C SPQERD.
C
C***FIRST EXECUTABLE STATEMENT  SSETV2
      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 SSETV2 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)1006,1001,1001
 1001 DO 1005 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))1005,1005,1002
 1002   M=M+1
        DO 1003 J=1,NUMP
          V(M,J)=PWR(I,J)
 1003     CONTINUE
        DO 1004 J=NMPP1,NMPPQ
          V(M,J)=-UPTBL(I)*PWR(I,J)
 1004     CONTINUE
        V(M,N)=0.0
 1005   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
 1006 IF(IOPT-(IOPT/100)*100-10)1012,1007,1007
 1007 DO 1011 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))1011,1011,1008
 1008   M=M+1
        DO 1009 J=1,NUMP
          V(M,J)=-PWR(I,J)
 1009     CONTINUE
        DO 1010 J=NMPP1,NMPPQ
          V(M,J)=ALTBL(I)*PWR(I,J)
 1010     CONTINUE
        V(M,N)=0.0
 1011   CONTINUE
C
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).
 1012 DO 1016 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))1016,1016,1013
 1013   M=M+2
        MM1=M-1
        DO 1014 J=1,NUMP
          V(MM1,J)=-WTBLE(I)*PWR(I,J)
          V(M,J)=-V(MM1,J)
 1014     CONTINUE
        DEL1=-DELTK+FTBWT(I)
        DEL2=-DELTK-FTBWT(I)
        DO 1015 J=NMPP1,NMPPQ
          V(MM1,J)=DEL1*PWR(I,J)
          V(M,J)=DEL2*PWR(I,J)
 1015     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)
 1016   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 1017 I=1,M
        V(I,NP1)=0.0
 1017   CONTINUE
C
C NOW SET THE CONSTRAINTS OF THE FORM -Q(J).LE.1.0
C AND Q(J).LE.1.0
      DO 1019 I=1,NUMQ
        M=M+2
        MM1=M-1
        DO 1018 J=1,N
          V(MM1,J)=0.0
          V(M,J)=0.0
 1018     CONTINUE
        NNN=NUMP+I
        V(MM1,NNN)=-1.0
        V(M,NNN)=1.0
        V(MM1,NP1)=1.0
        V(M,NP1)=1.0
 1019   CONTINUE
C
C IF THE ONES DIGIT OF IOPT EQUALS 1, THEN FORCE Q.GE.EPSQ.
      IF(IOPT-(IOPT/10)*10)10010,10010,10000
10000 DO 1022 I=1,NUMGR
        M=M+1
        DO 1020 J=1,N
          V(M,J)=0.0
 1020     CONTINUE
        DO 1021 J=NMPP1,NMPPQ
          V(M,J)=-PWR(I,J)
 1021     CONTINUE
        V(M,NP1)=-EPSQ
 1022   CONTINUE
C
C NOW SET THE BOTTOM ROW.  TO MINIMIZE W WE MAXIMIZE -W.
10010 MP1=M+1
      DO 10 J=1,NP1
        V(MP1,J)=0.0
   10   CONTINUE
      V(MP1,N)=1.0
      RETURN
      END
    1    2    3    1  101
      1010
                 0.0                 2.0
    1    7    8    1   21
     21010
                 0.03.141592653589793238
    1    7    8    1   41
    121010
                 0.03.141592653589793238
    2    3    2    5    5
      1000
                 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
    0    0    0    0    0