MINICERN/mathlib/gen/d/fumili.F

   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