C**********************************************************
C
C     Rosenbrock's method for unconstrained minimization.     
C
C**********************************************************
C     Neculai Andrei, 
C     Research Institute for Informatics
C     Advanced Modeling and Optimization Laboratory     
C**********************************************************  
C
	SUBROUTINE ROSE(N,NN,X,X1,X2,Z,H,HM,XD,P,Y)
C
	DOUBLE PRECISION X(N),X1(N),X2(N),Z(N),H(N),Y(N)
	DOUBLE PRECISION HM(NN),XD(NN),P(NN)
C
	DOUBLE PRECISION EPS,FX,FY,S,T,ALFA,BETA
C
	CHARACTER*1 FD
C
	COMMON /PREC/EPS
	COMMON /TYPE/NO
	COMMON /WRIT/NX
C
	NFUNC=0    
	JTYPE=0
	IR=0
C
	WRITE(NX,150)
150	FORMAT(10X,'Start Optimization Process by Rosenbrock',
     +             ' Method')
C
110	WRITE(*,111)'   Printing [F=File/D=Display] ? '
111	FORMAT(A,$)
	READ(5,112,ERR=113) FD
112	FORMAT(A1)
	IF(FD.NE.'F'.AND.FD.NE.'D') THEN
	WRITE(NX,114)
114	FORMAT(10X,'Error.   Please try again')
	GO TO 110
	ENDIF
	IF(FD.EQ.'F') NO=6
	IF(FD.EQ.'D') NO=5
	GO TO 115
113	WRITE(NX,114)
	GO TO 110
115	CONTINUE
C
	IF(NO.EQ.6) THEN
	OPEN(6,FILE='PRINT.DAT',ACCESS='SEQUENTIAL',
     +         FORM='FORMATTED')
	ENDIF
C
	ALFA=3.0D0
	BETA=5.0D-1
C
	DO 1 I=1,N
1	H(I)=8.0D-1
C
	ITER=1
C
	CALL FUNC(N,X,FX)
	NFUNC=NFUNC+1
C
	WRITE(NO,50)
50	FORMAT(10X,'***** Optimization Method of Rosenbrock *****')
	WRITE(NO,51)
51	FORMAT(10X,'Initial Point : ')
	WRITE(NO,52)(X(I),I=1,N)
52	FORMAT(10X,F18.6)
	WRITE(NO,53)EPS
53	FORMAT(10X,'Accuracy EPS = ',F18.6)
	WRITE(NO,54)FX
54	FORMAT(10X,'Val. Func. in Initial Point = ',F18.6)
C
	DO 3 I=1,N
	DO 2 J=1,N
	L=J+(I-1)*N
2	XD(L)=0.0
	L=I+(I-1)*N
3	XD(L)=1.0
C
	DO 4 I=1,N
	X1(I)=0.0
	X2(I)=0.0
4	Z(I)=0.0
C
5	CONTINUE
C
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,100) ITER
100           FORMAT(10X,'ITERATION = ',I3)
              WRITE(NO,101)
101           FORMAT(10X,'Current Point: ')
              WRITE(NO,102)(I,X(I),I=1,N)
102           FORMAT(10X,I4,5X,F18.6)
              WRITE(NO,103)FX
103           FORMAT(10X,'Func. Value in current point =',F18.6)
      ENDIF
C
	I=1
6	DO 7 J=1,N
	L=J+(I-1)*N
7	Y(J)=X(J)+H(I)*XD(L)
C
	CALL FUNC(N,Y,FY)
	NFUNC=NFUNC+1
C
	IF(DABS(FX-FY).LE.EPS) IR=IR+1
	IF(IR.GT.20) GO TO 30
C
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,104)
104           FORMAT(10X,'Point Y(I)')
              WRITE(NO,102)(I1,Y(I1),I1=1,N)
              WRITE(NO,105)FY
105           FORMAT(10X,'Val func in point Y=',F18.6)
      ENDIF
C
      IF(FY.LT.FX) THEN
C      
C--------------------------------------------------------Succes.
C
              IF(JTYPE.EQ.1) WRITE(NO,106)I
106           FORMAT(10X,'I=',I3,'   Success')
C
	Z(I)=Z(I)+H(I)
	H(I)=ALFA*H(I)
	X1(I)=1.0
	FX=FY
	DO 8 J=1,N
8	X(J)=Y(J)
	ITER=ITER+1
	IF(ITER.GT.1001) THEN
	WRITE(NX,40)
40	FORMAT(10X,'Too many iterations.',/,
     +       10X,'This function is too complicated for the',/,
     +       10X,'power of the Rosenbrock Optimization Method.',/,
     +       10X,'May be this function is a ridge or a valley.',//,
     +       10X,'Please consider another initial point, or',/,
     +       10X,'try another optimization method')
	RETURN
	ENDIF
	GO TO 9
C
	ELSE
C      
C-------------------------------------------------------Failure.
C
              IF(JTYPE.EQ.1) WRITE(NO,107) I
