Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/fmfp.ssp
There are 2 other files named fmfp.ssp in the archive. Click here to see a list.
C                                                                       FMFP  10
C     ..................................................................FMFP  20
C                                                                       FMFP  30
C        SUBROUTINE FMFP                                                FMFP  40
C                                                                       FMFP  50
C        PURPOSE                                                        FMFP  60
C           TO FIND A LOCAL MINIMUM OF A FUNCTION OF SEVERAL VARIABLES  FMFP  70
C           BY THE METHOD OF FLETCHER AND POWELL                        FMFP  80
C                                                                       FMFP  90
C        USAGE                                                          FMFP 100
C           CALL FMFP(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H)                FMFP 110
C                                                                       FMFP 120
C        DESCRIPTION OF PARAMETERS                                      FMFP 130
C           FUNCT  - USER-WRITTEN SUBROUTINE CONCERNING THE FUNCTION TO FMFP 140
C                    BE MINIMIZED. IT MUST BE OF THE FORM               FMFP 150
C                    SUBROUTINE FUNCT(N,ARG,VAL,GRAD)                   FMFP 160
C                    AND MUST SERVE THE FOLLOWING PURPOSE               FMFP 170
C                    FOR EACH N-DIMENSIONAL ARGUMENT VECTOR  ARG,       FMFP 180
C                    FUNCTION VALUE AND GRADIENT VECTOR MUST BE COMPUTEDFMFP 190
C                    AND, ON RETURN, STORED IN VAL AND GRAD RESPECTIVELYFMFP 200
C           N      - NUMBER OF VARIABLES                                FMFP 210
C           X      - VECTOR OF DIMENSION N CONTAINING THE INITIAL       FMFP 220
C                    ARGUMENT WHERE THE ITERATION STARTS. ON RETURN,    FMFP 230
C                    X HOLDS THE ARGUMENT CORRESPONDING TO THE          FMFP 240
C                    COMPUTED MINIMUM FUNCTION VALUE                    FMFP 250
C           F      - SINGLE VARIABLE CONTAINING THE MINIMUM FUNCTION    FMFP 260
C                    VALUE ON RETURN, I.E. F=F(X).                      FMFP 270
C           G      - VECTOR OF DIMENSION N CONTAINING THE GRADIENT      FMFP 280
C                    VECTOR CORRESPONDING TO THE MINIMUM ON RETURN,     FMFP 290
C                    I.E. G=G(X).                                       FMFP 300
C           EST    - IS AN ESTIMATE OF THE MINIMUM FUNCTION VALUE.      FMFP 310
C           EPS    - TESTVALUE REPRESENTING THE EXPECTED ABSOLUTE ERROR.FMFP 320
C                    A REASONABLE CHOICE IS 10**(-6), I.E.              FMFP 330
C                    SOMEWHAT GREATER THAN 10**(-D), WHERE D IS THE     FMFP 340
C                    NUMBER OF SIGNIFICANT DIGITS IN FLOATING POINT     FMFP 350
C                    REPRESENTATION.                                    FMFP 360
C           LIMIT  - MAXIMUM NUMBER OF ITERATIONS.                      FMFP 370
C           IER    - ERROR PARAMETER                                    FMFP 380
C                    IER = 0 MEANS CONVERGENCE WAS OBTAINED             FMFP 390
C                    IER = 1 MEANS NO CONVERGENCE IN LIMIT ITERATIONS   FMFP 400
C                    IER =-1 MEANS ERRORS IN GRADIENT CALCULATION       FMFP 410
C                    IER = 2 MEANS LINEAR SEARCH TECHNIQUE INDICATES    FMFP 420
C                    IT IS LIKELY THAT THERE EXISTS NO MINIMUM.         FMFP 430
C           H      - WORKING STORAGE OF DIMENSION N*(N+7)/2.            FMFP 440
C                                                                       FMFP 450
C        REMARKS                                                        FMFP 460
C            I) THE SUBROUTINE NAME REPLACING THE DUMMY ARGUMENT  FUNCT FMFP 470
C               MUST BE DECLARED AS EXTERNAL IN THE CALLING PROGRAM.    FMFP 480
C           II) IER IS SET TO 2 IF , STEPPING IN ONE OF THE COMPUTED    FMFP 490
C               DIRECTIONS, THE FUNCTION WILL NEVER INCREASE WITHIN     FMFP 500
C               A TOLERABLE RANGE OF ARGUMENT.                          FMFP 510
C               IER = 2 MAY OCCUR ALSO IF THE INTERVAL WHERE F          FMFP 520
C               INCREASES IS SMALL AND THE INITIAL ARGUMENT WAS         FMFP 530
C               RELATIVELY FAR AWAY FROM THE MINIMUM SUCH THAT THE      FMFP 540
C               MINIMUM WAS OVERLEAPED. THIS IS DUE TO THE SEARCH       FMFP 550
C               TECHNIQUE WHICH DOUBLES THE STEPSIZE UNTIL A POINT      FMFP 560
C               IS FOUND WHERE THE FUNCTION INCREASES.                  FMFP 570
C                                                                       FMFP 580
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  FMFP 590
C           FUNCT                                                       FMFP 600
C                                                                       FMFP 610
C        METHOD                                                         FMFP 620
C           THE METHOD IS DESCRIBED IN THE FOLLOWING ARTICLE            FMFP 630
C           R. FLETCHER AND M.J.D. POWELL, A RAPID DESCENT METHOD FOR   FMFP 640
C           MINIMIZATION,                                               FMFP 650
C           COMPUTER JOURNAL VOL.6, ISS. 2, 1963, PP.163-168.           FMFP 660
C                                                                       FMFP 670
C     ..................................................................FMFP 680
C                                                                       FMFP 690
      SUBROUTINE FMFP(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H)                FMFP 700
