C*******************************************************
C
C     Searching method by Hook-Jeeves.
C     
C*******************************************************
C     Neculai Andrei, 
C     Research Institute for Informatics
C     Advanced Modeling and Optimization Laboratory 
C*******************************************************
C
	SUBROUTINE HOOKJ(N,X1,X2,Y,Z,S,ITECH)
C
	DOUBLE PRECISION X1(N),X2(N),Y(N),Z(N),S(N)
	DOUBLE PRECISION F,FP,FM,BL,FY,F2,H,EPS,A
C
	CHARACTER*1 FD
C
      COMMON /PREC/EPS
      COMMON /TYPE/NO
      COMMON /WRIT/NX
	COMMON /L3A/ICOD
	COMMON /NOF/NBL
C
	WRITE(NX,150)
150	FORMAT(10X,'Start Optimization process by Hook-Jeeves Method')
C
	GO TO (151,152,153,154,155)ITECH
151	WRITE(NX,161)
161	FORMAT(10X,'Using Golden Searching Technique')
	GO TO 170
152	WRITE(NX,162)
162	FORMAT(10X,'Using Fibanacci Searching Technique')
	GO TO 170
153	WRITE(NX,163)
163	FORMAT(10X,'Using Quadratic Interpolation of Powell')
	GO TO 170
154	WRITE(NX,164)
164	FORMAT(10X,'Using Dichotomous Searching Technique')
	GO TO 170
155	WRITE(NX,165)
165	FORMAT(10X,'Using Lambda* = 1.0')
C
170	CONTINUE
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
	ICOD=0         
	JTYPE=0
	H=8.0D-1
	NFUNC=0
	NBL=0
C
	K=1
	CALL FUNC(N,X1,F)
	NFUNC=NFUNC+1
C
	IF(NO.EQ.6) THEN
	OPEN(6,FILE='PRINT.DAT',ACCESS='SEQUENTIAL',
     +         FORM='FORMATTED')
	ENDIF
C
	WRITE(NO,50)
50	FORMAT(10X,'*****Optimization Method of HOOK-JEEVES.*****')
	IF(ITECH.EQ.1) WRITE(NO,161)
	IF(ITECH.EQ.2) WRITE(NO,162)
	IF(ITECH.EQ.3) WRITE(NO,163)
	IF(ITECH.EQ.4) WRITE(NO,164)
	IF(ITECH.EQ.5) WRITE(NO,165)
C
	WRITE(NO,51)
51	FORMAT(10X,'Initial Point :')
	WRITE(NO,52)(X1(I),I=1,N)
52	FORMAT(10X,F12.6)
	WRITE(NO,53) H
53	FORMAT(10X,'Step H        = ',F12.6)
	WRITE(NO,54)EPS
54	FORMAT(10X,'Precision EPS = ',g30.15)
	WRITE(NO,55)F
55	FORMAT(10X,'Val. Func. in Initial point = ',F18.6,/)
C
              IF(JTYPE.EQ.1) WRITE(NO,56)K
56	        FORMAT(5X,'Iteration K =',I4)
20	CONTINUE
	        IF(JTYPE.EQ.1) WRITE(NO,72)
72	        FORMAT(10X,'--- Initialization ---')
	DO 1 J=1,N
1	Y(J)=X1(J)
C
10	DO 5 I=1,N
              IF(JTYPE.EQ.1) WRITE(NO,49)I
49	        FORMAT(8X,'I = ',I3)
	DO 2 J=1,N
2	Z(J)=Y(J)
C
	CALL FUNC(N,Z,F)
	NFUNC=NFUNC+1
              IF(JTYPE.EQ.1) WRITE(NO,57)
57	        FORMAT(10X,'Bazic Point')
              IF(JTYPE.EQ.1) WRITE(NO,52)(Z(J),J=1,N)
              IF(JTYPE.EQ.1) WRITE(NO,58)F
58	        FORMAT(10X,'Function value in the basic point = ',F12.6)
	Z(I)=Z(I)+H
C
	CALL FUNC(N,Z,FP)
	NFUNC=NFUNC+1
C
              IF(JTYPE.EQ.1) WRITE(NO,59)
59	        FORMAT(10X,'Pozitive Point')
              IF(JTYPE.EQ.1) WRITE(NO,52)(Z(J),J=1,N)
	        IF(JTYPE.EQ.1) WRITE(NO,60)FP
60	        FORMAT(10X,'Function value in pozitive point = ',F12.6)
	IF(FP.LT.F) GO TO 3
	Z(I)=Z(I)-2.*H
C
	CALL FUNC(N,Z,FM)
	NFUNC=NFUNC+1
C
              IF(JTYPE.EQ.1) WRITE(NO,61)
61	        FORMAT(10X,'Negative Point')
              IF(JTYPE.EQ.1) WRITE(NO,52)(Z(J),J=1,N)
              IF(JTYPE.EQ.1) WRITE(NO,62)FM
62	        FORMAT(10X,'Function value in negative point = ',F12.6)
	IF(FM.LT.F.AND.F.LT.FP) GO TO 3
	A=FP
	IF(FM.LT.A) A=FM
	IF(F.LT.A) GO TO 5
3	DO 4 J=1,N
4	Y(J)=Z(J)
5	CONTINUE
C
C
	IS=0
	DO 6 I=1,N
	IF(X1(I).NE.Y(I)) IS=1
6	CONTINUE
	IF(IS.EQ.1) THEN
	DO 7 I=1,N
7	X2(I)=Y(I)
	GO TO 8
	ELSE
	IF(H.LT.EPS) GO TO 100
	H=H/2.
              IF(JTYPE.EQ.1) WRITE(NO,71)H