107           FORMAT(10X,'I=',I3,'   Failure')
C
	H(I)=-BETA*H(I)
	X2(I)=1.0
	GO TO 9
C
	ENDIF
C
9	I=I+1
	IF(I.LE.N) GO TO 6
C
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,108)
108           FORMAT(10X,'Arrays X1 and X2')
              WRITE(NO,109)(I1,X1(I1),X2(I1),I1=1,N)
109           FORMAT(10X,I2,3X,F5.2,3X,F5.2)
      ENDIF
C
	IS=1
	DO 10 J=1,N
	IF(X1(J)*X2(J).NE.0.0) GO TO 10
	IS=0
	GO TO 11
10	CONTINUE
11	IF(IS.EQ.0) GO TO 5
C     
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,210)
210           FORMAT(10X,'Array Z(I)')
              WRITE(NO,102)(I1,Z(I1),I1=1,N)
      ENDIF                  
C
	IS=1
	DO 12 J=1,N
	IF(DABS(Z(J)).LE.EPS) GO TO 12
	IS=0
	GO TO 13
12	CONTINUE
13	IF(IS.EQ.1) GO TO 30
C
	DO 15 I=1,N
	DO 14 J=1,N
	L=J+(I-1)*N
14	HM(L)=Z(J)
	Z(I)=0.0
15	CONTINUE
C     
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,211)
211           FORMAT(10X,'Array HM:')
              WRITE(NO,102)(I1,HM(I1),I1=1,NN)
C
              WRITE(NO,213)
213           FORMAT(10X,'Array XD: ')
              WRITE(NO,102)(I1,XD(I1),I1=1,NN)
      ENDIF              
C
	DO 16 I=1,N
	DO 16 J=1,N
	L=I+(J-1)*N
	P(L)=0.0
	DO 16 K=1,N
	L1=I+(K-1)*N
	L2=K+(J-1)*N
16	P(L)=P(L)+XD(L1)*HM(L2)
C             
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,212)
212           FORMAT(10X,'Array P:')
              WRITE(NO,102)(I1,P(I1),I1=1,NN)
      ENDIF        
C
	S=0.0
	DO 17 J=1,N
17	S=S+P(J)*P(J)
	T=SQRT(S)
	DO 18 J=1,N
18	XD(J)=P(J)/T
C
	DO 25 I=2,N
	DO 19 K=1,N
	L=K+(I-1)*N
19	Y(K)=P(L)
C          
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,220)
220           FORMAT(15X,'Array Y: cycle 19')
              WRITE(NO,102)(I1,Y(I1),I1=1,N)
      ENDIF            
C
	DO 22 J=1,I-1
	S=0.0
	DO 20 K=1,N
	L1=K+(I-1)*N
	L2=K+(J-1)*N
20	S=S+P(L1)*XD(L2)
C
              IF(JTYPE.EQ.1) WRITE(NO,221)S
221           FORMAT(10X,'S = (cycle 20)',F18.6)
C
	DO 21 K=1,N
	L=K+(J-1)*N
21	Y(K)=Y(K)-S*XD(L)
C
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,222)
222           FORMAT(10X,'Array Y: cycle 21')
              WRITE(NO,102)(I1,Y(I1),I1=1,N)
      ENDIF
C            
22	CONTINUE
C
	S=0.0
	DO 23 K=1,N
23	S=S+Y(K)*Y(K)
C
              IF(JTYPE.EQ.1) WRITE(NO,223) S
223           FORMAT(10X,'S = (cycle 23)',F18.6)
C
	T=SQRT(S)
	DO 24 K=1,N
	L=K+(I-1)*N
24	XD(L)=Y(K)/T
25	CONTINUE
C
      IF(JTYPE.EQ.1) THEN
              WRITE(NO,216)
216           FORMAT(10X,'Array XD - New directions')
              WRITE(NO,102)(I1,XD(I1),I1=1,NN)
      ENDIF        
C
	DO 26 I=1,N
	X1(I)=0.0D0
	X2(I)=0.0D0
	Z(I)=0.0D0
26	H(I)=8.0D-1
C
	GO TO 5
C
30	CONTINUE
C
	WRITE(NO,31)
31	FORMAT(/5X,'------------------------------------------')
	WRITE(NO,32)ITER,NFUNC
32	FORMAT(10X,'Iterations = ',I4,/,
     +         10X,'No. Function Evaluations = ',I5)
	WRITE(NO,33)EPS
33	FORMAT(10X,'Accuracy =   ',F12.9)
	WRITE(NO,34)
34	FORMAT(5X,'********** Minimum Point **********')
	WRITE(NO,35)(X(J),J=1,N)
35	FORMAT(10X,F18.6)
	WRITE(NO,36)FX
36	FORMAT(10X,'Val. Func. in Minimum Point =',F18.6)
C
	RETURN
	END
C
C***************************Last card of ROSE.FOR