C                                                                       FMFP 710
C        DIMENSIONED DUMMY VARIABLES                                    FMFP 720
      DIMENSION H(1),X(1),G(1)                                          FMFP 730
C                                                                       FMFP 740
C        COMPUTE FUNCTION VALUE AND GRADIENT VECTOR FOR INITIAL ARGUMENTFMFP 750
      CALL FUNCT(N,X,F,G)                                               FMFP 760
C                                                                       FMFP 770
C        RESET ITERATION COUNTER AND GENERATE IDENTITY MATRIX           FMFP 780
      IER=0                                                             FMFP 790
      KOUNT=0                                                           FMFP 800
      N2=N+N                                                            FMFP 810
      N3=N2+N                                                           FMFP 820
      N31=N3+1                                                          FMFP 830
    1 K=N31                                                             FMFP 840
      DO 4 J=1,N                                                        FMFP 850
      H(K)=1.                                                           FMFP 860
      NJ=N-J                                                            FMFP 870
      IF(NJ)5,5,2                                                       FMFP 880
    2 DO 3 L=1,NJ                                                       FMFP 890
      KL=K+L                                                            FMFP 900
    3 H(KL)=0.                                                          FMFP 910
    4 K=KL+1                                                            FMFP 920
C                                                                       FMFP 930
C        START ITERATION LOOP                                           FMFP 940
    5 KOUNT=KOUNT +1                                                    FMFP 950
C                                                                       FMFP 960
C        SAVE FUNCTION VALUE, ARGUMENT VECTOR AND GRADIENT VECTOR       FMFP 970
      OLDF=F                                                            FMFP 980
      DO 9 J=1,N                                                        FMFP 990
      K=N+J                                                             FMFP1000
      H(K)=G(J)                                                         FMFP1010
      K=K+N                                                             FMFP1020
      H(K)=X(J)                                                         FMFP1030
C                                                                       FMFP1040
C        DETERMINE DIRECTION VECTOR H                                   FMFP1050
      K=J+N3                                                            FMFP1060
      T=0.                                                              FMFP1070
      DO 8 L=1,N                                                        FMFP1080
      T=T-G(L)*H(K)                                                     FMFP1090
      IF(L-J)6,7,7                                                      FMFP1100
    6 K=K+N-L                                                           FMFP1110
      GO TO 8                                                           FMFP1120
    7 K=K+1                                                             FMFP1130
    8 CONTINUE                                                          FMFP1140
    9 H(J)=T                                                            FMFP1150
C                                                                       FMFP1160
C        CHECK WHETHER FUNCTION WILL DECREASE STEPPING ALONG H.         FMFP1170
      DY=0.                                                             FMFP1180
      HNRM=0.                                                           FMFP1190
      GNRM=0.                                                           FMFP1200
