[CONTACT]

[ABOUT]

[POLICY]

DIFCOR differential correction algori

Found at: ftp.icm.edu.pl:70/packages/netlib/a/difcor

From: TAYLOR%CSUGREEN.BITNET@WISCVM.ARPA  (G. D. TAYLOR, MATH, COLORADO STATE UNIV.)
Subject:  DIFCOR differential correction algorithm 
 


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
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
        IF(INDF(I))10020,10020,10010
      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.
      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
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.
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
      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
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.
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.
        IF(V(2*I,NP1)-EPSQ+EPRES)10140,10120,10120
C HERE NO USER SUPPLIED CONSTRAINT WAS VIOLATED BY MORE THEN EPRES.
C IF WE ARE IN THE FINAL CHECK WE RETURN.
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.
      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
      GO TO 1026
C HERE ITLIM ITERATIONS HAVE BEEN PERFORMED AND WE WILL RETURN
C AFTER SETTING IFLAG AND INDLP (SEE BELOW).
      IF(INDIC)10170,10180,10180
      GO TO 10300
      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
      GO TO 10187
 1046 IADJS=2-IADJS/2
C RESTORE THE PREVIOUS IYRCT FOR THE NEXT ATTEMPT.
        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).
C SET UP TO TRY THE LOEB INITIALIZATION.
      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.
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.
      IF(DLOLD)10240,10260,10260
      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
      RETURN
      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.
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.
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.
      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
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.
C END OF SETTING MACHINE DEPENDENT PARAMETERS FOR SSETV1.
C
      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
C SET THE CONSTRAINTS OF THE FORM -Q(J) .LE. 1.0 AND Q(J)
C .LE. 1.0 FOR J=2,...,NUMQ.
        M=M+2
        MM1=M-1
        DO 10070 J=1,N1
          V(MM1,J)=0.0
          V(M,J)=0.0
        NNN1=NUMP+I-1
        V(MM1,NNN1)=-1.0
        V(M,NNN1)=1.0
        V(MM1,N1P1)=1.0
        V(M,N1P1)=1.0
C
C NOW SET THE BOTTOM ROW.  TO MINIMIZE W WE MAXIMIZE -W.
      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
      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.
      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
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.
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.
C RESET 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).
        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.
      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
        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.
      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


.

NEW PAGES:

[ODDNUGGET]

[GOPHER]