This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / d / d501l2.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:02:19  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10       SUBROUTINE D501L2(K,N,X,NX,MODE,EPS,MAXIT,IPRT,M,A,AL,AU,
11      1                  PHI1,DPHI,IAFR,MFR,STD,P,LAMU,DSCAL,B,W1,W2,
12      2                  AM,COV,SUB,NERROR)
13
14 #include "gen/imp64.inc"
15 #include "gen/def64.inc"
16      +    JP2,LAMBDA,LAMU,LK,MY
17       LOGICAL LFN,LID,LRP,LPR
18       DIMENSION X(*),A(*),AL(*),AU(*),DPHI(*),STD(*),P(*),LAMU(*)
19       DIMENSION DSCAL(*),B(*),W1(*),W2(*),AM(M,*),COV(M,*),IAFR(*)
20       PARAMETER (Z0 = 0, Z1 = 1, HALF = Z1/2, R3 = Z1/3, R10 = Z1/10)
21       PARAMETER (SIG1 = R10, SIG2 = 11*R10, COEF = R10**3, STEP = Z1)
22       PARAMETER (RHO1 = R10**4, RHO2 = Z1/4, RHO3 = 3*Z1/4)
23       EXTERNAL SUB
24
25 ************************************************************************
26 *   LEAMAX, VERSION: 15.03.1993
27 ************************************************************************
28 *
29 *   THIS ROUTINE IS ONE OF THE MAIN ROUTINES OF THE  LEAMAX  PACKAGE
30 *   FOR SOLVING THE PROBLEM OF MAXIMUM LIKELIHOOD ESTIMATION.
31 *   BOUNDS ON THE VARIABLES MAY BE SET.
32 *
33 ************************************************************************
34
35 ************************************************************************
36 *   COMPUTE AN APPROXIMATION  EPS0  TO THE RELATIVE MACHINE PRECISION
37 ************************************************************************
38
39       EPS0=Z1
40    10 EPS0=EPS0/10
41       IF (Z1+EPS0 .NE. Z1) GO TO 10
42       EPS0=10*EPS0
43
44 ************************************************************************
45 *   CHECK THE VALUES OF INPUT PARAMETERS
46 ************************************************************************
47
48       NERROR=0
49
50       CALL D501P1(K,N,M,X,NX,Y,SY,MODE,EPS0,EPS,MAXIT,IPRT,M,A,AL,AU,
51      +            NERROR,'DMAXLK')
52
53       IF(NERROR .NE. 0) RETURN
54
55 ************************************************************************
56 *   SET INITIAL VALUES
57 ************************************************************************
58
59       EPS1=10*EPS0
60
61       LFN=.FALSE.
62       LID=.FALSE.
63       LRP=.FALSE.
64       LPR=IPRT .NE. 0
65
66       ITER=0
67
68       CALL DVSET(M,Z1,DSCAL(1),DSCAL(2))
69       CALL DVSET(M,Z0,LAMU(1),LAMU(2))
70       CALL DVSET(M,Z0,STD(1),STD(2))
71
72 ************************************************************************
73 *   COMPUTE INITIAL VALUE  PHI1  OF OBJECTIVE FUNCTION
74 ************************************************************************
75
76       PHI1=0
77       IX=1
78       DO 20 I=1,N
79       CALL SUB(K,X(IX),M,A,F0,ZZ,0,NERROR)
80       IF(F0 .LE. 0  .OR.  NERROR .NE. 0) THEN
81        NERROR=3
82        RETURN
83       ENDIF
84       PHI1=PHI1-LOG(F0)
85    20 IX=IX+NX
86
87 ************************************************************************
88 *   COMPUTE J(TRANS) X J, B, DPHI, DSCAL, LAMU, MFR, IAFR
89 ************************************************************************
90
91       CALL D501N2(K,N,M,A,AL,AU,X,NX,W1,B,DPHI,DSCAL,LAMU,AM,COV,
92      +            IAFR,MFR,SUB,EPS0,EPS1,MODE,NERROR)
93       IF(NERROR .NE. 0) RETURN
94
95 ************************************************************************
96 *   IF  MFR = 0  MINIMUM IN A CORNER; STOP ITERATION
97 ************************************************************************
98
99       IF(MFR .EQ. 0) GO TO 190
100
101 ************************************************************************
102 *   ITERATION BEGINS
103 ************************************************************************
104
105       DELTA=0
106       LAMBDA=0
107
108 ************************************************************************
109 *   COMPUTE THE L2-NORM OF THE PROJECTED GRADIENT
110 ************************************************************************
111
112       DPHINO=SQRT(DVMPY(MFR,B(1),B(2),B(1),B(2)))
113
114       IF(LPR) CALL D501P2(LRP,M,A,DPHI,STD,LAMU,PHI1,DPHINO,ITER,LFN,
115      +                    MODE,'DMAXLK')
116
117       DA=0
118       DO 40 I=1,MFR
119    40 DA=DA+(DSCAL(I)*A(IAFR(I)))**2
120       DA=SQRT(DA)
121
122 ************************************************************************
123 *   ITERATION WITH GAUSS-NEWTON STEP
124 ************************************************************************
125
126    50 LAMBDA=0
127
128 ************************************************************************
129 *   SOLVE NORMAL EQUATIONS
130 ************************************************************************
131
132       CALL DVSCL(MFR,-Z1,B(1),B(2),W1(1),W1(2))
133       CALL DSINV(MFR,AM,M,NERROR)
134
135       IF(NERROR .NE. 0) THEN
136        NERROR=4
137        RETURN
138       ENDIF
139       CALL DMMPY(MFR,MFR,AM(1,1),AM(1,2),AM(2,1),W1(1),W1(2),P(1),P(2))
140
141
142 ************************************************************************
143 *   COMPUTE THE L2-NORM OF THE SCALED VECTOR  P
144 ************************************************************************
145
146       DP=0
147       DO 60 I=1,MFR
148    60 DP=DP+(DSCAL(I)*P(I))**2
149       DP=SQRT(DP)
150
151 ************************************************************************
152 *   COMPUTE THE STEP SIZE  ALFA
153 ************************************************************************
154
155       ALFA=1
156       DO 70 I=1,MFR
157       IF(P(I) .NE. 0) THEN
158        IF(P(I) .GT. 0) THEN
159         ALFA1=AU(IAFR(I))
160        ELSE
161         ALFA1=AL(IAFR(I))
162        ENDIF
163        ALFA1=(ALFA1-A(IAFR(I)))/P(I)
164        IF(ALFA1 .EQ. 0) THEN
165         P(I)=0
166        ELSE
167         ALFA=MIN(ALFA,ALFA1)
168        ENDIF
169       ENDIF
170    70 CONTINUE
171
172 ************************************************************************
173 *   COMPUTE INITIAL DELTA IF NECESSARY
174 ************************************************************************
175
176       IF(.NOT.LID) THEN
177        DELTA=STEP*MAX(DA,DP/SIG2)
178        LID=.TRUE.
179       ENDIF
180       IF(DELTA .LE. EPS*DA) GO TO 190
181
182 ************************************************************************
183 *   CONTINUATION WITH GAUSS-NEWTON OR SWITCHING TO LEVENBERG-MARQUARDT?
184 ************************************************************************
185
186       IF(DP .GT. SIG2*DELTA) THEN
187
188 ************************************************************************
189 *   DO THE LEVENBERG - MARQUARDT STEP, (HEBDEN'S METHOD).
190 *   - COMPUTE THE LM - PARAMETER LAMBDA
191 *   - COMPUTE THE CORRESPONDING STEP P, ITS DP AND ALFA.
192 ************************************************************************
193
194        CALL DVSCL(MFR,-Z1,B(1),B(2),STD(1),STD(2))
195        UK=0
196        DO 80 I=1,MFR
197    80  UK=UK+(DPHI(IAFR(I))/DSCAL(I))**2
198        UK=SQRT(UK)/DELTA
199
200 ************************************************************************
201 *   COMPUTE INITIAL LAMBDA
202 ************************************************************************
203
204        LAMBDA=COEF*UK
205
206        LK=0
207        ITERA=0
208
209    90  ITERA=ITERA+1
210        IF(ITERA .GE. 50) GO TO 190
211
212 ************************************************************************
213 *   RESET LAMBDA IF NECESSARY
214 ************************************************************************
215
216        IF(LK .GE. LAMBDA .OR. LAMBDA .GE. UK)
217      +    LAMBDA=MAX(COEF*UK,SQRT(LK*UK))
218
219 ************************************************************************
220 *   COMPUTE NEW P FOR NEW LAMBDA
221 ************************************************************************
222
223        CALL DMCPY(MFR,MFR,COV(1,1),COV(1,2),COV(2,1),
224      +                    AM(1,1),AM(1,2),AM(2,1))
225        DO 100 I=1,MFR
226   100  AM(I,I)=AM(I,I)+LAMBDA*DSCAL(I)**2
227
228 ************************************************************************
229 *   SOLVE NORMAL EQUATIONS
230 ************************************************************************
231
232        CALL DSINV(MFR,AM,M,NERROR)
233        IF(NERROR .NE. 0) THEN
234         NERROR=4
235         RETURN
236        ENDIF
237        CALL DMMPY(MFR,MFR,AM(1,1),AM(1,2),AM(2,1),STD(1),STD(2),
238      +                    P(1),P(2))
239
240 ************************************************************************
241 *   COMPUTE THE L2-NORM OF THE SCALED VECTOR  P
242 ************************************************************************
243
244        DP=0
245        DO 110 I=1,MFR
246   110  DP=DP+(DSCAL(I)*P(I))**2
247        DP=SQRT(DP)
248
249 ************************************************************************
250 *   STOP ITERATION IN THE CASE OF NORM EQUAL TO ZERO
251 ************************************************************************
252
253        IF (DP .LE. 0) GO TO 190
254
255        IF(SIG1*DELTA .GT. DP .OR. DP .GT. SIG2*DELTA) THEN
256
257 ************************************************************************
258 *   CONTINUE ITERATION FOR LAMBDA
259 ************************************************************************
260
261         P1=DP-DELTA
262         DO 120 I=1,MFR
263   120   W1(I)=DSCAL(I)**2*P(I)
264         P1P=-DMBIL(MFR,W1(1),W1(2),AM(1,1),AM(1,2),AM(2,1),
265      +                 W1(1),W1(2))/DP
266
267 ************************************************************************
268 *   UPDATE LK, UK, LAMBDA
269 ************************************************************************
270
271         IF(P1 .LT. 0) UK=LAMBDA
272         LK=MAX(LK,LAMBDA-P1/P1P)
273         IF(LK .GE. UK) UK=2*LK
274         LAMBDA=LAMBDA-(DP/DELTA)*(P1/P1P)
275         GO TO 90
276        ENDIF
277       ENDIF
278
279 ************************************************************************
280 *   END OF LEVENBERG - MARQUARDT STEP
281 ************************************************************************
282
283       ALFA=1
284       DO 130 I=1,MFR
285       IF(P(I) .NE. 0) THEN
286        IF(P(I) .GT. 0) THEN
287         ALFA1=AU(IAFR(I))
288        ELSE
289         ALFA1=AL(IAFR(I))
290        ENDIF
291        ALFA1=(ALFA1-A(IAFR(I)))/P(I)
292        IF(ALFA1 .EQ. 0) THEN
293         P(I)=0
294        ELSE
295         ALFA=MIN(ALFA,ALFA1)
296        ENDIF
297       ENDIF
298   130 CONTINUE
299
300 ************************************************************************
301 *   COMPUTE  A + ALPHA * P
302 ************************************************************************
303       CALL DVCPY(M,A(1),A(2),W1(1),W1(2))
304       DO 140 I=1,MFR
305   140 W1(IAFR(I))=A(IAFR(I))+ALFA*P(I)
306
307 ************************************************************************
308 *   COMPUTE VALUE  PHI2  OF OBJECTIVE FUNCTION
309 ************************************************************************
310
311       PHI2=0
312       IX=1
313       DO 150 I=1,N
314       CALL SUB(K,X(IX),M,W1,F0,ZZ,0,NERROR)
315       IF(F0 .LE. 0  .OR.  NERROR .NE. 0) THEN
316        NERROR=3
317        RETURN
318       ENDIF
319       PHI2=PHI2-LOG(F0)
320   150 IX=IX+NX
321
322       PHMAXI=1
323       AAU=MAX(ABS(PHI1),ABS(PHI2))
324       IF(AAU .GT. 0) PHMAXI=1/AAU
325
326       CALL DVSCL(MFR,PHMAXI,P(1),P(2),W2(1),W2(2))
327       JP2=DMBIL(MFR,W2(1),W2(2),COV(1,1),COV(1,2),COV(2,1),W2(1),W2(2))
328
329 ************************************************************************
330 *   COMPUTE THE APPROXIMATION MEASURE  RHO  AND THE UPDATING FACTOR  MY
331 *   FOR  DELTA
332 ************************************************************************
333
334        IF(PHI1 .LE. PHI2) THEN
335         RHO=0
336         MY=R10
337        ELSE
338         S2=2*(PHI1-PHI2)*PHMAXI**2
339         S3=LAMBDA*(DP*PHMAXI)**2
340         RHO=S2/(JP2+2*S3)
341         MY=-JP2-S3
342         S2=S2+2*MY
343         IF(S2 .EQ. 0) THEN
344          MY=R10
345         ELSE
346          MY=MIN(MAX(MY/S2,R10),HALF)
347         ENDIF
348        ENDIF
349
350 ************************************************************************
351 *   END OF COMPUTATTION OF RHO AND MY
352 ************************************************************************
353
354 ************************************************************************
355 *   IF RHO .LE. RHO1, REDUCE DELTA BY FACTOR MY AND MAKE NEW LEVENBERG-
356 *   MARQUARDT STEP, OTHERWISE ACCEPT P
357 ************************************************************************
358
359       IF(RHO .LE. RHO1) THEN
360        DELTA=MY*DELTA
361        DA=0
362        DO 160 I=1,MFR
363   160  DA=DA+(DSCAL(I)*A(IAFR(I)))**2
364        DA=SQRT(DA)
365        GO TO 50
366       ENDIF
367       CALL DVCPY(M,W1(1),W1(2),A(1),A(2))
368       DA=0
369       DO 170 I=1,MFR
370   170 DA=DA+(DSCAL(I)*A(IAFR(I)))**2
371       DA=SQRT(DA)
372       MFROLD=MFR
373
374 ************************************************************************
375 *   COMPUTE J(TRANS) X J, B, DPHI, DSCAL, LAMU, MFR, IAFR
376 ************************************************************************
377
378       CALL D501N2(K,N,M,A,AL,AU,X,NX,W1,B,DPHI,DSCAL,LAMU,AM,COV,
379      +            IAFR,MFR,SUB,EPS0,EPS1,MODE,NERROR)
380       IF(NERROR .NE. 0) RETURN
381
382 ************************************************************************
383 *   IF  MFR = 0  MINIMUM IN A CORNER; STOP ITERATION
384 ************************************************************************
385
386       IF(MFR .EQ. 0) THEN
387        ITER=ITER+1
388        GO TO 190
389       ENDIF
390
391 ************************************************************************
392 *   COMPUTE THE L2-NORM OF THE PROJECTED GRADIENT
393 ************************************************************************
394
395       DPHINO=SQRT(DVMPY(MFR,B(1),B(2),B(1),B(2)))
396
397 ************************************************************************
398 *   TERMINATION CRITERION
399 ************************************************************************
400
401       IF (     PHI2      .LE. PHI1
402      1   .AND. PHI1-PHI2 .LE. EPS*(1+ABS(PHI2))
403      2   .AND. DP        .LE. SQRT(EPS)*(1+DA)
404      3   .AND. DPHINO    .LE. EPS**R3*(1+ABS(PHI2)))    LFN=.TRUE.
405
406       ITER=ITER+1
407       PHI1=PHI2
408
409       IF(.NOT.LFN) THEN
410        IF(ITER .GE. MAXIT) THEN
411         IF(LPR) CALL D501P2(LRP,M,A,DPHI,STD,LAMU,PHI1,DPHINO,ITER,LFN,
412      +                      MODE,'DMAXLK')
413         NERROR=2
414         GO TO 190
415        ENDIF
416
417        IF(LPR) THEN
418           IF(MOD(ITER,IPRT) .EQ. 0)
419      1       CALL D501P2(LRP,M,A,DPHI,STD,LAMU,PHI1,DPHINO,ITER,LFN,
420      2                   MODE,'DMAXLK')
421        ENDIF
422
423 ************************************************************************
424 *   UPDATE DELTA AND GO BACK TO GAUSS-NEWTON STEP
425 ************************************************************************
426
427        IF(MFROLD .NE. MFR) LID=.FALSE.
428        IF(RHO .LE. RHO2) THEN
429         DELTA=MY*DELTA
430        ELSE IF(RHO .GE. RHO3 .OR. LAMBDA .EQ. 0) THEN
431         DELTA=2*DP
432        ENDIF
433
434        GO TO 50
435
436       ENDIF
437
438 ************************************************************************
439 *   END OF ITERATION
440 ************************************************************************
441
442   190 LFN=.TRUE.
443
444 ************************************************************************
445 *   COMPUTE J(TRANS) X J, B, DPHI, DSCAL, LAMU, MFR, IAFR
446 ************************************************************************
447
448       CALL D501N2(K,N,M,A,AL,AU,X,NX,W1,B,DPHI,DSCAL,LAMU,AM,COV,
449      +            IAFR,MFR,SUB,EPS0,EPS1,MODE,MERROR)
450       IF(MERROR .NE. 0) THEN
451        NERROR=MERROR
452        RETURN
453       ENDIF
454
455 ************************************************************************
456 *   COMPUTE THE L2-NORM OF THE PROJECTED GRADIENT
457 ************************************************************************
458
459       DPHINO=SQRT(DVMPY(MFR,B(1),B(2),B(1),B(2)))
460
461 ************************************************************************
462 *   COMPUTE THE VALUE  PHI1  OF THE OBJECTIVE FUNCTION
463 ************************************************************************
464
465       PHI1=0
466       IX=1
467       DO 200 I=1,N
468        CALL SUB(K,X(IX),M,A,F0,ZZ,0,MERROR)
469        IF(F0 .LE. 0  .OR.  MERROR .NE. 0) THEN
470         NERROR=3
471         RETURN
472        ENDIF
473        PHI1=PHI1-LOG(F0)
474   200 IX=IX+NX
475
476
477 ************************************************************************
478 *   PRINT LAST ITERATION RESULTS
479 ************************************************************************
480
481       IF(LPR) CALL D501P2(LRP,M,A,DPHI,STD,LAMU,PHI1,DPHINO,ITER,LFN,
482      +                    MODE,'DMAXLK')
483
484       RETURN
485
486       END
487
488
489