C                                                                       FMFP1210
C        CALCULATE DIRECTIONAL DERIVATIVE AND TESTVALUES FOR DIRECTION  FMFP1220
C        VECTOR H AND GRADIENT VECTOR G.                                FMFP1230
      DO 10 J=1,N                                                       FMFP1240
      HNRM=HNRM+ABS(H(J))                                               FMFP1250
      GNRM=GNRM+ABS(G(J))                                               FMFP1260
   10 DY=DY+H(J)*G(J)                                                   FMFP1270
C                                                                       FMFP1280
C        REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DIRECTIONAL  FMFP1290
C        DERIVATIVE APPEARS TO BE POSITIVE OR ZERO.                     FMFP1300
      IF(DY)11,51,51                                                    FMFP1310
C                                                                       FMFP1320
C        REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DIRECTION    FMFP1330
C        VECTOR H IS SMALL COMPARED TO GRADIENT VECTOR G.               FMFP1340
   11 IF(HNRM/GNRM-EPS)51,51,12                                         FMFP1350
C                                                                       FMFP1360
C        SEARCH MINIMUM ALONG DIRECTION H                               FMFP1370
C                                                                       FMFP1380
C        SEARCH ALONG H FOR POSITIVE DIRECTIONAL DERIVATIVE             FMFP1390
   12 FY=F                                                              FMFP1400
      ALFA=2.*(EST-F)/DY                                                FMFP1410
      AMBDA=1.                                                          FMFP1420
C                                                                       FMFP1430
C        USE ESTIMATE FOR STEPSIZE ONLY IF IT IS POSITIVE AND LESS THAN FMFP1440
C        1. OTHERWISE TAKE 1. AS STEPSIZE                               FMFP1450
      IF(ALFA)15,15,13                                                  FMFP1460
   13 IF(ALFA-AMBDA)14,15,15                                            FMFP1470
   14 AMBDA=ALFA                                                        FMFP1480
   15 ALFA=0.                                                           FMFP1490
C                                                                       FMFP1500
C        SAVE FUNCTION AND DERIVATIVE VALUES FOR OLD ARGUMENT           FMFP1510
   16 FX=FY                                                             FMFP1520
      DX=DY                                                             FMFP1530
C                                                                       FMFP1540
C        STEP ARGUMENT ALONG H                                          FMFP1550
      DO 17 I=1,N                                                       FMFP1560
   17 X(I)=X(I)+AMBDA*H(I)                                              FMFP1570
C                                                                       FMFP1580
C        COMPUTE FUNCTION VALUE AND GRADIENT FOR NEW ARGUMENT           FMFP1590
      CALL FUNCT(N,X,F,G)                                               FMFP1600
      FY=F                                                              FMFP1610
C                                                                       FMFP1620
C        COMPUTE DIRECTIONAL DERIVATIVE DY FOR NEW ARGUMENT.  TERMINATE FMFP1630
C        SEARCH, IF DY IS POSITIVE. IF DY IS ZERO THE MINIMUM IS FOUND  FMFP1640
      DY=0.                                                             FMFP1650
      DO 18 I=1,N                                                       FMFP1660
   18 DY=DY+G(I)*H(I)                                                   FMFP1670
      IF(DY)19,36,22                                                    FMFP1680
C                                                                       FMFP1690
C        TERMINATE SEARCH ALSO IF THE FUNCTION VALUE INDICATES THAT     FMFP1700
C        A MINIMUM HAS BEEN PASSED                                      FMFP1710
   19 IF(FY-FX)20,22,22                                                 FMFP1720
C                                                                       FMFP1730
C        REPEAT SEARCH AND DOUBLE STEPSIZE FOR FURTHER SEARCHES         FMFP1740
   20 AMBDA=AMBDA+ALFA                                                  FMFP1750
      ALFA=AMBDA                                                        FMFP1760
C        END OF SEARCH LOOP                                             FMFP1770
C                                                                       FMFP1780
C        TERMINATE IF THE CHANGE IN ARGUMENT GETS VERY LARGE            FMFP1790
      IF(HNRM*AMBDA-1.E10)16,16,21                                      FMFP1800
