[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / divon / locsch.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:03:27  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      SUBROUTINE LOCSCH (IIMIN,N,NFREE,IFREE,X,FX,XLOW,XHI,DELTA,SUFTOL,
     1NFCNT,Y,FY,Z,FZ,P)
      INTEGER IIMIN, N, NFREE, NFCNT
      INTEGER IFREE(N)
      DOUBLE PRECISION FX, DELTA, SUFTOL, FY, FZ
      DOUBLE PRECISION X(N), XLOW(N), XHI(N), Y(N), Z(N), P(N)
      INTEGER I, II, ILOC, ITEST
      DOUBLE PRECISION A, B, B1, D, DELX, DFUN, E, EPSMCH, ETA, FA,
     1                 FBEST, FSAV, FTEST, FU, FV, FW, GTEST1, GTEST2,
     2                 GTP, GU, OLDF, PE, PNORM, R, RR, RTEPS,
     3                 SCXBD, SFTBND, SMAX, SNMAX, SOPP, SPMAX, SS,
     4                 SSAV, STP, STPDN, STPNEG, STPPOS, STPUP, T,
     5                 TOL, U, XBEST, XLAMDA, XV, XW, YSAV
      DATA EPSMCH/2.22D-16/
      RTEPS=SQRT(EPSMCH)
      DELX=DELTA
      DO 10 I=1,N
      Y(I)=X(I)
      Z(I)=X(I)
 10   CONTINUE
      FZ=FX
      FY=FX
      SMAX=1.0D+30
      DO 40 I=1,NFREE
      II=IFREE(I)
      STPUP=XHI(II)-X(II)
      STPDN=X(II)-XLOW(II)
      IF(STPUP.LT.STPDN) GOTO 20
      P(I)=1.0D+0
      STP=STPUP
      GOTO 30
 20   P(I)=-1.0D+0
      STP=STPDN
 30   IF(STP.LT.SMAX) SMAX=STP
 40   CONTINUE
      SMAX=0.9D+0*SMAX
 50   IF(DELX.GT.SMAX) DELX=SMAX
      DO 60 I=1,NFREE
      II=IFREE(I)
      Y(II)=X(II)+DELX*P(I)
 60   CONTINUE
      FY=DFUN(N,Y)
      IF(IIMIN.EQ.2) FY=-FY
      NFCNT=NFCNT+1
      FTEST=SUFTOL*(1.0D+0+ABS(FX))
      IF(ABS(FX-FY).GT.FTEST.OR.DELX.EQ.SMAX) GOTO 70
      DELX=5.0D+0*DELX
      GOTO 50
 70   DELX=DELTA
      CALL ORTHVC(N,NFREE,IFREE,X,FX,XLOW,XHI,Y,FY,P,SPMAX,SNMAX)
      SMAX=0.9D+0*SPMAX
      SOPP=0.9D+0*SNMAX
      IF(SMAX.GE.SNMAX.OR.NFREE.EQ.1) GOTO 90
      DO 80 I=1,NFREE
      P(I)=-P(I)
 80   CONTINUE
      SSAV=SMAX
      SMAX=SOPP
      SOPP=SSAV
 90   IF(DELX.GT.SMAX) DELX=SMAX
 100  DO 110 I=1,NFREE
      II=IFREE(I)
      Z(II)=Y(II)+DELX*P(I)
 110  CONTINUE
      FZ=DFUN(N,Z)
      IF(IIMIN.EQ.2) FZ=-FZ
      NFCNT=NFCNT+1
      FTEST=SUFTOL*(1.0D+0+ABS(FY))
      IF(ABS(FY-FZ).GT.FTEST.OR.(5.0D+0*DELX).GT.SMAX) GOTO 120
      DELX=5.0D+0*DELX
      GOTO 100
 120  IF(FY.EQ.FZ) GOTO 210
      IF(FY.LT.FZ) GOTO 130
      GTP=(FZ-FY)/DELX
      XLAMDA=SMAX/0.9D+0
      U=MIN(2.0D+0*DELX,SMAX)
      GOTO 160
 130  DO 140 I=1,NFREE
      P(I)=-P(I)
 140  CONTINUE
      U=MIN(2.0D+0*DELX,SOPP)
      XLAMDA=SOPP/0.9D+0+DELX
      GTP=(FY-FZ)/DELX
      DO 150 I=1,N
      YSAV=Y(I)
      Y(I)=Z(I)
      Z(I)=YSAV
 150  CONTINUE
      FSAV=FY
      FY=FZ
      FZ=FSAV
 160  CALL RLEN(NFREE,P,PNORM)
      PE=PNORM+RTEPS
      ILOC=1
      FU=FY
      GU=GTP
      SFTBND=0.0D+0
      ETA=RTEPS
      T=RTEPS/PE
 170  CALL NEWPTQ(RTEPS,T,ETA,SFTBND,XLAMDA,U,FU,GU,XBEST,FBEST,XW,FW,X
     1V,FV,A,FA,B,OLDF,B1,SCXBD,E,D,RR,SS,GTEST1,GTEST2,TOL,ILOC,ITEST)
      IF(ITEST.NE.1) GOTO 190
      R=XBEST+U
      DO 180 I=1,NFREE
      II=IFREE(I)
      Z(II)=Y(II)+R*P(I)
 180  CONTINUE
      FU=DFUN(N,Z)
      IF(IIMIN.EQ.2) FU=-FU
      NFCNT=NFCNT+1
      FZ=FU
      GOTO 170
 190  IF(ITEST.NE.0) GOTO 210
      DO 200 I=1,NFREE
      II=IFREE(I)
      Z(II)=Y(II)+XBEST*P(I)
 200  CONTINUE
      FZ=FBEST
 210  IF(FZ.EQ.FX) RETURN
      DO 220 I=1,NFREE
      II=IFREE(I)
      P(I)=Z(II)-X(II)
 220  CONTINUE
      SPMAX=1.0D+30
      SNMAX=1.0D+30
      DO 230 I=1,NFREE
      II=IFREE(I)
      CALL MXSTEP(X(II),XLOW(II),XHI(II),P(I),STPPOS,STPNEG)
      IF(STPPOS.LT.SPMAX) SPMAX=STPPOS
      IF(STPNEG.LT.SNMAX) SNMAX=STPNEG
 230  CONTINUE
      CALL RLEN(NFREE,P,PNORM)
      GTP=FZ-FX
      U=MIN(2.0D+0,0.9D+0*SPMAX)
      XLAMDA=SPMAX
      IF(FZ.LT.FX) GOTO 270
      DELX=DELTA/(PNORM+RTEPS)
      IF(DELX.GT.0.9D+0*SNMAX) DELX=0.9D+0*SNMAX
      DO 240 I=1,NFREE
      II=IFREE(I)
      Y(II)=X(II)-DELX*P(I)
 240  CONTINUE
      FY=DFUN(N,Y)
      IF(IIMIN.EQ.2) FY=-FY
      NFCNT=NFCNT+1
      IF(FY.LT.FX) GOTO 250
      IF(DELX.LT.1.0D+0) GTP=(FX-FY)/DELX
      GOTO 270
 250  GTP=(FY-FX)/DELX
      DO 260 I=1,NFREE
      P(I)=-P(I)
 260  CONTINUE
      U=MIN(2.0D+0*DELX,0.9D+0*SNMAX)
      XLAMDA=SNMAX
 270  ILOC=1
      PE=PNORM+RTEPS
      FU=FX
      GU=GTP
      SFTBND=0.0D+0
      ETA=RTEPS
      T=RTEPS/PE
 280  CALL NEWPTQ(RTEPS,T,ETA,SFTBND,XLAMDA,U,FU,GU,XBEST,FBEST,XW,FW,X
     1V,FV,A,FA,B,OLDF,B1,SCXBD,E,D,RR,SS,GTEST1,GTEST2,TOL,ILOC,ITEST)
      IF(ITEST.NE.1) GOTO 300
      R=XBEST+U
      DO 290 I=1,NFREE
      II=IFREE(I)
      Y(II)=X(II)+R*P(I)
 290  CONTINUE
      FU=DFUN(N,Y)
      IF(IIMIN.EQ.2) FU=-FU
      NFCNT=NFCNT+1
      GOTO 280
 300  IF(ITEST.NE.0) RETURN
      DO 310 I=1,NFREE
      II=IFREE(I)
      Z(II)=X(II)+XBEST*P(I)
 310  CONTINUE
      FZ=FBEST
      RETURN
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1996/04/01 15:03:27 mclareni
	6	* Mathlib gen
	7	*
	8	*
	9	#include "gen/pilot.h"
	10	SUBROUTINE LOCSCH (IIMIN,N,NFREE,IFREE,X,FX,XLOW,XHI,DELTA,SUFTOL,
	11	1NFCNT,Y,FY,Z,FZ,P)
	12	INTEGER IIMIN, N, NFREE, NFCNT
	13	INTEGER IFREE(N)
	14	DOUBLE PRECISION FX, DELTA, SUFTOL, FY, FZ
	15	DOUBLE PRECISION X(N), XLOW(N), XHI(N), Y(N), Z(N), P(N)
	16	INTEGER I, II, ILOC, ITEST
	17	DOUBLE PRECISION A, B, B1, D, DELX, DFUN, E, EPSMCH, ETA, FA,
	18	1 FBEST, FSAV, FTEST, FU, FV, FW, GTEST1, GTEST2,
	19	2 GTP, GU, OLDF, PE, PNORM, R, RR, RTEPS,
	20	3 SCXBD, SFTBND, SMAX, SNMAX, SOPP, SPMAX, SS,
	21	4 SSAV, STP, STPDN, STPNEG, STPPOS, STPUP, T,
	22	5 TOL, U, XBEST, XLAMDA, XV, XW, YSAV
	23	DATA EPSMCH/2.22D-16/
	24	RTEPS=SQRT(EPSMCH)
	25	DELX=DELTA
	26	DO 10 I=1,N
	27	Y(I)=X(I)
	28	Z(I)=X(I)
	29	10 CONTINUE
	30	FZ=FX
	31	FY=FX
	32	SMAX=1.0D+30
	33	DO 40 I=1,NFREE
	34	II=IFREE(I)
	35	STPUP=XHI(II)-X(II)
	36	STPDN=X(II)-XLOW(II)
	37	IF(STPUP.LT.STPDN) GOTO 20
	38	P(I)=1.0D+0
	39	STP=STPUP
	40	GOTO 30
	41	20 P(I)=-1.0D+0
	42	STP=STPDN
	43	30 IF(STP.LT.SMAX) SMAX=STP
	44	40 CONTINUE
	45	SMAX=0.9D+0*SMAX
	46	50 IF(DELX.GT.SMAX) DELX=SMAX
	47	DO 60 I=1,NFREE
	48	II=IFREE(I)
	49	Y(II)=X(II)+DELX*P(I)
	50	60 CONTINUE
	51	FY=DFUN(N,Y)
	52	IF(IIMIN.EQ.2) FY=-FY
	53	NFCNT=NFCNT+1
	54	FTEST=SUFTOL*(1.0D+0+ABS(FX))
	55	IF(ABS(FX-FY).GT.FTEST.OR.DELX.EQ.SMAX) GOTO 70
	56	DELX=5.0D+0*DELX
	57	GOTO 50
	58	70 DELX=DELTA
	59	CALL ORTHVC(N,NFREE,IFREE,X,FX,XLOW,XHI,Y,FY,P,SPMAX,SNMAX)
	60	SMAX=0.9D+0*SPMAX
	61	SOPP=0.9D+0*SNMAX
	62	IF(SMAX.GE.SNMAX.OR.NFREE.EQ.1) GOTO 90
	63	DO 80 I=1,NFREE
	64	P(I)=-P(I)
65	80 CONTINUE
66	SSAV=SMAX
67	SMAX=SOPP
68	SOPP=SSAV
69	90 IF(DELX.GT.SMAX) DELX=SMAX
70	100 DO 110 I=1,NFREE
71	II=IFREE(I)
72	Z(II)=Y(II)+DELX*P(I)
73	110 CONTINUE
74	FZ=DFUN(N,Z)
75	IF(IIMIN.EQ.2) FZ=-FZ
76	NFCNT=NFCNT+1
77	FTEST=SUFTOL*(1.0D+0+ABS(FY))
78	IF(ABS(FY-FZ).GT.FTEST.OR.(5.0D+0*DELX).GT.SMAX) GOTO 120
79	DELX=5.0D+0*DELX
80	GOTO 100
81	120 IF(FY.EQ.FZ) GOTO 210
82	IF(FY.LT.FZ) GOTO 130
83	GTP=(FZ-FY)/DELX
84	XLAMDA=SMAX/0.9D+0
85	U=MIN(2.0D+0*DELX,SMAX)
86	GOTO 160
87	130 DO 140 I=1,NFREE
88	P(I)=-P(I)
89	140 CONTINUE
90	U=MIN(2.0D+0*DELX,SOPP)
91	XLAMDA=SOPP/0.9D+0+DELX
92	GTP=(FY-FZ)/DELX
93	DO 150 I=1,N
94	YSAV=Y(I)
95	Y(I)=Z(I)
96	Z(I)=YSAV
97	150 CONTINUE
98	FSAV=FY
99	FY=FZ
100	FZ=FSAV
101	160 CALL RLEN(NFREE,P,PNORM)
102	PE=PNORM+RTEPS
103	ILOC=1
104	FU=FY
105	GU=GTP
106	SFTBND=0.0D+0
107	ETA=RTEPS
108	T=RTEPS/PE
109	170 CALL NEWPTQ(RTEPS,T,ETA,SFTBND,XLAMDA,U,FU,GU,XBEST,FBEST,XW,FW,X
110	1V,FV,A,FA,B,OLDF,B1,SCXBD,E,D,RR,SS,GTEST1,GTEST2,TOL,ILOC,ITEST)
111	IF(ITEST.NE.1) GOTO 190
112	R=XBEST+U
113	DO 180 I=1,NFREE
114	II=IFREE(I)
115	Z(II)=Y(II)+R*P(I)
116	180 CONTINUE
117	FU=DFUN(N,Z)
118	IF(IIMIN.EQ.2) FU=-FU
119	NFCNT=NFCNT+1
120	FZ=FU
121	GOTO 170
122	190 IF(ITEST.NE.0) GOTO 210
123	DO 200 I=1,NFREE
124	II=IFREE(I)
125	Z(II)=Y(II)+XBEST*P(I)
126	200 CONTINUE
127	FZ=FBEST
128	210 IF(FZ.EQ.FX) RETURN
129	DO 220 I=1,NFREE
130	II=IFREE(I)
131	P(I)=Z(II)-X(II)
132	220 CONTINUE
133	SPMAX=1.0D+30
134	SNMAX=1.0D+30
135	DO 230 I=1,NFREE
136	II=IFREE(I)
137	CALL MXSTEP(X(II),XLOW(II),XHI(II),P(I),STPPOS,STPNEG)
138	IF(STPPOS.LT.SPMAX) SPMAX=STPPOS
139	IF(STPNEG.LT.SNMAX) SNMAX=STPNEG
140	230 CONTINUE
141	CALL RLEN(NFREE,P,PNORM)
142	GTP=FZ-FX
143	U=MIN(2.0D+0,0.9D+0*SPMAX)
144	XLAMDA=SPMAX
145	IF(FZ.LT.FX) GOTO 270
146	DELX=DELTA/(PNORM+RTEPS)
147	IF(DELX.GT.0.9D+0SNMAX) DELX=0.9D+0SNMAX
148	DO 240 I=1,NFREE
149	II=IFREE(I)
150	Y(II)=X(II)-DELX*P(I)
151	240 CONTINUE
152	FY=DFUN(N,Y)
153	IF(IIMIN.EQ.2) FY=-FY
154	NFCNT=NFCNT+1
155	IF(FY.LT.FX) GOTO 250
156	IF(DELX.LT.1.0D+0) GTP=(FX-FY)/DELX
157	GOTO 270
158	250 GTP=(FY-FX)/DELX
159	DO 260 I=1,NFREE
160	P(I)=-P(I)
161	260 CONTINUE
162	U=MIN(2.0D+0DELX,0.9D+0SNMAX)
163	XLAMDA=SNMAX
164	270 ILOC=1
165	PE=PNORM+RTEPS
166	FU=FX
167	GU=GTP
168	SFTBND=0.0D+0
169	ETA=RTEPS
170	T=RTEPS/PE
171	280 CALL NEWPTQ(RTEPS,T,ETA,SFTBND,XLAMDA,U,FU,GU,XBEST,FBEST,XW,FW,X
172	1V,FV,A,FA,B,OLDF,B1,SCXBD,E,D,RR,SS,GTEST1,GTEST2,TOL,ILOC,ITEST)
173	IF(ITEST.NE.1) GOTO 300
174	R=XBEST+U
175	DO 290 I=1,NFREE
176	II=IFREE(I)
177	Y(II)=X(II)+R*P(I)
178	290 CONTINUE
179	FU=DFUN(N,Y)
180	IF(IIMIN.EQ.2) FU=-FU
181	NFCNT=NFCNT+1
182	GOTO 280
183	300 IF(ITEST.NE.0) RETURN
184	DO 310 I=1,NFREE
185	II=IFREE(I)
186	Z(II)=X(II)+XBEST*P(I)
187	310 CONTINUE
188	FZ=FBEST
189	RETURN
190	END