* * $Id: hfumil.F,v 1.1.1.1 1996/01/16 17:07:38 mclareni Exp $ * * $Log: hfumil.F,v $ * Revision 1.1.1.1 1996/01/16 17:07:38 mclareni * First import * * #include "hbook/pilot.h" *CMZ : 4.19/00 19/04/93 10.27.31 by Rene Brun *-- Author : SUBROUTINE HFUMIL(TFUNC,S,M,IT,MC,Z, +Z0,G ,A ,DF ,PL0 ,SIGMA ,PL ,R ,DA,AMX,AMN,EXDA) *.==========> *. HFUMIL IS SPECIAL VERSION OF THE FUMILI *. (D510 CERN LIBRARY) *. THIS VERSION ALLOWS TO USE DYNAMIC MEMORY MANAGEMENT *. AND FITTING FUNCTIONS AS FOLLOWS: *. IFLFUN=1-EXTERNAL FUNCTION TO BE WRITTEN BY USER *. =2-INTERNAL HBOOK FUNCTION I.E. HGAUSS HPOLY HEXP. *. (THESE FUNCTIONS ARE USED INSTEAD OF THE HARITH) *..=========> ( I.Ivanchenko ) DIMENSION Z(1),Z0(1), +G(1),A(1),DF(1),PL0(1),SIGMA(1),PL(1),R(1),DA(1),AMX(1),AMN(1), +EXDA(1) #include "hbook/hcfit1.inc" #include "hbook/hcfit2.inc" #include "hbook/hcunit.inc" #include "hbook/hcflag.inc" #include "hbook/hcrlf.inc" * 10.*MAXIMUM RELATIVE PRECISION * DATA RP/1.E-14/ *.___________________________________________ N1=2 N2=1 N3=113 EPS=0.1 INDFLG(3)=0 IF(IT.GE.0)WRITE(LOUT,84)CRLF, NAMFUN,ID 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. PL(I)=PL0(I) 2 CONTINUE * * START NEW ITERATION * 3 NN1=1 T1=1. * * REPEAT ITERATION WITH SMALLER STEP * 4 S=0. N0=0 DO 7 I=1,N G(I)=0. IF (PL0(I).LE.0.)GO TO 7 N0=N0+1 IF (PL(I).GT.0.)THEN PL0(I)=PL(I) ENDIF 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 * * CALCULATE OBJECTIVE FUNCTION * CALL HSGZ (TFUNC,M,S,Z,G,A,DF,EXDA,PL0,AMX,AMN) SP=RP*ABS(S) IF (NN0.GE.1)THEN DO 10 I=1,NN0 Z0(I)=Z(I) 10 CONTINUE ENDIF IF (NN3.LE.0)GO TO 19 IF (NN1-N1.GT.0)GO TO 19 T=2.*(S-OLDS-GT) IF (INDFLG(1).NE.0)GO TO 16 IF (ABS(S-OLDS).LE.SP.AND.-GT.LE.SP) GO TO 19 IF (0.59*T+GT.LT.0.)GO TO 19 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 * * 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).LE.0.)GO TO 30 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).LE.0)GO TO 28 IF (PL(J).GT.0.)THEN Z(K2)=Z0(K1) K2=K2+1 ENDIF K1=K1+1 28 CONTINUE 29 I1=I1+1 30 CONTINUE * * INVERT Z * I1=1 L=I1 DO 32 I=1,N IF (PL(I).LE.0.)GO TO 32 R(I)=Z(L) I1=I1+1 L=L+I1 32 CONTINUE N0=I1-1 CALL HMCONV (N0,Z,PL,R) IF (INDFLG(1).NE.0)THEN INDFLG(1)=0 INDFLG(2)=1 GO TO 49 ENDIF * * CALCULATE THEORETICAL STEP TO MINIMUM I1=1 DO 41 I=1,N DA(I)=0. IF (PL(I).LE.0.)GO TO 41 L1=1 DO 40 L=1,N IF (PL(L).LE.0.)GO TO 40 IF (I1-L1.LE.0)THEN K=L1*(L1-1)/2+I1 ELSE K=I1*(I1-1)/2+L1 ENDIF DA(I)=DA(I)-G(L)*Z(K) L1=L1+1 40 CONTINUE I1=I1+1 41 CONTINUE * * CHECK FOR PARAMETERS ON BOUNDARY * AFIX=0. IFIX=0 I1=1 L=I1 DO 47 I=1,N IF (PL(I).LE.0.)GO TO 47 SIGI=SQRT(ABS(Z(L))) R(I)=R(I)*Z(L) IF (EPS.GT.0.)THEN SIGMA(I)=SIGI ENDIF IF ((A(I).LT.AMX(I).OR.DA(I).LE.0.).AND. + (A(I).GT.AMN(I).OR.DA(I).GE.0.)) GO TO 46 AKAP=ABS(DA(I)/SIGI) IF (AKAP-AFIX.GT.0.)THEN AFIX=AKAP IFIX=I IFIX1=I ENDIF 46 I1=I1+1 L=L+I1 47 CONTINUE IF (IFIX.EQ.0)GO TO 50 PL(IFIX)=-1. 49 FIXFLG=FIXFLG+1. FI=0. * * REPEAT CALCULATION OF THEORETICAL STEP * AFTER FIXING EACH PARAMETER * GO TO 19 * * CALCULATE STEP CORRECTION FACTOR * 50 ALAMBD=1. AKAPPA=0. IMAX=0 DO 60 I=1,N IF (PL(I).LE.0.)GO TO 60 BM=AMX(I)-A(I) ABI=A(I)+PL(I) ABM=AMX(I) IF (DA(I).LE.0.)THEN BM=A(I)-AMN(I) ABI=A(I)-PL(I) ABM=AMN(I) ENDIF BI=PL(I) IF (BI-BM.GT.0.)THEN BI=BM ABI=ABM ENDIF IF (ABS(DA(I))-BI.GT.0.)THEN AL=ABS(BI/DA(I)) IF (ALAMBD-AL.GT.0.)THEN IMAX=I AIMAX=ABI ALAMBD=AL ENDIF ENDIF AKAP=ABS(DA(I)/SIGMA(I)) IF (AKAP-AKAPPA.LE.0.)GO TO 60 AKAPPA=AKAP 60 CONTINUE * * CALCULATE NEW CORRECTED STEP * GT=0. AMB=1.E+18 IF (ALAMBD.GT.0)THEN AMB=0.25/ALAMBD ENDIF DO 67 I=1,N IF (PL(I).LE.0.)GO TO 67 IF (NN2-N2.GT.0)THEN IF (ABS(DA(I)/PL(I))-AMB.GE.0.)THEN PL(I)=4.*PL(I) T1=4. ENDIF ENDIF DA(I)=DA(I)*ALAMBD GT=GT+DA(I)*G(I) 67 CONTINUE * * 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.LT.0)GO TO 85 IF (AKAPPA-ABS(EPS).GE.0.)GO TO 75 IF (FIXFLG.EQ.0)THEN ENDFLG=1. GO TO 85 ENDIF IF (ENDFLG) 85,77,73 73 IF (IFIX1) 85,85,76 74 IF (FI-FIXFLG) 76,76,77 75 IF (FIXFLG.NE.0)GO TO 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 HMONIT(S,M,NN3,IT,GT,AKAPPA,ALAMBD,A,SIGMA,R,PL,PL0) IF (ENDFLG) 83,79,83 * * 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 * * NEXT ITERATION * 79 ENDFLG=0. DO 80 I=1,N A(I)=A(I)+DA(I) 80 CONTINUE IF (IMAX.GT.0)THEN A(IMAX)=AIMAX ENDIF OLDS=S NN2=NN2+1 NN3=NN3+1 GO TO 3 83 MC=ENDFLG * 84 FORMAT(A, +/' ********************************************************', +/' * *',/ +5X,'* FUNCTION MINIMIZATION BY SUBROUTINE HFIT',A2,' ID=',I6,'*' +/' * *', +/' ********************************************************', +/'0 S = VALUE OF OBJECTIVE FUNCTION ', +/' 2S = CHISQUARE ', +/' EC = EXPECTED CHANGE IN S DURING THE NEXT ITERATION ', +/' KAPPA = ESTIMATED DISTANCE TO MINIMUM ', +/' LAMBDA= STEP LENGTH MODIFIER '//) END