C                                                                       FMFP1810
C        LINEAR SEARCH TECHNIQUE INDICATES THAT NO MINIMUM EXISTS       FMFP1820
   21 IER=2                                                             FMFP1830
      RETURN                                                            FMFP1840
C                                                                       FMFP1850
C        INTERPOLATE CUBICALLY IN THE INTERVAL DEFINED BY THE SEARCH    FMFP1860
C        ABOVE AND COMPUTE THE ARGUMENT X FOR WHICH THE INTERPOLATION   FMFP1870
C        POLYNOMIAL IS MINIMIZED                                        FMFP1880
   22 T=0.                                                              FMFP1890
   23 IF(AMBDA)24,36,24                                                 FMFP1900
   24 Z=3.*(FX-FY)/AMBDA+DX+DY                                          FMFP1910
      ALFA=AMAX1(ABS(Z),ABS(DX),ABS(DY))                                FMFP1920
      DALFA=Z/ALFA                                                      FMFP1930
      DALFA=DALFA*DALFA-DX/ALFA*DY/ALFA                                 FMFP1940
      IF(DALFA)51,25,25                                                 FMFP1950
   25 W=ALFA*SQRT(DALFA)                                                FMFP1960
      ALFA=DY-DX+W+W                                                    FMFP1970
      IF(ALFA) 250,251,250                                              FMFP1971
  250 ALFA=(DY-Z+W)/ALFA                                                FMFP1972
      GO TO 252                                                         FMFP1973
  251 ALFA=(Z+DY-W)/(Z+DX+Z+DY)                                         FMFP1974
  252 ALFA=ALFA*AMBDA                                                   FMFP1975
      DO 26 I=1,N                                                       FMFP1980
   26 X(I)=X(I)+(T-ALFA)*H(I)                                           FMFP1990
C                                                                       FMFP2000
C        TERMINATE, IF THE VALUE OF THE ACTUAL FUNCTION AT X IS LESS    FMFP2010
C        THAN THE FUNCTION VALUES AT THE INTERVAL ENDS. OTHERWISE REDUCEFMFP2020
C        THE INTERVAL BY CHOOSING ONE END-POINT EQUAL TO X AND REPEAT   FMFP2030
C        THE INTERPOLATION.  WHICH END-POINT IS CHOOSEN DEPENDS ON THE  FMFP2040
C        VALUE OF THE FUNCTION AND ITS GRADIENT AT X                    FMFP2050
C                                                                       FMFP2060
      CALL FUNCT(N,X,F,G)                                               FMFP2070
      IF(F-FX)27,27,28                                                  FMFP2080
   27 IF(F-FY)36,36,28                                                  FMFP2090
   28 DALFA=0.                                                          FMFP2100
      DO 29 I=1,N                                                       FMFP2110
   29 DALFA=DALFA+G(I)*H(I)                                             FMFP2120
      IF(DALFA)30,33,33                                                 FMFP2130
   30 IF(F-FX)32,31,33                                                  FMFP2140
   31 IF(DX-DALFA)32,36,32                                              FMFP2150
   32 FX=F                                                              FMFP2160
      DX=DALFA                                                          FMFP2170
      T=ALFA                                                            FMFP2180
      AMBDA=ALFA                                                        FMFP2190
      GO TO 23                                                          FMFP2200
   33 IF(FY-F)35,34,35                                                  FMFP2210
   34 IF(DY-DALFA)35,36,35                                              FMFP2220
   35 FY=F                                                              FMFP2230
      DY=DALFA                                                          FMFP2240
      AMBDA=AMBDA-ALFA                                                  FMFP2250
      GO TO 22                                                          FMFP2260
C                                                                       FMFP2270
C        TERMINATE, IF FUNCTION HAS NOT DECREASED DURING LAST ITERATION FMFP2280
   36 IF(OLDF-F+EPS)51,38,38                                            FMFP2290
C                                                                       FMFP2300
C        COMPUTE DIFFERENCE VECTORS OF ARGUMENT AND GRADIENT FROM       FMFP2310
C        TWO CONSECUTIVE ITERATIONS                                     FMFP2320
   38 DO 37 J=1,N                                                       FMFP2330
      K=N+J                                                             FMFP2340
      H(K)=G(J)-H(K)                                                    FMFP2350
      K=N+K                                                             FMFP2360
   37 H(K)=X(J)-H(K)                                                    FMFP2370
