[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / d / fumili.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:21  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      SUBROUTINE FUMILI (S,M,N1,N2,N3,EPS,AKAPPA,ALAMBD,IT,MC)
C-----ENTRY FOR CHISQ MINIMISATION
#include "d510pl.inc"
#include "d510si.inc"
#include "d510ui.inc"
#include "d510uo.inc"
C-----10.*MAXIMUM RELATIVE PRECISION ON CDC 6000
      DATA RP/1.E-14/
      INDFLG(3)=0
 1    IF (IT.GE.0) WRITE(6,84)
#if defined(CERNLIB_IBM)||defined(CERNLIB_CDC)
      CALL VFILL(AMX,100,1.0E75)
#endif
#if (!defined(CERNLIB_CDC))&&(!defined(CERNLIB_IBM))
      CALL VFILL(AMX,100,1.0E37)
#endif
      CALL VCOPYN(AMX,AMN,100)
      NN2=0
      N=M
      FIXFLG=0.
      ENDFLG=0.
      INDFLG(2)=0
      IFIX1=0.
      FI=0.
      NN3=0
      DO 2 I=1,N
      R(I)=0.
      IF (EPS.GT.0.) SIGMA(I)=0.
 2    PL(I)=PL0(I)
C-----START NEW ITERATION
 3    NN1=1
      T1=1.
C-----REPEAT ITERATION WITH SMALLER STEP
 4    S=0.
      N0=0
      DO 7 I=1,N
      G(I)=0.
      IF (PL0(I)) 7,7,5
 5    N0=N0+1
      IF (PL(I)) 7,7,6
 6    PL0(I)=PL(I)
 7    CONTINUE
      NN0=N0*(N0+1)/2
      IF (NN0.LT.1) GO TO 9
      DO 8 I=1,NN0
      Z(I)=0.
 8    CONTINUE
 9    NA=M
      INDFLG(1)=0
C-----CALCULATE OBJECTIVE FUNCTION
      CALL SGZ (M,S)
      SP=RP*ABS(S)
      IF (NN0.LT.1) GO TO 11
      DO 10 I=1,NN0
      Z0(I)=Z(I)
 10   CONTINUE
 11   IF (NN3) 19,19,12
 12   IF (NN1-N1) 13,13,19
 13   T=2.*(S-OLDS-GT)
      IF (INDFLG(1)) 16,14,16
 14   IF (ABS(S-OLDS).LE.SP.AND.-GT.LE.SP) GO TO 19
      IF (0.59*T+GT) 19,15,15
 15   T=-GT/T
      IF (T-0.25) 16,17,17
 16   T=0.25
 17   GT=GT*T
      T1=T1*T
      NN2=0
      DO 18 I=1,N
      IF (PL(I).LE.0.) GO TO 18
      A(I)=A(I)-DA(I)
      PL(I)=PL(I)*T
      DA(I)=DA(I)*T
      A(I)=A(I)+DA(I)
 18   CONTINUE
      NN1=NN1+1
      GO TO 4
C-----REMOVE CONTRIBUTION OF FIXED PARAMETERS FROM Z
   19 IF(INDFLG(1).EQ.0) GO TO 20
      ENDFLG=-4.
      GO TO 85
 20   K1=1
      K2=1
      I1=1
      DO 30 I=1,N
      IF (PL0(I)) 30,30,21
 21   IF (PL(I).EQ.0.) PL(I)=PL0(I)
      IF (PL(I)) 23,23,24
 22   PL(I)=0.
 23   K1=K1+I1
      GO TO 29
 24   IF (A(I).GE.AMX(I).AND.G(I).LT.0.) GO TO 22
      IF (A(I).LE.AMN(I).AND.G(I).GT.0.) GO TO 22
      DO 28 J=1,I
      IF (PL0(J)) 28,28,25
 25   IF (PL(J)) 27,27,26
 26   Z(K2)=Z0(K1)
      K2=K2+1
 27   K1=K1+1
 28   CONTINUE
 29   I1=I1+1
 30   CONTINUE
C-----INVERT Z
      I1=1
      L=I1
      DO 32 I=1,N
      IF (PL(I)) 32,32,31
 31   R(I)=Z(L)
      I1=I1+1
      L=L+I1
 32   CONTINUE
      N0=I1-1
      CALL MCONV (N0)
      IF (INDFLG(1)) 33,34,33
 33   INDFLG(1)=0
      INDFLG(2)=1
      GO TO 49
 34   CONTINUE
C-----CALCULATE THEORETICAL STEP TO MINIMUM
      I1=1
      DO 41 I=1,N
      DA(I)=0.
      IF (PL(I)) 41,41,35
 35   L1=1
      DO 40 L=1,N
      IF (PL(L)) 40,40,36
 36   IF (I1-L1) 37,37,38
 37   K=L1*(L1-1)/2+I1
      GO TO 39
 38   K=I1*(I1-1)/2+L1
 39   DA(I)=DA(I)-G(L)*Z(K)
      L1=L1+1
 40   CONTINUE
      I1=I1+1
 41   CONTINUE
C-----CHECK FOR PARAMETERS ON BOUNDARY
      AFIX=0.
      IFIX=0
      I1=1
      L=I1
      DO 47 I=1,N
      IF (PL(I)) 47,47,42
 42   SIGI=SQRT(ABS(Z(L)))
      R(I)=R(I)*Z(L)
      IF (EPS) 44,44,43
 43   SIGMA(I)=SIGI
 44   IF ((A(I).LT.AMX(I).OR.DA(I).LE.0.).AND.(A(I).GT.AMN(I).OR.DA(I).G
     1E.0.)) GO TO 46
      AKAP=ABS(DA(I)/SIGI)
      IF (AKAP-AFIX) 46,46,45
 45   AFIX=AKAP
      IFIX=I
      IFIX1=I
 46   I1=I1+1
      L=L+I1
 47   CONTINUE
      IF (IFIX) 48,50,48
 48   PL(IFIX)=-1.
 49   FIXFLG=FIXFLG+1.
      FI=0.
C-----REPEAT CALCULATION OF THEORETICAL STEP AFTER FIXING EACH PARAMETER
      GO TO 19
C-----CALCULATE STEP CORRECTION FACTOR
 50   ALAMBD=1.
      AKAPPA=0.
      IMAX=0
      DO 60 I=1,N
      IF (PL(I)) 60,60,51
 51   BM=AMX(I)-A(I)
      ABI=A(I)+PL(I)
      ABM=AMX(I)
      IF (DA(I)) 52,52,53
 52   BM=A(I)-AMN(I)
      ABI=A(I)-PL(I)
      ABM=AMN(I)
 53   BI=PL(I)
      IF (BI-BM) 55,55,54
 54   BI=BM
      ABI=ABM
 55   IF (ABS(DA(I))-BI) 58,58,56
 56   AL=ABS(BI/DA(I))
      IF (ALAMBD-AL) 58,58,57
 57   IMAX=I
      AIMAX=ABI
      ALAMBD=AL
 58   AKAP=ABS(DA(I)/SIGMA(I))
      IF (AKAP-AKAPPA) 60,60,59
 59   AKAPPA=AKAP
 60   CONTINUE
C-----CALCULATE NEW CORRECTED STEP
      GT=0.
      AMB=1.E18
      IF (ALAMBD) 62,62,61
 61   AMB=0.25/ALAMBD
 62   CONTINUE
      DO 67 I=1,N
      IF (PL(I)) 67,67,63
 63   IF (NN2-N2) 66,66,64
 64   IF (ABS(DA(I)/PL(I))-AMB) 66,65,65
 65   PL(I)=4.*PL(I)
      T1=4.
 66   DA(I)=DA(I)*ALAMBD
      GT=GT+DA(I)*G(I)
 67   CONTINUE
C-----CHECK IF MINIMUM ATTAINED AND SET EXIT MODE
      IF (-GT.GT.SP.OR.T1.GE.1..OR.ALAMBD.GE.1.) GO TO 68
      ENDFLG=-1.
 68   IF (ENDFLG) 85,69,69
 69   IF (AKAPPA-ABS(EPS)) 70,75,75
 70   IF (FIXFLG) 72,71,72
 71   ENDFLG=1.
      GO TO 85
 72   IF (ENDFLG) 85,77,73
 73   IF (IFIX1) 85,85,76
 74   IF (FI-FIXFLG) 76,76,77
 75   IF (FIXFLG) 74,76,74
 76   FI=FI+1.
      ENDFLG=0.
 85   IF(ENDFLG.EQ.0..AND.NN3.GE.N3) ENDFLG=-3.
      IF(ENDFLG.GT.0..AND.INDFLG(2).GT.0) ENDFLG=-2.
      CALL MONITO (S,M,NN3,IT,EPS,GT,AKAPPA,ALAMBD)
      IF (ENDFLG) 83,79,83
C-----CHECK IF FIXING ON BOUND IS CORRECT
 77   ENDFLG=1.
      FIXFLG=0.
      IFIX1=0
      DO 78 I=1,M
 78   PL(I)=PL0(I)
      INDFLG(2)=0
      GO TO 19
C-----NEXT ITERATION
 79   ENDFLG=0.
      DO 80 I=1,N
      A(I)=A(I)+DA(I)
 80   CONTINUE
      IF (IMAX) 82,82,81
 81   A(IMAX)=AIMAX
 82   OLDS=S
      NN2=NN2+1
      NN3=NN3+1
      GO TO 3
 83   MC=ENDFLG
      RETURN
C-----ENTRY FOR MAXIMUM LIKLEHOOD
#if (defined(CERNLIB_CDC))&&(defined(CERNLIB_F4))
      ENTRY LIKELM
#endif
#if !defined(CERNLIB_CDC)||!defined(CERNLIB_F4)
      ENTRY LIKELM(S,M,N1,N2,N3,EPS,AKAPPA,ALAMBD,IT,MC)
#endif
      INDFLG(3)=1
      GO TO 1
C
 84   FORMAT('1',43X,'FUNCTION MINIMISATION BY SUBROUTINE FUMILI/LIKE',
     +'LM'/'0',55X,'IN THE FOLLOWING PRINT-OUT'/
     +     '0',27X,'S = VALUE OF OBJECTIVE FUNCTION,',
     + 'EC = EXPECTED CHANGE IN S DURING NEXT ITERATION'/
     +   '0',34X,'KAPPA = ESTIMATED DISTANCE TO MINIMUM,  LAMBDA =',
     + 'STEP LENGTH MODIFIER'///)
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1996/04/01 15:02:21 mclareni
	6	* Mathlib gen
	7	*
	8	*
	9	#include "gen/pilot.h"
	10	SUBROUTINE FUMILI (S,M,N1,N2,N3,EPS,AKAPPA,ALAMBD,IT,MC)
	11	C-----ENTRY FOR CHISQ MINIMISATION
	12	#include "d510pl.inc"
	13	#include "d510si.inc"
	14	#include "d510ui.inc"
	15	#include "d510uo.inc"
	16	C-----10.*MAXIMUM RELATIVE PRECISION ON CDC 6000
	17	DATA RP/1.E-14/
	18	INDFLG(3)=0
	19	1 IF (IT.GE.0) WRITE(6,84)
	20	#if defined(CERNLIB_IBM)\|\|defined(CERNLIB_CDC)
	21	CALL VFILL(AMX,100,1.0E75)
	22	#endif
	23	#if (!defined(CERNLIB_CDC))&&(!defined(CERNLIB_IBM))
	24	CALL VFILL(AMX,100,1.0E37)
	25	#endif
	26	CALL VCOPYN(AMX,AMN,100)
	27	NN2=0
	28	N=M
	29	FIXFLG=0.
	30	ENDFLG=0.
	31	INDFLG(2)=0
	32	IFIX1=0.
	33	FI=0.
	34	NN3=0
	35	DO 2 I=1,N
	36	R(I)=0.
	37	IF (EPS.GT.0.) SIGMA(I)=0.
	38	2 PL(I)=PL0(I)
	39	C-----START NEW ITERATION
	40	3 NN1=1
	41	T1=1.
	42	C-----REPEAT ITERATION WITH SMALLER STEP
	43	4 S=0.
	44	N0=0
	45	DO 7 I=1,N
	46	G(I)=0.
	47	IF (PL0(I)) 7,7,5
	48	5 N0=N0+1
	49	IF (PL(I)) 7,7,6
	50	6 PL0(I)=PL(I)
	51	7 CONTINUE
	52	NN0=N0*(N0+1)/2
	53	IF (NN0.LT.1) GO TO 9
	54	DO 8 I=1,NN0
	55	Z(I)=0.
	56	8 CONTINUE
	57	9 NA=M
	58	INDFLG(1)=0
	59	C-----CALCULATE OBJECTIVE FUNCTION
	60	CALL SGZ (M,S)
	61	SP=RP*ABS(S)
	62	IF (NN0.LT.1) GO TO 11
	63	DO 10 I=1,NN0
	64	Z0(I)=Z(I)
65	10 CONTINUE
66	11 IF (NN3) 19,19,12
67	12 IF (NN1-N1) 13,13,19
68	13 T=2.*(S-OLDS-GT)
69	IF (INDFLG(1)) 16,14,16
70	14 IF (ABS(S-OLDS).LE.SP.AND.-GT.LE.SP) GO TO 19
71	IF (0.59*T+GT) 19,15,15
72	15 T=-GT/T
73	IF (T-0.25) 16,17,17
74	16 T=0.25
75	17 GT=GT*T
76	T1=T1*T
77	NN2=0
78	DO 18 I=1,N
79	IF (PL(I).LE.0.) GO TO 18
80	A(I)=A(I)-DA(I)
81	PL(I)=PL(I)*T
82	DA(I)=DA(I)*T
83	A(I)=A(I)+DA(I)
84	18 CONTINUE
85	NN1=NN1+1
86	GO TO 4
87	C-----REMOVE CONTRIBUTION OF FIXED PARAMETERS FROM Z
88	19 IF(INDFLG(1).EQ.0) GO TO 20
89	ENDFLG=-4.
90	GO TO 85
91	20 K1=1
92	K2=1
93	I1=1
94	DO 30 I=1,N
95	IF (PL0(I)) 30,30,21
96	21 IF (PL(I).EQ.0.) PL(I)=PL0(I)
97	IF (PL(I)) 23,23,24
98	22 PL(I)=0.
99	23 K1=K1+I1
100	GO TO 29
101	24 IF (A(I).GE.AMX(I).AND.G(I).LT.0.) GO TO 22
102	IF (A(I).LE.AMN(I).AND.G(I).GT.0.) GO TO 22
103	DO 28 J=1,I
104	IF (PL0(J)) 28,28,25
105	25 IF (PL(J)) 27,27,26
106	26 Z(K2)=Z0(K1)
107	K2=K2+1
108	27 K1=K1+1
109	28 CONTINUE
110	29 I1=I1+1
111	30 CONTINUE
112	C-----INVERT Z
113	I1=1
114	L=I1
115	DO 32 I=1,N
116	IF (PL(I)) 32,32,31
117	31 R(I)=Z(L)
118	I1=I1+1
119	L=L+I1
120	32 CONTINUE
121	N0=I1-1
122	CALL MCONV (N0)
123	IF (INDFLG(1)) 33,34,33
124	33 INDFLG(1)=0
125	INDFLG(2)=1
126	GO TO 49
127	34 CONTINUE
128	C-----CALCULATE THEORETICAL STEP TO MINIMUM
129	I1=1
130	DO 41 I=1,N
131	DA(I)=0.
132	IF (PL(I)) 41,41,35
133	35 L1=1
134	DO 40 L=1,N
135	IF (PL(L)) 40,40,36
136	36 IF (I1-L1) 37,37,38
137	37 K=L1*(L1-1)/2+I1
138	GO TO 39
139	38 K=I1*(I1-1)/2+L1
140	39 DA(I)=DA(I)-G(L)*Z(K)
141	L1=L1+1
142	40 CONTINUE
143	I1=I1+1
144	41 CONTINUE
145	C-----CHECK FOR PARAMETERS ON BOUNDARY
146	AFIX=0.
147	IFIX=0
148	I1=1
149	L=I1
150	DO 47 I=1,N
151	IF (PL(I)) 47,47,42
152	42 SIGI=SQRT(ABS(Z(L)))
153	R(I)=R(I)*Z(L)
154	IF (EPS) 44,44,43
155	43 SIGMA(I)=SIGI
156	44 IF ((A(I).LT.AMX(I).OR.DA(I).LE.0.).AND.(A(I).GT.AMN(I).OR.DA(I).G
157	1E.0.)) GO TO 46
158	AKAP=ABS(DA(I)/SIGI)
159	IF (AKAP-AFIX) 46,46,45
160	45 AFIX=AKAP
161	IFIX=I
162	IFIX1=I
163	46 I1=I1+1
164	L=L+I1
165	47 CONTINUE
166	IF (IFIX) 48,50,48
167	48 PL(IFIX)=-1.
168	49 FIXFLG=FIXFLG+1.
169	FI=0.
170	C-----REPEAT CALCULATION OF THEORETICAL STEP AFTER FIXING EACH PARAMETER
171	GO TO 19
172	C-----CALCULATE STEP CORRECTION FACTOR
173	50 ALAMBD=1.
174	AKAPPA=0.
175	IMAX=0
176	DO 60 I=1,N
177	IF (PL(I)) 60,60,51
178	51 BM=AMX(I)-A(I)
179	ABI=A(I)+PL(I)
180	ABM=AMX(I)
181	IF (DA(I)) 52,52,53
182	52 BM=A(I)-AMN(I)
183	ABI=A(I)-PL(I)
184	ABM=AMN(I)
185	53 BI=PL(I)
186	IF (BI-BM) 55,55,54
187	54 BI=BM
188	ABI=ABM
189	55 IF (ABS(DA(I))-BI) 58,58,56
190	56 AL=ABS(BI/DA(I))
191	IF (ALAMBD-AL) 58,58,57
192	57 IMAX=I
193	AIMAX=ABI
194	ALAMBD=AL
195	58 AKAP=ABS(DA(I)/SIGMA(I))
196	IF (AKAP-AKAPPA) 60,60,59
197	59 AKAPPA=AKAP
198	60 CONTINUE
199	C-----CALCULATE NEW CORRECTED STEP
200	GT=0.
201	AMB=1.E18
202	IF (ALAMBD) 62,62,61
203	61 AMB=0.25/ALAMBD
204	62 CONTINUE
205	DO 67 I=1,N
206	IF (PL(I)) 67,67,63
207	63 IF (NN2-N2) 66,66,64
208	64 IF (ABS(DA(I)/PL(I))-AMB) 66,65,65
209	65 PL(I)=4.*PL(I)
210	T1=4.
211	66 DA(I)=DA(I)*ALAMBD
212	GT=GT+DA(I)*G(I)
213	67 CONTINUE
214	C-----CHECK IF MINIMUM ATTAINED AND SET EXIT MODE
215	IF (-GT.GT.SP.OR.T1.GE.1..OR.ALAMBD.GE.1.) GO TO 68
216	ENDFLG=-1.
217	68 IF (ENDFLG) 85,69,69
218	69 IF (AKAPPA-ABS(EPS)) 70,75,75
219	70 IF (FIXFLG) 72,71,72
220	71 ENDFLG=1.
221	GO TO 85
222	72 IF (ENDFLG) 85,77,73
223	73 IF (IFIX1) 85,85,76
224	74 IF (FI-FIXFLG) 76,76,77
225	75 IF (FIXFLG) 74,76,74
226	76 FI=FI+1.
227	ENDFLG=0.
228	85 IF(ENDFLG.EQ.0..AND.NN3.GE.N3) ENDFLG=-3.
229	IF(ENDFLG.GT.0..AND.INDFLG(2).GT.0) ENDFLG=-2.
230	CALL MONITO (S,M,NN3,IT,EPS,GT,AKAPPA,ALAMBD)
231	IF (ENDFLG) 83,79,83
232	C-----CHECK IF FIXING ON BOUND IS CORRECT
233	77 ENDFLG=1.
234	FIXFLG=0.
235	IFIX1=0
236	DO 78 I=1,M
237	78 PL(I)=PL0(I)
238	INDFLG(2)=0
239	GO TO 19
240	C-----NEXT ITERATION
241	79 ENDFLG=0.
242	DO 80 I=1,N
243	A(I)=A(I)+DA(I)
244	80 CONTINUE
245	IF (IMAX) 82,82,81
246	81 A(IMAX)=AIMAX
247	82 OLDS=S
248	NN2=NN2+1
249	NN3=NN3+1
250	GO TO 3
251	83 MC=ENDFLG
252	RETURN
253	C-----ENTRY FOR MAXIMUM LIKLEHOOD
254	#if (defined(CERNLIB_CDC))&&(defined(CERNLIB_F4))
255	ENTRY LIKELM
256	#endif
257	#if !defined(CERNLIB_CDC)\|\|!defined(CERNLIB_F4)
258	ENTRY LIKELM(S,M,N1,N2,N3,EPS,AKAPPA,ALAMBD,IT,MC)
259	#endif
260	INDFLG(3)=1
261	GO TO 1
262	C
263	84 FORMAT('1',43X,'FUNCTION MINIMISATION BY SUBROUTINE FUMILI/LIKE',
264	+'LM'/'0',55X,'IN THE FOLLOWING PRINT-OUT'/
265	+ '0',27X,'S = VALUE OF OBJECTIVE FUNCTION,',
266	+ 'EC = EXPECTED CHANGE IN S DURING NEXT ITERATION'/
267	+ '0',34X,'KAPPA = ESTIMATED DISTANCE TO MINIMUM, LAMBDA =',
268	+ 'STEP LENGTH MODIFIER'///)
269	END