# 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
```

```.
```