+++ /dev/null
-*
-* $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