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