C                                                                       FMFP2380
C        TEST LENGTH OF ARGUMENT DIFFERENCE VECTOR AND DIRECTION VECTOR FMFP2390
C        IF AT LEAST N ITERATIONS HAVE BEEN EXECUTED. TERMINATE, IF     FMFP2400
C        BOTH ARE LESS THAN  EPS                                        FMFP2410
      IER=0                                                             FMFP2420
      IF(KOUNT-N)42,39,39                                               FMFP2430
   39 T=0.                                                              FMFP2440
      Z=0.                                                              FMFP2450
      DO 40 J=1,N                                                       FMFP2460
      K=N+J                                                             FMFP2470
      W=H(K)                                                            FMFP2480
      K=K+N                                                             FMFP2490
      T=T+ABS(H(K))                                                     FMFP2500
   40 Z=Z+W*H(K)                                                        FMFP2510
      IF(HNRM-EPS)41,41,42                                              FMFP2520
   41 IF(T-EPS)56,56,42                                                 FMFP2530
C                                                                       FMFP2540
C        TERMINATE, IF NUMBER OF ITERATIONS WOULD EXCEED  LIMIT         FMFP2550
   42 IF(KOUNT-LIMIT)43,50,50                                           FMFP2560
C                                                                       FMFP2570
C        PREPARE UPDATING OF MATRIX H                                   FMFP2580
   43 ALFA=0.                                                           FMFP2590
      DO 47 J=1,N                                                       FMFP2600
      K=J+N3                                                            FMFP2610
      W=0.                                                              FMFP2620
      DO 46 L=1,N                                                       FMFP2630
      KL=N+L                                                            FMFP2640
      W=W+H(KL)*H(K)                                                    FMFP2650
      IF(L-J)44,45,45                                                   FMFP2660
   44 K=K+N-L                                                           FMFP2670
      GO TO 46                                                          FMFP2680
   45 K=K+1                                                             FMFP2690
   46 CONTINUE                                                          FMFP2700
      K=N+J                                                             FMFP2710
      ALFA=ALFA+W*H(K)                                                  FMFP2720
   47 H(J)=W                                                            FMFP2730
C                                                                       FMFP2740
C        REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF RESULTS      FMFP2750
C        ARE NOT SATISFACTORY                                           FMFP2760
      IF(Z*ALFA)48,1,48                                                 FMFP2770
C                                                                       FMFP2780
C        UPDATE MATRIX H                                                FMFP2790
   48 K=N31                                                             FMFP2800
      DO 49 L=1,N                                                       FMFP2810
      KL=N2+L                                                           FMFP2820
      DO 49 J=L,N                                                       FMFP2830
      NJ=N2+J                                                           FMFP2840
      H(K)=H(K)+H(KL)*H(NJ)/Z-H(L)*H(J)/ALFA                            FMFP2850
   49 K=K+1                                                             FMFP2860
      GO TO 5                                                           FMFP2870
C        END OF ITERATION LOOP                                          FMFP2880
C                                                                       FMFP2890
C        NO CONVERGENCE AFTER  LIMIT  ITERATIONS                        FMFP2900
   50 IER=1                                                             FMFP2910
      RETURN                                                            FMFP2920
C                                                                       FMFP2930
C        RESTORE OLD VALUES OF FUNCTION AND ARGUMENTS                   FMFP2940
   51 DO 52 J=1,N                                                       FMFP2950
      K=N2+J                                                            FMFP2960
   52 X(J)=H(K)                                                         FMFP2970
      CALL FUNCT(N,X,F,G)                                               FMFP2980
C                                                                       FMFP2990
C        REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DERIVATIVE   FMFP3000
C        FAILS TO BE SUFFICIENTLY SMALL                                 FMFP3010
      IF(GNRM-EPS)55,55,53                                              FMFP3020
C                                                                       FMFP3030
C        TEST FOR REPEATED FAILURE OF ITERATION                         FMFP3040
   53 IF(IER)56,54,54                                                   FMFP3050
   54 IER=-1                                                            FMFP3060
      GOTO 1                                                            FMFP3070
   55 IER=0                                                             FMFP3080
   56 RETURN                                                            FMFP3090
      END                                                               FMFP3100