71	        FORMAT(10X,'*********** Value of H = ',F12.6)
	GO TO 10
	ENDIF
C
C
8	DO 9 I=1,N
9	S(I)=X2(I)-X1(I)
              IF(JTYPE.EQ.1) WRITE(NO,63)
63	        FORMAT(10X,'Direction S :')
              IF(JTYPE.EQ.1) WRITE(NO,52)(S(J),J=1,N)
C
C----------------------------------------------------------------
	GO TO (171,172,173,174,175) ITECH
171	CALL L1(N,X2,S,BL)
	GO TO 180
172	CALL L2(N,X2,S,BL)
	IF(EPS.EQ.0.0) RETURN
	GO TO 180
173	CALL L3(N,X2,S,BL)
	IF(ICOD.EQ.1) RETURN
	GO TO 180
174	CALL L4(N,X2,S,BL)
	GO TO 180
175	BL=1.0D0
	IF(EPS.EQ.0.0) EPS=0.0000001
C
180	CONTINUE
C---------------------------------------------------------------
C
              IF(JTYPE.EQ.1) WRITE(NO,64)BL
64	        FORMAT(12X,'Steplength =',F12.6)
C
	NFUNC=NFUNC+NBL
C
	DO 11 I=1,N
11	Y(I)=X2(I)+BL*S(I)
              IF(JTYPE.EQ.1)WRITE(NO,65)
65	        FORMAT(10X,'New Point: Y = X + lambda*S')
              IF(JTYPE.EQ.1) WRITE(NO,52)(Y(J),J=1,N)

	K=K+1
	IF(K.GT.1000) THEN
	WRITE(NX,201)
201	FORMAT(10X,'Too many iterations.',/,
     +         10X,'This function is too complicated for the power',/,
     +         10X,'of Hook-Jeeves 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
C
              IF(JTYPE.EQ.1) WRITE(NO,66) K
66            FORMAT(5X,'Iteration K = ',I4)
C
	CALL FUNC(N,Y,F)
	NFUNC=NFUNC+1
C
              IF(JTYPE.EQ.1) WRITE(NO,67)F
67            FORMAT(10X,'Function value in the new point = ',F12.6)
	DO 15 I=1,N
              IF(JTYPE.EQ.1) WRITE(NO,68)I
68            FORMAT(8X,'I = ',I3)
	DO 12 J=1,N
12	Z(J)=Y(J)
C
	CALL FUNC(N,Z,F)
	NFUNC=NFUNC+1
C
              IF(JTYPE.EQ.1) WRITE(NO,57)
              IF(JTYPE.EQ.1) WRITE(NO,52)(Z(J),J=1,N)
              IF(JTYPE.EQ.1) WRITE(NO,58)F
	Z(I)=Z(I)+H
C
	CALL FUNC(N,Z,FP)
	NFUNC=NFUNC+1
C
              IF(JTYPE.EQ.1) WRITE(NO,59)
              IF(JTYPE.EQ.1) WRITE(NO,52)(Z(J),J=1,N)
              IF(JTYPE.EQ.1) WRITE(NO,60)FP
	IF(FP.LT.F) GO TO  13
	Z(I)=Z(I)-2.*H
C
	CALL FUNC(N,Z,FM)
	NFUNC=NFUNC+1
C
              IF(JTYPE.EQ.1) WRITE(NO,61)
              IF(JTYPE.EQ.1) WRITE(NO,52)(Z(J),J=1,N)
              IF(JTYPE.EQ.1) WRITE(NO,62)FM
	IF(FM.LT.F.AND.F.LT.FP) GO TO 13
	A=FP
	IF(FM.LT.A) A=FM
	IF(F.LT.A) GO TO 15
13	DO 14 J=1,N
14	Y(J)=Z(J)
15	CONTINUE
C
C
	CALL FUNC(N,X2,F2)
	CALL FUNC(N,Y,FY)
	NFUNC=NFUNC+2
C
              IF(JTYPE.EQ.1) WRITE(NO,69)F2
69            FORMAT(12X,'Function value in X =',F12.6)
              IF(JTYPE.EQ.1) WRITE(NO,70)FY
70            FORMAT(12X,'Function value in Y =',F12.6)
C
C
	IF(FY.LT.F2) THEN
	DO 16 I=1,N
16	X1(I)=X2(I)
	DO 17 I=1,N
17	X2(I)=Y(I)
C
	IF(DABS(FY-F2).LE.EPS) THEN
	IF(H.LE.EPS) GO TO 100
	H=H/2.
	ENDIF
C
	GO TO 8
	ELSE
	IF(H.LE.EPS) GO TO 100
	H=H/2.
              IF(JTYPE.EQ.1) WRITE(NO,71)H
	DO 18 I=1,N
18	X1(I)=X2(I)
	GO TO 20
	ENDIF
C
C
100	CONTINUE
C
	WRITE(NO,79)
79	FORMAT(/5X,'----------------------------------------')
	WRITE(NO,80)K,NFUNC
80	FORMAT(10X,'Iterations = ',I4,/,
     +         10X,'No. Function Evaluations = ',I5)
	WRITE(NO,81)EPS
81	FORMAT(10X,'Precision = ',F12.9)
	WRITE(NO,82)
82	FORMAT(5X,'********** Minimum Point *********')
	WRITE(NO,83)(Z(J),J=1,N)
83	FORMAT(10X,F12.6)
	WRITE(NO,84)F
84	FORMAT(10X,'Val. Func. in Minimum Point =',F12.6)
C
C
	RETURN
	END
C
C******************************* Last card  of HOOKJ.FOR
