c From CMUVM.BITNET!32ZEJ7N Fri Mar 03 08:25:31 EDT 1989 IMPLICIT REAL*8 (A-H,O-Z) C 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(9),ALC(101,16) C C THE APPROXIMATION PROBLEMS RUN WITH THIS DEMONSTRATION DRIVER C ARE AS FOLLOWS. NOTE THAT P(1) MEANS THE FIRST NUMERATOR C COEFFICIENT, WHILE P(1.0) (RESP. P(Z)) MEANS THE VALUE OF THE C NUMERATOR AT 1.0 (RESP. AT Z), AND SIMILARLY FOR Q. 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, (P/Q)(Z) .GE. 0.0 FOR C 1.0 .LT. Z .LE. 2.0, AND Q(Z) .GE. 1000.0*B**(-T) FOR ALL Z, C WHERE T IS THE NUMBER OF BASE B DIGITS IN THE MANTISSA OF C FLOATING POINT NUMBERS. (THIS LAST REQUIREMENT IS USED IN THE C FIRST FOUR JOBS TO AVOID THE POSSIBILITY THAT A BAD DENOMINATOR C VALUE AT A POINT WHERE F IS UNDEFINED, SO THE DIFFERENTIAL C CORRECTION INEQUALITIES DO NOT PROVIDE ANY CONTROL ON THE C DENOMINATOR, MIGHT PREMATURELY TERMINATE THE PROGRAM.) 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 C .LE. PI. (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. C WE REQUIRE (P/Q)(Z) .GE. 0.0 FOR PI/2 .LT. Z .LE. PI AND C Q(Z) .GE. 1000.0*B**(-T) FOR ALL Z. 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). WE REQUIRE C Q(Z1,Z2) .GE. 1000.0*B**(-T) FOR ALL (Z1,Z2). 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 IN JOB 8 THE GRID IS A 101 POINT EQUALLY SPACED SUBDIVISION OF C THE CLOSED INTERVAL (0.0,2.0), F(Z) IS EXP(Z) FOR ALL Z, AND C (P/Q)(Z) IS (P(1)+P(2)*Z+P(3)*Z**2)/Q(1). WE REQUIRE Q(1)=1.0, C (P/Q)(1.0) = F(1.0), AND 2.0*P(3) .GE. 3.0. THIS EXAMPLE C ILLUSTRATES THE USE OF EQUALITY CONSTRAINTS (INCLUDING AN C INTERPOLATORY CONSTRAINT) AND AN ADDITIONAL INEQUALITY C CONSTRAINT FROM OPTION 4. NOTE THAT IT IS NOT NECESSARY TO C EXPLICITLY IMPOSE THE CONSTRAINT Q(1)=1.0 FOR (GENERALIZED) C POLYNOMIAL APPROXIMATION AS LONG AS NO CONSTRAINTS FROM C OPTION 4 HAVE NONZERO RIGHT SIDES, BUT IT MAKES LESS WORK FOR C THE LINEAR PROGRAMMING SUBROUTINE IF THE CONSTRAINT IS IMPOSED. C C JOB 9 IS THE SAME AS JOB 1 EXCEPT THAT THE CONSTRAINTS Q(Z) .GE. C 1000.0*B**(-T) FOR ALL Z ARE REPLACED BY THE CONSTRAINTS Q(Z) .GE. C 0.5 FOR ALL Z, AND THE CONSTRAINTS (P/Q)(Z) .GE. 0.9999 FOR C 0.0 .LE. Z .LE. 0.3 ARE REPLACED BY THE CONSTRAINT P(1) = C (Q(1)+Q(2))/2.0 + 1.0. 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 1011 C 0.0 2.0 C 1 7 8 1 21 C 21011 C 0.03.141592653589793238 C 1 7 8 1 41 C 121011 C 0.03.141592653589793238 C 2 3 2 5 5 C 1001 C 0.0 0.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 1 3 1 1 101 C 2 C 2 1 C 0.0 2.0 C 1 2 3 1 101 C 1013 C 1 0 C 0.0 2.0 C 0 0 0 0 0 C C SET MACHINE DEPENDENT PARAMETERS FOR THE DEMONSTRATION DRIVER. C SET INPUT AND OUTPUT UNIT NUMBERS. NREAD=I1MACH(1) NWRIT=I1MACH(2) C SET PRECISION DEPENDENT CONSTANTS. ONE=1.0D0 ZERO=ONE-ONE TWO=ONE+ONE THREE=ONE+TWO FOUR=TWO+TWO TEN=FOUR+FOUR+TWO SIXTN=TEN+THREE+THREE 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=D1MACH(3) C SET THE ERROR NORM TEST TOLERANCE ERTOL=MAX(10.0**(-11), C 1000.0*SPCMN). ERTOL=SPCMN*TEN**3 IF(ERTOL-TEN**(-11))10,20,20 10 ERTOL=TEN**(-11) C END OF SETTING MACHINE AND PRECISION DEPENDENT CONSTANTS FOR THE C DEMONSTRATION DRIVER. C 20 JOBNM=1 C SET THE ERROR NORMS COMPUTED WITH THE IBM 3090 AT CENTRAL MICHIGAN C UNIVERSITY FOR COMPARISON PURPOSES. THIS MACHINE HAS 14 BASE 16 C DIGITS IN DOUBLE PRECISION. ERTST(1)=0.369559497889575D0 ERTST(2)=0.8792063120*TEN**(-5) ERTST(3)=0.54401686050*TEN**(-4) ERTST(4)=0.447875655532295D0 ERTST(5)=0.22*TEN**(-13) ERTST(6)=0.999998000*TEN**(-6) ERTST(7)=0.200000000000000D0 ERTST(8)=0.148781030293754D0 ERTST(9)=0.380060239141926D0 C 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) C NMPP1=NUMP+1 NMPPQ=NUMP+NUMQ 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 C READ THE OPTION SWITCH IOPT. INSTRUCTIONS FOR SETTING IT C ARE IN THE INITIAL COMMENTS OF CDFCOR. READ(NREAD, 80)IOPT 80 FORMAT(I10) C IONE=IOPT-(IOPT/10)*10 C C IF THE ONES DIGIT OF IOPT IS .GT. 1, READ THE NUMBER MEQ OF C EQUALITY CONSTRAINTS AND THE NUMBER MINEQ OF EXTRA INEQUALITY C CONSTRAINTS (NOT INCLUDING THE DENOMINATOR OR RESTRICTED C RANGE CONSTRAINTS SET ELSEWHERE). IF(IONE-1)87,87,83 83 READ(NREAD,85)MEQ,MINEQ 85 FORMAT(2I5) C 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. 87 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, *ANEP2,T) C 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 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, 160, 160),JOBNM C C INPUT OF BASIS FUNCTIONS FOR JOBS 1, 5, 6, 7, 8, AND 9. 160 PWR(I,1)=ONE 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)=ONE 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, 533, 310),JOBNM C C INPUT OF BASIS FUNCTIONS FOR JOBS 2 AND 3. 230 PWR(I,1)=ONE 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)=ONE 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 C INPUT OF BASIS FUNCTIONS FOR JOB 4. 300 PWR(I,1)=ONE PWR(I,2)=T(IDBM)**2 PWR(I,3)=T(IDB)**2 PWR(I,4)=ONE PWR(I,5)=PWR(I,2)*PWR(I,3) GO TO 430 C C INPUT OF F, INDF, DENOMINATOR LOWER BOUND, AND LOWER CONSTRAINTS C FOR JOBS 1 AND 9, AND ADDITIONAL EQUALITY CONSTRAINT FOR JOB 9. 310 NUM1=3*(NUMGR-1)/20+1 NUM2=(NUMGR-1)/2+1 NUM3=3*(NUMGR-1)/5+1 C IF I=1 SET THE DENOMINATOR LOWER BOUND FOR JOBS 1 AND 9 AND THE C ADDITIONAL EQUALITY CONSTRAINT FOR JOB 9. NOTE THAT THESE THINGS C NEED BE SET ONLY ONCE IN ANY ONE JOB. IF(I-1)311,311,317 311 IF(JOBNM-1)312,312,315 312 EPSQ=SPCMN*TEN**3 GO TO 317 315 EPSQ=ONE/TWO ALC(1,1)=ONE ALC(1,2)=ZERO ALC(1,3)=-ONE/TWO ALC(1,4)=-ONE/TWO ALC(1,5)=ZERO ALC(1,6)=ONE C SET RESTRICTED RANGE CONSTRAINTS FOR JOBS 1 AND 9. 317 IF(I-NUM1) 318, 318, 330 318 IF(JOBNM-1)320,320,340 320 INDF(I)=1 FTBLE(I)=ONE INDLW(I)=1 ALTBL(I)=ONE-TEN**(-4) GO TO 140 330 IF(I-NUM2) 340, 340, 350 340 INDF(I)=1 FTBLE(I)=ONE INDLW(I)=0 GO TO 140 350 IF(I-NUM3) 360, 370, 370 360 INDF(I)=0 INDLW(I)=1 ALTBL(I)=ZERO GO TO 140 370 INDF(I)=1 FTBLE(I)=ZERO INDLW(I)=1 ALTBL(I)=ZERO GO TO 140 C C INPUT OF F, INDF, DENOMINATOR LOWER BOUND, WEIGHT FUNCTION, AND C LOWER CONSTRAINTS FOR JOBS 2 AND 3. 380 NUM1=(NUMGR-1)/2+1 NUM2=3*(NUMGR-1)/5+1 EPSQ=SPCMN*TEN**3 IF(I-NUM1) 390, 390, 400 390 INDF(I)=1 FTBLE(I)=ONE WTBLE(I)=ONE/TWO INDLW(I)=0 GO TO 140 400 IF(I-NUM2) 410, 420, 420 410 INDF(I)=0 INDLW(I)=1 ALTBL(I)=ZERO GO TO 140 420 INDF(I)=1 FTBLE(I)=ZERO WTBLE(I)=ONE INDLW(I)=1 ALTBL(I)=ZERO GO TO 140 C C INPUT OF F, INDF, AND DENOMINATOR LOWER BOUND FOR JOB 4. 430 EPSQ=SPCMN*TEN**3 C WE USE 16.0001 INSTEAD OF 16.0 IN THE NEXT TEST TO AVOID TROUBLE C DUE TO ROUND OFF ERROR. IF(T(IDBM)**2+T(IDB)**2-(SIXTN+TEN**(-4))) 440, 440, 450 440 INDF(I)=1 IF(SIXTN-T(IDBM)**2-T(IDB)**2)447,447,443 443 FTBLE(I)=SQRT(SIXTN-T(IDBM)**2-T(IDB)**2) GO TO 140 447 FTBLE(I)=ZERO GO TO 140 450 INDF(I)=0 GO TO 140 C 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 IF(I-1)470,470,480 470 FTBLE(I)=ONE EPSQ=TEN**(-7) GO TO 140 480 FTBLE(I)=ZERO GO TO 140 C 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)=THREE/(ONE+T(I))-TWO*T(I)+TWO/TEN GO TO 140 510 IF(I-NUM2) 520, 520, 530 520 FTBLE(I)=THREE/(ONE+T(I))+T(I)/TWO-THREE/TEN GO TO 140 530 FTBLE(I)=THREE/(ONE+T(I))-FOUR*T(I)/TEN+(THREE+THREE)/TEN GO TO 140 C C INPUT OF F AND EQUALITY CONSTRAINTS FOR JOB 8. WE NEED SET C ALC ONLY ONCE. 533 FTBLE(I)=EXP(T(I)) IF(I-1)535,535,140 535 DO 537 J=1,NUMP ALC(1,J)=ZERO ALC(2,J)=ONE ALC(3,J)=ZERO 537 CONTINUE ALC(1,NMPP1)=ONE ALC(1,NMPP1+1)=ONE ALC(2,NMPP1)=ZERO ALC(2,NMPP1+1)=EXP(ONE) ALC(3,NUMP)=-TWO ALC(3,NMPP1)=ZERO ALC(3,NMPP1+1)=-THREE C END COMPUTING OF BASIS FUNCTIONS AND F. GO TO 140 C 540 CALL CDFCOR(NUMP,NUMQ,NUMGR,FTBLE,PWR,IOPT,EPSQ,INDLW,ALTBL, *INDUP,UPTBL,INDF,WTBLE,MEQ,MINEQ,ALC,MEQA,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.) C 620 WRITE(NWRIT, 630)NUMDM 630 FORMAT(/33H WE ARE DEALING WITH FUNCTIONS OF,I4, *12H VARIABLES) C WRITE(NWRIT, 640)NUMP,NUMQ,NUMGR,NRWGR,NCLGR 640 FORMAT(/6H P HAS,I4,24H COEFFICIENTS Q HAS,I4, *15H COEFFICIENTS//13H THE GRID HAS,I5, *18H POINTS ARRANGED,I5,5H BY,I5) C WRITE(NWRIT, 650)IOPT 650 FORMAT(/8H IOPT IS,I10) C C IF THE ONES DIGIT OF IOPT EXCEEDS 1, WRITE THE NUMBERS OF C EQUALITY AND EXTRA INEQUALITY CONSTRAINTS, AND THE NUMBER OF C REDUNDANT EQUALITY CONSTRAINTS. IF(IONE-1)657,657,653 653 MEQDF=MEQ-MEQA WRITE(NWRIT,655)MEQ,MEQDF,MINEQ 655 FORMAT(/10H THERE ARE,I5,31H EQUALITY CONSTRAINTS OF WHICH, *I5,15H ARE REDUNDANT//10H THERE ARE,I5, *30H EXTRA INEQUALITY CONSTRAINTS) C 657 IF(NUMDM-2) 660, 680, 660 660 WRITE(NWRIT, 670)ASWP1,ASWP2 670 FORMAT(/22H THE LEFT END POINT IS,E25.15, *//23H THE RIGHT END POINT IS,E25.15) GO TO 700 680 WRITE(NWRIT, 690)ASWP1,ASWP2,ANEP1,ANEP2 690 FORMAT(/40H THE COORDINATES OF THE SW AND NE POINTS, *4H ARE/2E25.15,6H AND/2E25.15) C C WRITE IFLAG AND INDLP. 700 WRITE(NWRIT, 710)IFLAG,INDLP 710 FORMAT(/9H IFLAG IS,I5,13H INDLP IS,I4) C 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/23H ACHIEVED AT GRID POINT,I5) C C WRITE ITER AND DELTK. DELTK=ERROR(NUMGR+1) WRITE(NWRIT, 730)ITER,DELTK 730 FORMAT(/21H THE ERROR NORM AFTER,I5,13H ITERATIONS, *3H IS,E25.15) C 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/(/3E25.15)) C WRITE(NWRIT, 750)(X(I),I=NMPP1,NMPPQ) 750 FORMAT(/26H THE COEFFICIENTS OF Q ARE/(/3E25.15)) C C WRITE THE ERRORS. 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(/4E20.10) 780 CONTINUE C JOBNM=JOBNM+1 GO TO 30 END SUBROUTINE CDFCOR(NUMP,NUMQ,NUMGR,FTBLE,PWR,IOPT,EPSQ,INDLW, *ALTBL,INDUP,UPTBL,INDF,WTBLE,MEQ,MINEQ,ALC,MEQA,X,ERROR,QMIN, *IQMIN,ITER,INDLP,IFLAG) C***BEGIN PROLOGUE CDFCOR C***DATE WRITTEN 720120 (YYMMDD) C***REVISION DATE 880226 (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, AND ADDITIONAL ARBITRARY LINEAR CONSTRAINTS C ON THE COEFFICIENTS OF P 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, MEQ, MINEQ, ALC, C MEQA (OPTIONAL INPUT OR OUTPUT) C THESE ELEVEN VARIABLES ARE USED FOR INPUT OR OUTPUT FOR THE C EXTRA OPTIONS DESCRIBED IN THE LONG DESCRIPTION BELOW. NONE C OF 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. FOR EXAMPLE, INDLP=44 C USUALLY INDICATES INCORRECT INPUT TO THE PROGRAM. 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 CDFCOR 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), ALC(101,16). 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 UP TO 101 GRID POINTS (SO NUMGR .LE. 101), AND UP TO 101 C ADDITIONAL LINEAR CONSTRAINTS FROM OPTION 4 BELOW (SO C MEQ+MINEQ .LE. 101), EVEN IF ALL THE EXTRA OPTIONS DESCRIBED C IN THE LONG DESCRIPTION BELOW ARE ACTIVATED. SEE THE LONG C DESCRIPTION FOR THE DIMENSION PATTERN REQUIRED IN THE DRIVER C PROGRAM AND SUBROUTINES IF MORE BASIS FUNCTIONS, GRID POINTS, C OR ADDITIONAL LINEAR CONSTRAINTS 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 SEVEN 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 (IF OPTION 4 IS NOT USED) OR TO 3 (IF OPTION 4 IS C USED) 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 CONSTRAINTS (ADDITIONAL LINEAR CONSTRAINTS) C C TO ACTIVATE THIS OPTION, SET THE ONES DIGIT OF C IOPT TO 2 (IF OPTION 1 IS NOT USED) OR TO 3 C (IF OPTION 1 IS USED) AND SET MEQ, MINEQ, AND C ALC AS FOLLOWS. C C MEQ (INPUT) THIS IS THE NUMBER OF ADDITIONAL EQUALITY C CONSTRAINTS. C C MINEQ (INPUT) THIS IS THE NUMBER OF ADDITIOINAL INEQUALITY C CONSTRAINTS. C C ALC (INPUT) FOR I=1,...,MEQ, THE ITH ROW OF ALC SPECIFIES C THE EQUALITY CONSTRAINT C ALC(I,1)*X(1)+...+ALC(I,NUMP+NUMQ)*X(NUMP+NUMQ) = C ALC(I,NUMP+NUMQ+1) C AND FOR I=MEQ+1,...,MEQ+MINEQ, THE ITH ROW OF ALC C SPECIFIES THE INEQUALITY CONSTRAINT C ALC(I,1)*X(1)+...+ALC(I,NUMP+NUMQ)*X(NUMP+NUMQ) .LE. C ALC(I,NUMP+NUMQ+1). C C OPTION 5 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 PROGRAMMING HINT... IF OPTION 5 IS USED WITH A NONTRIVIAL C DENOMINATOR, IT IS A GOOD IDEA TO ALSO USE OPTION 1 OR 4 C TO IMPOSE A POSITIVE LOWER BOUND ON THE DENOMINATOR AT C LEAST AT THE POINTS WHERE F IS IGNORED, SINCE THE C DIFFERENTIAL CORRECTION INQUALITIES ARE NOT PRESENT TO C CONTROL THE DENOMINATOR AT SUCH POINTS, AND A DENOMINATOR C VALUE WHICH IS SMALL IN ABSOLUTE VALUE OR NEGATIVE COULD C TERMINATE THE PROGRAM BEFORE A BEST APPROXIMATION IS FOUND. C A REASONABLE LOWER BOUND IS 1000.0*B**(-T) (WHERE T IS THE C NUMBER OF BASE B DIGITS IN THE MANTISSA). NOTE THAT AT C POINTS WHERE F IS NOT IGNORED, IF C (ABS(F*Q-WT*P)-DELTK*Q)/QK .LE. W .LT. 0.0, THEN C Q .GE. (ABS(F*Q-WT*P)+ABS(W)*QK)/DELTK .GT. 0.0 C C OPTION 6 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 7 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 THREE C POSSIBILITIES. C C IFLAG=7777 MEANS OPTION 4 ABOVE WAS USED, AND THE C LINEAR INEQUALITY CONSTRAINTS APPEARED TO BE C INCONSISTENT. 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, OR ELSE THE APPROXIMATION VIOLATED SOME CONSTRAINT C BY MORE THAN EPRES = 10000.0*B**(-T). NO APPROXIMATION C WAS RETURNED, ALL THE COEFFICIENTS OF P AND Q WERE SET C 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 IF IT IS DESIRED TO CHANGE FROM DOUBLE PRECISION TO SINGLE PRECISION C (FOR EXAMPLE), ONLY THREE THINGS NEED TO BE DONE. C C (1) REPLACE ALL THE STATEMENTS ONE=1.0D0 BY ONE=1.0 C C (2) REPLACE ALL THE STATEMENTS IMPLICIT REAL*8 (A-H,O-Z) BY C C IMPLICIT REAL*8 (A-H,O-Z) C C (3 REPLACE ALL OCCURRENCES OF D1MACH BY R1MACH C (OR MODIFY THIS PROGRAM OR THE SUBPROGRAM D1MACH SUPPLIED BY C NETLIB SO THAT WHEREVER THE STATEMENT SPCMN=D1MACH(3) C APPEARS IN THE CURRENT VERSION, SPCMN WILL BE COMPUTED AS C B**(-T), WHERE B IS THE BASE FOR FLOATING POINT NUMBERS AND C T IS THE NUMBER OF BASE B DIGITS IN THE MANTISSA) C (NOTE THAT THE ONLY VALUES WHICH THIS VERSION OF THIS PACKAGE C USES FROM THE NETLIB SUPPLIED SUBPROGRAMS I1MACH, R1MACH, AND C D1MACH ARE D1MACH(3), I1MACH(1), AND I1MACH(2), WITH THE LATTER C TWO VALUES USED ONLY IN THE SAMPLE DRIVER PROGRAM.) 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), MGR GRID POINTS (THUS C NUMGR .LE. MGR), AND MALC ADDITIONAL LINEAR CONSTRAINTS IN C OPTION 4 ABOVE (THUS MEQ+MINEQ .LE. MALC) 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), ALC(MALC,MFN+1), XTEMP(MFN+1), FTBWT(MGR), IYRCT(MAXCN), C IYSAV(MAXCN), V(MAXCN+1,MFN+2), IYCCT(MFN+1), IXRCT(MAXCN), Y(MAXCN), C EQ(MFN,MFN+1), NSCOL(MFN), IRESL(2*MFN), C C WHERE MAXCN=5*MGR+2*MFN+MALC-2. NOTE THAT THE LAST 11 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, AND MAXCN COULD BE C REDUCED BY MALC IF OPTION 4 IS NOT USED. 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 (AND THE OTHER INEQUALITY CONSTRAINTS) 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 THE EQUALITY CONSTRAINTS (IF ANY, SEE OPTION 4 ABOVE) ARE C HANDLED BY SOLVING THEM FOR MEQA OF THE COEFFICIENTS IN C TERMS OF THE REMAINING COEFFICIENTS (USING A TOTAL C PIVOTING STRATEGY) AND IN EFFECT SUBSTITUTING TO REDUCE C THE NUMBER OF VARIABLES IN THE REST OF THE PROBLEM BY MEQA. 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. THIS IDEA IS DUE C TO ALI SAZEGARI. 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 VIOLATED SOME CONSTRAINT BY MORE THAN C EPRES, OR WHICH HAD A DENOMINATOR WHICH WAS VERY SMALL IN C ABSOLUTE VALUE OR NEGATIVE AT SOME POINT (AND THE MAIN C ALGORITHM WAS UNABLE TO CORRECT THIS PROBLEM). C THE THREE METHODS ARE C C METHOD 1 THE USER SUPPLIES AN INITIAL APPROXIMATION (SEE C OPTION 7 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, ABS(Q(J)) .LE. 1.0 C FOR J=1,...,NUMQ (AND THE OTHER C INEQUALITY CONSTRAINTS). 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. C IF SLNPRO FAILS TO PRODUCE AN APPROXIMATION P/Q WE C TRY WITH BOTH LOOSENED AND TIGHTENED TOLERANCES. 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***REFER TO (NONE) C***ROUTINES CALLED SLNPRO,SPQERD,SSETV,EQELM,SCRCH,EXPND, C ICNCK,NRMLZ C***END PROLOGUE CDFCOR C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION FTBLE(101),PWR(101,15),V(635,17),WTBLE(101), * UPTBL(101),ALTBL(101),INDUP(101),INDLW(101), * ERROR(102),IYRCT(634),IYSAV(634),XTEMP(16),X(16), * INDF(101),FTBWT(101),ALC(101,16),EQ(15,16),NSCOL(15), * IRESL(30) C C SET MACHINE DEPENDENT PARAMETERS FOR SUBROUTINE CDFCOR. 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 CDFCOR SPCMN=D1MACH(3) C SET PRECISION DEPENDENT CONSTANTS FOR CDFCOR. ONE=1.0D0 ZERO=ONE-ONE TWO=ONE+ONE FOUR=TWO+TWO TEN=FOUR+FOUR+TWO TENTH=ONE/TEN 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=SPCMN*TEN**4 C SET TOLEL=100.0*SPCMN. IN SUBROUTINE EQELM, NUMBERS WITH C ABSOLUTE VALUE .LE. TOLEL WILL BE TREATED AS ZEROES. C THIS WILL ALSO BE THE CASE IN SUBROUTINE NRMLZ. TOLEL=TEN*TEN*SPCMN C END OF SETTING MACHINE AND PRECISION DEPENDENT PARAMETERS FOR CDFCOR. C ITLIM=100 ITER=0 NMPP1=NUMP+1 NMPPQ=NUMP+NUMQ N=NMPPQ+1 NP1=N+1 NLP=N NLP1=NP1 IADJS=0 INDIC=0 IRET1=0 ILBFL=0 ILST1=0 QMIN=ZERO IQMIN=0 INITL=0 IFLAG=0 IPBIT=1 IONE=IOPT-(IOPT/10)*10 C IF THE ONES DIGIT OF IOPT IS .LE. 1, SET MEQ AND MINEQ TO ZERO. IF(IONE-1)100,100,200 100 MEQ=0 MINEQ=0 200 MEQA=0 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)300,500,500 300 DO 400 I=1,NUMGR INDF(I)=1 400 CONTINUE C C IF WE ARE DEALING WITH UNWEIGHTED APPROXIMATION SET THE C WEIGHT FUNCTION IDENTICALLY EQUAL TO 1. 500 IF(IOPT-(IOPT/100000)*100000-10000)600,800,800 600 DO 700 I=1,NUMGR WTBLE(I)=ONE 700 CONTINUE C 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. 800 IF(IOPT-(IOPT/100000)*100000-20000)900,1200,1200 900 DO 1100 I=1,NUMGR IF(INDF(I))1100,1100,1000 1000 FTBWT(I)=FTBLE(I) 1100 CONTINUE GO TO 1500 1200 DO 1400 I=1,NUMGR IF(INDF(I))1400,1400,1300 1300 FTBWT(I)=WTBLE(I)*FTBLE(I) 1400 CONTINUE C C IF THERE ARE EQUALITY CONSTRAINTS, ATTEMPT TO SOLVE FOR MEQA C OF THE COEFFICIENTS IN TERMS OF THE REMAINING NUMP+NUMQ-MEQA C COEFFICIENTS, WHERE MEQA IS THE NUMBER OF INDEPENDENT C EQUALITY CONSTRAINTS. 1500 IF(MEQ)2000,2000,1600 C HERE THERE ARE MEQ EQUALITY CONSTRAINTS, AND WE COPY THE FIRST C MEQ ROWS OF ALC INTO EQ. 1600 DO 1800 I=1,MEQ DO 1700 J=1,N EQ(I,J)=ALC(I,J) 1700 CONTINUE 1800 CONTINUE CALL EQELM(NMPPQ,MEQ,TOLEL,EQ,NSCOL,IRESL,MEQA,IFLAG) C HERE IFLAG WILL BE EITHER 0 OR 7777. IN THE LATTER CASE THE C EQUALITY CONSTRAINTS WERE INCONSISTENT AND WE RETURN. IF(IFLAG)1900,1900,9600 C RESET N AND NP1 FOR USE IN SLNPRO. 1900 NLP=N-MEQA NLP1=NP1-MEQA 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. 2000 IF(IOPT/100000)2200,2200,2100 C SET INITL=-1 TO INDICATE WE ARE USING A USER SUPPLIED INITIAL C APPROXIMATION. 2100 INITL=-1 GO TO 3400 C C HERE WE ATTEMPT TO COMPUTE AN INITIAL APPROXIMATION USING C LOEBS ALGORITHM. 2200 INITL=0 C WE TRY TO FIND P, Q SUCH THAT Q IS POSITIVE (AND NOT TOO C SMALL) AT EVERY GRID POINT, P AND Q SATISFY THE CONSTRAINTS C (INCLUDING THE CONSTRAINTS MAX(ABS(Q(J))) .LE. 1.0), AND C MAX(ABS(F*Q-WT*P)) IS MINIMIZED AT THOSE GRID POINTS WHERE C F IS TO BE APPROXIMATED. 2300 CALL SSETV(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,INDLW,ALTBL, *INDUP,UPTBL,INDF,WTBLE,MEQ,MINEQ,ALC,0,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 REDUCE V IF MEQA .GT. 0 IF(MEQA)2500,2500,2400 2400 CALL SCRCH(M,NMPPQ,MEQA,NSCOL,IRESL,EQ,V) 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. 2500 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 NLP=-NLP AS A SIGNAL TO TIGHTEN THE TOLERANCES IN SLNPRO C (I. E. DECREASE REA AND REA1) AND TRY ONCE AGAIN. SLNPRO WILL C RESET NLP TO ITS ORIGINAL VALUE. NLP=NLP*(1+IADJS-IADJS**2) CALL SLNPRO(V,M,NLP,IYRCT,X,INDIC) C SAVE INDIC IN CASE THE APPROXIMATION JUST COMPUTED IS THE ONE C RETURNED. IRET1=INDIC C IF THIS WAS THE FIRST LOEB ATTEMPT WE SAVE INDIC IN ILBFL FOR C OUTPUT PURPOSES IN CASE BOTH LOEB AND THE FIRST MAIN C ITERATION FAIL. IF(IADJS)2600,2600,2700 2600 ILBFL=INDIC C IF INDIC IS POSITIVE SOMETHING WENT WRONG IN SLNPRO. 2700 IF(INDIC)3200,3200,2800 2800 IF(IADJS-1)2900,2900,3000 C BEFORE CONCEDING DEFEAT WE WILL TRY AGAIN WITH CHANGED C TOLERANCES. 2900 IADJS=IADJS+1 GO TO 2300 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. 3000 INITL=1 DO 3100 J=1,NMPPQ X(J)=ZERO 3100 CONTINUE X(NMPP1)=ONE GO TO 3500 C C HERE INDIC WAS NONPOSITIVE, SO SLNPRO PRODUCED AN APPROXIMATION. C IF MEQA .GT. 0, EXPAND X TO INCLUDE THE MEQA VARIABLES C SOLVED FOR IN SUBROUTINE EQELM. 3200 IF(MEQA)3400,3400,3300 3300 CALL EXPND(NMPPQ,MEQA,NSCOL,IRESL,EQ,X) C NORMALIZE P AND Q IF POSSIBLE. 3400 CALL NRMLZ(NUMP,NUMQ,MEQ,MINEQ,NSCOL,IRESL,ALC,TOLEL,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. 3500 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X, *V,QMIN,IQMIN) DELTK=V(2*NUMGR+1,NP1) DLOLD=DELTK C RESET IADJS. IADJS=0 C PUT THE ERRORS AND ERROR NORM INTO ERROR. DO 3600 I=1,NUMGR ERROR(I)=V(I,N) 3600 CONTINUE ERROR(NUMGR+1)=DELTK C 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 ACCEPTABLE 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)3700,3800,3800 3700 V(2*NUMGR+1,NP1)=-V(2*NUMGR+1,NP1)-TWO C SET ICON0=1 TO INDICATE THIS IS NOT AN ACCEPTABLE C APPROXIMATION. ICON0=1 GO TO 3900 C C HERE DELTK .GE. 0.0. MAKE ICON0=0 IF THE CONSTRAINTS ARE C SATISFIED WITHIN EPRES AND ICON0=1 OTHERWISE. 3800 ICON0=ICNCK(IOPT,NUMP,NUMQ,NUMGR,INDLW,ALTBL,INDUP, *UPTBL,EPSQ,X,V,MEQ,MINEQ,ALC,EPRES) C C END OF INITIALIZATION PHASE. C C C THE NEXT GROUP OF STATEMENTS (THROUGH STATEMENT 5400) ARE C FOR THE FIRST MAIN DIFFERENTIAL CORRECTION ITERATION. C IF THIS ITERATION FAILS TO PRODUCE AN ACCEPTABLE C APPROXIMATION AND THE INITIAL APPROXIMATION IS ALSO C UNACCEPTABLE WE QUIT, OTHERWISE WE TRY THE NEXT C INITIALIZATION STAGE. C C SET UP THE LINEAR PROGRAMMING PROBLEM. 3900 CALL SSETV(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,INDLW, *ALTBL,INDUP,UPTBL,INDF,WTBLE,MEQ,MINEQ,ALC,1,V,M) C SET IYRCT(1)=-1 AS A SIGNAL TO DO THE INITIAL EXCHANGES IN C SLNPRO STRICTLY ACCORDING TO A PIVOTING STRATEGY RATHER C THAN ATTEMPTING TO USE INFORMATION ABOUT A PREVIOUS VERTEX. IYRCT(1)=-1 C REDUCE V IF MEQA .GT. 0 IF(MEQA)4100,4100,4000 4000 CALL SCRCH(M,NMPPQ,MEQA,NSCOL,IRESL,EQ,V) C SET M=-M AS A SIGNAL TO LOOSEN THE TOLERANCES IN SLNPRO IF C IADJS=2 AND LEAVE M UNCHANGED IF IADJS=0 OR 1. SLNPRO WILL C RESET M TO ITS ORIGINAL VALUE IF IT IS CHANGED HERE. 4100 M=M*(1+IADJS-IADJS**2) C SET NLP=-NLP AS A SIGNAL TO TIGHTEN THE TOLERANCES IN SLNPRO IF C IADJS=1 AND LEAVE NLP UNCHANGED IF IADJS=0 OR 2. SLNPRO WILL C RESET NLP TO ITS ORIGINAL VALUE IF IT IS CHANGED HERE. NLP=NLP*(1-4*IADJS+2*IADJS**2) CALL SLNPRO(V,M,NLP,IYRCT,XTEMP,INDIC) C IF THIS WAS THE FIRST SLNPRO CALL IN THIS ITERATION, SO THE C TOLERANCES IN SLNPRO WERE NOT ADJUSTED BY CDFCOR, FOR C POSSIBLE OUTPUT PURPOSES WE SAVE INDIC IN ILST1. IF(IADJS)4200,4200,4300 4200 ILST1=INDIC C IF INDIC IS POSITIVE SOMETHING WENT WRONG IN SLNPRO SO WE C TRY AGAIN WITH BOTH TIGHTER AND LOOSER TOLERANCES IF THIS C HAS NOT ALREADY BEEN TRIED. 4300 IF(INDIC)4900,4900,4400 C C HERE WE DO NOT HAVE AN ACCEPTABLE APPROXIMATION. 4400 IF(IADJS-1)4500,4500,4600 4500 IADJS=IADJS+1 C CALL SPQERD WITH THE COEFFICIENTS FROM THE INITIALIZATION C (WHICH ARE IN X) TO RECOMPUTE THE VALUES OF P AND Q AND C PUT THEM IN V FOR USE BY SSETV. CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X,V, *QMIN,IQMIN) GO TO 3900 C C HERE WE HAVE TRIED AND FAILED WITH NORMAL, TIGHTENED, AND C LOOSENED TOLEARANCES SO THE FIRST MAIN DIFFERENTIAL C CORRECTION ITERATION HAS FAILED TO PRODUCE AN ACCEPTABLE C APPROXIMATION. IF THE APPROXIMATION FROM THE INITIALIZATION C PHASE WAS ACCEPTABLE, WE SET IFLAG AND INDLP, AND RETURN. 4600 IF(ICON0)8500,8500,4700 C HERE THE INITIAL APPROXIMATION WAS UNACCEPTABLE. IF IT WAS C USER SUPPLIED WE RESET IADJS AND TRY LOEB, IF IT CAME C FROM LOEB WE TRY P/Q = 0.0/1.0, AND IF IT WAS P/Q = 0.0/1.0 C WE SET IFLAG AND INDLP, AND RETURN. 4700 IF(INITL)4800,3000,8500 4800 IADJS=0 GO TO 2200 C END OF BLOCK WITHIN FIRST MAIN DIFFERENTIAL CORRECTION C ITERATION TO DEAL WITH AN UNACCEPTABLE APPROXIMATION. C C HERE INDIC WAS NONPOSITIVE, SO SLNPRO PLACED AN APPROXIMATION C IN XTEMP, WHICH WE NOW CHECK FOR ACCEPTABLENESS. FIRST C EXPAND XTEMP IF MEQA .GT. O 4900 IF(MEQA)5100,5100,5000 5000 CALL EXPND(NMPPQ,MEQA,NSCOL,IRESL,EQ,XTEMP) C CALL SPQERD TO COMPUTE P, Q, THE ERRORS, AND THE C UNIFORM ERROR NORM (AND TO PLACE THESE VALUES IN V FOR C USE BY OTHER ROUTINES). 5100 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE, *XTEMP,V,QMIN,IQMIN) DELTK=V(2*NUMGR+1,NP1) C IF DELTK IS NEGATIVE (INDICATING A VERY SMALL OR NEGATIVE C DENOMINATOR VALUE AT SOME POINT) WE REJECT THIS C APPROXIMATION AND TRY AGAIN. IF(DELTK)4400,5200,5200 C HERE DELTK IS NONNEGATIVE, AND WE NORMALIZE XTEMP AND CHECK C THE CONSTRAINTS. 5200 CALL NRMLZ(NUMP,NUMQ,MEQ,MINEQ,NSCOL,IRESL,ALC,TOLEL,XTEMP) ICON1=ICNCK(IOPT,NUMP,NUMQ,NUMGR,INDLW,ALTBL,INDUP, *UPTBL,EPSQ,XTEMP,V,MEQ,MINEQ,ALC,EPRES) C IF THE APPROXIMATION IS UNACCEPTABLE WE TRY AGAIN. IF(ICON1)5300,5300,4400 C HERE THE APPROXIMATION COMPUTED BY THE FIRST MAIN CDFCOR C ITERATION SATISFIES THE CONSTRAINTS WITHIN EPRES, AND WE WILL C ACCEPT IT IF THE APPROXIMATION FROM THE INITIALIZATION WAS C UNACCEPTABLE OR DOES NOT HAVE A SMALLER ERROR NORM. C OTHERWISE WE WILL REJECT IT AND TRY AGAIN SINCE SOMETHING C IS FISHY. 5300 IF(ICON0)5400,5400,5500 5400 IF(DELTK-DLOLD)5500,5500,4400 C C END OF GROUP OF STATEMENTS FOR FIRST MAIN DIFFERENTIAL C CORRECTION ITERATION. FROM NOW ON ONLY STATEMENTS BELOW C THIS POINT IN CDFCOR WILL BE EXECUTED. C C C HERE THE PRESENT APPROXIMATION WILL BE KEPT, AND WE PREPARE C TO COMPUTE ANOTHER ONE. C COPY XTEMP INTO X. 5500 DO 5600 I=1,NMPPQ X(I)=XTEMP(I) 5600 CONTINUE C PUT THE ERRORS AND ERROR NORM INTO ERROR. DO 5700 I=1,NUMGR ERROR(I)=V(I,N) 5700 CONTINUE ERROR(NUMGR+1)=DELTK ITER=ITER+1 IF(ITER-ITLIM)5800,6000,6000 5800 DLOLD=DELTK C SAVE INDIC IN CASE THE APPROXIMATION JUST COMPUTED IS THE ONE C RETURNED. IRET1=INDIC C RESET IADJS. 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)6300,6300,5900 5900 V(2*NUMGR+1,NP1)=V(2*NUMGR+1,NP1)/TWO GO TO 6300 C C HERE ITLIM ITERATIONS HAVE BEEN PERFORMED AND WE WILL RETURN C AFTER SETTING IFLAG AND INDLP (SEE BELOW). 6000 INDLP=5-INDIC IF(INDIC)6100,6200,6200 6100 IFLAG=-12 GO TO 9400 6200 IFLAG=-2 GO TO 9400 C C SAVE IYRCT IN IYSAV SO WE CAN LATER USE INFORMATION FROM THIS C SUCCESSFUL SLNPRO ITERATION TO RESTART SLNPRO. 6300 DO 6400 I=1,M IYSAV(I)=IYRCT(I) 6400 CONTINUE C CALL SSETV TO SET UP FOR THE NEXT ITERATION. 6500 CALL SSETV(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ,INDLW, *ALTBL,INDUP,UPTBL,INDF,WTBLE,MEQ,MINEQ,ALC,1,V,M) C REDUCE V IF MEQA .GT. 0 IF(MEQA)6700,6700,6600 6600 CALL SCRCH(M,NMPPQ,MEQA,NSCOL,IRESL,EQ,V) C IF IADJS=1, SET NLP=-NLP (AND LEAVE M ALONE) AS A SIGNAL TO C TIGHTEN THE TOLERANCES IN SLNPRO, AND IF IADJS=2 SET M=-M C (AND LEAVE NLP ALONE) AS A SIGNAL TO LOOSEN THE TOLERANCES C IN SLNPRO. IF IADJS=0, LEAVE THE TOLERANCES ALONE. M AND C NLP WILL BE RESET TO THEIR ORIGINAL VALUES BY SLNPRO IF THEY C ARE CHANGED HERE. 6700 M=M*(1+IADJS-IADJS**2) NLP=NLP*(1-4*IADJS+2*IADJS**2) CALL SLNPRO(V,M,NLP,IYRCT,XTEMP,INDIC) C IF THIS WAS THE FIRST SLNPRO CALL IN THIS ITERATION, SO THE C TOLERANCES IN SLNPRO WERE NOT ADJUSTED BY CDFCOR, FOR C POSSIBLE OUTPUT PURPOSES WE SAVE INDIC IN ILST1. IF(IADJS)6800,6800,6900 6800 ILST1=INDIC C IF INDIC IS POSITIVE SOMETHING WENT WRONG IN SLNPRO SO WE GO C INTO THE BLOCK OF STATEMENTS DESIGNED TO HANDLE FAILURE TO C PRODUCE AN ACCEPTABLE APPROXIMATION. 6900 IF(INDIC)7600,7600,7000 C C HERE WE DID NOT GET AN ACCEPTABLE APPROXIMATION. IF WE ARE C STILL DOING PARAMETRIC BISECTION, TO SAVE TIME WE DO NOT TRY C SLNPRO AGAIN WITH ADJUSTED TOLERANCES SINCE THE PARAMETRIC C BISECTION PROBABLY CAUSED THE FAILURE, BUT RATHER WE SWITCH C PERMANENTLY INTO ORDINARY DIFFERENTIAL CORRECTION. 7000 IF(IPBIT)7200,7200,7100 7100 IPBIT=0 GO TO 7400 C IF WE WERE NOT DOING PARAMETRIC BISECTION WE TRY AGAIN WITH C TIGHTER AND LOOSER TOLERANCES, IF THIS HAS NOT ALREADY BEEN C TRIED. IF IT HAS BEEN TRIED, WE SET IFLAG AND INDLP AND RETURN. 7200 IF(IADJS-1)7300,7300,8500 7300 IADJS=IADJS+1 C CALL SPQERD WITH THE COEFFICIENTS FROM THE PREVIOUS ITERATION C (WHICH ARE IN X) TO RECOMPUTE THE VALUES OF P AND Q AND PUT C THEM IN V FOR USE BY SSETV. 7400 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X,V, *QMIN,IQMIN) C RESTORE IYRCT TO ITS VALUE AT THE PREVIOUS ITERATION. DO 7500 I=1,M IYRCT(I)=IYSAV(I) 7500 CONTINUE GO TO 6500 C END OF BLOCK OF STATEMENTS WITHIN THE SECOND OR LATER MAIN C DIFFERENTIAL CORRECTION ITERATION TO HANDLE FAILURE TO C PRODUCE AN ACCEPTABLE APPROXIMATION. C C HERE INDIC WAS NONPOSITIVE, SO SLNPRO PRODUCED AN C APPROXIMATION IN XTEMP WHICH WE NOW CHECK FOR C ACCEPTABILITY. FIRST EXPAND XTEMP IF MEQA .GT. 0 7600 IF(MEQA)7800,7800,7700 7700 CALL EXPND(NMPPQ,MEQA,NSCOL,IRESL,EQ,XTEMP) C CALL SPQERD TO COMPUTE P, Q, THE ERRORS, AND THE UNIFORM C ERROR NORM (AND PLACE THEM IN V FOR USE BY OTHER ROUTINES). 7800 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE, *XTEMP,V,QMIN,IQMIN) DELTK=V(2*NUMGR+1,NP1) C IF DELTK IS NEGATIVE (INDICATING A VERY SMALL OR NEGATIVE C DENOMINATOR VALUE AT SOME POINT) WE REJECT THIS C APPROXIMATION AND TRY AGAIN. IF(DELTK)7000,7900,7900 C IF INDIC IS NEGATIVE HERE IT WAS NECESSARY C TO ADJUST TOLERANCES IN SLNPRO, SO THE APPROXIMATION C COMPUTED IS MORE LIKELY TO HAVE BEEN AFFECTED C BY ROUNDOFF ERROR. THUS IN THIS CASE WE CHECK THE C CONSTRAINTS, AND REJECT THE APPROXIMATION IF ANY OF THEM ARE C VIOLATED BY MORE THAN EPRES. NOTE THAT INDIC CANNOT BE C POSITIVE HERE SINCE SLNPRO PRODUCED AN APPROXIMATION. 7900 IF(INDIC)8000,8100,8100 8000 ICON=ICNCK(IOPT,NUMP,NUMQ,NUMGR,INDLW,ALTBL,INDUP, *UPTBL,EPSQ,XTEMP,V,MEQ,MINEQ,ALC,EPRES) IF(ICON)8100,8100,7000 C REJECT THIS APPROXIMATION AND TRY AGAIN IF THE ERROR NORM HAS C NOT DECREASED. 8100 IF(DELTK-DLOLD)8200,7000,7000 C RARELY AN APPROXIMATION WITH ALL COEFFICIENTS VERY SMALL C WILL BE PRODUCED. TO BE SAFE, IF ALL DENOMINATOR C COEFFICIENTS ARE .LE. 0.1 IN ABSOLUTE VALUE WE ATTEMPT TO C NORMALIZE THE APPROXIMATION AND CHECK THE CONSTRAINTS. IF C THEY ARE SATISFIED WITHIN EPRES WE ACCEPT THIS C APPROXIMATION, AND OTHERWISE WE TRY AGAIN. 8200 DO 8400 J=NMPP1,NMPPQ IF(XTEMP(J)-TENTH)8300,8300,5500 8300 IF(XTEMP(J)+TENTH)5500,8400,8400 8400 CONTINUE CALL NRMLZ(NUMP,NUMQ,MEQ,MINEQ,NSCOL,IRESL,ALC,TOLEL, *XTEMP) ICON=ICNCK(IOPT,NUMP,NUMQ,NUMGR,INDLW,ALTBL,INDUP, *UPTBL,EPSQ,XTEMP,V,MEQ,MINEQ,ALC,EPRES) IF(ICON)5500,5500,7000 C END OF COMPUTATION OF APPROXIMATIONS. C 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 CDFCOR 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 CDFCOR 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. 8500 INDLP=10*(5-ILST1)+5-IRET1 C C WE NOW SET UP IFLAG FOR THE INFORMATION OF THE USER. IF(ITER)8600,8600,9200 C C HERE ITER=0, SO THE INITIAL APPROXIMATION WAS NOT IMPROVED. 8600 IF(INITL)9100,9100,8700 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, OR IF IT VIOLATED A CONSTRAINT (EITHER OF C WHICH WOULD BE INDICATED BY ICON0=1), THEN SET IFLAG=6666, C ITER=-1, AND ZERO OUT THE ERRORS AND COEFFICIENTS AND RETURN. C 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. 8700 INDLP=10*(5-ILST1)+5-ILBFL IF(ICON0)9000,9000,8800 8800 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)=ZERO DO 8900 I=1,NUMGR ERROR(I)=ZERO 8900 CONTINUE RETURN 9000 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. 9100 IFLAG=-1 C C IF IRET1 .LT. 0 SET THE TENS DIGIT OF IFLAG TO INDICATE C THAT THE TOLERANCES IN SLNPRO HAD TO BE ADJUSTED DURING THE C COMPUTATION OF THE APPROXIMATION RETURNED. 9200 IF(IRET1)9300,9400,9400 9300 IFLAG=IFLAG-10 C NOW WE CHECK THE USER SUPPLIED CONSTRAINTS (IF ANY). C CALL SPQERD TO RECOMPUTE THE VALUES OF P AND Q. QMIN AND C IQMIN ARE ALSO RECOMPUTED BY SPQERD. 9400 CALL SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X, *V,QMIN,IQMIN) ICON=ICNCK(IOPT,NUMP,NUMQ,NUMGR,INDLW,ALTBL,INDUP, *UPTBL,EPSQ,X,V,MEQ,MINEQ,ALC,EPRES) IF(ICON)9600,9600,9500 9500 IFLAG=IFLAG-100 9600 RETURN END SUBROUTINE STCOMP(NUMDM,NRWGR,NCLGR,ASWP1,ASWP2, *ANEP1,ANEP2,T) C***BEGIN PROLOGUE STCOMP C***DATE WRITTEN 720120 (YYMMDD) C***REVISION DATE 880226 (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 C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION T(202) C C SET PRECISION DEPENDENT CONSTANTS FOR STCOMP. C***FIRST EXECUTABLE STATEMENT STCOMP ONE=1.0D0 ZERO=ONE-ONE C END SETTING PRECISION DEPENDENT CONSTANTS FOR STCOMP. C 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=ZERO 6 IF(NRWGR-1)7,8,7 7 VSPAC=(ANEP2-ASWP2)/FLOAT(NRWGR-1) GO TO 9 8 VSPAC=ZERO 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 SPQERD(NUMP,NUMQ,NUMGR,FTBWT,PWR,INDF,WTBLE,X, *V,QMIN,IQMIN) C***BEGIN PROLOGUE SPQERD C***REFER TO CDFCOR 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 C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION FTBWT(101),PWR(101,15),V(635,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=D1MACH(3) C SET PRECISION DEPENDENT CONSTANTS FOR SPQERD. ONE=1.0D0 ZERO=ONE-ONE TWO=ONE+ONE FOUR=TWO+TWO TEN=FOUR+FOUR+TWO 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=TEN*TEN*SPCMN C END OF SETTING MACHINE AND PRECISION 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)=ZERO ABAD=ONE 1 I=I+1 IF(I-NUMGR)3,3,1001 C 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 C 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)=ZERO DO 4 J=1,NUMP V(IDBM,NP1)=V(IDBM,NP1)+X(J)*PWR(I,J) 4 CONTINUE C C COMPUTE THE VALUE OF Q AT GRID POINT I AND PUT IT IN C V(2*I,N+1). V(IDB,NP1)=ZERO 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 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)=ZERO 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=-ONE V(I,N)=ZERO 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)=TEN*EPS2 GO TO 1 END SUBROUTINE SLNPRO(V,M,N,IYRCT,X,INDIC) C***BEGIN PROLOGUE SLNPRO C***REFER TO CDFCOR 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 C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION V(635,17),IYRCT(634),X(16),Y(634), * IXRCT(634),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 DEPENDENT 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=D1MACH(3) C SET PRECISION DEPENDENT CONSTANTS FOR SLNPRO. ONE=1.0D0 ZERO=ONE-ONE TWO=ONE+ONE FOUR=TWO+TWO TEN=FOUR+TWO+TWO 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=TEN*TEN*SPCMN IF(REA-TEN**(-8))10000,10010,10010 10000 REA=TEN**(-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=TEN*SPCMN C END OF INITIAL SETTING OF MACHINE DEPENDENT PARAMETERS FOR C SLNPRO. THESE PARAMETERS MAY BE ADJUSTED BY A COMMAND FROM C CDFCOR. 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/(TEN*TEN) 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)=ZERO 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=ONE AMPRV=ONE AMPR2=ONE 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=ZERO 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=ONE 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=ONE 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=ONE 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-ONE/TWO)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=ONE 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=-ONE 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+ONE/TWO)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)=ZERO 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 C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION V(635,17) C C SET PRECISION DEPENDENT CONSTANTS FOR SJELIM. C***FIRST EXECUTABLE STATEMENT SJELIM ONE=1.0D0 C EMD OF SETTING PRECISION DEPENDENT CONSTANTS FOR SJELIM. C C DIVIDE THE ENTRIES IN THE RESOLVENT ROW (EXCEPT FOR THE C RESOLVENT) BY THE RESOLVENT. 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)=ONE/RESOL RETURN END SUBROUTINE SSETV(NUMP,NUMQ,NUMGR,FTBWT,PWR,IOPT,EPSQ, *INDLW,ALTBL,INDUP,UPTBL,INDF,WTBLE,MEQ,MINEQ,ALC,MODE,V,M) C***BEGIN PROLOGUE SSETV C***REFER TO CDFCOR C***ROUTINES CALLED (NONE) C***PURPOSE THIS SUBROUTINE SETS UP THE LINEAR PROGRAMMING PROBLEM C FOR THE INITIAL (IF MODE=0) OR MAIN (IF MODE=1) C ITERATIONS OF THE DIFFERENTIAL CORRECTION ALGORITHM. C***END PROLOGUE SSETV2 C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION FTBWT(101),PWR(101,15),V(635,17),WTBLE(101), * UPTBL(101),ALTBL(101),INDUP(101),INDLW(101), * INDF(101),ALC(101,16) C C IN THE INITIALIZATION WE MINIMIZE W WITH THE RESTRICTIONS C ABS(F*Q-WT*P) .LE. W AT ALL GRID POINTS WHERE F IS TO BE C APPROXIMATED, C Q .GE. EPS AT ALL GRID POINTS, C ABS(Q(J)) .LE. 1.0 FOR ALL J, C PLUS THE OTHER CONSTRAINTS (IF ANY). C IN THE MAIN ITERTATIONS WE MINIMIZE W WITH THE RESTRICTIONS C (ABS(F*Q-WT*P)-DELTK*Q)/QK .LE. W AT ALL GRID POINTS WHERE C F IS TO BE APPROXIMATED, C ABS(Q(J)) .LE. 1.0 FOR ALL J, C PLUS THE OTHER CONSTRAINTS (IF ANY). C EACH ABSOLUTE VALUE CONSTRAINT IS BROKEN INTO TWO CONSTRAINTS C TO ACCOUNT FOR THE ABSOLUTE VALUE. IN THE MAIN ITERATIONS C THE SUBROUTINE USES THE FACT THAT THE VALUE OF QK AT GRID C POINT I IS STORED IN V(2*I,N+1) AND DELTK IS IN C V(2*NUMGR+1,N+1) FROM A PREVIOUS CALL TO SPQERD. C C***FIRST EXECUTABLE STATEMENT SSETV2 IONE=IOPT-(IOPT/10)*10 C C SET MACHINE AND PRECISION DEPENDENT CONSTANTS FOR C SUBROUTINE SSETV. ONE=1.0D0 ZERO=ONE-ONE TWO=ONE+ONE FOUR=TWO+TWO TEN=FOUR+FOUR+TWO C IF WE ARE IN THE MAIN PART OF DIFFERENTIAL CORRECTION AND C IONE IS ODD WE SET EPS=EPSQ. IN THIS CASE WE WILL REQUIRE C THAT THE DENOMINATOR BE .GE. EPS AT EVERY GRID POINT. IF(MODE)300,300,100 100 IF(IONE-(IONE/2)*2)800,800,200 200 EPS=EPSQ GO TO 800 C HERE WE ARE IN THE INITIALIZATION AND WE AT LEAST TEMPORARILY C COMPUTE EPS=MAX(0.0001,1000.0*SPCMN) 300 SPCMN=D1MACH(3) EPS=SPCMN*TEN**3 IF(EPS-TEN**(-4))400,500,500 400 EPS=TEN**(-4) C IF IONE IS ODD HERE WE INCREASE EPS TO EPSQ IF NECESSARY. 500 IF(IONE-(IONE/2)*2)800,800,600 600 IF(EPS-EPSQ)700,800,800 700 EPS=EPSQ C END OF SETTING MACHINE AND PRECISION DEPENDENT CONSTANTS C FOR SSETV. C 800 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 SSETV 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)1400,900,900 900 DO 1300 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))1300,1300,1000 1000 M=M+1 DO 1100 J=1,NUMP V(M,J)=PWR(I,J) 1100 CONTINUE DO 1200 J=NMPP1,NMPPQ V(M,J)=-UPTBL(I)*PWR(I,J) 1200 CONTINUE V(M,N)=ZERO 1300 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 1400 IF(IOPT-(IOPT/100)*100-10)2000,1500,1500 1500 DO 1900 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))1900,1900,1600 1600 M=M+1 DO 1700 J=1,NUMP V(M,J)=-PWR(I,J) 1700 CONTINUE DO 1800 J=NMPP1,NMPPQ V(M,J)=ALTBL(I)*PWR(I,J) 1800 CONTINUE V(M,N)=ZERO 1900 CONTINUE C 2000 IF(MODE)2100,2100,2700 C C HERE WE ARE IN THE INITIALIZATION PHASE, AND FOR EACH GRID C POINT AT WHICH F IS TO BE APPROXIMATED WE SET THE CONSTRAINTS C -WT*P+F*Q-W .LE. 0.0 AND WT*P-F*Q-W .LE. 0.0 2100 DO 2500 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))2500,2500,2200 2200 M=M+2 MM1=M-1 DO 2300 J=1,NUMP V(MM1,J)=-WTBLE(I)*PWR(I,J) V(M,J)=-V(MM1,J) 2300 CONTINUE DO 2400 J=NMPP1,NMPPQ V(MM1,J)=FTBWT(I)*PWR(I,J) V(M,J)=-V(MM1,J) 2400 CONTINUE V(MM1,N)=-ONE V(M,N)=-ONE 2500 CONTINUE C SINCE WE ARE IN THE INITIALIZATION, COLUMN N+1 IS NOT NEEDED C FOR STORAGE OF QK, SO ZERO IT OUT FOR THOSE ROWS ALREADY SET. DO 2600 I=1,M V(I,NP1)=ZERO 2600 CONTINUE GO TO 3300 C C HERE WE ARE IN THE MAIN DIFFERENTIAL CORRECTION PHASE, AND 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). 2700 DO 3100 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))3100,3100,2800 2800 M=M+2 MM1=M-1 DO 2900 J=1,NUMP V(MM1,J)=-WTBLE(I)*PWR(I,J) V(M,J)=-V(MM1,J) 2900 CONTINUE DEL1=-DELTK+FTBWT(I) DEL2=-DELTK-FTBWT(I) DO 3000 J=NMPP1,NMPPQ V(MM1,J)=DEL1*PWR(I,J) V(M,J)=DEL2*PWR(I,J) 3000 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) 3100 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 3200 I=1,M V(I,NP1)=ZERO 3200 CONTINUE C IF THE ONES DIGIT OF IOPT IS ODD HERE WE FORCE Q .GE. EPS. IF(IONE-(IONE/2)*2)3700,3700,3300 C C HERE WE FORCE Q .GE. EPS. 3300 DO 3600 I=1,NUMGR M=M+1 DO 3400 J=1,N V(M,J)=ZERO 3400 CONTINUE DO 3500 J=NMPP1,NMPPQ V(M,J)=-PWR(I,J) 3500 CONTINUE V(M,NP1)=-EPS 3600 CONTINUE C C NOW INCLUDE THE ADDITIONAL INEQUALITY CONSTRAINTS FROM ALC, C IF THERE ARE ANY. 3700 IF(MINEQ)4100,4100,3800 3800 IALC=MEQ DO 4000 I=1,MINEQ M=M+1 IALC=IALC+1 DO 3900 J=1,NMPPQ V(M,J)=ALC(IALC,J) 3900 CONTINUE V(M,N)=ZERO V(M,NP1)=ALC(IALC,N) 4000 CONTINUE C C NOW SET THE CONSTRAINTS OF THE FORM -Q(J).LE.1.0 C AND Q(J).LE.1.0 4100 DO 4300 I=1,NUMQ M=M+2 MM1=M-1 DO 4200 J=1,N V(MM1,J)=ZERO V(M,J)=ZERO 4200 CONTINUE NNN=NUMP+I V(MM1,NNN)=-ONE V(M,NNN)=ONE V(MM1,NP1)=ONE V(M,NP1)=ONE 4300 CONTINUE C C NOW SET THE BOTTOM ROW. TO MINIMIZE W WE MAXIMIZE -W. MP1=M+1 DO 4400 J=1,NP1 V(MP1,J)=ZERO 4400 CONTINUE V(MP1,N)=ONE RETURN END SUBROUTINE EQELM(NVAR,MEQ,TOLEL,EQ,NSCOL,IRESL,MEQA,IFLAG) C***BEGIN PROLOGUE EQELM C***REFER TO CDFCOR C***ROUTINES CALLED (NONE) C***PURPOSE THIS SUBROUTINE ATTEMPTS TO SOLVE FOR SOME OF THE C VARIABLES IN THE MEQ EQUALITY CONSTRAINTS IN TERMS C OF THE OTHER VARIABLES. C***END PROLOGUE EQELM C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION EQ(15,16),NSCOL(15),IRESL(30) C C GIVEN THE EQUATIONS C EQ(I,1)*X(1)+...+EQ(I,NVAR)*X(NVAR) = EQ(I,NVAR+1) C FOR I=1,...,MEQ, WE CHOOSE RESOLVENTS BY TOTAL PIVOTING IN THE C LEFT MEQ BY NVAR PORTION OF EQ, STORE THE LTH RESOLVENT C POSITION IN (IRESL(2*L-1),IRESL(2*L)), AND APPLY ROW C OPERATIONS TO MAKE 1.0 IN THE RESOLVENT POSITION AND 0.0 C ELSEWHERE IN THE RESOLVENT COLUMN. ONLY RESOLVENTS WITH ABSOLUTE C VALUE .GT. TOLEL ARE ACCEPTED. MEQA (WHICH STANDS FOR MEQ ACTUAL) C WILL BE THE NUMBER OF ACCEPTABLE RESOLVENTS FOUND. IF C MEQA .LT. MEQ, AND IF SOME ROW NOT CONTAINING A RESOLVENT C HAS LAST ELEMENT .GT. TOLEL IN ABSOLUTE VALUE, WE SET C IFLAG=7777 AS AN INCONSISTENCY WARNING AND RETURN. IF C MEQA .LT. MEQ BUT THERE IS NO INCONSISTENCY THEN THE MEQ-MEQA C ROWS OF EQ NOT CONTAINING RERSOLVENTS WILL BE IGNORED (THIS IS C THE REDUNDANT CONSTRAINTS SITUATION). C NSCOL (WHICH STANDS FOR NON SOLVED FOR COLUMN) IS DEFINED BY C NSCOL(K) = THE SUBSCRIPT OF THE KTH UNSOLVED FOR VARIABLE C (THAT IS, VARIABLE WITH NO RESOLVENT IN ITS COLUMN) FOR K=1,..., C NVAR-MEQA. THUS FOR K=1,...,MEQA WE WILL HAVE C X(IRESL(2*K)) = -EQ(IRESL(2*K-1),NSCOL(1))*X(NSCOL(1))-...- C EQ(IRESL(2*K-1),NSCOL(NVAR-MEQA))*X(NSCOL(NVAR-MEQA))+ C EQ(IRESL(2*K-1),NVAR+1). C C SET PRECISION DEPENDENT CONSTANTS FOR EQELM. C***FIRST EXECUTABLE STATEMENT EQELM ONE=1.0D0 ZERO=ONE-ONE C END OF SETTING PRECISION DEPENDENT CONSTANTS FOR EQELM. C N=NVAR+1 NRESL=0 MEQA=MEQ C IF MEQ IS NOT POSITIVE THERE IS NOTHING TO DO HERE AND WE RETURN. IF(MEQ)30,30,70 30 RETURN C ZERO OUT NSCOL. LATER EACH ELEMENT OF NSCOL CORRESPONDING TO A C COLUMN CONTAINING A RESOLVENT WILL BE SET TO 1, AND NSCOL WILL C BE READJUSTED TO ITS FINAL STATE AT THE END OF THIS SUBROUTINE C IF THE EQUATIONS IN EQ ARE CONSISTENT. 70 DO 100 J=1,NVAR NSCOL(J)=0 100 CONTINUE C C WE ATTEMPT TO FIND MEQ RESOLVENTS USING TOTAL PIVOTING. DO 1400 L=1,MEQ AMAX=ZERO DO 800 I=1,MEQ DO 700 J=1,NVAR IF(NRESL)500,500,200 C WHILE SEARCHING, IGNORE ROWS AND COLUMNS CONTAINING AN C EARLIER RESOLVENT. 200 DO 400 K=1,NRESL IF(I-IRESL(2*K-1))300,800,300 300 IF(J-IRESL(2*K))400,700,400 400 CONTINUE C HERE LOCATION (I,J) IS NOT IN A ROW OR COLUMN WITH AN C EARLIER RESOLVENT. 500 IF(ABS(EQ(I,J))-AMAX)700,700,600 600 AMAX=ABS(EQ(I,J)) IMAX=I JMAX=J 700 CONTINUE 800 CONTINUE IF(AMAX-TOLEL)1500,1500,900 C C HERE WE HAVE FOUND AN ACCEPTABLE RESOLVENT. 900 NRESL=NRESL+1 IRESL(2*NRESL-1)=IMAX IRESL(2*NRESL)=JMAX NSCOL(JMAX)=1 C NOW DO ROW OPERATIONS TO MAKE 1.0 IN THE (IMAX,JMAX) POSITION C AND 0.0 ELSEWHERE IN COLUMN JMAX (THUS SOLVING FOR X(JMAX) C IN TERMS OF THE OTHER VARIABLES). C FIRST DIVIDE ROW IMAX BY THE RESOLVENT. RESOL=EQ(IMAX,JMAX) DO 1000 J=1,N EQ(IMAX,J)=EQ(IMAX,J)/RESOL 1000 CONTINUE DO 1300 I=1,MEQ IF(I-IMAX)1100,1300,1100 1100 FACT=EQ(I,JMAX) DO 1200 J=1,N EQ(I,J)=EQ(I,J)-FACT*EQ(IMAX,J) 1200 CONTINUE 1300 CONTINUE 1400 CONTINUE GO TO 2100 C C HERE WE WERE UNABLE TO FIND AN ACCEPTABLE RESOLVENT. C RESET MEQA TO THE NUMBER OF RESOLVENTS FOUND AND CHECK TO C SEE IF ANY OF THE RIGHT SIDES OF ROWS NOT CONTAINING A C RESOLVENT ARE .GT. TOLEL IN ABSOLUTE VALUE. IF SO, THE C EQUATIONS IN EQ WERE APPARENTLY INCONSISTENT, AND IF NOT, C THEY WERE APPARENTLY REDUNDANT. 1500 MEQA=NRESL DO 1900 I=1,MEQ IF(NRESL)1800,1800,1600 1600 DO 1700 K=1,NRESL IF(I-IRESL(2*K-1))1700,1900,1700 1700 CONTINUE C HERE ROW I DOES NOT CONTAIN A RESOLVENT. 1800 IF(ABS(EQ(I,N))-TOLEL)1900,1900,2000 1900 CONTINUE GO TO 2100 C C HERE THE CONSTRAINTS WERE INCONSISTENT. 2000 IFLAG=7777 RETURN C C HERE THE CONSTRAINTS WERE CONSISTENT, AND WE SET NSCOL AND C RETURN. 2100 IARG=0 DO 2300 J=1,NVAR IF(NSCOL(J))2200,2200,2300 2200 IARG=IARG+1 NSCOL(IARG)=J 2300 CONTINUE RETURN END SUBROUTINE SCRCH(M,NVAR,MEQA,NSCOL,IRESL,EQ,V) C***BEGIN PROLOGUE SCRCH C***REFER TO CDFCOR C***ROUTINES CALLED (NONE) C***PURPOSE THIS SUBROUTINE USES THE NONREDUNDANT EQUALITY C CONSTRAINTS TO ELIMINATE MEQA VARIABLES, THUS C SCRUNCHING V FROM MP1 BY NP1 TO MP1 BY NP1-MEQA. C***END PROLOGUE SCRCH C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION NSCOL(15),IRESL(30),EQ(15,16),V(635,17) C C IF MEQA=0 THERE IS NO SCRUNCHING TO BE DONE, AND WE RETURN. C***FIRST EXECUTABLE STATEMENT SCRCH IF(MEQA)100,100,200 100 RETURN C 200 N=NVAR+1 NP1=N+1 MP1=M+1 NS=NVAR-MEQA C DO 800 L=1,MEQA C THE RESOLVENT IS AT POSITION (LROW,LCOL) OF EQ. LROW=IRESL(2*L-1) LCOL=IRESL(2*L) C IF ALL THE VARIABLES X(1),...,X(NVAR) ARE SOLVED FOR WE GO C TO COLUMN NP1. IF(NS)600,600,300 300 DO 500 J=1,NS JCOL=NSCOL(J) FACT=EQ(LROW,JCOL) C WE REPLACE COLUMN JCOL OF V BY ITSELF MINUS EQ(LROW,JCOL) C TIMES COLUMN LCOL OF V. DO 400 I=1,MP1 V(I,JCOL)=V(I,JCOL)-FACT*V(I,LCOL) 400 CONTINUE 500 CONTINUE C C NOW DO COLUMN NP1 OF V. 600 FACT=EQ(LROW,N) DO 700 I=1,MP1 V(I,NP1)=V(I,NP1)-FACT*V(I,LCOL) 700 CONTINUE 800 CONTINUE C C NOW COUMNS IRESL(2),...,IRESL(2*MEQA) OF V ARE NO LONGER C NEEDED, AND WE COMPACTIFY V. IF(NS)1300,1300,900 900 DO 1200 J=1,NS JCOL=NSCOL(J) C COLUMN JCOL OF V IS NOW MOVED TO THE COLUMN J POSITION IF IT C IS NOT ALREADY THERE. IF(JCOL-J)1000,1200,1000 1000 DO 1100 I=1,MP1 V(I,J)=V(I,JCOL) 1100 CONTINUE 1200 CONTINUE C NOW MOVE THE LAST TWO COLUMNS OF V TO THE LEFT MEQA SPACES. C NOTE THAT THE VARIABLE W CORRESPONDING TO COLUMN N = NVAR+1 C OF V IS AUTOMATICALLY AN UNSOLVED FOR VARIABLE SINCE IT IS C NOT INCLUDED IN EQ. 1300 NMEQ=N-MEQA NPMEQ=NMEQ+1 DO 1400 I=1,MP1 V(I,NMEQ)=V(I,N) V(I,NPMEQ)=V(I,NP1) 1400 CONTINUE RETURN END SUBROUTINE EXPND(NVAR,MEQA,NSCOL,IRESL,EQ,X) C***BEGIN PROLOGUE EXPND C***REFER TO CDFCOR C***ROUTINES CALLED (NONE) C***PURPOSE THIS SUBROUTINE COMPUTES THE VALUES OF THE MEQA C SOLVED FOR VARIABLES AND SPLICES THEM INTO X. C***END PROLOGUE EXPND C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION NSCOL(15),IRESL(30),EQ(15,16),X(16) C C IF MEQA=0 NO EXPANDING NEEDS TO BE DONE AND WE RETURN. IF(MEQA)100,100,200 100 RETURN C 200 N=NVAR+1 C SHIFT THE LAST ELEMENT IN X INTO ITS PROPER PLACE. NOTE C THAT THIS ELEMENT CORRESPONDS TO THE VARIABLE W, WHICH IS C UNSOLVED FOR SINCE IT IS NOT INCLUDED IN EQ. NMEQ=N-MEQA X(N)=X(NMEQ) NS=NVAR-MEQA C SHIFT THE NS UNSOLVED FOR VARIABLES (NOT COUNTING W) INTO THE C PROPER PLACES IN X, IF THERE ARE ANY. IF(NS)500,500,300 300 DO 400 J=1,NS JJ=NS-J+1 X(NSCOL(JJ))=X(JJ) 400 CONTINUE C C NOW COMPUTE THE MEQA SOLVED FOR VARIABLES. 500 DO 800 L=1,MEQA C THE LTH RESOLVENT POSITION IN EQ IS (LROW,LCOL). WE C COMPUTE X(LCOL). LROW=IRESL(2*L-1) LCOL=IRESL(2*L) X(LCOL)=EQ(LROW,N) IF(NS)800,800,600 600 DO 700 J=1,NS JCOL=NSCOL(J) X(LCOL)=X(LCOL)-EQ(LROW,JCOL)*X(JCOL) 700 CONTINUE 800 CONTINUE RETURN END SUBROUTINE NRMLZ(NUMP,NUMQ,MEQ,MINEQ,NSCOL,IRESL,ALC, *TOLNR,X) C***BEGIN PROLOGUE NRMLZ C***REFER TO CDFCOR C***ROUTINES CALLED (NONE) C***PURPOSE THIS SUBROUTINE ATTEMPTS TO NORMALIZE X SO THAT C MAX(ABS(Q(J)))=1.0, OR SOME INEQUALITY CONSTRAINT C IN ALC WITH LEFT AND RIGHT SIDES .GT. TOLNR HOLDS C WITH EQUALITY. THE NORMALIZATION IS SKIPPED IF IT C WOULD REQUIRE MULTIPLYING BY A FACTOR LESS THAN 1.0, C OR IF SOME EQUALITY CONSTRAINT HAS RIGHT SIDE WITH C ABSOLUTE VALUE .GT. TOLNR. C***END PROLOGUE NRMLZ C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION NSCOL(15),IRESL(30),ALC(101,16),X(16) C C SET PRECISION DEPENDENT CONSTANTS FOR NRMLZ. C***FIRST EXECUTABLE STATEMENT NRMLZ ONE=1.0D0 C END OF SETTING PRECISION DEPENDENT CONSTANTS FOR NRMLZ. C NMPPQ=NUMP+NUMQ N=NMPPQ+1 IF(MEQ)300,300,100 C IF THE RIGHT SIDE OF SOME EQUALITY CONSTRAINT IS .GT. TOLNR, C WE SKIP THE NORMALIZATION. 100 DO 200 I=1,MEQ IF(ABS(ALC(I,N))-TOLNR)200,200,1500 200 CONTINUE C C COMPUTE AMAX = MAX(TOLNR,ABS(Q(1)),...,ABS(Q(NUMQ))). 300 AMAX=TOLNR DO 500 J=1,NUMQ NPJ=NUMP+J IF(ABS(X(NPJ))-AMAX)500,500,400 400 AMAX=ABS(X(NPJ)) 500 CONTINUE FACT=ONE/AMAX IF(MINEQ)1200,1200,600 C C CHECK THE INEQUALITY CONSTRAINTS IN ALC WITH SIGNIFICANTLY C POSITIVE LEFT AND RIGHT SIDES. MULTIPLYING ALL COEFFICIENTS C BY A FACTOR .GT. 1.0 WILL NOT AFFECT INEQUALITY CONSTRAINTS C WITH LEFT AND RIGHT SIDES OF OPPOSITE SIGN OR ZERO RIGHT SIDES C AND CAN ONLY HELP INEQUALITY CONSTRAINTS WITH NEGATIVE LEFT C AND RIGHT SIDES. 600 DO 1100 I=1,MINEQ MPI=MEQ+I IF(ALC(MPI,N)-TOLNR)1100,1100,700 700 ALFT=ALC(MPI,1)*X(1) DO 800 J=2,NMPPQ ALFT=ALFT+ALC(MPI,J)*X(J) 800 CONTINUE IF(ALFT-TOLNR)1100,1100,900 900 IF(ALC(MPI,N)/ALFT-FACT)1000,1100,1100 1000 FACT=ALC(MPI,N)/ALFT 1100 CONTINUE C C IF FACT .GT. 1.0, MULTIPLY ALL COEFFICIENTS BY FACT. 1200 IF(FACT-ONE)1500,1500,1300 1300 DO 1400 J=1,NMPPQ X(J)=FACT*X(J) 1400 CONTINUE 1500 RETURN END FUNCTION ICNCK(IOPT,NUMP,NUMQ,NUMGR,INDLW,ALTBL,INDUP, *UPTBL,EPSQ,X,V,MEQ,MINEQ,ALC,EPRES) C***BEGIN PROLOGUE ICNCK C***REFER TO CDFCOR C***ROUTINES CALLED (NONE) C***PURPOSE THIS SUBPROGRAM CHECKS ALL CONSTRAINTS, RETURNING C WITH VALUE 0 IF ALL ARE OK AND VALUE 1 IF ANY ARE C NOT SATISFIED WITHIN EPRES. C***END PROLOGUE ICNCK C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION INDLW(101),ALTBL(101),INDUP(101),UPTBL(101), *X(16),V(635,17),ALC(101,16) C C SET PRECISION DEPENDENT CONSTANTS FOR FUNCTION ICNCK. C***FIRST EXECUTABLE STATEMENT ICNCK ONE=1.0D0 C END OF SETTING PRECISION DEPENDENT CONSTANTS FOR ICNCK. C NMPPQ=NUMP+NUMQ N=NMPPQ+1 NP1=N+1 IONE=IOPT-(IOPT/10)*10 C C IF V(2*NUMGR+1,N+1) IS NEGATIVE, THEN SPQERD DETECTED A POINT C WHERE THE DENOMINATOR WAS VERY SMALL OR NEGATIVE, SO WE SET C ICNCK=1 AND RETURN. IF(V(2*NUMGR+1,NP1))2200,100,100 C C CHECK LOWER RESTRICTED RANGE CONSTRAINTS, IF THERE ARE ANY. 100 IF(IOPT-(IOPT/100)*100-10)500,200,200 200 DO 400 I=1,NUMGR IF(INDLW(I))400,400,300 300 IF(V(2*I-1,NP1)/V(2*I,NP1)-ALTBL(I)+EPRES)2200,400,400 400 CONTINUE C C CHECK UPPER RESTRICTED RANGE CONSTRAINTS, IF THERE ARE ANY. 500 IF(IOPT-(IOPT/1000)*1000-100)900,600,600 600 DO 800 I=1,NUMGR IF(INDUP(I))800,800,700 700 IF(V(2*I-1,NP1)/V(2*I,NP1)-UPTBL(I)-EPRES)800,800,2200 800 CONTINUE C C CHECK DENOMINATOR LOWER BOUND CONSTRAINTS, IF THERE ARE ANY. 900 IF(IONE-(IONE/2)*2)1200,1200,1000 1000 DO 1100 I=1,NUMGR IF(V(2*I,NP1)-EPSQ+EPRES)2200,1100,1100 1100 CONTINUE C C CHECK ADDITIONAL EQUALITY CONSTRAINTS, IF THERE ARE ANY. 1200 IF(MEQ)1600,1600,1300 1300 DO 1500 L=1,MEQ CHK=-ALC(L,N) DO 1400 J=1,NMPPQ CHK=CHK+ALC(L,J)*X(J) 1400 CONTINUE IF(ABS(CHK)-EPRES)1500,1500,2200 1500 CONTINUE C C CHECK ADDITIONAL INEQUALITY CONSTRAINTS, IF THERE ARE ANY. 1600 IF(MINEQ)2000,2000,1700 1700 DO 1900 L=1,MINEQ MEQL=MEQ+L CHK=-ALC(MEQL,N) DO 1800 J=1,NMPPQ CHK=CHK+ALC(MEQL,J)*X(J) 1800 CONTINUE IF(CHK-EPRES)1900,1900,2200 1900 CONTINUE C C CHECK THE DENOMINATOR COEFFICIENT BOUNDS. 2000 OPE=ONE+EPRES DO 2100 J=1,NUMQ NPJ=NUMP+J IF(ABS(X(NPJ))-OPE)2100,2100,2200 2100 CONTINUE C C IF WE REACH THIS SPOT, THEN ALL CONSTRAINTS ARE OK WITHIN EPRES. ICNCK=0 RETURN C C HERE AT LEAST ONE CONSTRAINT WAS NOT SATISFIED WITHIN EPRES. 2200 ICNCK=1 RETURN END 1 2 3 1 101 1011 0.0 2.0 1 7 8 1 21 21011 0.03.141592653589793238 1 7 8 1 41 121011 0.03.141592653589793238 2 3 2 5 5 1001 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 1 3 1 1 101 2 2 1 0.0 2.0 1 2 3 1 101 1013 1 0 0.0 2.0 0 0 0 0 0