From TAYLOR%CSUGREEN.BITNET@WISCVM.ARPA Thu Feb 19 14:53:43 1987 Return-Path: Received: from anl-mcs.ARPA by dasher.sun.com (3.2/SMI-3.2) id AA23328; Thu, 19 Feb 87 14:53:27 CST Received: from wiscvm.wisc.edu (wiscvm.wisc.edu.ARPA) by anl-mcs.ARPA (4.12/4.9) id AA24730; Thu, 19 Feb 87 14:47:31 cst Received: from CSUGREEN.BITNET by wiscvm.wisc.edu on 02/19/87 at 12:58:13 CST Message-Id: <870219114541.00000ECD.MKNK.AB@CSUGREEN> (UMass-Mailer 4.03) Date: Thu, 19 Feb 87 11:45:41 MST From: TAYLOR%CSUGREEN.BITNET@WISCVM.ARPA (G. D. TAYLOR, MATH, COLORADO STATE UNIV.) Subject: DIFCOR To: dongarra%dasher@anl-mcs.arpa In-Reply-To: Re: differential correction algorithm <8702190444.AA22025@dasher.su Status: RO Jack, attached is the code SDIFCOR complete with driver and input at end of code. I ran it on our CYBER 830 before sending it to be safe and it ran. (This the most recent version of the code and I had just received it from Kaufman via a tape.) Let me if I can be of any help. Jerry PROGRAM DRIVR(INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT) DIMENSION X(16),PWR(101,15),FTBLE(101),T(202), * WTBLE(101),UPTBL(101),INDUP(101),ALTBL(101), * INDLW(101),INDF(101),ERROR(102),ERTST(7) C C THE APPROXIMATION PROBLEMS RUN WITH THIS DEMONSTRATION DRIVER C ARE AS FOLLOWS. C C IN JOB 1, THE GRID IS A 101 POINT EQUALLY SPACED SUBDIVISION OF C THE CLOSED INTERVAL (0.0,2.0). F(Z) IS 1.0 FOR 0.0 .LE. Z .LE. 1.0, C UNDEFINED FOR 1.0 .LT. Z .LT. 1.2, AND 0.0 FOR 1.2 .LE. Z .LE. 2.0. C (P/Q)(Z) IS (P(1)+P(2)*Z)/(Q(1)+Q(2)*Z+Q(3)*Z**2). WE REQUIRE C (P/Q)(Z) .GE. 0.9999 FOR 0.0 .LE. Z .LE. 0.3 AND (P/Q)(Z) .GE. 0.0 FOR C 1.0 .LT. Z .LE. 2.0. C C IN JOB 2, THE GRID IS A 21 POINT EQUALLY SPACED SUBDIVISION OF C THE CLOSED INTERVAL (0.0,PI). F(Z) IS 1.0 FOR 0.0 .LE. Z .LE. PI/2, C UNDEFINED FOR PI/2 .LT. Z .LT. 3*PI/5, AND 0.0 FOR 3*PI/5 .LE. Z .LE. PI. C (P/Q)(Z) IS (P(1)+P(2)*COSZ+...+P(7)*COS6Z)/(Q(1)+Q(2)*Z+...+ C Q(8)*COS7Z). THE WEIGHT FUNCTION (FOR WEIGHTED ERROR CURVE C APPROXIMATION) IS 0.5 FOR 0.0 .LE. Z .LE. PI/2, UNDEFINED FOR C PI/2 .LT. Z .LT. 3*PI/5, AND 1.0 FOR 3*PI/5 .LE. Z .LE. PI. WE REQUIRE C (P/Q)(Z) .GE. 0.0 FOR PI/2 .LT. Z .LE. PI. C C JOB 3 IS THE SAME AS JOB 2 EXCEPT THE GRID IS A 41 POINT EQUALLY C SPACED SUBDIVISION OF THE CLOSED INTERVAL (0.0,PI), AND WE USE C THE P/Q PRODUCED IN JOB 2 FOR INITIALIZATION. C C IN JOB 4 THE GRID IS THE SET OF ORDERED PAIRS (Z1,Z2) WHERE Z1 C AND Z2 ARE NONNEGATIVE INTEGERS WITH Z1**2+Z2**2 .LE. 16.0. C F(Z1,Z2) IS SQRT(16.0-Z1**2-Z2**2). (P/Q)(Z1,Z2) IS (P(1)+ C P(2)*Z1**2+P(3)*Z2**2)/(Q(1)+Q(2)*(Z1*Z2)**2). C C IN JOB 5 THE GRID IS A 21 POINT EQUALLY SPACED SUBDIVISION OF THE C CLOSED INTERVAL (0.0,2.0). F(Z) IS 1.0 FOR Z=0.0 AND 0.0 FOR C 0.0 .LT. Z .LE. 2.0. (P/Q)(Z) IS P(1)/(Q(1)+Q(2)*Z). THIS IS AN C EXAMPLE WHERE NO BEST APPROXIMATION EXISTS. C C JOB 6 IS THE SAME AS JOB 5 EXCEPT WE REQUIRE Q(Z) .GE. 10.0**(-7) ON C THE GRID. HERE A BEST APPROXIMATION DOES EXIST. C C IN JOB 7 THE GRID IS A 101 POINT EQUALLY SPACED SUBDIVISION OF THE C CLOSED INTERVAL (0.0,2.0). F(Z) IS 3.0/(1.0+Z) + G(Z), WHERE C G HAS VALUES 0.2, -0.2, 0.2, AND -0.2 AT Z = 0.0, 0.2, 1.0, C AND 2.0 RESPECTIVELY, AND G IS LINEAR BETWEEN THESE POINTS. C (P/Q)(Z) IS (P(1)+P(2)*Z)/(Q(1)+Q(2)*Z+Q(3)*Z**2). THIS IS AN C EXAMPLE WHERE THE BEST APPROXIMATION IS DEGENERATE. C C THE DATA CARDS NEEDED FOR THESE JOBS ARE AS FOLLOWS (WITH THE C CS IN THE FIRST COLUMN REPLACED BY BLANKS). C 1 2 3 1 101 C 1010 C 0.0 2.0 C 1 7 8 1 21 C 21010 C 0.03.141592653589793238 C 1 7 8 1 41 C 121010 C 0.03.141592653589793238 C 2 3 2 5 5 C 1000 C 0.0 0.0 4.0 4.0 C 1 1 2 1 21 C 0 C 0.0 2.0 C 1 1 2 1 21 C 1 C 0.0 2.0 C 1 2 3 1 101 C 0 C 0.0 2.0 C 0 0 0 0 0 C SET MACHINE DEPENDENT PARAMETERS FOR THE DEMONSTRATION DRIVER. C SET INPUT AND OUTPUT UNIT NUMBERS. NREAD=I1MACH(1) NWRIT=I1MACH(2) C SET SPCMN=B**(-T), WHERE B IS THE BASE AND T IS THE NUMBER C OF BASE B DIGITS IN THE MANTISSA. SPCMN IS THE MINIMUM C RELATIVE SPACING ABS((X1-X2)/X2) BETWEEN TWO SUCCESSIVE C FLOATING POINT NUMBERS, SO IT IS THE SPACING BETWEEN TWO C SUCCESSIVE FLOATING POINT NUMBERS IN THE CLOSED INTERVAL C (0.1,1.0). WE ALSO HAVE SPCMN=10.0**(-T*(LOG10)(B))= C 10.0**(-TNMAN), WHERE TNMAN IS THE BASE 10 EQUIVALENT OF C THE LENGTH OF THE MANTISSA. SPCMN=R1MACH(3) C SET THE ERROR NORM TEST TOLERANCE ERTOL=MAX(10.0**(-12), C 100.0*SPCMN). ERTOL=100.0*SPCMN IF(ERTOL-1.0E-12)10,20,20 10 ERTOL=1.0E-12 20 JOBNM=1 C SET THE ERROR NORMS COMPUTED WITH THE CDC CYBER 172 AT C CENTRAL MICHIGAN UNIVERSITY FOR COMPARISON PURPOSES. ERTST(1)=0.36955949788958 ERTST(2)=0.879206364E-5 ERTST(3)=0.5440168621E-4 ERTST(4)=0.44787565553231 ERTST(5)=0.463E-11 ERTST(6)=0.99999800E-6 ERTST(7)=0.20000000000001 C READ NUMDM, NUMP=NUMBER OF COEFFICIENTS IN THE NUMERATOR, C NUMQ=NUMBER OF COEFFICIENTS IN THE DENOMINATOR, NRWGR= C NUMBER OF ROWS IN THE GRID, NCLGR=NUMBER OF C COLUMNS IN THE GRID. 30 READ(NREAD, 40 )NUMDM,NUMP,NUMQ,NRWGR,NCLGR 40 FORMAT(5I5) NMPP1=NUMP+1 NMPPQ=NUMP+NUMQ NUMG1=NMPPQ IF(NUMDM) 60 , 50 , 60 C A DATA CARD WITH 0 FOR NUMDM IS USED AS A SIGNAL TO STOP. 50 STOP 60 WRITE(NWRIT, 70 )JOBNM 70 FORMAT(///30H ***********THIS IS JOB NUMBER,I4) NUMGR=NRWGR*NCLGR C READ THE OPTION SWITCH IOPT. INSTRUCTIONS FOR SETTING IT C ARE IN THE INITIAL COMMENTS OF SDFCOR. READ(NREAD, 80)IOPT 80 FORMAT(I10) C READ THE ENDPOINTS INTO ASWP1, ASWP2 IN THE ONE DIMENSIONAL CASE. C IN THE TWO DIMENSIONAL CASE ASWP1, ASWP2 AND ANEP1, ANEP2 ARE THE C COORDINATES OF THE LOWER LEFT AND UPPER RIGHT CORNERS OF C THE GRID RESPECTIVELY. IF(NUMDM-1)90,90,110 90 READ(NREAD,100)ASWP1,ASWP2 100 FORMAT(2F20.10) GO TO 130 110 READ(NREAD, 120)ASWP1,ASWP2,ANEP1,ANEP2 120 FORMAT(4F20.10) C PUT THE GRID POINTS INTO T. 130 CALL STCOMP(NUMDM,NRWGR,NCLGR,ASWP1,ASWP2,ANEP1, 1ANEP2,T) C SET UP PWR AND FTBLE. I=0 IDBM=-1 IDB=0 140 I=I+1 IF(I-NUMGR) 150, 150, 540 150 IDBM=IDBM+2 IDB=IDB+2 C COMPUTE THE BASIS FUNCTIONS AND F AND, IF NECESSARY, C WTBLE, INDUP, UPTBL, INDLW, ALTBL, AND INDF. C IN THE ONE DIMENSIONAL CASE THE COORDINATE WILL BE T(I), C WHILE IN THE TWO DIMENSIONAL CASE THE COORDINATES WILL BE C (T(IDBM),T(IDB)). GO TO ( 160, 230, 230, 300, 160, 160, 160),JOBNM C INPUT OF BASIS FUNCTIONS FOR JOBS 1, 5, 6, AND 7. 160 PWR(I,1)=1.0 IF(NUMP-1) 190, 190, 170 170 DO 180 J=2,NUMP PWR(I,J)=(T(I))**(J-1) 180 CONTINUE 190 PWR(I,NMPP1)=1.0 IF(NUMQ-1) 220, 220, 200 200 DO 210 J=2,NUMQ JJ=NUMP+J PWR(I,JJ)=(T(I))**(J-1) 210 CONTINUE 220 GO TO ( 310, 310, 310, 310, 460, 460, 490),JOBNM C INPUT OF BASIS FUNCTIONS FOR JOBS 2 AND 3. 230 PWR(I,1)=1.0 IF(NUMP-1) 260, 260, 240 240 DO 250 J=2,NUMP PWR(I,J)=COS((J-1)*T(I)) 250 CONTINUE 260 PWR(I,NMPP1)=1.0 IF(NUMQ-1) 290, 290, 270 270 DO 280 J=2,NUMQ JJ=NUMP+J PWR(I,JJ)=COS((J-1)*T(I)) 280 CONTINUE 290 CONTINUE GO TO 380 C INPUT OF BASIS FUNCTIONS FOR JOB 4. 300 PWR(I,1)=1.0 PWR(I,2)=T(IDBM)**2 PWR(I,3)=T(IDB)**2 PWR(I,4)=1.0 PWR(I,5)=PWR(I,2)*PWR(I,3) GO TO 430 C INPUT OF F, INDF, AND LOWER CONSTRAINTS FOR JOB 1. 310 NUM1=3*(NUMGR-1)/20+1 NUM2=(NUMGR-1)/2+1 NUM3=3*(NUMGR-1)/5+1 IF(I-NUM1) 320, 320, 330 320 INDF(I)=1 FTBLE(I)=1.0 INDLW(I)=1 ALTBL(I)=0.9999 GO TO 140 330 IF(I-NUM2) 340, 340, 350 340 INDF(I)=1 FTBLE(I)=1.0 INDLW(I)=0 GO TO 140 350 IF(I-NUM3) 360, 370, 370 360 INDF(I)=0 INDLW(I)=1 ALTBL(I)=0.0 GO TO 140 370 INDF(I)=1 FTBLE(I)=0.0 INDLW(I)=1 ALTBL(I)=0.0 GO TO 140 C INPUT OF F, INDF, WEIGHT FUNCTION, AND LOWER CONSTRAINTS C FOR JOBS 2 AND 3. 380 NUM1=(NUMGR-1)/2+1 NUM2=3*(NUMGR-1)/5+1 IF(I-NUM1) 390, 390, 400 390 INDF(I)=1 FTBLE(I)=1.0 WTBLE(I)=0.5 INDLW(I)=0 GO TO 140 400 IF(I-NUM2) 410, 420, 420 410 INDF(I)=0 INDLW(I)=1 ALTBL(I)=0.0 GO TO 140 420 INDF(I)=1 FTBLE(I)=0.0 WTBLE(I)=1.0 INDLW(I)=1 ALTBL(I)=0.0 GO TO 140 C INPUT OF F AND INDF FOR JOB 4. C WE USE 16.0001 INSTEAD OF 16.0 IN THE NEXT TEST TO AVOID TROUBLE C DUE TO ROUND OFF ERROR. 430 IF(T(IDBM)**2+T(IDB)**2-16.0001) 440, 440, 450 440 INDF(I)=1 FTBLE(I)=SQRT(16.0-T(IDBM)**2-T(IDB)**2) GO TO 140 450 INDF(I)=0 GO TO 140 C INPUT OF F FOR JOBS 5 AND 6. WE ALSO INPUT EPSQ HERE, ALTHOUGH C IT IS NOT ACTUALLY USED IN JOB 5. 460 EPSQ=1.0E-7 IF(I-1)470,470,480 470 FTBLE(I)=1.0 GO TO 140 480 FTBLE(I)=0.0 GO TO 140 C INPUT OF F FOR JOB 7. 490 NUM1=(NUMGR-1)/10+1 NUM2=(NUMGR-1)/2+1 IF(I-NUM1) 500, 500, 510 500 FTBLE(I)=3.0/(1.0+T(I))-2.0*T(I)+0.2 GO TO 140 510 IF(I-NUM2) 520, 520, 530 520 FTBLE(I)=3.0/(1.0+T(I))+0.5*T(I)-0.3 GO TO 140 530 FTBLE(I)=3.0/(1.0+T(I))-0.4*T(I)+0.6 C END COMPUTING OF BASIS FUNCTIONS AND F. GO TO 140 540 CALL SDFCOR(NUMP,NUMQ,NUMGR,FTBLE,PWR,IOPT,EPSQ,INDLW,ALTBL, *INDUP,UPTBL,INDF,WTBLE,X,ERROR,QMIN,IQMIN,ITER,INDLP,IFLAG) C CHECK THE OUTPUT. FIRST MAKE SURE IFLAG IS NONPOSITIVE. IF(IFLAG)570,570,550 550 WRITE(NWRIT,560) 560 FORMAT(/42H THE PROGRAM FAILED IN THE INITIALIZATION,, *28H SO WE GIVE A FULL PRINTOUT.) GO TO 620 C CHECK TO SEE IF THE COMPUTED ERROR NORM IS CORRECT WITHIN ERTOL. 570 IF(ABS(ERROR(NUMGR+1)-ERTST(JOBNM))-ERTOL)580,580,600 580 WRITE(NWRIT,590) 590 FORMAT(/33H THE PROGRAM HAS WORKED NORMALLY.) C C**********MODIFICATION IF A FULL PRINTOUT IS DESIRED EVEN IF THE C PROGRAM SUCCEEDS, THE NEXT STATEMENT SHOULD BE GO TO 620 C OTHERWISE THE NEXT STATEMENT SHOULD BE GO TO 790 GO TO 620 C**********END MODIFICATION C 600 WRITE(NWRIT,610)ERTOL 610 FORMAT(/36H THE ERROR NORM WAS OFF BY MORE THAN,E12.2, *30H SO WE GIVE A FULL PRINTOUT.) 620 WRITE(NWRIT, 630)NUMDM 630 FORMAT(/33H WE ARE DEALING WITH FUNCTIONS OF,I4, 112H VARIABLES) WRITE(NWRIT, 640)NUMP,NUMQ,NUMGR,NRWGR,NCLGR 640 FORMAT(/6H P HAS,I4,24H COEFFICIENTS Q HAS,I4, 131H COEFFICIENTS THE GRID HAS,I5, 218H POINTS ARRANGED,I5,5H BY,I5) WRITE(NWRIT, 650)IOPT 650 FORMAT(/8H IOPT IS,I10) IF(NUMDM-2) 660, 680, 660 660 WRITE(NWRIT, 670)ASWP1,ASWP2 670 FORMAT(/22H THE LEFT END POINT IS,E24.14, 125H THE RIGHT END POINT IS,E24.14) GO TO 700 680 WRITE(NWRIT, 690)ASWP1,ASWP2,ANEP1,ANEP2 690 FORMAT(/40H THE COORDINATES OF THE SW AND NE POINTS, 14H ARE,2E15.5,6H AND,2E15.5) C WRITE IFLAG AND INDLP. 700 WRITE(NWRIT, 710)IFLAG,INDLP 710 FORMAT(/9H IFLAG IS,I5,13H INDLP IS,I4) C WRITE THE DENOMINATOR MINIMUM AND THE INDEX OF THE POINT C WHERE IT OCCURRED. WRITE(NWRIT,720)QMIN,IQMIN 720 FORMAT(/40H THE MINIMUM VALUE OF THE DENOMINATOR IS, *E25.15,25H, ACHIEVED AT GRID POINT,I5) C WRITE ITER AND DELTK. DELTK=ERROR(NUMGR+1) WRITE(NWRIT, 730)ITER,DELTK 730 FORMAT(/21H THE ERROR NORM AFTER,I5,13H ITERATIONS, 13H IS,E25.15) C WRITE THE COEFFICIENTS OF THE NUMERATOR AND DENOMINATOR WRITE(NWRIT, 740)(X(I),I=1,NUMP) 740 FORMAT(/26H THE COEFFICIENTS OF P ARE/(/5E24.14)) WRITE(NWRIT, 750)(X(I),I=NMPP1,NMPPQ) 750 FORMAT(/26H THE COEFFICIENTS OF Q ARE/(/5E24.14)) WRITE(NWRIT, 760) 760 FORMAT(/15H THE ERRORS ARE) DO 780 I=1,NRWGR IND1=1+(I-1)*NCLGR IND2=I*NCLGR WRITE(NWRIT, 770)(ERROR(IT),IT=IND1,IND2) 770 FORMAT(/6E20.10) 780 CONTINUE 790 JOBNM=JOBNM+1 GO TO 30 END FUNCTION I1MACH(I) C C THIS SUBPROGRAM IS NOT TO BE INCLUDED IN THE SLATEC LIBRARY. C IT APPEARS HERE FOR TESTING PURPOSES TO MIMIC THE BEHAVIOR OF C A FUNCTION ALREADY IN THE LIBRARY FOR ASSIGNING MACHINE C DEPENDENT PARAMETERS, IN THE CASE OF THE CDC CYBER 172 AT C CENTRAL MICHIGAN UNIVERSITY. C C I1MACH(1) IS THE INPUT UNIT NUMBER, AND I1MACH(2) IS THE OUTPUT C UNIT NUMBER. IF(I-1)1,1,2 1 I1MACH=5 RETURN 2 I1MACH=6 RETURN END FUNCTION R1MACH(I) C C THIS SUBPROGRAM IS NOT TO BE INCLUDED IN THE SLATEC LIBRARY. C IT APPEARS HERE FOR TESTING PURPOSES TO MIMIC THE BEHAVIOR C OF A FUNCTION ALREADY IN THE LIBRARY FOR ASSIGNING MACHINE C DEPENDENT PARAMETERS, IN THE CASE OF THE CDC CYBER 172 AT C CENTRAL MICHIGAN UNIVERSITY. C C R1MACH(3) IS B**(-T), WHERE B IS THE BASE FOR FLOATING POINT NUMBERS C AND T IS THE NUMBER OF BASE B DIGITS IN THE MANTISSA. R1MACH=2.0**(-48) RETURN END SUBROUTINE SDFCOR(NUMP,NUMQ,NUMGR,FTBLE,PWR,IOPT,EPSQ,INDLW, *ALTBL,INDUP,UPTBL,INDF,WTBLE,X,ERROR,QMIN,IQMIN,ITER,INDLP,IFLAG) C***BEGIN PROLOGUE SDFCOR C***DATE WRITTEN 720120 (YYMMDD) C***REVISION DATE 860328 (YYMMDD) C***CATEGORY NO. E2C1 C***KEYWORDS APPROXIMATION,CHEBYSHEV,RATIONAL,UNIFORM,BEST, C GENERALIZED,DIFFERENTIAL CORRECTION,MULTIVARIATE, C CONSTRAINTS,RESTRICTED RANGE,WEIGHTED C***AUTHOR KAUFMAN, EDWIN H. JR., C DEPARTMENT OF MATHEMATICS C CENTRAL MICHIGAN UNIVERSITY C MOUNT PLEASANT, MICHIGAN 48859 C TAYLOR, GERALD D., C DEPARTMENT OF MATHEMATICS C COLORADO STATE UNIVERSITY C FORT COLLINS, COLORADO 80523 C***PURPOSE THIS SUBROUTINE COMPUTES A BEST UNIFORM GENERALIZED C RATIONAL APPROXIMATION P/Q TO A FUNCTION F DEFINED C ON A FINITE SET OF GRID POINTS, WITH THE POSSIBILITY C OF INCLUDING A WEIGHT FUNCTION AND CONSTRAINTS ON C P/Q AND Q. C***DESCRIPTION C C THE VARIABLES IN THE CALLING SEQUENCE ARE AS FOLLOWS. NONE OF C THE STRICTLY INPUT VARIABLES IS CHANGED BY THIS PROGRAM. C C NUMP (INPUT) THIS IS THE NUMBER OF BASIS FUNCTIONS FOR THE C NUMERATOR. C C NUMQ (INPUT) THIS IS THE NUMBER OF BASIS FUNCTIONS FOR THE C DENOMINATOR. C C NUMGR (INPUT) THIS IS THE NUMBER OF GRID POINTS. C C FTBLE (INPUT) FTBLE(I) IS THE VALUE OF THE FUNCTION F AT GRID C POINT I (I=1,...,NUMGR). C C PWR (INPUT) PWR(I,J) IS THE VALUE OF THE JTH BASIS FUNCTION C (WHERE THE NUMERATOR BASIS FUNCTIONS PRECEDE THE C DENOMINATOR BASIS FUNCTIONS) AT GRID POINT I (J=1,..., C NUMP+NUMQ, I=1,...,NUMGR). C C IOPT (INPUT) THIS IS THE OPTION SWITCH. SETTING IOPT=0 TURNS C OFF ALL EXTRA OPTIONS. THE PROPER SETTING OF IOPT TO C ACTIVATE ONE OR MORE OF THE EXTRA OPTIONS IS DISCUSSED IN C THE LONG DESCRIPTION BELOW. C C EPSQ, INDLW, ALTBL, INDUP, UPTBL, INDF, WTBLE (OPTIONAL INPUT) C THESE SEVEN VARIABLES ARE USED FOR INPUT FOR THE EXTRA C OPTIONS DESCRIBED IN THE LONG DESCRIPTION BELOW. NONE OF C THEM NEED BE ASSIGNED VALUES IF IOPT=0. C C X (OUTPUT AND OPTIONAL INPUT) THE PROGRAM WILL PLACE THE C COMPUTED COEFFICIENT OF THE JTH BASIS FUNCTION IN X(J) C (J=1,...,NUMP+NUMQ). X IS ALSO USED IF THE USER WISHES C TO SUPPLY AN INITIAL APPROXIMATION, IN WHICH CASE IOPT C MUST BE SET APPROPRIATELY (SEE LONG DESCRIPTION BELOW). C C ERROR (OUTPUT) THE PROGRAM WILL PLACE THE ERROR AT GRID POINT C I, WHICH WILL BE (F-P/Q)(GRID POINT I) IF NEITHER WEIGHT C FUNCTION OPTION IS USED (SEE LONG DESCRIPTION BELOW FOR C THE WEIGHT FUNCTION OPTIONS), IN ERROR(I) (I=1,...,NUMGR). C THE UNIFORM ERROR NORM WILL BE PLACED IN ERROR(NUMGR+1). C C QMIN (OUTPUT) THE PROGRAM WILL PLACE THE MINIMUM VALUE OF THE C DENOMINATOR ON THE GRID IN QMIN. C C IQMIN (OUTPUT) THE PROGRAM WILL PLACE THE INDEX OF THE GRID POINT C AT WHICH THE MINIMUM VALUE OF THE DENOMINATOR IS ACHIEVED C IN IQMIN. C C ITER (OUTPUT) THE PROGRAM WILL PLACE THE NUMBER OF ITERATIONS C (NOT COUNTING INITIALIZATION) REQUIRED IN ITER. C C INDLP (OUTPUT) THE PROGRAM WILL PLACE INFORMATION ABOUT THE C PERFORMANCE OF THE LINEAR PROGRAMMING SUBROUTINE SLNPRO C IN THE TWO DIGIT VARIABLE INDLP (SEE COMMENTS AT THE END OF C THIS SUBROUTINE AND THE BEGINNING OF SUBROUTINE SLNPRO). C THE USER NORMALLY NEED NOT CONSIDER INDLP UNLESS IFLAG (SEE C BELOW) OR SOMETHING ELSE INDICATES SOMETHING IS WRONG. IN C THAT CASE A CAREFUL STUDY OF THE INDICATED COMMENTS MAY C PROVIDE A CLUE AS TO WHAT HAPPENED. C C IFLAG (OUTPUT) IFLAG IS THE ERROR FLAG FOR THE PROGRAM. IFLAG=0 C MEANS THE PROGRAM APPEARED TO WORK NORMALLY, WITH NONE OF C THE CONDITIONS DESCRIBED BELOW FOR IFLAG .NE. O OCCURRING. C IFLAG POSITIVE IS A WARNING OF FAILURE IN THE INITIALIZATION. C IFLAG NEGATIVE INDICATES CAUTION SHOULD BE USED, ALTHOUGH C THE RETURNED APPROXIMATION MAY WELL BE ACCEPTABLE. FOR C DETAILS CONCERNING SPECIFIC NONZERO VALUES OF IFLAG, SEE C THE LONG DESCRIPTION BELOW. C C THE DIMENSIONS OF THE ARRAYS IN THE CALLING SEQUENCE FOR C SDFCOR ARE C C FTBLE(101), PWR(101,15), INDLW(101), ALTBL(101), INDUP(101), C UPTBL(101), INDF(101), WTBLE(101), X(16), ERROR(102). C C IN ADDITION, IF SUBROUTINE STCOMP IS TO BE USED, C DIMENSION T(202) MUST BE ADDED. C C THESE DIMENSIONS ARE SET SO THAT THE PROGRAM CAN BE APPLIED TO C ANY PROBLEM WITH UP TO 15 BASIS FUNCTIONS (SO NUMP+NUMQ .LE. 15) C AND UP TO 101 GRID POINTS (SO NUMGR .LE. 101), EVEN IF ALL THE C EXTRA OPTIONS DESCRIBED IN THE LONG DESCRIPTION BELOW ARE C ACTIVATED. SEE THE LONG DESCRIPTION FOR THE DIMENSION C PATTERN REQUIRED IN THE DRIVER PROGRAM AND SUBROUTINES IF C MORE BASIS FUNCTIONS OR GRID POINTS ARE DESIRED, AND FOR THE C DIMENSIONS OF THE INTERNAL WORK ARRAYS. C C***LONG DESCRIPTION C C THE USER CAN, BY SETTING IOPT AND THE CORRESPONDING OPTIONAL C VARIABLES, USE ANY COMBINATION OF THE SIX EXTRA OPTIONS C DESCRIBED BELOW. AN OPTIONAL VARIABLE NEED BE ASSIGNED A VALUE C BY THE USER ONLY IF ITS OPTION IS TO BE USED. SETTING IOPT=0 C TURNS OFF ALL THE EXTRA OPTIONS. C C OPTION 1 CONSTRAINTS (DENOMINATOR) C C TO ACTIVATE THIS OPTION, SET THE ONES DIGIT OT IOPT C TO 1 AND SET EPSQ AS FOLLOWS. C C EPSQ (INPUT) THIS IS THE CONSTANT LOWER BOUND ON THE VALUES C OF THE DENOMINATOR. THUS THE PROGRAM WILL SET UP C EXTRA CONSTRAINTS IN SUBROUTINE SLNPRO TO FORCE C Q .GE. EPSQ EVERYWHERE ON THE GRID. THIS OPTION IS C USEFUL BECAUSE THE COMPUTED APPROXIMATION P/Q MAY BE OF C MORE VALUE IF ITS DENOMINATOR IS KEPT FARTHER AWAY FROM C ZERO, AND ALSO BECAUSE IN DIFFICULT CASES THIS OPTION C MAY ALLEVIATE NUMERICAL DIFFICULTIES IN THE PROGRAM. C C OPTION 2 CONSTRAINTS (LOWER RESTRICTED RANGE) C C TO ACTIVATE THIS OPTION, SET THE TENS DIGIT OF IOPT C TO 1 AND SET INDLW AND ALTBL AS FOLLOWS. C C INDLW (INPUT) INDLW(I)=1 IF THERE IS A LOWER RESTRICTED C RANGE CONSTRAINT AT GRID POINT I, AND INDLW(I)=0 C OTHERWISE (I=1,...,NUMGR). C C ALTBL (INPUT) ALTBL(I) CONTAINS THE LOWER BOUND AT GRID POINT I C (THUS WE ARE FORCING (P/Q)(GRID POINT I) .GE. ALTBL(I)) C IF THERE IS ONE, AND ALTBL(I) NEED NOT BE ASSIGNED A C VALUE OTHERWISE (I=1,...,NUMGR). C C OPTION 3 CONSTRAINTS (UPPER RESTRICTED RANGE) C C TO ACTIVATE THIS OPTION, SET THE HUNDREDS DIGIT OF C IOPT TO 1 AND SET INDUP AND UPTBL AS FOLLOWS. C C INDUP (INPUT) INDUP(I)=1 IF THERE IS AN UPPER RESTRICTED C RANGE CONSTRAINT AT GRID POINT I, AND INDUP(I)=0 C OTHERWISE (I=1,...,NUMGR). C C UPTBL (INPUT) UPTBL(I) CONTAINS THE UPPER BOUND AT GRID POINT I C (THUS WE ARE FORCING (P/Q)(GRID POINT I) .LE. UPTBL(I)) C IF THERE IS ONE, AND UPTBL(I) NEED NOT BE ASSIGNED A VALUE C OTHERWISE (I=1,...,NUMGR). C C OPTION 4 IGNORE F AT SPECIFIED POINTS C C TO ACTIVATE THIS OPTION, SET THE THOUSANDS DIGIT OF C IOPT TO 1 AND SET INDF AS FOLLOWS. C C INDF (INPUT) INDF(I)=1 IF F IS TO BE APPROXIMATED AT GRID C POINT I, AND INDF(I)=0 IF F IS TO BE IGNORED AT GRID C POINT I (ALTHOUGH OPTIONS 1, 2, AND 3 CAN BE USED AT C SUCH POINTS) (I=1,...,NUMGR). NEITHER FTBLE(I) NOR C WTBLE(I) NEED BE ASSIGNED A VALUE IF INDF(I)=0. C C OPTION 5 WEIGHTED APPROXIMATION C C TO ACTIVATE STRAIGHT WEIGHTED APPROXIMATION, THAT IS, C TO MINIMIZE MAX(ABS(F-WT*P/Q)), SET THE TEN THOUSANDS C DIGIT OF IOPT TO 1 AND SET WTBLE AS BELOW. C C TO ACTIVATE WEIGHTED ERROR CURVE APPROXIMATION, THAT IS, C TO MINIMIZE MAX(ABS(WT*(F-P/Q))), SET THE TEN THOUSANDS C DIGIT OF IOPT TO 2 AND SET WTBLE AS BELOW. C C WTBLE (INPUT) WTBLE(I) IS THE VALUE OF THE WEIGHT FUNCTION WT C (WHICH NEED NOT BE POSITIVE) AT GRID POINT I (I=1,...,NUMGR). C C OPTION 6 USER SUPPLIED INITIAL APPROXIMATION C C TO ACTIVATE THIS OPTION, SET THE HUNDRED THOUSANDS DIGIT C OF IOPT TO 1 AND PLACE THE NUMP+NUMQ INITIAL C COEFFICIENTS IN X (WHERE X IS DEFINED IN THE C DESCRIPTION SECTION ABOVE). SOMETIMES ROUND OFF ERROR C CAN BE ALLEVIATED, AND COMPUTER TIME CAN BE SAVED, BY C FIRST RUNNING A PROBLEM WITH A COARSE GRID AND THEN C USING THE RESULT TO INITIALIZE THE PROBLEM ON THE C DESIRED GRID. C C THE MEANINGS OF SPECIFIC NONZERO VALUES OF IFLAG ARE AS FOLLOWS. C C IFLAG POSITIVE (INITIALIZATION FAILURE) RARELY OCCURS IF THE C DATA FOR THE PROGRAM IS ENTERED CORRECTLY. THERE ARE TWO C POSSIBILITIES. C C IFLAG=3333 MEANS THE APPROXIMATION RETURNED TO THE USER IS C THE INITIAL APPROXIMATION P/Q=0.0/PWR(.,NUMP+1) (THAT IS, C ALL COEFFICIENTS ARE 0.0 EXCEPT THE FIRST DENOMINATOR C COEFFICIENT, WHICH IS 1.0). THIS APPROXIMATION IS LIKELY C TO BE OF LITTLE VALUE. C C IFLAG=6666 MEANS EVEN THE APPROXIMATION P/Q=0.0/PWR(.,NUMP+1) C WAS REJECTED BY SUBROUTINE SPQERD BECAUSE ITS DENOMINATOR C WAS VERY SMALL IN ABSOLUTE VALUE OR NEGATIVE AT SOME GRID C POINT. NO APPROXIMATION WAS RETURNED, ALL THE COEFFICIENTS C OF P AND Q WERE SET TO 0.0, AND ITER WAS SET TO -1. C C IN CASE IFLAG IS NEGATIVE (INDICATING CAUTION), IT MAY BE C ADVISABLE TO EXAMINE THE ERRORS AT THE GRID POINTS, AS C NORMALLY THE NUMBER OF POINTS WHERE THE ERROR NORM IS C ACHIEVED OR A USER SUPPLIED CONSTRAINT IS SATISFIED WITH C EQUALITY, PLUS THE NUMBER OF DENOMINATOR COEFFICIENTS C WITH ABSOLUTE VALUE 1.0, WILL BE AT LEAST NUMP+NUMQ+1. C THE SPECIFIC POSSIBILITIES FOR IFLAG NEGATIVE, MORE THAN C ONE OF WHICH CAN OCCUR, ARE AS FOLLOWS. C C IF THE ONES DIGIT OF IFLAG IS 1, THE INITIAL APPROXIMATION C WAS NOT IMPROVED, AND WAS RETURNED TO THE USER. THIS C FREQUENTLY HAPPENS NATURALLY (AND WITH NO ILL C CONSEQUENCES) WHEN ONE IS DOING GENERALIZED POLYNOMIAL C APPROXIMATION, THAT IS, WHEN ONE HAS SET NUMQ=1 AND C PWR(I,NUMP+1)=1.0 FOR I=1,...,NUMGR. C C IF THE ONES DIGIT OF IFLAG IS 2, THE PROGRAM WAS TERMINATED C BECAUSE ITLIM ITERATIONS WERE PERFORMED, WHERE CURRENTLY C ITLIM=100. ONE MAY WISH TO RESTART THE PROGRAM, C INITIALIZING WITH THE CURRENT APPROXIMATION, TO SEE IF THIS C APPROXIMATION CAN BE IMPROVED. C C IF THE TENS DIGIT OF IFLAG IS 1, THEN DURING THE COMPUTATION C OF THE RETURNED APPROXIMATION IT WAS NECESSARY TO ADJUST C THE TOLERANCES IN SUBROUTINE SLNPRO TO AVOID FAILURE. C THUS THE DANGER THAT THE RETURNED APPROXIMATION HAS BEEN C SERIOUSLY AFFECTED BY ROUND OFF ERROR IS GREATER. C C IF THE HUNDREDS DIGIT OF IFLAG IS 1, A USER SUPPLIED C CONSTRAINT WAS VIOLATED BY MORE THAN EPRES, WHERE C EPRES=10000.0*B**(-T). (B IS THE BASE FOR FLOATING POINT C NUMBERS FOR THE MACHINE AND T IS THE NUMBER OF BASE B C DIGITS IN THE MANTISSA.) C C THE PATTERN FOR DIMENSION DECLARATIONS WHICH WILL SET ASIDE C SUFFICIENT SPACE TO SOLVE ANY PROBLEM WITH UP TO MFN BASIS C FUNCTIONS (THUS NUMP+NUMQ .LE. MFN) AND MGR GRID POINTS (THUS C NUMGR .LE. MGR) IS C C FTBLE(MGR), PWR(MGR,MFN), INDLW(MGR), ALTBL(MGR), INDUP(MGR), C UPTBL(MGR), INDF(MGR), WTBLE(MGR), X(MFN+1), ERROR(MGR+1), C T(2*MGR), XTEMP(MFN+1), FTBWT(MGR), IYRCT(MAXCN), IYSAV(MAXCN), C V(MAXCN+1,MFN+2), IYCCT(MFN+1), IXRCT(MAXCN), Y(MAXCN), C C WHERE MAXCN=5*MGR+2*MFN-2. NOTE THAT THE LAST 8 ARRAYS LISTED C ARE INTERNAL WORK ARRAYS AND NEED NOT BE DIMENSIONED IN THE C USERS DRIVER PROGRAM. WE NOTE THAT FOR EACH OF OPTIONS 2 OR 3 C WHICH IS NOT USED, MAXCN COULD BE REDUCED BY MGR. C C THE METHOD USED BY THIS PROGRAM IS THE DIFFERENTIAL CORRECTION C ALGORITHM, WHICH REQUIRES AT EACH ITERATION THE SOLUTION OF THE C LINEAR PROGRAMMING PROBLEM C MINIMIZE MAX(ABS(F*Q-WT*P)-DELTK*Q)/QK, C SUBJECT TO ABS(Q(J)) .LE. 1.0 FOR J=1,...,NUMQ, C WHERE QK AND DELTK ARE THE DENOMINATOR AND THE UNIFORM ERROR C NORM FOR THE PREVIOUS APPROXIMATION, WT IS THE WEIGHT FUNCTION C (NOTE THAT WT WILL BE 1.0 FOR UNWEIGHTED APPROXIMATION, AND F C WILL HAVE BEEN MULTIPLIED BY WT FOR WEIGHTED ERROR CURVE C APPROXIMATION), AND Q(J) IS THE JTH COEFFICIENT OF THE NEW Q. C HERE MAX MEANS THE MAXIMUM COMPUTED OVER THE GRID. C IN THEORY, THE ERROR NORMS OF THE APPROXIMATIONS PRODUCED BY C THE ALGORITHM WILL CONVERGE TO THE INFIMUM OF ALL POSSIBLE C ERROR NORMS. C C IN ORDER TO SPEED UP CONVERGENCE FAR FROM THE SOLUTION, FOR C ITERATIONS BEYOND THE INITIALIZATION AND FIRST MAIN ITERATION C THE DELTK IN THE ABOVE EXPRESSION IS REPLACED BY 0.5*DELTK C UNTIL THE PROGRAM FAILS TO PRODUCE A BETTER APPROXIMATION, C AFTER WHICH THE PROGRAM REVERTS TO ORDINARY DIFFERENTIAL C CORRECTION FROM THEN ON. IN THEORY THIS PROCESS, KNOWN AS C PARAMETRIC BISECTION, WILL RESULT IN MORE THAN HALVING THE C ERROR NORM AT EACH ITERATION IN WHICH THE ERROR NORM IS C MORE THAN TWICE THE OPTIMAL ERROR NORM. C C IN THE EVENT OF FAILURE OF THE LINEAR PROGRAMMING SUBROUTINE C SLNPRO TO PRODUCE AN APPROXIMATION P/Q BETTER THAN THE PREVIOUS C ONE (IN THE ERROR NORM SENSE), THE PROGRAM TRIES TWICE MORE, ONCE C WITH LOOSENED AND ONCE WITH TIGHTENED TOLERANCES IN SLNPRO. AN C APPROXIMATION PRODUCED BY SLNPRO IS TESTED IN SUBROUTINE SPQERD TO C SEE IF ITS DENOMINATOR IS NEGATIVE OR VERY SMALL IN ABSOLUTE VALUE C (THAT IS, LESS THAN 100.0*B**-(T)) AT ANY GRID POINT. IF SO, THE C APPROXIMATION IS REJECTED (EXCEPT IN THE INITIALIZATION, WHERE C THE PROGRAM RESETS THE DENOMINATOR AT THE BAD POINTS AND C ATTEMPTS TO CONTINUE, AND IF A BETTER APPROXIMATION IS NOT C PRODUCED, THEN THIS APPROXIMATION IS REJECTED). C TO SAVE TIME, THE TWO EXTRA LINEAR PROGRAMMING ATTEMPTS ARE C NOT DONE IF THE FAILURE OCCURS DURING THE PARAMETRIC BISECTION C PHASE SINCE IN THAT PHASE THE PARAMETRIC BISECTION ATTEMPT WAS C THE PROBABLE CAUSE OF THE FAILURE, AND THE PROGRAM WILL HAVE C ANOTHER CHANCE USING ORDINARY DIFFERENTIAL CORRECTION. C C THE INITIAL APPROXIMATION FOR THE ALGORITHM IS OBTAINED BY ONE C OF THREE METHODS, WITH THE USER BEING ABLE TO SELECT THE FIRST C OR (BY DEFAULT) THE SECOND. IF A METHOD FAILS, THE PROGRAM C WILL TRY THE NEXT METHOD IN THE SEQUENCE. FAILURE OF A METHOD C MEANS EITHER NO RATIONAL APPROXIMATION P/Q COULD BE PRODUCED, C OR ONE WAS PRODUCED WHICH HAD A DENOMINATOR WHICH WAS VERY C SMALL IN ABSOLUTE VALUE OR NEGATIVE AT SOME POINT (AND THE MAIN C ALGORITHM WAS UNABLE TO CORRECT THIS PROBLEM). THE THREE C METHODS ARE C C METHOD 1 THE USER SUPPLIES AN INITIAL APPROXIMATION (SEE C OPTION 6 ABOVE). C C METHOD 2 ONE ITERATION OF LOEBS ALGORITHM IS PERFORMED, THAT C IS, SUBROUTINE SLNPRO IS USED TO SOLVE THE PROBLEM C MINIMIZE MAX(ABS(F*Q-WT*P)) C SUBJECT TO Q .GE. MAX(0.00001,10000.0*B**(-T)) C EVERYWHERE ON THE GRID, AND Q(1)=1.0 C IF WE ARE PLACING A LOWER BOUND ON THE DENOMINATOR C IN THE MAIN ITERATIONS, THE LOWER BOUND HERE IS C INCREASED (IF NECESSARY) TO MATCH THAT BOUND, AND C THE ABSOLUTE VALUES OF THE FREE DENOMINATOR C COEFFICIENTS ARE BOUNDED BY 1.0 C IF SLNPRO FAILS TO PRODUCE AN APPROXIMATION P/Q WE C TRY WITH BOTH LOOSENED AND TIGHTENED TOLERANCES, C AND ALSO WITH Q(1)=-1.0 C C METHOD 3 THE INITIAL APPROXIMATION IS TAKEN TO BE P/Q= C 0.0/PWR(.,NUMP+1). C C THE ALGORITHM IS TERMINATED WHEN A BETTER APPROXIMATION CANNOT C BE PRODUCED OR ITLIM ITERATIONS HAVE BEEN PERFORMED, AND THE C BEST APPROXIMATON FOUND SO FAR IS THEN RETURNED TO THE USER. C IF IT WAS NECESSARY TO ADJUST THE TOLERANCES IN SLNPRO TO C OBTAIN AN APPROXIMAITON P/Q, THEN VIOLATION OF A USER SUPPLIED C CONSTRAINT BY MORE THAN 10000.0*B**(-T) WILL CAUSE THAT C APPROXIMATION TO BE REJECTED SINCE THE ADJUSTED TOLERANCES C MAY HAVE ALLOWED REDUCTION OF THE UNIFORM ERROR NORM AT THE C EXPENSE OF A SERIOUS VIOLATION OF THE USER SUPPLIED CONSTRAINTS. C C***REFERENCES CHENEY, E. W. AND LOEB, H. L., TWO NEW ALGORITHMS C FOR RATIONAL APPROXIMATION, NUMER. MATH. 4, C PP. 124-127, 1962. C BARRODALE, I., POWELL, M. J. D., AND ROBERTS, C F. D. K., THE DIFFERENTIAL CORRECTION ALGORITHM C FOR RATIONAL L-INFINITY APPROXIMATION, SIAM J. C NUMER. ANAL. 9, PP. 493-504, 1972. C KAUFMAN, EDWIN H. JR., LEEMING, D. J., AND TAYLOR, C G. D., UNIFORM RATIONAL APPROXIMATION BY C DIFFERENTIAL CORRECTION AND REMES-DIFFERENTIAL C CORRECTION, INT. J. FOR NUMER. METH. IN ENGRG. C 17, PP. 1273-1278, 1981. C***ROUTINES CALLED SINTPQ,SLNPRO,SPQERD,SSETV1,SSETV2 C***END PROLOGUE SDFCOR DIMENSION FTBLE(101),PWR(101,15),V(534,17),WTBLE(101), * UPTBL(101),ALTBL(101),INDUP(101),INDLW(101), * ERROR(102),IYRCT(533),IYSAV(533),XTEMP(16),X(16), * INDF(101),FTBWT(101) C C THE BASIC PATTERN OF STATEMENT NUMBERS IN THIS PROGRAM IS THAT C STATEMENTS NUMBERED 1, 2,... ARE FROM A VERSION OF THE PROGRAM C PUBLISHED PRIOR TO THE THIRD REFERENCE CITED ABOVE, STATEMENT C NUMBERS 1001, 1002,... ARE FROM THE VERSION OF THAT PAPER, AND C STATEMENT NUMBERS 10000 OR HIGHER REPRESENT LATER MODIFICATIONS. C C SET MACHINE DEPENDENT PARAMETERS FOR SUBROUTINE SDFCOR. C SET SPCMN=B**(-T), WHERE B IS THE BASE AND T IS THE NUMBER C OF BASE B DIGITS IN THE MANTISSA. SPCMN IS THE MINIMUM C RELATIVE SPACING ABS((X1-X2)/X2) BETWEEN TWO SUCCESSIVE C FLOATING POINT NUMBERS, SO IT IS THE SPACING BETWEEN TWO C SUCCESSIVE FLOATING POINT NUMBERS IN THE CLOSED INTERVAL C (0.1,1.0). WE ALSO HAVE SPCMN=10.0**(-T*(LOG10)(B))= C 10.0**(-TNMAN), WHERE TNMAN IS THE BASE 10 EQUIVALENT OF C THE LENGTH OF THE MANTISSA. C C***FIRST EXECUTABLE STATEMENT SDFCOR SPCMN=R1MACH(3) C SET EPRES=10000.0*SPCMN. THUS EPRES=10.0**(-(TNMAN-4)). C AN APPROXIMATION WILL BE REJECTED IF A USER SUPPLIED C CONSTRAINT IS VIOLATED BY MORE THAN EPRES AND IT WAS NECESSARY C TO ADJUST THE TOLERANCES IN SUBROUTINE SLNPRO WHILE COMPUTING C THE APPROXIMATION. EPRES=10000.0*SPCMN C END OF SETTING MACHINE DEPENDENT PARAMETERS FOR SDFCOR. C ITLIM=100 ITER=0 NMPP1=NUMP+1 NMPPQ=NUMP+NUMQ N=NMPPQ+1 NP1=N+1 IADJS=0 INDIC=0 IRET1=0 ILBFL=0 ILST1=0 QMIN=0.0 IQMIN=0 INITL=0 ICHK=0 IFLAG=0 IPBIT=1 C C IF F IS NOT TO BE IGNORED AT ANY OF THE GRID POINTS, C SET INDF IDENTICALLY EQUAL TO 1. IF(IOPT-(IOPT/10000)*10000-1000)1001,1003,1003 1001 DO 1002 I=1,NUMGR INDF(I)=1 1002 CONTINUE C C IF WE ARE DEALING WITH UNWEIGHTED APPROXIMATION SET THE C WEIGHT FUNCTION IDENTICALLY EQUAL TO 1. 1003 IF(IOPT-(IOPT/100000)*100000-10000)1004,1006,1006 1004 DO 1005 I=1,NUMGR WTBLE(I)=1.0 1005 CONTINUE C COPY FTBLE INTO FTBWT FOR USE BY THE REST OF THIS PROGRAM C UNLESS WEIGHTED ERROR CURVE APPROXIMATION IS DESIRED. IN C THAT CASE SET FTBWT(I)=WTBLE(I)*FTBLE(I) AT THOSE POINTS C WHERE F IS TO BE APPROXIMATED, SO WE WILL HAVE C FTBWT(I)-WTBLE(I)*P/Q=WTBLE(I)*(FTBLE(I)-P/Q). THE VECTOR C FTBWT IS INTRODUCED SO THAT FTBLE NEED NOT BE CHANGED BY C THE PROGRAM. 1006 IF(IOPT-(IOPT/100000)*100000-20000)10000,1007,1007 10000 DO 10020 I=1,NUMGR IF(INDF(I))10020,10020,10010 10010 FTBWT(I)=FTBLE(I) 10020 CONTINUE GO TO 1010 1007 DO 1009 I=1,NUMGR IF(INDF(I))1009,1009,1008 1008 FTBWT(I)=WTBLE(I)*FTBLE(I) 1009 CONTINUE C C IF THE 100000 DIGIT OF IOPT IS 1, THEN THE C INITIALIZATION BY LOEBS ALGORITHM WILL BE SKIPPED, AND C THE COEFFICIENTS PLACED IN X BY THE USER WILL BE TAKEN FOR C THE INITIAL APPROXIMATION. 1010 IF(IOPT/100000)1012,1012,10030 C SET INITL=-1 TO INDICATE WE ARE USING A USER SUPPLIED INITIAL C APPROXIMATION. 10030 INITL=-1 GO TO 1024 C C DURING THE LOEB INITIALIZATION X(N) WILL BE USED TO STORE C THE VALUE ASSIGNED TO Q(1) (EITHER 1.0 OR -1.0). 1012 X(N)=1.0 C WE TRY TO FIND P,Q SUCH THAT Q(1)=1.0, Q IS POSITIVE (AND C NOT TOO SMALL) AT EVERY GRID POINT, P/Q SATISFIES THE C RESTRICTED RANGE CONSTRAINTS, AND MAX ABS(F*Q-WT*P) IS C MINIMIZED AT THOSE GRID POINTS WHERE F IS TO BE C APPROXIMATED. 1013 CALL SSETV1(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,INDLW,ALTBL, *INDUP,UPTBL,INDF,WTBLE,X,V,M) C SET IYRCT(1)=-1 AS A SIGNAL TO DO THE INITIAL EXCHANGES C IN SLNPRO STRICTLY ACCORDING TO A PIVOTING STRATEGY. IYRCT(1)=-1 C IF SLNPRO FAILED DURING THE INITIALIZATION WITH NORMAL TOLERANCES, C AS INDICATED BY IADJS=1, WE SET M=-M AS A SIGNAL TO LOOSEN THE C TOLERANCES IN SLNPRO (I. E. INCREASE REA AND REA1) AND TRY AGAIN. C SLNPRO WILL RESET M TO ITS ORIGINAL VALUE. ANY LOSS OF ACCURACY C DUE TO THIS LOOSENING WILL PROBABLY BE CORRECTED IN LATER ITERATIONS. M=M*(1-4*IADJS+2*IADJS**2) C IF SLNPRO FAILED DURING THE INITIALIZATION WITH BOTH NORMAL C TOLERANCES AND LOOSENED TOLERANCES, AS INDICATED BY IADJS=2, C WE SET NMPPQ=-NMPPQ AS A SIGNAL TO TIGHTEN THE TOLERANCES IN SLNPRO C (I. E. DECREASE REA AND REA1) AND TRY ONCE AGAIN. SLNPRO WILL C RESET NMPPQ TO ITS ORIGINAL VALUE. NMPPQ=NMPPQ*(1+IADJS-IADJS**2) CALL SLNPRO(V,M,NMPPQ,IYRCT,X,INDIC) IRET1=INDIC C IF THIS WAS THE FIRST SLNPRO CALL IN THE LOEB INITIALIZATION, C SO WE HAD X(N)=1.0 AND THE TOLERANCES IN SLNPRO WERE NOT C ADJUSTED BY SDFCOR, FOR POSSIBLE OUTPUT PURPOSES WE SAVE C INDIC IN ILBFL. IF(X(N))10060,10060,10040 10040 IF(IADJS)10050,10050,10060 10050 ILBFL=INDIC C IF INDIC IS POSITIVE SOMETHING WENT WRONG IN SLNPRO. C OTHERWISE WE PUT THE COEFFICIENTS OF Q INTO THEIR PROPER C PLACES IN X AND NORMALIZE P AND Q. 10060 IF(INDIC)1023,1023,1016 C IF WE CANNOT FIND P/Q (WITH Q POSITIVE) SATISFYING THE C CONSTRAINTS WITH Q(1)=1.0, WE TRY WITH Q(1)=-1.0. 1016 IF(X(N))1018,1018,1017 1017 X(N)=-1.0 GO TO 1013 1018 IF(IADJS-1)1019,1019,1020 C BEFORE CONCEDING DEFEAT WE WILL TRY AGAIN WITH CHANGED C TOLERANCES. 1019 IADJS=IADJS+1 GO TO 1012 C C HERE THE LOEB INITIALIZATION HAS FAILED, AND WE TAKE OUR INITIAL C APPROXIMATION TO BE P/Q=0.0/PWR(.,NUMP+1). WE ALSO SET INITL=1 C AS A WARNING THAT THIS HAS BEEN DONE. 1020 INITL=1 DO 10070 J=1,NMPPQ X(J)=0.0 10070 CONTINUE X(NMPP1)=1.0 C RESET PARAMETERS. IADJS=0 INDIC=0 IRET1=0 QMIN=0.0 IQMIN=0 GO TO 1024 C C HERE INDIC WAS NONPOSITIVE, SO SLNPRO PRODUCED AN APPROXIMATION. C PUT THE COEFFICIENTS OF P AND Q IN THE PROPER PLACES IN X AND C NORMALIZE Q. 1023 CALL SINTPQ(NUMP,NUMQ,X) C C COMPUTE THE VALUES OF P AND Q AT THE GRID POINTS, THE C ERRORS AT THE GRID POINTS, AND THE UNIFORM ERROR NORM. 1024 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X, *V,QMIN,IQMIN) DELTK=V(2*NUMGR+1,NP1) DLOLD=DELTK IADJS=0 C PUT THE ERRORS AND ERROR NORM INTO ERROR. DO 6 I=1,NUMGR ERROR(I)=V(I,N) 6 CONTINUE ERROR(NUMGR+1)=DELTK C IF THE DENOMINATOR WAS VERY SMALL IN ABSOLUTE VALUE OR C NEGATIVE AT SOME POINT (WHICH IS INDICATED BY C V(2*NUMGR+1,NP1) BEING NEGATIVE), WE RESET C V(2*NUMGR+1,NP1) TO THE ERROR NORM COMPUTED WITH THE BAD C POINTS IGNORED FOR USE IN SSETV2, AND ATTEMPT TO CONTINUE. C NOTE THAT DELTK, DLOLD, AND ERROR(NUMGR+1) WILL REMAIN C NEGATIVE AS A WARNING UNLESS AN APPROXIMATION WITH C SIGNIFICANTLY POSITIVE DENOMINATOR IS FOUND BELOW. C THIS SITUATION SOMETIMES OCCURS WHEN AN APPROXIMATION C COMPUTED ON A COARSE GRID IS USED AS A STARTING C APPROXIMATION ON A FINE GRID. IF(DELTK)1025,1026,1026 1025 V(2*NUMGR+1,NP1)=-V(2*NUMGR+1,NP1)-2.0 C END OF INITIALIZATION PHASE. C C C SET UP AND SOLVE THE LINEAR PROGRAMMING PROBLEM TO GET THE C NEXT APPROXIMATION. 1026 CALL SSETV2(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,INDLW, *ALTBL,INDUP,UPTBL,INDF,WTBLE,V,M) C IF WE ARE IN THE FIRST STEP OF THE MAIN ALGORITHM, SET C IYRCT(1)=-1 AS A SIGNAL TO DO THE INITIAL EXCHANGES IN C SLNPRO STRICTLY ACCORDING TO A PIVOTING STRATEGY. C HENCEFORTH WHEN SLNPRO IS CALLED THE OLD VERTEX WILL BE C TAKEN AS THE STARTING POINT. IF(ITER)1027,1027,1028 1027 IYRCT(1)=-1 C SAVE IYRCT IN IYSAV. 1028 DO 1029 I=1,M IYSAV(I)=IYRCT(I) 1029 CONTINUE C SET M=-M AS A SIGNAL TO LOOSEN THE TOLERANCES IN SLNPRO (I.E. C INCREASE REA AND REA1) IF IADJS=1, AND LEAVE M UNCHANGED IF C IADJS=0 OR IADJS=2. SLNPRO WILL AUTOMATICALLY RESET M IF IT C IS CHANGED HERE. M=M*(1-4*IADJS+2*IADJS**2) C SET N=-N AS A SIGNAL TO TIGHTEN THE TOLERANCES IN SLNPRO (I.E. C DECREASE REA AND REA1) IF IADJS=2, AND LEAVE N UNCHANGED IF C IADJS=0 OR IADJS=1. SLNPRO WILL AUTOMATICALLY RESET N IF IT C IS CHANGED HERE. N=N*(1+IADJS-IADJS**2) CALL SLNPRO(V,M,N,IYRCT,XTEMP,INDIC) C IF THIS WAS THE FIRST SLNPRO CALL IN THIS ITERATION, SO THE C TOLERANCES IN SLNPRO WERE NOT ADJUSTED BY SDFCOR, FOR POSSIBLE C OUTPUT PURPOSES WE SAVE INDIC IN ILST1. IF(IADJS)10080,10080,10090 10080 ILST1=INDIC C IF INDIC IS POSITIVE SOMETHING WENT WRONG IN SLNPRO, SO C WE RETURN WITH THE X AND ERROR COMPUTED AT THE PREVIOUS C ITERATION AFTER TRYING AGAIN WITH BOTH TIGHTER AND LOOSER C TOLERANCES, IF THIS HAS NOT ALREADY BEEN TRIED. 10090 IF(INDIC)1030,1030,1045 C C HERE INDIC WAS NONPOSITIVE, SO SLNPRO PRODUCED AN APPROXIMATION. C COMPUTE THE VALUES OF P AND Q AT THE GRID POINTS, THE C ERRORS AT THE GRID POINTS, AND THE UNIFORM ERROR NORM. 1030 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,XTEMP, *V,QMIN,IQMIN) DELTK=V(2*NUMGR+1,NP1) C REJECT THIS APPROXIMATION IF THE DENOMINATOR WAS VERY C SMALL IN ABSOLUTE VALUE OR NEGATIVE AT SOME POINT. C THIS IS INDICATED BY DELTK BEING NEGATIVE. IF(DELTK)1045,1031,1031 C IF INDIC IS NEGATIVE HERE IT WAS NECESSARY TO ADJUST THE C TOLERANCES IN SLNPRO TO GET AN ACCEPTABLE APPROXIMATION, SO THE C APPROXIMATION COMPUTED IS MORE LIKELY TO HAVE BEEN C AFFECTED BY ROUND OFF ERROR. THIS PRESENTS NO DANGER IF NEITHER C RESTRICTED RANGE NOR DENOMINATOR CONSTRAINTS WERE USED, SINCE AN C APPROXIMATION WITH ERROR NORM WORSE THAN THE PREVIOUS ONE C WILL BE REJECTED BELOW ANYWAY, BUT IF THERE ARE RESTRICTED C RANGE OR DENOMINATOR CONSTRAINTS THE ERROR NORM COULD HAVE C BEEN REDUCED AT THE EXPENSE OF VIOLATING THE CONSTRAINTS. C THUS WE REJECT THE APPROXIMATION IF ANY CONSTRAINT IS VIOLATED C BY MORE THAN EPRES. 1031 IF(INDIC)1032,1040,1040 C C CHECK UPPER CONSTRAINTS, IF THERE ARE ANY. 1032 IF(IOPT-(IOPT/1000)*1000-100)1036,1033,1033 1033 DO 1035 I=1,NUMGR IF(INDUP(I))1035,1035,1034 1034 IF(V(2*I-1,NP1)/V(2*I,NP1)-UPTBL(I)-EPRES)1035,1035,10140 1035 CONTINUE C CHECK LOWER CONSTRAINTS, IF THERE ARE ANY. 1036 IF(IOPT-(IOPT/100)*100-10)10100,1037,1037 1037 DO 1039 I=1,NUMGR IF(INDLW(I))1039,1039,1038 1038 IF(V(2*I-1,NP1)/V(2*I,NP1)-ALTBL(I)+EPRES)10140,1039,1039 1039 CONTINUE C CHECK DENOMINATOR CONSTRAINTS, IF THERE ARE ANY. 10100 IF(IOPT-(IOPT/10)*10)10130,10130,10110 10110 DO 10120 I=1,NUMGR IF(V(2*I,NP1)-EPSQ+EPRES)10140,10120,10120 10120 CONTINUE C HERE NO USER SUPPLIED CONSTRAINT WAS VIOLATED BY MORE THEN EPRES. C IF WE ARE IN THE FINAL CHECK WE RETURN. 10130 IF(ICHK)1040,1040,15 C HERE AT LEAST ONE USER SUPPLIED CONSTRAINT WAS VIOLATED BY MORE C THAN EPRES. IF WE ARE IN THE FINAL CHECK WE ADJUST IFLAG AND C RETURN. 10140 IF(ICHK)1045,1045,10150 10150 IFLAG=IFLAG-100 RETURN C C IF DLOLD IS NEGATIVE, THEN WE WILL HAVE ITER=0, AND THE C INITIAL APPROXIMATION WILL HAVE A DENOMINATOR WHICH IS C VERY SMALL IN ABSOLUTE VALUE OR NEGATIVE AT SOME POINT. C THUS THE NEW APPROXIMATION, WHICH HAS A SIGNIFICANTLY C POSITIVE DENOMINATOR EVERYWHERE ON THE GRID, WILL BE C PREFERRED. 1040 IF(DLOLD)10,1041,1041 C IF THE ERROR NORM HAS NOT DECREASED, RETURN WITH THE X C AND ERROR COMPUTED AT THE PREVIOUS ITERATION AFTER TRYING C AGAIN WITH CHANGED TOLERANCES, IF BOTH TIGHTER AND LOOSER C TOLERANCES HAVE NOT ALREADY BEEN TRIED. 1041 IF(DELTK-DLOLD)1042,1045,1045 C IN RARE INSTANCES AN APPROXIMATION WITH ALL COEFFICIENTS C SMALL BUT WITH A GOOD COMPUTED ERROR NORM WILL BE C PRODUCED. TO AVOID THIS WE REJECT AN APPROXIMATION IF C MAX(ABS(Q(J))).LE.0.1 1042 DO 1043 J=NMPP1,NMPPQ IF(ABS(XTEMP(J))-0.1)1043,1043,10 1043 CONTINUE GO TO 1045 C C HERE THE PRESENT APPROXIMATION WILL BE KEPT, AND WE PREPARE C TO COMPUTE ANOTHER ONE. C COPY XTEMP INTO X. 10 DO 11 I=1,NMPPQ X(I)=XTEMP(I) 11 CONTINUE C PUT THE ERRORS AND ERROR NORM INTO ERROR. DO 12 I=1,NUMGR ERROR(I)=V(I,N) 12 CONTINUE ERROR(NUMGR+1)=DELTK ITER=ITER+1 IF(ITER-ITLIM)14,10160,10160 14 DLOLD=DELTK IRET1=INDIC IADJS=0 C IF IPBIT=1 WE ARE STILL IN THE PARAMETRIC BISECTION PHASE, C AND WE CUT V(2*NUMGR+1,N+1) IN HALF IN ORDER TO TRY FOR AT C LEAST A HALVING OF THE ERROR NORM IN THE NEXT STEP. IF(IPBIT)1026,1026,10155 10155 V(2*NUMGR+1,NP1)=0.5*V(2*NUMGR+1,NP1) GO TO 1026 C HERE ITLIM ITERATIONS HAVE BEEN PERFORMED AND WE WILL RETURN C AFTER SETTING IFLAG AND INDLP (SEE BELOW). 10160 INDLP=5-INDIC IF(INDIC)10170,10180,10180 10170 IFLAG=-12 GO TO 10300 10180 IFLAG=-2 GO TO 10300 C C HERE SLNPRO FAILED TO PRODUCED AN APPROXIMATION WHICH IS BETTER C THAN THE PREVIOUS APPROXIMATION. IF WE HAVE NOT ALREADY C DONE SO, WE TRY AGAIN WITH TIGHTER TOLERANCES. IF THIS C HAS BEEN TRIED, WE TRY ONCE MORE WITH LOOSER TOLERANCES. C THE ONLY EXCEPTION TO THIS IS THAT IF WE ARE STILL IN THE C PARAMETRIC BISECTION PHASE AND ITER .GT. 0, TO SAVE TIME WE C DO NOT TRY WITH CHANGED TOLERANCES, SINCE THE FAILURE WAS C PROBABLY DUE TO THE PARAMETRIC BISECTION. INSTEAD WE TRY C WITH ORDINARY DIFFERENTIAL CORRECTION, AND USE ORDINARY C DIFFERENTIAL CORRECTION FROM NOW ON. 1045 IF(ITER*IPBIT)10185,10185,10183 10183 IPBIT=0 GO TO 10187 10185 IF(IADJS-1)1046,10190,1046 1046 IADJS=2-IADJS/2 C RESTORE THE PREVIOUS IYRCT FOR THE NEXT ATTEMPT. 10187 DO 1047 I=1,M IYRCT(I)=IYSAV(I) 1047 CONTINUE C CALL SPQERD WITH THE PREVIOUS APPROXIMATION (WHICH IS IN X) C TO RECOMUTE THE PREVIOUS VALUES OF P AND Q AND PUT THEM C IN V FOR USE BY SSETV2. CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X, *V,QMIN,IQMIN) GO TO 1026 C C HERE THE APPROXIMATION HAS NOT BEEN IMPROVED EVEN BY ADJUSTING C TOLERANCES TWICE IN SLNPRO. IF IT IS AN INITIAL APPROXIMATION C WHICH WAS FOUND BY SPQERD TO HAVE A DENOMINATOR WHICH IS VERY C SMALL IN ABSOLUTE VALUE OR NEGATIVE AT SOME POINT (INDCIATED BY C DLOLD NEGATIVE), AND IF IT IS A USER SUPPLIED APPROXIMATION C (INDICATED BY INITL=-1), THEN WE ATTEMPT A LOEB REINITIALIZATION. C IF IT WAS A LOEB APPROXIMATION (INDICATED BY INITL=0) WITH C UNACCEPTABLE DENOMINATOR, THEN WE TRY P/Q=0.0/PWR(.,NUMP+1). 10190 IF(DLOLD)10200,1048,1048 10200 IF(INITL)10210,1020,1048 C SET UP TO TRY THE LOEB INITIALIZATION. 10210 INITL=0 IADJS=0 INDIC=0 IRET1=0 QMIN=0.0 IQMIN=0 GO TO 1012 C C HERE WE ARE FINISHED COMPUTING APPROXIMATIONS, SO WE SET INDLP, C IFLAG, QMIN, AND IQMIN, AND WE RETURN. C WE CARRY INFORMATION BACK TO THE CALLING PROGRAM ABOUT TWO C SLNPRO CALLS BY SETTING THE TWO DIGIT VARIABLE INDLP WITH C TENS DIGIT 5-ILAST AND ONES DIGIT 5-IRET, WHERE IN MOST C CASES ILAST AND IRET GIVE INFORMATION ON THE LAST (REJECTED) C SLNPRO ATTEMPT AND THE RETURNED APPROXIMATION RESPECTIVELY. C TO BE EXACT, IF THE PROGRAM WAS TERMINATED BECAUSE ITLIM C ITERATIONS WERE PERFORMED SUCCESSFULLY, ILAST WILL BE 5 (AND C SO 5-ILAST WILL BE 0). OTHERWISE, ILAST WILL BE THE VALUE OF C INDIC PRODUCED BY SLNPRO AT THE LAST (FAILED) MAIN ITERATION, C WITH THE TOLERANCES IN SLNPRO NOT HAVING BEEN ADJUSTED BY C SDFCOR YET. C IF THE PROGRAM FAILED IN THE INITIALIZATION (SO IFLAG WAS C POSITIVE), IRET WILL BE THE VALUE OF INDIC PRODUCED BY C SLNPRO DURING THE (FAILED) LOEB INITIALIZATION, WITH THE C TOLERANCES IN SLNPRO NOT HAVING BEEN ADJUSTED BY SDFCOR C YET. IF A USER SUPPLIED INITIAL APPROXIMATION WAS RETURNED C TO THE USER, IRET WILL BE O. OTHERWISE, IRET WILL BE THE C VALUE OF INDIC PRODUCED BY SLNPRO DURING THE COMPUTATION C OF THE APPROXIMATION WHICH WAS ACTUALLY RETURNED TO THE USER. C THUS IF THE APPROXIMATION RETURNED WAS COMPUTED BY SLNPRO C (WHICH WILL BE THE CASE IF EITHER THE ONES DIGIT OF IFLAG IS C 0 OR 2, OR IT IS 1 AND THE RETURNED APPROXIMATION IS NOT A USER C SUPPLIED INITIAL APPROXIMATION) AND THE ONES DIGIT OF INDLP IS 5, C THEN SLNPRO APPEARED TO WORK NORMALLY DURING THE COMPUTATION OF C THE RETURNED APPROXIMATION. 1048 INDLP=10*(5-ILST1)+5-IRET1 C C WE NOW SET UP IFLAG FOR THE INFORMATION OF THE USER. IF(ITER)10220,10220,10280 C HERE ITER=0, SO THE INITIAL APPROXIMATION WAS NOT IMPROVED. 10220 IF(INITL)10270,10270,10230 C HERE THE UNIMPROVED INITIAL APPROXIMATION WAS P/Q= C 0.0/PWR(.,NUMP+1). IF IT WAS FOUND BY SPQERD TO HAVE A C DENOMINATOR WHICH IS VERY SMALL IN ABSOLUTE VALUE OR C NEGATIVE AT SOME POINT, THEN SET IFLAG=6666, ITER=-1, AND C ZERO OUT THE ERRORS (WHICH WERE INVALID ANYWAY) AND C COEFFICIENTS AND RETURN. OTHERWISE SET IFLAG=3333 AND RETURN. C FIRST RESET INDLP SO THAT 5-(ONES DIGIT OF INDIC) IS THE VALUE C OF INDIC THE FIRST TIME SLNPRO WAS CALLED IN THE (FAILED) LOEB C INITIALIZATION. 10230 INDLP=10*(5-ILST1)+5-ILBFL IF(DLOLD)10240,10260,10260 10240 IFLAG=6666 ITER=-1 C NOTE THAT HERE X(I) ALREADY EQUALS 0.0 FOR I=1,...,NMPPQ, C I.NE.NUMP+1. X(NMPP1)=0.0 DO 10250 I=1,NUMGR ERROR(I)=0.0 10250 CONTINUE RETURN 10260 IFLAG=3333 RETURN C C HERE THE UNIMPROVED INITIAL APPROXIMATION WAS USER SUPPLIED OR C COMPUTED BY THE LOEB INITIALIZATION. SET THE ONES DIGIT OF C IFLAG TO INDICATE THE INITIAL APPROXIMATION WAS NOT IMPROVED. 10270 IFLAG=-1 C IF IRET1.LT.0 SET THE TENS DIGIT OF IFLAG TO INDICATE THAT C THE TOLERANCES IN SLNPRO HAD TO BE ADJUSTED DURING THE C COMPUTATION OF THE APPROXIMATION RETURNED. 10280 IF(IRET1)10290,10300,10300 10290 IFLAG=IFLAG-10 C NOW WE GO BACK AND CHECK THE USER SUPPLIED CONSTRAINTS (IF ANY). C SET ICHK=1 AS A SIGNAL AND CALL SPQERD TO RECOMPUTE THE VALUES C OF P AND Q. QMIN AND IQMIN ARE ALSO RECOMPUTED BY SPQERD. 10300 ICHK=1 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X, *V,QMIN,IQMIN) GO TO 1032 15 RETURN END SUBROUTINE STCOMP(NUMDM,NRWGR,NCLGR,ASWP1,ASWP2, *ANEP1,ANEP2,T) C***BEGIN PROLOGUE STCOMP C***DATE WRITTEN 720120 (YYMMDD) C***REVISION DATE 840301 (YYMMDD) C***CATEGORY NO. M2,Q C***KEYWORDS GRID CONSTRUCTION C***AUTHOR KAUFMAN, EDWIN H. JR., C DEPARTMENT OF MATHEMATICS C CENTRAL MICHIGAN UNIVERSITY C MOUNT PLEASANT, MICHIGAN 48859 C TAYLOR, GERALD D., C DEPARTMENT OF MATHEMATICS C COLORADO STATE UNIVERSITY C FORT COLLINS, COLORADO 80523 C***PURPOSE THIS SUBROUTINE, WHICH IS CALLED ONLY BY THE C USER, CAN BE USED TO COMPUTE THE COORDINATES C OF THE GRID POINTS IF THEY ARE AN EQUALLY C SPACED (IN EACH DIRECTION) SUBSET OF AN INTERVAL C OR RECTANGLE. C***DESCRIPTION C C NUMDM (INPUT) THIS IS THE NUMBER OF DIMENSIONS (1 OR 2). C C NRWGR (INPUT) THIS NEED NOT BE ASSIGNED A VALUE IN THE C 1 DIMENSIONAL CASE. IN THE 2 DIMENSIONAL CASE IT C IS THE NUMBER OF ROWS OF THE GRID. C C NCLGR (INPUT) THIS IS THE NUMBER OF GRID POINTS IN THE 1 C DIMENSIONAL CASE. IN THE 2 DIMENSIONAL CASE IT C IS THE NUMBER OF COLUMNS OF THE GRID. C C ASWP1 (INPUT) THIS IS THE LEFT END POINT IN THE 1 DIMENSIONAL C CASE. IN THE 2 DIMENSIONAL CASE IT IS THE FIRST C COORDINATE OF THE LOWER LEFT CORNER POINT. C C ASWP2 (INPUT) THIS IS THE RIGHT END POINT IN THE 1 DIMENSIONAL C CASE. IN THE 2 DIMENSIONAL CASE IT IS THE SECOND C COORDINATE OF THE LOWER LEFT CORNER POINT. C C ANEP1 (INPUT) THIS NEED NOT BE ASSIGNED A VALUE IN THE 1 C DIMENSIONAL CASE. IN THE 2 DIMENSIONAL CASE IT IS C THE FIRST COORDINATE OF THE UPPER RIGHT CORNER POINT. C C ANEP2 (INPUT) THIS NEED NOT BE ASSIGNED A VALUE IN THE 1 C DIMENSIONAL CASE. IN THE 2 DIMENSIONAL CASE IT IS C THE SECOND COORDINATE OF THE UPPER RIGHT CORNER POINT. C C T (OUTPUT) IN THE 1 DIMENSIONAL CASE THE COORDINATE C OF THE ITH GRID POINT WILL BE PLACED IN T(I). IN C THE 2 DIMENSIONAL CASE THE COORDINATES OF THE ITH C GRID POINT (READING ROW BY ROW, TOP TO BOTTOM) WILL C BE PLACED IN (T(2*I-1),T(2*I)). C***REFERENCES (NONE) C***ROUTINES CALLED (NONE) C***END PROLOGUE STCOMP DIMENSION T(501) C C***FIRST EXECUTABLE STATEMENT STCOMP IF(NUMDM-2)1,3,1 C C THIS IS THE 1 DIMENSIONAL CASE. 1 HSPAC=(ASWP2-ASWP1)/FLOAT(NCLGR-1) DO 2 I=1,NCLGR T(I)=ASWP1+HSPAC*FLOAT(I-1) 2 CONTINUE RETURN C C THIS IS THE 2 DIMENSIONAL CASE. 3 NCGDB=2*NCLGR IF(NCLGR-1)4,5,4 4 HSPAC=(ANEP1-ASWP1)/FLOAT(NCLGR-1) GO TO 6 5 HSPAC=0.0 6 IF(NRWGR-1)7,8,7 7 VSPAC=(ANEP2-ASWP2)/FLOAT(NRWGR-1) GO TO 9 8 VSPAC=0.0 9 DO 11 INDEX=1,NRWGR YCOR=ANEP2-VSPAC*FLOAT(INDEX-1) IND1=2+NCGDB*(INDEX-1) IND2=IND1+NCGDB-2 C FILL IN THE Y COORDINATES IN ROW INDEX. DO 10 J=IND1,IND2,2 T(J)=YCOR 10 CONTINUE 11 CONTINUE DO 13 INDEX=1,NCLGR XCOR=ASWP1+HSPAC*FLOAT(INDEX-1) IND1=2*INDEX-1 IND2=IND1+NCGDB*(NRWGR-1) C FILL IN THE X COORDINATES IN COLUMN INDEX. DO 12 J=IND1,IND2,NCGDB T(J)=XCOR 12 CONTINUE 13 CONTINUE RETURN END SUBROUTINE SSETV1(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,INDLW, *ALTBL,INDUP,UPTBL,INDF,WTBLE,X,V,M) C***BEGIN PROLOGUE SSETV1 C***REFER TO SDFCOR C***ROUTINES CALLED (NONE) C***PURPOSE THIS SUBROUTINE SETS UP THE LINEAR PROGRAMMING PROBLEM C FOR THE LOEB INITIALIZATION OF THE DIFFERENTIAL C CORRECTION ALGORITHM. C***END PROLOGUE SSETV1 DIMENSION FTBWT(101),PWR(101,15),V(534,17),WTBLE(101), * UPTBL(101),ALTBL(101),INDUP(101),INDLW(101),X(16), * INDF(101) C C THIS SUBROUTINE SETS UP V FOR THE LINEAR PROGRAMMING C PROBLEM OF MINIMIZING W WITH THE RESTRICTIONS C ABS(F*Q-WT*P).LE.W, WITH Q(1) FIXED AS X(NMPPQ+1). C C SET MACHINE DEPENDENT PARAMETERS FOR SUBROUTINE SSETV1. C SET SPCMN=B**(-T), WHERE B IS THE BASE AND T IS THE NUMBER C OF BASE B DIGITS IN THE MANTISSA. SPCMN IS THE MINIMUM C RELATIVE SPACING ABS((X1-X2)/X2) BETWEEN TWO SUCCESSIVE C FLOATING POINT NUMBERS, SO IT IS THE SPACING BETWEEN TWO C SUCCESSIVE FLOATING POINT NUMBERS IN THE CLOSED INTERVAL C (0.1,1.0). WE ALSO HAVE SPCMN=10.0**(-T*(LOG10)(B))= C 10.0**(-TNMAN), WHERE TNMAN IS THE BASE 10 EQUIVALENT OF C THE LENGTH OF THE MANTISSA. C C***FIRST EXECUTABLE STATEMENT SSETV1 SPCMN=R1MACH(3) C SET EPS=MAX(0.0001,10000.0*SPCMN). THUS EPS=10.0**(-EXEPS), C WHERE EXEPS=MIN(4,TNMAN-4). C WE WILL REQUIRE Q.GE.EPS AT EVERY GRID POINT. EPS=10000.0*SPCMN IF(EPS-0.0001)10000,10010,10010 10000 EPS=0.0001 C IF WE ARE PLACING A LOWER BOUND ON THE DENOMINATOR IN THE MAIN C ITERATIONS (INDICATED BY THE ONES DIGIT OF IOPT BEING 1) WE C INCREASE EPS, IF NECESSARY, TO MATCH THIS LOWER BOUND. 10010 IF(IOPT-(IOPT/10)*10)10040,10040,10020 10020 IF(EPSQ-EPS)10040,10040,10030 10030 EPS=EPSQ C END OF SETTING MACHINE DEPENDENT PARAMETERS FOR SSETV1. C 10040 NMPP1=NUMP+1 N1=NUMP+NUMQ NMPQ1=N1-1 N1P1=N1+1 C M WILL KEEP TRACK OF WHICH ROW OF V IS BEING SET. AT THE C CONCLUSION OF SSETV1 M WILL BE THE TOTAL NUMBER OF C CONSTRAINTS. M=0 C C SET THE UPPER CONSTRAINTS, IF THERE ARE ANY. THEY WILL BE C OF THE FORM P-Q*UP.LE.0.0, WITH THE CONSTANT Q(1)*UP TERM C TRANSPOSED. IF(IOPT-(IOPT/1000)*1000-100)1008,1001,1001 1001 DO 1007 I=1,NUMGR C INDUP(I)=1 IF THERE IS AN UPPER CONSTRAINT AT GRID POINT C I, AND INDUP(I)=0 OTHERWISE. IF(INDUP(I))1007,1007,1002 1002 M=M+1 DO 1003 J=1,NUMP V(M,J)=PWR(I,J) 1003 CONTINUE IF(NUMQ-1)1006,1006,1004 1004 DO 1005 J=NMPP1,NMPQ1 JP1=J+1 V(M,J)=-UPTBL(I)*PWR(I,JP1) 1005 CONTINUE 1006 V(M,N1)=0.0 V(M,N1P1)=X(N1P1)*UPTBL(I)*PWR(I,NMPP1) 1007 CONTINUE C C SET THE LOWER CONSTRAINTS, IF THERE ARE ANY. THEY WILL BE C OF THE FORM -P+Q*LW.LE.0.0, WITH THE CONSTANT Q(1)*LW TERM C TRANSPOSED. 1008 IF(IOPT-(IOPT/100)*100-10)1016,1009,1009 1009 DO 1015 I=1,NUMGR C INDLW(I)=1 IF THERE IS A LOWER CONSTRAINT AT GRID POINT I, C AND INDLW(I)=0 OTHERWISE. IF(INDLW(I))1015,1015,1010 1010 M=M+1 DO 1011 J=1,NUMP V(M,J)=-PWR(I,J) 1011 CONTINUE IF(NUMQ-1)1014,1014,1012 1012 DO 1013 J=NMPP1,NMPQ1 JP1=J+1 V(M,J)=ALTBL(I)*PWR(I,JP1) 1013 CONTINUE 1014 V(M,N1)=0.0 V(M,N1P1)=-X(N1P1)*ALTBL(I)*PWR(I,NMPP1) 1015 CONTINUE C C FOR EACH GRID POINT AT WHICH F IS TO BE APPROXIMATED SET C THE CONSTRAINTS -WT*P+F*Q-W.LE.0.0 AND WT*P-F*Q-W.LE.0.0, C WITH THE CONSTANT Q(1)*F TERM TRANSPOSED. 1016 DO 1022 I=1,NUMGR C INDF(I)=1 IF F IS TO BE APPROXIMATED AT GRID POINT I, AND C INDF(I)=0 IF F IS TO BE IGNORED AT GRID POINT I. IF(INDF(I))1022,1022,1017 1017 M=M+2 MM1=M-1 DO 1018 J=1,NUMP V(MM1,J)=-WTBLE(I)*PWR(I,J) V(M,J)=-V(MM1,J) 1018 CONTINUE IF(NUMQ-1)1021,1021,1019 1019 DO 1020 J=NMPP1,NMPQ1 JP1=J+1 V(MM1,J)=FTBWT(I)*PWR(I,JP1) V(M,J)=-V(MM1,J) 1020 CONTINUE 1021 V(MM1,N1)=-1.0 V(M,N1)=-1.0 V(MM1,N1P1)=-X(N1P1)*FTBWT(I)*PWR(I,NMPP1) V(M,N1P1)=-V(MM1,N1P1) 1022 CONTINUE C C NOW SET THE CONSTRAINTS OF THE FORM -Q.LE.-EPS, WITH THE C CONSTANT Q(1) TERM TRANSPOSED. DO 1027 I=1,NUMGR M=M+1 DO 1023 J=1,NUMP V(M,J)=0.0 1023 CONTINUE IF(NUMQ-1)1026,1026,1024 1024 DO 1025 J=NMPP1,NMPQ1 JP1=J+1 V(M,J)=-PWR(I,JP1) 1025 CONTINUE 1026 V(M,N1)=0.0 V(M,N1P1)=X(N1P1)*PWR(I,NMPP1)-EPS 1027 CONTINUE C C IF WE ARE PLACING A LOWER BOUND ON THE DENOMINATOR (INDICATED C BY THE ONES COEFFICIENT OF IOPT BEING 1), THEN WE BOUND THE C FREE DENOMINATOR COEFFICIENTS (IF THERE ARE ANY) SO THAT WHEN C THE DENOMINATOR IS NORMALIZED IN SINTPQ, THE DENOMINATOR BOUND C WILL NOT BE VIOLATED. IF(IOPT-(IOPT/10)*10)10090,10090,10050 10050 IF(NUMQ-1)10090,10090,10060 C SET THE CONSTRAINTS OF THE FORM -Q(J) .LE. 1.0 AND Q(J) C .LE. 1.0 FOR J=2,...,NUMQ. 10060 DO 10080 I=2,NUMQ M=M+2 MM1=M-1 DO 10070 J=1,N1 V(MM1,J)=0.0 V(M,J)=0.0 10070 CONTINUE NNN1=NUMP+I-1 V(MM1,NNN1)=-1.0 V(M,NNN1)=1.0 V(MM1,N1P1)=1.0 V(M,N1P1)=1.0 10080 CONTINUE C C NOW SET THE BOTTOM ROW. TO MINIMIZE W WE MAXIMIZE -W. 10090 MP1=M+1 DO 6 J=1,N1P1 V(MP1,J)=0.0 6 CONTINUE V(MP1,N1)=1.0 RETURN END SUBROUTINE SINTPQ(NUMP,NUMQ,X) C***BEGIN PROLOGUE SINTPQ C***REFER TO SDFCOR C***ROUTINES CALLED (NONE) C***PURPOSE THIS SUBROUTINE PUTS THE DENOMINATOR COEFFICIENTS C INTO THE PROPER LOCATIONS IN X AND NORMALIZES P/Q C DURING THE LOEB INITIALIZATION OF THE DIFFERENTIAL C CORRECTION ALGORITHM. C***END PROLOGUE SINTPQ DIMENSION X(16) C C***FIRST EXECUTABLE STATEMENT SINTPQ NMPP1=NUMP+1 NMPPQ=NUMP+NUMQ NMQ1=NUMQ-1 C PUT THE COEFFICIENTS OF Q IN THEIR PROPER PLACES IN X. IF(NUMQ-1)1002,1002,1001 1001 DO 1 J=1,NMQ1 K=NMPPQ+1-J X(K)=X(K-1) 1 CONTINUE 1002 X(NMPP1)=X(NMPPQ+1) C C WE DIVIDE BY THE LARGEST ABSOLUTE VALUE OF THE C COEFFICIENTS OF Q TO NORMALIZE P AND Q. AMAX=0.0 DO 12 J=NMPP1,NMPPQ IF(ABS(X(J))-AMAX)12,12,11 11 AMAX=ABS(X(J)) 12 CONTINUE DO 13 J=1,NMPPQ X(J)=X(J)/AMAX 13 CONTINUE RETURN END SUBROUTINE SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X, *V,QMIN,IQMIN) C***BEGIN PROLOGUE SPQERD C***REFER TO SDFCOR C***ROUTINES CALLED (NONE) C***PURPOSE THIS SUBROUTINE COMPUTES THE VALUES OF P, Q, AND C THE ERROR F-WT*P/Q, AS WELL AS THE UNIFORM ERROR NORM C AND THE MINIMUM OF THE DENOMINATOR, FOR THE C DIFFERENTIAL CORRECTION ALGORITHM. C***END PROLOGUE SPQERD DIMENSION FTBWT(101),PWR(101,15),V(534,17),WTBLE(101), * X(16),INDF(101) C C FOR EACH GRID POINT I THIS SUBROUTINE COMPUTES P,Q, AND C THE ERROR F-WT*P/Q. THE UNIFORM ERROR NORM C MAX ABS(F-WT*P/Q) IS ALSO COMPUTED. RECALL THAT IF WE ARE C DOING WEIGHTED ERROR CURVE APPROXIMATION, F WILL ACTUALLY C BE WT*F. C C SET MACHINE DEPENDENT PARAMETERS FOR SUBROUTINE SPQERD. C SET SPCMN=B**(-T), WHERE B IS THE BASE AND T IS THE NUMBER C OF BASE B DIGITS IN THE MANTISSA. SPCMN IS THE MINIMUM C RELATIVE SPACING ABS((X1-X2)/X2) BETWEEN TWO SUCCESSIVE C FLOATING POINT NUMBERS, SO IT IS THE SPACING BETWEEN TWO C SUCCESSIVE FLOATING POINT NUMBERS IN THE CLOSED INTERVAL C (0.1,1.0). WE ALSO HAVE SPCMN=10.0**(-T*(LOG10)(B))= C 10.0**(-TNMAN), WHERE TNMAN IS THE BASE 10 EQUIVALENT OF C THE LENGTH OF THE MANTISSA. C C***FIRST EXECUTABLE STATEMENT SPQERD SPCMN=R1MACH(3) C SET EPS2=100.0*SPCMN. THUS EPS2=10.0**(-(TNMAN-2)). C THE DENOMINATOR WILL NOT BE ALLOWED TO BE SMALLER THAN C EPS2 AT ANY GRID POINT. THIS WILL AVOID A DIVIDE FAULT C AND ALSO REJECT SOME BAD APPROXIMATIONS. C EPS2=100.0*SPCMN C END OF SETTING MACHINE DEPENDENT PARAMETERS FOR SPQERD. C NMPP1=NUMP+1 NMPPQ=NUMP+NUMQ N=NMPPQ+1 NP1=N+1 NN=2*NUMGR+1 I=0 IDBM=-1 IDB=0 V(NN,NP1)=0.0 ABAD=1.0 1 I=I+1 IF(I-NUMGR)3,3,1001 C REPLACE V(NN,NP1) BY -V(NN,NP1)-2.0 IF THE DENOMINATOR AT C SOME POINT WAS VERY SMALL IN ABSOLUTE VALUE OR NEGATIVE. 1001 V(NN,NP1)=ABAD*V(NN,NP1)+ABAD-ABS(ABAD) RETURN 3 IDBM=IDBM+2 IDB=IDB+2 C C COMPUTE THE VALUE OF P AT GRID POINT I AND PUT IT IN C V(2*I-1,N+1). V(IDBM,NP1)=0.0 DO 4 J=1,NUMP V(IDBM,NP1)=V(IDBM,NP1)+X(J)*PWR(I,J) 4 CONTINUE C COMPUTE THE VALUE OF Q AT GRID POINT I AND PUT IT IN C V(2*I,N+1). V(IDB,NP1)=0.0 DO 5 J=NMPP1,NMPPQ V(IDB,NP1)=V(IDB,NP1)+X(J)*PWR(I,J) 5 CONTINUE C KEEP TRACK OF THE MINIMUM OF THE DENOMINATOR IN QMIN AND THE C INDEX OF THE GRID POINT WHERE IT OCCURS IN IQMIN. IF(I-1)10010,10010,10000 10000 IF(V(IDB,NP1)-QMIN)10010,1003,1003 10010 QMIN=V(IDB,NP1) IQMIN=I 1003 IF(V(IDB,NP1)-EPS2)1005,10020,10020 C NOW COMPUTE THE ERROR AT GRID POINT I AND PUT IT IN C V(I,N). KEEP TRACK OF THE ERROR NORM IN V(2*NUMGR+1,N+1). C IF INDF(I)=0, THEN F IS TO BE IGNORED AT GRID POINT I, C SO SET V(I,N)=0.0 AND LOOK AT THE NEXT GRID POINT. 10020 IF(INDF(I))10030,10030,1004 10030 V(I,N)=0.0 GO TO 1 1004 V(I,N)=FTBWT(I)-WTBLE(I)*V(IDBM,NP1)/V(IDB,NP1) ABERR=ABS(V(I,N)) IF(ABERR-V(NN,NP1))1,1,6 6 V(NN,NP1)=ABERR GO TO 1 C C HERE V(2*I,N+1) IS VERY SMALL IN ABSOLUTE VALUE OR NEGATIVE. C WE SET THE ERROR EQUAL TO ZERO AND REPLACE THE FINAL ERROR C NORM BY ITS NEGATIVE, MINUS 2.0, AS A WARNING. THIS C SITUATION SOMETIMES OCCURS IF THE PREVIOUS APPROXIMATION C WAS NEAR BEST, OR IF AN APPROXIMATION COMPUTED ON A COARSE C GRID WAS USED TO INITIALIZE THE MAIN ALGORITHM ON A FINE C GRID. 1005 ABAD=-1.0 V(I,N)=0.0 C TO INCREASE THE PROBABILITY THAT THE DENOMINATOR OF THE C NEXT COMPUTED APPROXIMATION WILL BE .GE. EPS2 IF THE MAIN C PROGRAM CONTINUES ON, WE RESET THE DENOMINATOR OF THE C PRESENT APPROXIMATION AT THE BAD POINTS. V(IDB,NP1)=10.0*EPS2 GO TO 1 END SUBROUTINE SLNPRO(V,M,N,IYRCT,X,INDIC) C***BEGIN PROLOGUE SLNPRO C***REFER TO SDFCOR C***ROUTINES CALLED SJELIM C***PURPOSE THIS SUBROUTINE SOLVES THE LINEAR PROGRAMMING PROBLEM C MAXIMIZE Z = -V(M+1,1)*X(1)-...-V(M+1,N)*X(N) C WHERE X(1),...,X(N) ARE FREE, SUBJECT TO C V(I,1)*X(1)+...+V(I,N)*X(N) .LE. V(I,N+1), FOR I=1,..,N. C (INFORMATION CONCERNING TOLERANCES AND BASIC VARIABLES C IS ALSO TRANSMITTED USING M, N, AND IYRCT.) C***REFERENCES AVDEYEVA, L. I. AND ZUKHOVITSKIY, S. I., C LINEAR AND CONVEX PROGRAMMING, C SAUNDERS, PHILADELPHIA, 1966. C***END PROLOGUE SLNPRO DIMENSION V(534,17),IYRCT(533),X(16),Y(533), * IXRCT(533),IYCCT(16) C C GIVEN INTEGERS M AND N (WITH M.GE.N) AND A MATRIX V, C THIS SUBROUTINE SOLVES THE LINEAR PROGRAMMING PROBLEM C MAXIMIZE Z=-V(M+1,1)X(1)-...-V(M+1,N)X(N)+V(M+1,N+1) C SUBJECT TO THE CONSTRAINTS C V(I,1)X(1)+...+V(I,N)X(N).LE.V(I,N+1), I=1,...,M C USING ESSENTIALLY THE METHOD IN THE BOOK BY AVDEYEVA AND C ZUKHOVITSKIY. Y(I)=-V(I,1)X(1)-...-V(I,N)X(N)+V(I,N+1), C I=1,...,M ARE SLACK VARIABLES. THE METHOD HAS 4 PHASES. C C FIRST, XS ARE EXCHANGED FOR YS TO GET A PROBLEM C INVOLVING ONLY NONNEGATIVE VARIABLES, THE YS BEING C SELECTED IN THE ORDER DETERMINED BY IYRCT AND A PIVOTING C STRATEGY. AT THE BEGINNING OF THIS ROUTINE WE MUST HAVE C IYRCT(1) NONPOSITIVE, OR IYRCT MUST CONTAIN SOME C PERMUTATION OF THE INTEGERS 1,...,M (SEE BELOW). C SECOND, THE SLACK CONSTANTS OF THE DUAL PROBLEM ARE MADE C (SIGNIFICANTLY) NONNEGATIVE. C THIRD, THE COST COEFFICIENTS OF THE DUAL PROBLEM ARE MADE C (SIGNIFICANTLY) NONNEGATIVE. C FINALLY, THE SOLUTION VECTOR IS COMPUTED. C C THE VARIABLE INDIC WILL BE GIVEN VALUE C 0, IF A SOLUTION WAS FOUND NORMALLY C 1, IF THERE WAS TROUBLE IN PHASE 1 C 2, IF THERE WAS TROUBLE IN PHASE 2 (EITHER ROUND OFF C ERROR, OR THE ORIGINAL PROBLEM WAS INCONSISTENT OR C UNBOUNDED) C 3, IF THERE WAS TROUBLE IN PHASE 3 (EITHER ROUND OFF C ERROR, OR THE ORIGINAL CONSTRAINTS WERE INCONSISTENT) C 4, IF 300 MODIFIED JORDAN ELIMINATIONS WERE USED IN C PHASES 2 AND 3 C -1, IF A SOLUTION WAS FOUND BUT IN ORDER TO OVERCOME C TROUBLE IN PHASE 2 OR 3 IT WAS NECESSARY TO TEMPORARILY C RELAX THE RESTRICTION ON DIVISION BY NUMBERS WITH SMALL C ABSOLUTE VALUE. THUS THERE IS AN INCREASED CHANCE OF C SERIOUS ROUNDOFF ERROR IN THE RESULTS. C -2, IF A SOLUTION WAS FOUND NORMALLY, EXCEPT THAT C THE PARAMETERS REA AND REA1 WERE INCREASED BY A SIGNAL C FROM THE CALLING PROGRAM (NAMELY, M=-M). THE INCREASED C TOLERANCES MAY HAVE ALLOWED MORE ERROR THAN USUAL. C -3, IF IN ORDER TO COMPLETE PHASE 1 IT WAS NECESSARY TO C TEMPORARILY RELAX THE RESTRICTION ON DIVISION BY NUMBERS C WITH SMALL ABSOLUTE VALUE. THUS THERE IS AN INCREASED C CHANCE OF SERIOUS ROUNDOFF ERROR IN THE RESULTS. C -4, IF A SOLUTION WAS FOUND NORMALLY, EXCEPT THAT REA AND REA1 C WERE DECREASED BY A SIGNAL FROM THE CALLING PROGRAM (NAMELY C N=-N) IN ORDER TO TRY FOR MORE ACCURACY. THIS INCREASES THE C CHANCES OF SERIOUS ROUNDOFF ERROR IN THE RESULTS. C NOTE THAT INDIC=-3 WILL OVERWRITE (AND THUS CONCEAL) INDIC=-2 C OR INDIC=-4, AND INDIC=-1 WILL OVERWRITE INDIC=-2, -3, OR -4 C C SET MACHINE DE(ENDENT PARAMETERS FOR SUBROUTINE SLNPRO. C SET SPCMN=B**(-T), WHERE B IS THE BASE AND T IS THE NUMBER C OF BASE B DIGITS IN THE MANTISSA. SPCMN IS THE MINIMUM C RELATIVE SPACING ABS((X1-X2)/X2) BETWEEN TWO SUCCESSIVE C FLOATING POINT NUMBERS, SO IT IS THE SPACING BETWEEN TWO C SUCCESSIVE FLOATING POINT NUMBERS IN THE CLOSED INTERVAL C (0.1,1.0). WE ALSO HAVE SPCMN=10.0**(-T*(LOG10)(B))= C 10.0**(-TNMAN), WHERE TNMAN IS THE BASE 10 EQUIVALENT OF C THE LENGTH OF THE MANTISSA. C C***FIRST EXECUTABLE STATEMENT SLNPRO SPCMN=R1MACH(3) C SET REA (ROUND OFF ERROR ADJUSTMENT) = C MAX(10.0**(-8),100.0*SPCMN). THUS REA=10.0**(-EXREA), C WHERE EXREA=MIN(8,TNMAN-2). C DIVISION BY NUMBERS .LE. REA IN ABSOLUTE VALUE WILL NOT BE C PERMITTED. REA=100.0*SPCMN IF(REA-1.0E-8)10000,10010,10010 10000 REA=1.0E-8 C SET REA1=10.0*SPCMN. THUS REA1=10.0**(-(TNMAN-1)). C NUMBERS IN ROW M+1 OR COLUMN N+1 WHICH ARE .LE. REA1 IN C ABSOLUTE VALUE WILL BE TREATED AS ZEROES. SLNPRO ASSUMES C THAT 0.0.LT.REA1.LE.REA. 10010 REA1=10.0*SPCMN C END OF INITIAL SETTING OF MACHINE DEPENDENT PARAMETERS FOR C SLNPRO. THESE PARAMETERS MAY BE ADJUSTED BY A COMMAND FROM C SDFCOR. C INDIC=0 C M BEING NEGATIVE IS A SIGNAL TO INCREASE REA AND REA1, C THUS TREATING MORE NUMBERS WITH SMALL ABSOLUTE VALUES AS C ZEROES. THIS MAY GIVE THIS ROUTINE A BETTER CHANCE TO C SUCCEED, BUT MAY ALSO CAUSE LARGER ERRORS. IF(M)1001,10001,10001 C RESET M. 1001 M=-M REA=SQRT(REA) REA1=SQRT(REA1) INDIC=-2 C N BEING NEGATIVE IS A SIGNAL TO DECREASE REA AND REA1 TO TRY C FOR MORE ACCURACY. AMONG OTHER THINGS, THIS MAKES IT MORE C LIKELY THAT THE PREVIOUS VERTEX WILL BE RETAINED IN PHASE 1 C BELOW, BUT IT ALSO COULD INCREASE ROUND OFF ERROR. 10001 IF(N)10002,1002,1002 C RESET N. 10002 N=-N REA=REA1 REA1=REA1/100.0 INDIC=-4 C PRESERVE REA IN CASE IT MUST BE TEMPORARILY RELAXED. C IRLAX=0 INDICATES REA IS NOT RELAXED AT THIS STAGE. 1002 REAKP=REA IRLAX=0 C IN COLUMN N+1, NUMBERS .LE. REA2 IN ABSOLUTE VALUE WILL BE C TREATED AS ZEROES. REA2=REA1 NP1=N+1 MP1=M+1 KTJOR=0 IBACK=0 C SET V(MP1,NP1)=0.0 SO THE DESCRIPTIONS IN AND FOLLOWING THE C PROLOGUE WILL AGREE. V(MP1,NP1)=0.0 C THE ONLY REASON FOR THE FOLLOWING THREE STATEMENTS IS TO C AVOID THE ERROR MESSAGE ON SOME MACHINES THAT THESE C VARIABLES HAVE NOT BEEN ASSIGNED A VALUE. THEY WILL BE C REASSIGNED A VALUE BEFORE THE PROGRAM REACHES A SPOT WHERE C THEY WILL ACTUALLY BE USED. DIST=1.0 AMPRV=1.0 AMPR2=1.0 C SET IXRCT. IXRCT(I)=0 MEANS SOME Y IS IN ROW I, WHILE C IXRCT(I)=K.NE.0 MEANS X(K) IS IN ROW I. DO 1 I=1,M IXRCT(I)=0 1 CONTINUE C C EXCHANGE THE XS AT THE TOP OF THE TABLE FOR YS. C IF IYRCT(1) IS NONPOSITIVE, WE SET IYRCT AND CHOOSE THE C LARGEST POSSIBLE RESOLVENTS FOR THE EXCHANGES. IF C IYRCT(1) IS POSITIVE, IYRCT WILL HAVE BEEN PREVIOUSLY SET C AND WE TRY TO EXCHANGE IN ROWS IYRCT(1),...,IYRCT(N), C STILL EMPLOYING A PIVOTING STRATEGY, BUT IF WE CANNOT, WE C EXCHANGE IN ROWS IYRCT(N+1),...,IYRCT(M). IF(IYRCT(1))1003,1003,1005 1003 I10=1 I20=M C IF WE HAVE NO INFORMATION FROM A PREVIOUS VERTEX, WE GIVE C UP A LITTLE ACCURACY IN COLUMN N+1 TO HAVE A BETTER CHANCE C OF SUCCESS. REA2=REA C THIS ROUTINE HAS A BACKTRACKING OPTION WHICH SOMETIMES C INCREASES ACCURACY BUT SOMETIMES LEADS TO FAILURE DUE TO C CYCLING. IT IS SUGGESTED THAT THIS OPTION BE EMPLOYED IF C INFORMATION ABOUT A STARTING VERTEX IS AVAILABLE, AND C OTHERWISE BE DISABLED BY SETTING IBACK=1. IBACK=1 DO 1004 I=1,M IYRCT(I)=I 1004 CONTINUE GO TO 1006 1005 I10=1 I20=N 1006 J=0 C SET THE LOWER BOUND ON THE ABSOLUTE VALUE OF A RESOLVENT IN C PHASE 1. ALSO SET IFAIL=0 TO INDICATE THE RESOLVENT SEARCH C IN THIS COLUMN HAS NOT FAILED. REA3=REA IFAIL=0 2 J=J+1 IF(J-N)1007,1007,9 C SET I1, I2 ACCORDING TO THE STRATEGY WE ARE USING. 1007 I1=I10 I2=I20 AMAX=0.0 C SEARCH FOR A RESOLVENT IN ROWS IYRCT(I1),...,IYRCT(I2). 10003 DO 1012 I=I1,I2 IYTMP=IYRCT(I) IF(IXRCT(IYTMP))1012,1009,1012 1009 ABSV=ABS(V(IYTMP,J)) IF(ABSV-AMAX)1012,1012,1011 1011 IYRI=IYTMP AMAX=ABSV 1012 CONTINUE C CHECK TO SEE IF THE PROSPECTIVE RESOLVENT IS LARGE ENOUGH C IN ABSOLUTE VALUE. IF(AMAX-REA3)1013,1013,7 C EXCHANGE X(J) FOR Y(IYRI). 7 CALL SJELIM(MP1,1,NP1,IYRI,J,V) IXRCT(IYRI)=J IYCCT(J)=IYRI C IYCCT(J)=IYRI MEANS Y(IYRI) IS IN COLUMN J. C RESET REA3 AND IFAIL SINCE WE SUCCESSFULLY FOUND A RESOLVENT IN C THIS COLUMN, AND THE RESOLVENT SEARCH IN THE NEXT COLUMN HAS C NOT FAILED. REA3=REA IFAIL=0 GO TO 2 C WE HAVE NOT FOUND A SUITABLE RESOLVENT IN ROWS IYRCT(I1), C ...IYRCT(I2). IF I2.LT.M WE SEARCH THE REST OF COLUMN J. 1013 IF(I2-M)1014,10004,10004 1014 I1=I2+1 I2=M GO TO 10003 C HERE WE FAILED TO FIND A RESOLVENT IN COLUMN J WITH ABSOLUTE C VALUE .GT. REA3. IF IFAIL=0 WE SET INDIC=-3 AND TRY AGAIN C WITH REA3 REDUCED. IF THIS HAS ALREADY BEEN TRIED WE SET C INDIC=1 AND RETURN. 10004 IF(IFAIL)10005,10005,8 10005 IFAIL=1 INDIC=-3 REA3=REA1 GO TO 1007 8 INDIC=1 RETURN C C REARRANGE THE ROWS OF V SO THAT X(1),...,X(N) COME FIRST C IN THAT ORDER. REDEFINE IYRCT SO THAT AFTER THE C REARRANGEMENT IS DONE, IYRCT(I)=K WILL MEAN Y(K) IS IN C ROW I (FOR I GREATER THAN N). 9 DO 10 I=1,M IYRCT(I)=I 10 CONTINUE IROW=0 11 IROW=IROW+1 IF(IROW-M)12,12,20 12 IF(IXRCT(IROW))13,11,13 13 IF(IXRCT(IROW)-IROW)14,11,14 C NOW X(L) IS IN ROW IROW, BUT WE WANT IT IN ROW L. 14 L=IXRCT(IROW) LL=IXRCT(L) IF(LL)15,16,15 C X(L) IS IN ROW IROW, WHILE X(LL) IS IN ROW L. 15 IXRCT(IROW)=LL IXRCT(L)=L GO TO 17 C X(L) IS IN ROW IROW, WHILE Y(IYRCT(L)) IS IN ROW L. 16 IXRCT(IROW)=0 IYRCT(IROW)=IYRCT(L) IXRCT(L)=L C NOW EXCHANGE THE CONTENTS OF ROWS IROW AND L. 17 DO 18 J=1,NP1 TEMP=V(IROW,J) V(IROW,J)=V(L,J) V(L,J)=TEMP 18 CONTINUE IF(IXRCT(IROW))19,11,19 19 IF(IXRCT(IROW)-IROW)14,11,14 C NOW IXRCT IS NO LONGER NEEDED, SO STORE THE PRESENT IYCCT C IN IT. 20 DO 21 I=1,N IXRCT(I)=IYCCT(I) 21 CONTINUE C END OF PHASE 1. C C THE FIRST N ROWS OF V GIVE THE XS IN TERMS OF CERTAIN C YS. THESE ROWS WILL NOT BE CHANGED BY LATER OPERATIONS. C C WE NOW ATTACK THE MAXIMIZATION PROBLEM BY LOOKING AT THE C DUAL PROBLEM OF MINIMIZING A FORM GIVEN BY THE C COEFFICIENTS IN V(N+1,N+1) THROUGH V(M,N+1) WITH SLACK C TERMS IN THE BOTTOM ROW OF V. C SEARCH ROW M+1 FOR A SIGNIFICANTLY NEGATIVE ELEMENT. IF C THERE ARE NONE, PROCEED TO THE ACTUAL MINIMIZATION C PROBLEM. STICK WITH COLUMN JJ UNTIL V(M+1,JJ).GE.-REA1. JJ=0 22 JJ=JJ+1 IF(JJ-N)1015,1015,1035 1015 IF(V(MP1,JJ)+REA1)24,22,22 C C WE HAVE V(M+1,JJ) SIGNIFICANTLY NEGATIVE. SEARCH COLUMN C JJ FOR A POSITIVE ELEMENT, TREATING A VERY SMALL V(I,J) C AS A ZERO. IF THERE ARE NO POSITIVE ELEMENTS THE DUAL C CONSTRAINTS WERE INCONSISTENT, SO THE ORIGINAL PROBLEM WAS C INCONSISTENT OR UNBOUNDED. 24 I=N INAMP=0 25 I=I+1 IF(I-M)1016,1016,1020 1016 IF(V(I,JJ)-REA)25,25,1017 C C NOW V(I,JJ).GT.REA. WE SEARCH ROW I FOR INDICES K SUCH C THAT V(M+1,K).GE.0.0.OR.K.LT.JJ, AND V(I,K).LT.-REA, AND C FIND THE MAXIMUM RATIO (I.E. THE RATIO WITH SMALLEST C ABSOLUTE VALUE, IF V(M+1,K).GE.0.0) V(M+1,K)/V(I,K). IF C THERE IS NO SUCH K WE LOOK AT POSITIVE V(I,K) BELOW. 1017 INDST=0 DO 32 J=1,N IF(V(MP1,J))1018,28,28 1018 IF(J-JJ)28,32,32 28 IF(V(I,J)+REA)29,32,32 29 DIST1=V(MP1,J)/V(I,J) IF(INDST)31,31,30 30 IF(DIST1-DIST)32,32,31 31 DIST=DIST1 INDST=1 K=J 32 CONTINUE IF(INDST)35,35,33 C C WE NOW COMPUTE V(I,JJ)*DIST AND GO ON TO LOOK AT OTHER C ROWS TO MINIMIZE THIS QUANTITY (I.E. TO MAXIMIZE ITS C ABSOLUTE VALUE, IF V(M+1,K).GE.0.0). THIS IS THE NEGATIVE C OF THE CHANGE IN V(M+1,JJ). 33 BMPRV=V(I,JJ)*DIST IF(INAMP)34,34,1019 1019 IF(BMPRV-AMPRV)34,25,25 34 AMPRV=BMPRV INAMP=1 KPMP1=I KPMP2=K C (KPMP1,KPMP2) GIVES THE POSITION OF THE BEST PROSPECTIVE C RESOLVENT FOUND SO FAR. GO TO 25 C C IF THERE WAS NO INDEX K SUCH THAT V(M+1,K).GE.0.0.OR.K.LT. C JJ, AND V(I,K).LT.-REA, WE LOOK FOR THE SMALLEST (I.E. C LARGEST IN ABSOLUTE VALUE) RATIO V(M+1,K)/V(I,K) FOR C V(I,K).GT.REA AND V(M+1,K).LT.0.0, AND PERFORM AN C ELIMINATION WITH RESOLVENT V(I,K). THERE IS AT LEAST ONE C SUCH K, NAMELY JJ. C THIS WILL FINISH PHASE 2 UNLESS BACKTRACKING IS NECESSARY. 35 DIST=1.0 DO 39 J=1,N IF(V(MP1,J))36,39,39 36 IF(V(I,J)-REA)39,39,37 37 DIST1=V(MP1,J)/V(I,J) IF(DIST1-DIST)38,39,39 38 DIST=DIST1 K=J 39 CONTINUE GO TO 49 1020 IF(INAMP)1021,1021,1023 C AT THIS POINT INAMP IS POSITIVE IFF THERE WAS AT LEAST ONE C ELEMENT .GT. REA IN COLUMN JJ. IF THERE WERE NONE, WE C TEMPORARILY RELAX REA AND TRY AGAIN. 1021 IF(IRLAX)1022,1022,41 1022 IRLAX=1 INDIC=-1 REA=REA1 GO TO 24 41 INDIC=2 RETURN C CHECK TO SEE IF V(MP1,KPMP2) IS VERY SMALL IN ABSOLUTE C VALUE OR NEGATIVE. THIS INDICATES DEGENERACY. 1023 IF(V(MP1,KPMP2)-REA)1024,1024,43 C DO AN ELIMINATION WITH RESOLVENT V(KPMP1,KPMP2). 43 I=KPMP1 K=KPMP2 GO TO 49 C C WE ARE NOW STUCK IN DEGENERATE COLUMN KPMP2. WE SEARCH C EACH DEGENERATE COLUMN IN WHICH WE ARE STUCK FOR A C RESOLVENT WHICH WILL KEEP US FROM GETTING STUCK IN THIS C COLUMN NEXT TIME, AND TO REDUCE THE ROUND-OFF ERROR WE C TAKE THE SMALLEST OF THESE (I.E. LARGEST IN ABSOLUTE C VALUE) AS OUR ACTUAL RESOLVENT. 1024 AMIN=1.0 DO 1034 J=1,N C COLUMN J MAY BE DEGENERATE IF 0.0.LE.V(M+1,J).LE.REA, C OR V(M+1,J).LT.0.0.AND.J.LT.JJ. IF(V(MP1,J))1025,1026,1026 1025 IF(J-JJ)1027,1034,1034 1026 IF(V(MP1,J)-REA)1027,1027,1034 C WE WILL BE STUCK IN COLUMN J IFF THERE IS AN INDEX ID FOR C WHICH V(ID,JJ).GT.REA AND V(ID,J).LT.-REA. IF THIS IS THE C CASE, CHOOSING SUCH AN ID SO THAT V(ID,JJ)/V(ID,J) IS C MINIMIZED (I.E. MAXIMIZED IN ABSOLUTE VALUE) AND TAKING C V(ID,J) AS THE RESOLVENT WILL INSURE THAT WE DONT GET C STUCK IN COLUMN J NEXT TIME. 1027 DIST=1.0 DO 1031 I=NP1,M IF(V(I,JJ)-REA)1031,1031,1028 1028 IF(V(I,J)+REA)1029,1031,1031 1029 DIST1=V(I,JJ)/V(I,J) IF(DIST1-DIST)1030,1031,1031 1030 DIST=DIST1 ID=I 1031 CONTINUE IF(DIST-0.5)1032,1034,1034 C WE HAVE NOW DETERMINED THAT WE ARE STUCK IN COLUMN J. C IF V(ID,J).LT.AMIN THEN V(ID,J) IS THE BEST RESOLVENT C FOUND SO FAR. 1032 IF(V(ID,J)-AMIN)1033,1034,1034 1033 AMIN=V(ID,J) KPMP1=ID KPMP2=J 1034 CONTINUE C THE BEST RESOLVENT IS V(KPMP1,KPMP2), SO WE DO AN C ELIMINATION. GO TO 43 49 KTJOR=KTJOR+1 IF(KTJOR-300)50,50,73 50 CALL SJELIM(MP1,NP1,NP1,I,K,V) ITEMP=IYRCT(I) IYRCT(I)=IYCCT(K) IYCCT(K)=ITEMP C RESET REA AND IRLAX. REA=REAKP IRLAX=0 C IF NOW V(M+1,JJ) HAS BEEN MADE NOT SIGNIFICANTLY NEGATIVE, C WE GO TO THE NEXT COLUMN. OTHERWISE WE TRY AGAIN IN C COLUMN JJ. IF(V(MP1,JJ)+REA1)24,22,22 C IN THE UNLIKELY EVENT THAT SOME V(M+1,J) IS STILL VERY C SIGNIFICANTLY NEGATIVE WE BACKTRACK TO COLUMN J. THIS C COULD NOT HAPPEN IF THERE WERE NO ROUNDOFF ERROR AND WE C COULD ALLOW DIVISION BY NUMBERS WITH VERY SMALL ABSOLUTE C VALUE. OMIT BACKTRACKING IF IBACK=1. 1035 IF(IBACK)1036,1036,51 1036 DO 1038 J=1,N IF(V(MP1,J)+REA)1037,1037,1038 1037 JJ=J GO TO 24 1038 CONTINUE C END OF PHASE 2. C 51 I=N KKK=0 C C SEARCH FOR A SIGNIFICANTLY NEGATIVE ELEMENT BETWEEN C V(N+1,N+1) AND V(N+1,M). IF THERE ARE NONE WE HAVE THE C MINIMAL POINT OF THE DUAL PROBLEM (AND THUS THE MAXIMAL C POINT OF THE DIRECT PROBLEM) ALREADY. 52 I=I+1 IF(I-M)1039,1039,1043 1039 IF(V(I,NP1)+REA2)1040,52,52 C C SEARCH FOR A NEGATIVE ELEMENT IN ROW I, TREATING A NUMBER C WHICH IS VERY SMALL IN ABSOLUTE VALUE AS A ZERO. IF THERE C ARE NO NEGATIVE ELEMENTS THE DUAL PROBLEM WAS UNBOUNDED C BELOW, SO THE ORIGINAL CONSTRAINTS WERE INCONSISTENT. C FIND THE INDEX K OF THE LARGEST (I.E. SMALLEST IN ABSOLUTE C VALUE, IF V(M+1,K).GE.0.0) RATIO V(M+1,K)/V(I,K) WITH C V(I,K).LT.-REA. 1040 INDST=0 DO 58 J=1,N IF(V(I,J)+REA)55,58,58 55 DIST1=V(MP1,J)/V(I,J) IF(INDST)57,57,56 56 IF(DIST1-DIST)58,58,57 57 K=J INDST=1 DIST=DIST1 58 CONTINUE IF(INDST)1041,1041,60 C RELAX REA AND LOOK FOR NEGATIVE ELEMENTS WITH SMALLER C ABSOLUTE VALUE. 1041 IF(IRLAX)1042,1042,59 1042 IRLAX=1 INDIC=-1 REA=REA1 GO TO 1040 59 INDIC=3 RETURN C C COMPUTE THE IMPROVEMENT DIST*V(I,N+1) IN THE VALUE OF THE C FORM USING V(I,K) AS THE RESOLVENT. SET KKK=1 TO INDICATE C A SIGNIFICANTLY NEGATIVE V(I,N+1) WAS FOUND, AND LOOK AT C THE OTHER ROWS TO FIND THE RESOLVENT GIVING THE LARGEST C IMPROVEMENT. 60 BMPR2=DIST*V(I,NP1) C RESET IRLAX SO THAT THE NEXT ROW WHICH NEEDS RELAXING DOES C NOT TERMINATE THE ROUTINE. REA WILL REMAIN RELAXED UNTIL C AFTER THE NEXT ELIMINATION. IRLAX=0 IF(KKK)62,62,61 61 IF(BMPR2-AMPR2)52,52,62 62 KKK=1 KEEP=I KEEP1=K AMPR2=BMPR2 GO TO 52 C KKK=0 HERE IFF NONE OF THE COST COEFFICIENTS ARE C SIGNIFICANTLY NEGATIVE. 1043 IF(KKK)1048,1044,1048 C CHECK TO SEE IF ANY OF THE NUMBERS IN THE BOTTOM ROW HAVE C BECOME VERY SIGNIFICANTLY NEGATIVE. IF SO, WE MUST C BACKTRACK TO PHASE 2 (SEE COMMENT ABOVE STATEMENT 1035). C OMIT BACKTRACKING IF IBACK=1. 1044 IF(IBACK)1045,1045,74 1045 DO 1047 J=1,N IF(V(MP1,J)+REA)1046,1046,1047 1046 JJ=J GO TO 24 1047 CONTINUE GO TO 74 C CHECK TO SEE IF V(MP1,KEEP1) IS VERY SMALL IN ABSOLUTE C VALUE OR NEGATIVE. THIS INDICATES DEGENERACY. 1048 IF(V(MP1,KEEP1)-REA)1049,1049,65 C DO AN ELIMINATION WITH RESOLVENT V(KEEP,KEEP1). 65 I=KEEP K=KEEP1 GO TO 71 C C WE ARE NOW STUCK IN DEGENERATE COLUMN KEEP1. WE SEARCH C EACH DEGENERATE COLUMN IN WHICH WE ARE STUCK FOR A C RESOLVENT WHICH WILL KEEP US FROM GETTING STUCK IN THIS C COLUMN NEXT TIME. IF WE ARE NOT USING THE OPTION C DESCRIBED IN THE COMMENTS PRECEDING STATEMENT 1055, WE C TAKE THE SMALLEST OF THESE (I.E. THE LARGEST IN ABSOLUTE C VALUE) AS OUR ACTUAL RESOLVENT IN ORDER TO REDUCE THE C GROWTH OF ROUND-OFF ERROR. 1049 AMIN=1.0 MXRKN=NP1 DO 1072 J=1,N C COLUMN J MAY BE DEGENERATE IF V(M+1,J).LE.REA. IF(V(MP1,J)-REA)1050,1050,1072 C WE WILL BE STUCK IN COLUMN J IFF THERE IS AN INDEX ID FOR C WHICH V(ID,N+1).LT.-REA2 AND V(ID,J).LT.-REA. IF THIS C IS THE CASE, CHOOSING SUCH AN ID SO THAT V(ID,N+1)/V(ID,J) C IS MAXIMIZED AND TAKING V(ID,J) AS THE RESOLVENT WILL C INSURE THAT WE DONT GET STUCK IN COLUMN J NEXT TIME. 1050 DIST=-1.0 DO 1054 I=NP1,M IF(V(I,NP1)+REA2)1051,1054,1054 1051 IF(V(I,J)+REA)1052,1054,1054 1052 DIST1=V(I,NP1)/V(I,J) IF(DIST1-DIST)1054,1054,1053 1053 DIST=DIST1 ID=I 1054 CONTINUE IF(DIST+0.5)1072,1072,1055 C WE HAVE NOW DETERMINED THAT WE ARE STUCK IN COLUMN J. C THE FOLLOWING STATEMENTS ATTEMPT TO BREAK DEGENERACY C FASTER BY LOOKING ONE ITERATION INTO THE FUTURE, THAT IS, C BY CHOOSING FROM THE PROSPECTIVE RESOLVENTS FOUND ABOVE C THAT ONE WHICH MINIMIZES THE MINIMUM NUMBER OF STICKING C PLACES IN ANY ROW AT THE NEXT STAGE. C BECAUSE OF COMPUTER TIME AND THE POSSIBLE LOSS OF ACCURACY C DUE TO LESSENED PIVOTING (EVEN THOUGH TIES ARE ALWAYS C BROKEN IN FAVOR OF THE RESOLVENT WITH GREATEST ABSOLUTE C VALUE), IT IS SUGGESTED THAT THIS OPTION BE OMITTED IF C INFORMATION WAS AVAILABLE FROM A PREVIOUS VERTEX. THIS C WILL BE THE CASE IFF THE BACKTRACKING OPTION WAS USED, C THAT IS, IFF IBACK=0. 1055 IF(IBACK)1070,1070,1056 C COMPUTE WHAT THE NEW BOTTOM ROW WOULD BE (EXCEPT FOR C POSITION J) IF V(ID,J) WERE USED AS THE RESOLVENT, AND C PUT THE RESULTS INTO Y. 1056 ROWQ=V(MP1,J)/V(ID,J) DO 1058 L=1,N IF(L-J)1057,1058,1057 1057 Y(L)=V(MP1,L)-V(ID,L)*ROWQ 1058 CONTINUE LRKNT=-1 C WE LOOK FOR A ROW WHICH WILL HAVE A SIGNIFICANTLY NEGATIVE C LAST ELEMENT BUT A MINIMUM NUMBER OF PLACES WHERE WE WILL C BE STUCK IN DEGENERATE COLUMNS. LRKNT=-1 MEANS WE HAVE C NOT YET FOUND A ROW WHICH WILL HAVE A SIGNIFICANTLY C NEGATIVE LAST ELEMENT. DO 1068 II=NP1,M IF(II-ID)1059,1068,1059 1059 ROWQ=V(II,J)/V(ID,J) RTCOL=V(II,NP1)-V(ID,NP1)*ROWQ IF(RTCOL+REA2)1060,1068,1068 C IF WE HAVE ALREADY LOCATED A RESOLVENT WHICH WILL FINISH C THE ROUTINE, BUT THE PRESENT PROSPECTIVE RESOLVENT WOULD C GIVE A ROW WITH A SIGNIFICANTLY NEGATIVE LAST ELEMENT, WE C LOOK AT THE NEXT PROSPECTIVE RESOLVENT FOR PIVOTING C PURPOSES. 1060 IF(MXRKN+1)1061,1072,1061 1061 LRKNT=0 C NOW COUNT THE NUMBER (LRKNT) OF STICKING PLACES IN ROW II C AT THE NEXT ITERATION. DO 1065 JJ=1,N IF(JJ-J)1062,1065,1062 1062 IF(Y(JJ)-REA)1063,1063,1065 1063 IF(V(II,JJ)-V(ID,JJ)*ROWQ+REA)1064,1065,1065 1064 LRKNT=LRKNT+1 IF(LRKNT-MXRKN)1065,1065,1068 1065 CONTINUE IF(LRKNT-MXRKN)1067,1066,1068 1066 IF(V(ID,J)-AMIN)1067,1068,1068 1067 MXRKN=LRKNT AMIN=V(ID,J) KEEP=ID KEEP1=J 1068 CONTINUE C LRKNT=-1 HERE WOULD MEAN THIS RESOLVENT WOULD FINISH THE C ROUTINE. IF LRKNT.GE.0 THEN MXRKN.GE.0 ALSO, SO WE WILL C NOT HAVE EARLIER FOUND A RESOLVENT WHICH WILL FINISH THE C ROUTINE. IF(LRKNT+1)1072,1069,1072 1069 IF(MXRKN+1)1071,1070,1071 1070 IF(V(ID,J)-AMIN)1071,1072,1072 1071 MXRKN=-1 AMIN=V(ID,J) KEEP=ID KEEP1=J 1072 CONTINUE C THE BEST RESOLVENT IS V(KEEP,KEEP1), SO WE DO AN C ELIMINATION. GO TO 65 71 KTJOR=KTJOR+1 IF(KTJOR-300)72,72,73 72 CALL SJELIM(MP1,NP1,NP1,I,K,V) ITEMP=IYRCT(I) IYRCT(I)=IYCCT(K) IYCCT(K)=ITEMP C RESET REA AND IRLAX. REA=REAKP IRLAX=0 GO TO 51 73 INDIC=4 RETURN C END OF PHASE 3. WE NOW HAVE THE VERTEX WE ARE SEEKING. C C READ OFF THE Y VALUES FOR THIS VERTEX. 74 DO 75 J=1,N IYCJ=IYCCT(J) Y(IYCJ)=0.0 75 CONTINUE DO 76 I=NP1,M IYRI=IYRCT(I) Y(IYRI)=V(I,NP1) 76 CONTINUE C COMPUTE THE XS FROM THE YS. RECALL THAT IXRCT CONTAINS C THE FORMER IYCCT. DO 78 I=1,N X(I)=V(I,NP1) DO 77 J=1,N IXRJ=IXRCT(J) X(I)=X(I)-V(I,J)*Y(IXRJ) 77 CONTINUE 78 CONTINUE C C NOW PUT THE VALUES IN IYCCT INTO THE FIRST N POSITIONS OF C IYRCT IN DECREASING ORDER SO THAT WHEN SLNPRO IS CALLED C AGAIN THE YS AT THE PRESENT MINIMIZING VERTEX WILL BE C SCANNED FIRST, BEGINNING WITH THOSE CORRESPONDING TO C -1.0.LE.Q(J).LE.1.0. THUS THESE WILL PROBABLY APPEAR IN C THE INITIAL VERTEX AFTER EXCHANGE OF XS AND YS. C TO ACCOMPLISH THIS, MAKE IXRCT(I)=-1 IF IYCCT(J)=I FOR C SOME J, THEN SCAN IXRCT BACKWARDS. DO 79 J=1,N IYCJ=IYCCT(J) IXRCT(IYCJ)=-1 79 CONTINUE K=1 I=MP1 80 I=I-1 IF(I)83,83,81 81 IF(IXRCT(I)+1)80,82,80 82 IYRCT(K)=I K=K+1 GO TO 80 C NOW FILL IN THE REST OF IYRCT BY SCANNING IXRCT AGAIN. 83 I=MP1 84 I=I-1 IF(I)87,87,85 85 IF(IXRCT(I))84,86,86 86 IYRCT(K)=I K=K+1 GO TO 84 87 RETURN END SUBROUTINE SJELIM(L,LL,K,IR,IS,V) C***BEGIN PROLOGUE SJELIM C***REFER TO SLNPRO C***ROUTINES CALLED (NONE) C***PURPOSE THIS SUBROUTINE PERFORMS A MODIFIED JORDAN C ELIMINATION ON THE L-LL+1 BY K MATRIX C CONSISTING OF ROWS LL THROUGH L OF V AND C COLUMNS 1 THROUGH K OF V. THE RESOLVENT C IS V(IR,IS). C***END PROLOGUE SJELIM DIMENSION V(534,17) C C DIVIDE THE ENTRIES IN THE RESOLVENT ROW (EXCEPT FOR THE C RESOLVENT) BY THE RESOLVENT. C C***FIRST EXECUTABLE STATEMENT SJELIM RESOL=V(IR,IS) DO 2 J=1,K IF(J-IS)1001,2,1001 1001 V(IR,J)=V(IR,J)/RESOL 2 CONTINUE C SWEEP OUT IN ALL BUT ROW IR AND COLUMN IS. DO 6 I=LL,L IF(I-IR)1002,6,1002 1002 FACT=-V(I,IS) DO 5 J=1,K IF(J-IS)1003,5,1003 1003 V(I,J)=V(I,J)+V(IR,J)*FACT 5 CONTINUE 6 CONTINUE C DIVIDE THE ENTRIES IN THE RESOLVENT COLUMN (EXCEPT FOR THE C RESOLVENT) BY THE NEGATIVE OF THE RESOLVENT. DO 8 I=LL,L IF(I-IR)1004,8,1004 1004 V(I,IS)=-V(I,IS)/RESOL 8 CONTINUE C REPLACE THE RESOLVENT BY ITS RECIPROCAL. V(IR,IS)=1.0/RESOL RETURN END SUBROUTINE SSETV2(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ, *INDLW,ALTBL,INDUP,UPTBL,INDF,WTBLE,V,M) C***BEGIN PROLOGUE SSETV2 C***REFER TO SDFCOR C***ROUTINES CALLED (NONE) C***PURPOSE THIS SUBROUTINE SETS UP THE LINEAR PROGRAMMING PROBLEM C FOR THE MAIN ITERATIONS OF THE DIFFERENTIAL CORRECTION C ALGORITHM. C***END PROLOGUE SSETV2 DIMENSION FTBWT(101),PWR(101,15),V(534,17),WTBLE(101), * UPTBL(101),ALTBL(101),INDUP(101),INDLW(101), * INDF(101) C C THIS SUBROUTINE SETS UP V FOR THE LINEAR PROGRAMMING C PROBLEM OF MINIMIZING W WITH THE RESTRICTIONS C (ABS(F*Q-WT*P)-DELTK*Q)/QK.LE.W AND C ABS(Q(J)).LE.1.0 FOR ALL J. C EACH RESTRICTION IS BROKEN INTO TWO CONSTRAINTS TO ACCOUNT C FOR THE ABSOLUTE VALUE. THE SUBROUTINE USES THE FACT THAT C THE VALUE OF QK AT GRID POINT I IS STORED IN V(2*I,N+1) C AND DELTK IS IN V(2*NUMGR+1,N+1) FROM A PREVIOUS CALL TO C SPQERD. C C***FIRST EXECUTABLE STATEMENT SSETV2 NMPP1=NUMP+1 NMPPQ=NUMP+NUMQ N=NMPPQ+1 NP1=N+1 DELTK=V(2*NUMGR+1,NP1) C M WILL KEEP TRACK OF WHICH ROW OF V IS BEING SET. AT THE C CONCLUSION OF SSETV2 M WILL BE THE TOTAL NUMBER OF C CONSTRAINTS. M=0 C C SET THE UPPER CONSTRAINTS, IF THERE ARE ANY, EXCEPT FOR C THE LAST COLUMN. THEY WILL BE OF THE FORM P-Q*UP.LE.0.0 IF(IOPT-(IOPT/1000)*1000-100)1006,1001,1001 1001 DO 1005 I=1,NUMGR C INDUP(I)=1 IF THERE IS AN UPPER CONSTRAINT AT GRID POINT C I, AND INDUP(I)=0 OTHERWISE. IF(INDUP(I))1005,1005,1002 1002 M=M+1 DO 1003 J=1,NUMP V(M,J)=PWR(I,J) 1003 CONTINUE DO 1004 J=NMPP1,NMPPQ V(M,J)=-UPTBL(I)*PWR(I,J) 1004 CONTINUE V(M,N)=0.0 1005 CONTINUE C C SET THE LOWER CONSTRAINTS, IF THERE ARE ANY, EXCEPT FOR C THE LAST COLUMN. THEY WILL BE OF THE FORM -P+Q*LW.LE.0.0 1006 IF(IOPT-(IOPT/100)*100-10)1012,1007,1007 1007 DO 1011 I=1,NUMGR C INDLW(I)=1 IF THERE IS A LOWER CONSTRAINT AT GRID POINT I, C AND INDLW(I)=0 OTHERWISE. IF(INDLW(I))1011,1011,1008 1008 M=M+1 DO 1009 J=1,NUMP V(M,J)=-PWR(I,J) 1009 CONTINUE DO 1010 J=NMPP1,NMPPQ V(M,J)=ALTBL(I)*PWR(I,J) 1010 CONTINUE V(M,N)=0.0 1011 CONTINUE C C FOR EACH GRID POINT AT WHICH F IS TO BE APPROXIMATED SET C THE CONSTRAINTS -WT*P+(-DELTK+F)*Q-QK*W.LE.0.0 C AND WT*P+(-DELTK-F)*Q-QK*W.LE.0.0 C (EXCEPT FOR THE LAST COLUMN). 1012 DO 1016 I=1,NUMGR C INDF(I)=1 IF F IS TO BE APPROXIMATED AT GRID POINT I, C AND INDF(I)=0 IF F IS TO BE IGNORED AT GRID POINT I. IF(INDF(I))1016,1016,1013 1013 M=M+2 MM1=M-1 DO 1014 J=1,NUMP V(MM1,J)=-WTBLE(I)*PWR(I,J) V(M,J)=-V(MM1,J) 1014 CONTINUE DEL1=-DELTK+FTBWT(I) DEL2=-DELTK-FTBWT(I) DO 1015 J=NMPP1,NMPPQ V(MM1,J)=DEL1*PWR(I,J) V(M,J)=DEL2*PWR(I,J) 1015 CONTINUE C QK AT GRID POINT I IS STORED IN V(2*I,N+1) FROM AN C EARLIER CALL TO SPQERD. IDB=2*I V(MM1,N)=-V(IDB,NP1) V(M,N)=-V(IDB,NP1) 1016 CONTINUE C NOW THE LAST COLUMN IS NO LONGER NEEDED FOR THE STORAGE OF C QK, SO ZERO IT OUT FOR THOSE ROWS ALREADY SET. DO 1017 I=1,M V(I,NP1)=0.0 1017 CONTINUE C C NOW SET THE CONSTRAINTS OF THE FORM -Q(J).LE.1.0 C AND Q(J).LE.1.0 DO 1019 I=1,NUMQ M=M+2 MM1=M-1 DO 1018 J=1,N V(MM1,J)=0.0 V(M,J)=0.0 1018 CONTINUE NNN=NUMP+I V(MM1,NNN)=-1.0 V(M,NNN)=1.0 V(MM1,NP1)=1.0 V(M,NP1)=1.0 1019 CONTINUE C C IF THE ONES DIGIT OF IOPT EQUALS 1, THEN FORCE Q.GE.EPSQ. IF(IOPT-(IOPT/10)*10)10010,10010,10000 10000 DO 1022 I=1,NUMGR M=M+1 DO 1020 J=1,N V(M,J)=0.0 1020 CONTINUE DO 1021 J=NMPP1,NMPPQ V(M,J)=-PWR(I,J) 1021 CONTINUE V(M,NP1)=-EPSQ 1022 CONTINUE C C NOW SET THE BOTTOM ROW. TO MINIMIZE W WE MAXIMIZE -W. 10010 MP1=M+1 DO 10 J=1,NP1 V(MP1,J)=0.0 10 CONTINUE V(MP1,N)=1.0 RETURN END 1 2 3 1 101 1010 0.0 2.0 1 7 8 1 21 21010 0.03.141592653589793238 1 7 8 1 41 121010 0.03.141592653589793238 2 3 2 5 5 1000 0.0 0.0 4.0 4.0 1 1 2 1 21 0 0.0 2.0 1 1 2 1 21 1 0.0 2.0 1 2 3 1 101 0 0.0 2.0 0 0 0 0 0