C ALGORITHM 450, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN COMMUNICATIONS OF THE ACM C VOL. 16, NO. 8, August, 1973, PP.482--483. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Fortran/ # Fortran/Sp/ # Fortran/Sp/Drivers/ # Fortran/Sp/Drivers/Makefile # Fortran/Sp/Drivers/driver.f # Fortran/Sp/Drivers/res # Fortran/Sp/Src/ # Fortran/Sp/Src/src.f # This archive created: Thu Dec 15 13:28:13 2005 export PATH; PATH=/bin:$PATH if test ! -d 'Fortran' then mkdir 'Fortran' fi cd 'Fortran' if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << "SHAR_EOF" > 'Makefile' all: Res src.o: src.f $(F77) $(F77OPTS) -c src.f driver.o: driver.f $(F77) $(F77OPTS) -c driver.f DRIVERS= driver RESULTS= Res Objs1= driver.o src.o driver: $(Objs1) $(F77) $(F77OPTS) -o driver $(Objs1) $(SRCLIBS) Res: driver ./driver >Res diffres:Res res echo "Differences in results from driver" $(DIFF) Res res clean: rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS) SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << "SHAR_EOF" > 'driver.f' program main c*********************************************************************** c cc TOMS450_PRB tests ROMIN. c implicit none integer n, ncyc parameter ( n = 2 ) real f real f2 external funct integer i external monitr real step real x(n) real x2(n) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOMS450_PRB' write ( *, '(a)' ) ' Test TOMS algorithm 450, Rosenbrock' write ( *, '(a)' ) ' function minimization.' write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' Dimension of problem N = ', n x(1) = -1.2E+00 x(2) = 1.0E+00 ncyc = 3 step = 0.1E+00 call funct ( n, x, f ) write ( *, '(a)' ) ' ' write ( *, '(a,2g14.6)' ) ' Initial X = ', ( x(i), i = 1, n ) write ( *, '(a,g14.6)' ) ' Initial F(X) = ', f call romin ( n, x, funct, step, ncyc, monitr ) call funct ( n, x, f ) write ( *, '(a)' ) ' ' write ( *, '(a,2g14.6)' ) ' Final X = ', ( x(i), i = 1, n ) write ( *, '(a,g14.6)' ) ' Final F(X) = ', f x2(1) = 0.99513E+00 x2(2) = 0.99053E+00 call funct ( n, x2, f2 ) write ( *, '(a)' ) ' ' write ( *, '(a,2g14.6)' ) ' Expected X = ', ( x2(i), i = 1, n ) write ( *, '(a,g14.6)' ) ' Expected F(X) = ', f2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOMS452_PRB' write ( *, '(a)' ) ' Normal end of execution.' stop end subroutine monitr ( n, x, f, r, b, con, nr ) c*********************************************************************** c cc MONITR monitors the progress of the iteration. c c Discussion: c c This version of the monitor routine is designed simply c to stop the iteration after 200 steps have been taken. c c Parameters: c c R IS THE ACTUAL NUMBER OF FUNCTION EVALUATIONS (FOR THE c INITIAL ESTIMATE R = 0.) c c B IS THE VALUE OF THE EUCLIDEAN NORM OF THE VECTOR c REPRESENTING THE TOTAL PROGRESS MADE SINCE THE c AXES WERE LAST ROTATED. c c CON IS A LOGICAL VARIABLE. AT THE START OF THE c SUBROUTINE ROMIN CON=.FALSE. IF THE CONVERGENCE c CRITERIA OF THE ROUTINE MONITOR ARE SATISFIED c CON MUST BE SET TO .TRUE. TO STOP THE PROCESS. c c NR IS THE MONITOR INDEX. c implicit none integer n real b logical con real f integer nr integer r real x(n) if ( 200 .le. r ) then con = .true. end if return end subroutine funct ( n, x, f ) c*********************************************************************** c cc FUNCT evaluates the scalar function of N variables to be minimized. c implicit none integer n real f real x(n) f = 100.0E+00 * ( x(2) - x(1) * x(1) )**2 + ( 1.0E+00 - x(1) )**2 return end SHAR_EOF fi # end of overwriting check if test -f 'res' then echo shar: will not over-write existing file "'res'" else cat << "SHAR_EOF" > 'res' TOMS450_PRB Test TOMS algorithm 450, Rosenbrock function minimization. Dimension of problem N = 2 Initial X = -1.20000 1.00000 Initial F(X) = 24.2000 Final X = 1.00000 1.00000 Final F(X) = 0.568434E-13 Expected X = 0.995130 0.990530 Expected F(X) = 0.297826E-04 TOMS452_PRB Normal end of execution. SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << "SHAR_EOF" > 'src.f' SUBROUTINE ROMIN(N, X, FUNCT, STEP, NCYC, MONITR ) C Incorporates changes suggested in the remark by Bultheel in C cacm 17(8) August 1974 p 470 C Incorporates changes suggested in the remark by Davies in C TOMS 2(3) Sept. 1976, pp300--301 INTEGER N, IP REAL STEP, V, X, E, D, A, ALPHA DIMENSION X(N) LOGICAL CON INTEGER I, J, K, L, P, R, M, JCYC, NCYC REAL F0, F1, B, BETY C **** The following statements have been changed to incorporate C **** the remark by Davies TOMS 2(3) Sept. 1976, pp300--301 C DIMENSION A(30), D(30), V(30,30), ALPHA(30,30), BETA(30), C & E(30), AV(30) DIMENSION A(30), D(30), V(30,30), ALPHA(30,30), & E(30) C THIS SUBROUTINE MINIMIZES A FUNCTION OF N VARIABLES C USING THE METHOD OF ROSENBROCK. THE PARAMETERS ARE C DESCRIBED AS FOLLOWS: C N IS THE NUMBER OF INDEPENDENT VARIABLES. C X(N) IS AN ESTIMATE OF THE SOLUTION ( ON ENTRY - C AN INITIAL ESTIMATE, ON EXIT - THE BEST ESTIMATE C OF THE SOLUTION FOUND.) C FUNCT(N,X,F) IS A ROUTINE PROVIDED BY THE USER TO C CALCULATE THE VALUE F OF THE MINIMIZED FUNCTION C AT ANY POINT X. C STEP IS AN INITIAL STEP LENGTH FOR ALL COORDINATE C DIRECTIONS AT THE START OF THE PROCESS. C MONITR(N,X,F,R,S,CON,NR) IS A ROUTINE PROVIDED BY C THE USER FOR DIAGNOSTIC AND CONVERGENCE PURPOSES. C R IS THE ACTUAL NUMBER OF FUNCTION EVALUATIONS (FOR THE C INITIAL ESTIMATE R = 0.) C B IS THE VALUE OF THE EUCLIDEAN NORM OF THE VECTOR C REPRESENTING THE TOTAL PROGRESS MADE SINCE THE C AXES WERE LAST ROTATED. C CON IS A LOGICAL VARIABLE. AT THE START OF THE C SUBROUTINE ROMIN CON=.FALSE. IF THE CONVERGENCE C CRITERIA OF THE ROUTINE MONITOR ARE SATISFIED C CON MUST BE SET TO .TRUE. TO STOP THE PROCESS. C NR IS THE MONITOR INDEX. C INITIALIZE CON, D(I) AND R. C E(I) IS A SET OF STEPS TO BE TAKEN IN THE CORRESPONDING C COORDINATE DIRECTIONS. CON = .FALSE. DO 10 I=1,N E(I) = STEP 10 CONTINUE R = 0 C V(I,J) IS AN NXN MATRIX DEFINING A SET OF N MUTUALLY C ORTHOGONAL COORDINATE DIRECTIONS. V(I,J) IS THE UNIT C MATRIX AT THE START OF THE PROCESS. DO 30 I=1,N DO 20 J=1,N V(I,J) = 0.0 IF (I.EQ.J) V(I,J) = 1.0 20 CONTINUE 30 CONTINUE CALL FUNCT(N, X, F0) CALL MONITR(N, X, F0, R, 1.0E10, CON, 0) C START OF THE ITERATION LOOP 40 DO 50 I=1,N A(I) = 2.0 D(I) = 0.0 50 CONTINUE C EVALUATE F AT THE NEW POINT X. 60 DO 130 I=1,N DO 70 J=1,N X(J) = X(J) + E(I)*V(I,J) 70 CONTINUE R = R + 1 CALL FUNCT(N, X, F1 ) CALL MONITR(N, X, F1, R, 0.0, CON, 1) IF ( CON ) GO TO 290 IF (F1-F0) 80, 90, 90 C THE NEW VALUE OF THE FUNCTION IS LESS THAN THE OLD ONE. 80 D(I) = D(I) + E(I) E(I) = 3.0 * E(I) F0 = F1 IF (A(I).GT.1.5) A(I) = 1.0 GO TO 110 C THE NEW VALUE OF THE FUNCTION IS GREATER THAN OR EQUAL C TO THE OLD ONE. 90 DO 100 J=1,N X(J) = X(J) - E(I)*V(I,J) 100 CONTINUE E(I) = -0.5 * E(I) IF (A(I).LT.1.5) A(I) = 0.0 110 DO 120 J=1,N IF (A(J).GE.0.5) GO TO 130 120 CONTINUE GO TO 140 130 CONTINUE GO TO 60 C GRAM-SCHMIDT ORTHOGONALIZATION PROCESS. 140 DO 160 K=1,N DO 150 L=1,N ALPHA(K,L) = 0.0 150 CONTINUE 160 CONTINUE DO 190 I=1,N DO 180 J=1,N DO 170 L=I,N ALPHA(I,J) = ALPHA(I,J) + D(L) * V(L,J) 170 CONTINUE 180 CONTINUE 190 CONTINUE B = 0.0 DO 200 J=1,N B = B + ALPHA(1,J)**2 200 CONTINUE B = SQRT(B) C CALCULATE THE NEW SET OF ORTHONORMAL COORDINATE C DIRECTIONS ( THE NEW MATRIX V(I,J) ). DO 210 J=1,N V(1,J) = ALPHA(1,J) / B 210 CONTINUE C **** The following statements have been changed to incorporate C **** the remark by Davies TOMS 2(3) Sept. 1976, pp300--301 C DO 280 P=2,N C BETY = 0.0 C IP = P - 1 C DO 220 M=1,N C BETA(M) = 0.0 C220 CONTINUE C DO 250 J=1,N C DO 240 K=1,IP C AV(K) = 0.0 C DO 230 L=1,N C AV(K) = AV(K) + ALPHA(P,L) * V(K,L) C230 CONTINUE C BETA(J) = BETA(J) - AV(K) * V(K,J) C240 CONTINUE C250 CONTINUE C DO 260 J=1,N C BETA(J) = BETA(J) + ALPHA(P,J) C BETY = BETY + BETA(J)**2 C260 CONTINUE C BETY = SQRT(BETY) C DO 270 J=1,N C V(P,J) = BETA(J) / BETY C270 CONTINUE C280 CONTINUE C **** End of replaced code DO 303 JCYC = 1, NCYC DO 250 P = 2, N IP = P-1 DO 230 M = P, N BETY = 0.0 DO 220 K = 1, N BETY = BETY - ALPHA(M, K)*V(IP, K) 220 CONTINUE DO 225 J = 1, N ALPHA(M, J) = ALPHA(M, J) + BETY*V(IP,J) 225 CONTINUE 230 CONTINUE BETY = 0.0 DO 240 K = 1, N BETY = BETY + ALPHA(P,K)**2 240 CONTINUE BETY = SQRT(BETY) DO 245 K = 1, N V(P,K) = ALPHA(P,K)/BETY 245 CONTINUE 250 CONTINUE IF (JCYC .EQ. NCYC) GO TO 303 DO 302 I = 2, N DO 301 J = 1, N ALPHA(I,J) = V(I,J) 301 CONTINUE 302 CONTINUE 303 CONTINUE C **** End of replacement code C END OF GRAM-SCHMIDT PROCESS. CALL MONITR(N, X, F0, R, B, CON, 2) IF ( CON ) GO TO 290 C GO TO THE NEXT ITERATION GO TO 40 290 RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. # End of shell archive exit 0