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
.
