PROGRAM TNELMEAD C Use nelmead to minimize the Griewank function: C this is a driving program for the Nelder Mead simplex algorithm C as implemented in minim.f90 by Alan Miller et al, see his homepage C http://users.bigpond.net.au/amiller C this f77-distribution prepared by Hans Mittelmann (mittelmann@asu.edu) C f(p) = sum(p(i)^2)/4000 - product(cos(p(i))/sqrt(i)) + 1 C global minimum 0 at p = 0 DOUBLE PRECISION OBJECT, P(20), SIMP, STEP(20), STOPCR, VAR(20) INTEGER IER, IPRINT, IQUAD, MAXF, NLOOP, NOP, J LOGICAL FIRST EXTERNAL OBJFUN C NOP is number of variables. NOP <= 20 (see subroutine NELMEAD) NOP = 6 C Set up starting values & step sizes. DO J = 1, NOP P(J) = 1.D0 STEP(J) = 0.5D0 END DO C Set max. no. of function evaluations = 500, print every 20. MAXF = 500 IPRINT = 20 C Set value for stopping criterion. Stopping occurs when the C standard deviation of the values of the objective function at C the points of the current simplex < stopcr. STOPCR = 1.D-5 NLOOP = 6 C Fit a quadratic surface to be sure a minimum has been found. IQUAD = 1 C As function value is being evaluated in DOUBLE PRECISION, it C should be accurate to about 15 decimals. If we set simp = 1.d-6, C we should get about 9 dec. digits accuracy in fitting the surface. SIMP = 1.D-6 C Now call NELMEAD to do the work. FIRST = .TRUE. 1003 CONTINUE IF (FIRST) THEN CALL NELMEAD (P, STEP, NOP, OBJECT, MAXF, IPRINT, STOPCR, 1 NLOOP, IQUAD, SIMP, VAR, OBJFUN, IER) C If ier > 0, try a few more function evaluations. IF (IER .EQ. 0) GO TO 1002 IF (.NOT.FIRST) STOP FIRST = .FALSE. MAXF = 100 GO TO 1003 ENDIF 1002 CONTINUE C Successful termination. WRITE (*, 900) OBJECT 900 FORMAT(' Success C '/' Objective function = ',F12.6/' at: ') WRITE(*,*) (P(J), J=1,NOP) WRITE (*, 910) 910 FORMAT(' Elements of var = ') WRITE(*,*) (VAR(J), J=1,NOP) STOP END SUBROUTINE OBJFUN(NOP, P, FUNC) C This is the subroutine which the user must write. DOUBLE PRECISION P(NOP), FUNC, FAC INTEGER NOP,J FUNC=0.D0 FAC=1.D0 DO J=1,NOP FUNC=FUNC+P(J)*P(J) FAC=FAC*COS(P(J)/SQRT(DFLOAT(J))) ENDDO FUNC=1.D0 + FUNC/4.D3 - FAC RETURN END SUBROUTINE NELMEAD(P,STEP,NOP,FUNC,MAX,IPRINT,STOPCR,NLOOP,IQUAD, 1 SIMP,VAR,FUNCTN,IFAULT) C C A PROGRAM FOR FUNCTION MINIMIZATION USING THE SIMPLEX METHOD. C C FOR DETAILS, SEE NELDER & MEAD, THE COMPUTER JOURNAL, JANUARY 1965 C C PROGRAMMED BY D.E.SHAW, C CSIRO, DIVISION OF MATHEMATICS & STATISTICS C P.O. BOX 218, LINDFIELD, N.S.W. 2070 C C WITH AMENDMENTS BY R.W.M.WEDDERBURN C ROTHAMSTED EXPERIMENTAL STATION C HARPENDEN, HERTFORDSHIRE, ENGLAND C C Further amended by Alan Miller C CSIRO Division of Mathematics & Statistics C Private Bag 10, CLAYTON, VIC. 3168 C C ARGUMENTS:- C P() = INPUT, STARTING VALUES OF PARAMETERS C OUTPUT, FINAL VALUES OF PARAMETERS C STEP() = INPUT, INITIAL STEP SIZES C NOP = INPUT, NO. OF PARAMETERS, INCL. ANY TO BE HELD FIXED C FUNC = OUTPUT, THE FUNCTION VALUE CORRESPONDING TO THE FINAL C PARAMETER VALUES. C MAX = INPUT, THE MAXIMUM NO. OF FUNCTION EVALUATIONS ALLOWED. C Say, 20 times the number of parameters, NOP. C IPRINT = INPUT, PRINT CONTROL PARAMETER C < 0 NO PRINTING C = 0 PRINTING OF PARAMETER VALUES AND THE FUNCTION C VALUE AFTER INITIAL EVIDENCE OF CONVERGENCE. C > 0 AS FOR IPRINT = 0 PLUS PROGRESS REPORTS AFTER C EVERY IPRINT EVALUATIONS, PLUS PRINTING FOR THE C INITIAL SIMPLEX. C STOPCR = INPUT, STOPPING CRITERION. C The criterion is applied to the standard deviation of C the values of FUNC at the points of the simplex. C NLOOP = INPUT, THE STOPPING RULE IS APPLIED AFTER EVERY NLOOP C FUNCTION EVALUATIONS. Normally NLOOP should be slightly C greater than NOP, say NLOOP = 2*NOP. C IQUAD = INPUT, = 1 IF FITTING OF A QUADRATIC SURFACE IS REQUIRED C = 0 IF NOT C N.B. The fitting of a quadratic surface is strongly C recommended, provided that the fitted function is C continuous in the vicinity of the minimum. It is often C a good indicator of whether a premature termination of C the search has occurred. C SIMP = INPUT, CRITERION FOR EXPANDING THE SIMPLEX TO OVERCOME C ROUNDING ERRORS BEFORE FITTING THE QUADRATIC SURFACE. C The simplex is expanded so that the function values at C the points of the simplex exceed those at the supposed C minimum by at least an amount SIMP. C VAR() = OUTPUT, CONTAINS THE DIAGONAL ELEMENTS OF THE INVERSE OF C THE INFORMATION MATRIX. C FUNCTN = INPUT, NAME OF THE USER'S SUBROUTINE - ARGUMENTS C (NOP,P,FUNC) WHICH RETURNS THE FUNCTION VALUE FOR A GIVEN C SET OF PARAMETER VALUES IN ARRAY P. C**** FUNCTN MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM. C IFAULT = OUTPUT, = 0 FOR SUCCESSFUL TERMINATION C = 1 IF MAXIMUM NO. OF FUNCTION EVALUATIONS EXCEEDED C = 2 IF INFORMATION MATRIX IS NOT +VE SEMI-DEFINITE C = 3 IF NOP < 1 C = 4 IF NLOOP < 1 C C N.B. P, STEP AND VAR (IF IQUAD = 1) MUST HAVE DIMENSION AT LEAST NOP C IN THE CALLING PROGRAM. C THE DIMENSIONS BELOW ARE FOR A MAXIMUM OF 20 PARAMETERS. C C LATEST REVISION - 6 April 1985 C C***************************************************************************** C IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION P(NOP),STEP(NOP),VAR(NOP) DIMENSION G(21,20),H(21),PBAR(20),PSTAR(20),PSTST(20),AVAL(20), 1 BMAT(210),PMIN(20),VC(210),TEMP(20) EXTERNAL FUNCTN C C A = REFLECTION COEFFICIENT, B = CONTRACTION COEFFICIENT, AND C C = EXPANSION COEFFICIENT. C DATA A,B,C/1.0D0, 0.5D0, 2.0D0/ C C SET LOUT = LOGICAL UNIT NO. FOR OUTPUT C DATA LOUT/6/ C C IF PROGRESS REPORTS HAVE BEEN REQUESTED, PRINT HEADING C IF(IPRINT.GT.0) WRITE(LOUT,1000) IPRINT 1000 FORMAT(' Progress Report every',I4,' function evaluations'/, 1 ' EVAL. FUNC.VALUE.',10X,'PARAMETER VALUES') C C CHECK INPUT ARGUMENTS C IFAULT = 0 IF(NOP.LE.0) IFAULT = 3 IF(NLOOP.LE.0) IFAULT = 4 IF(IFAULT.NE.0) RETURN C C SET NAP = NO. OF PARAMETERS TO BE VARIED, I.E. WITH STEP.NE.0 C NAP = 0 NEVAL = 0 LOOP = 0 IFLAG = 0 DO 10 I = 1,NOP IF(STEP(I).NE.0.D0) NAP = NAP+1 10 CONTINUE C C IF NAP = 0 EVALUATE FUNCTION AT THE STARTING POINT AND RETURN C IF(NAP.GT.0) GO TO 30 CALL FUNCTN(NOP, P, FUNC) RETURN C C SET UP THE INITIAL SIMPLEX C 30 DO 40 I = 1,NOP 40 G(1,I) = P(I) IROW = 2 DO 60 I = 1,NOP IF(STEP(I).EQ.0.D0) GO TO 60 DO 50 J = 1,NOP 50 G(IROW,J) = P(J) G(IROW,I) = P(I)+STEP(I) IROW = IROW+1 60 CONTINUE C NP1 = NAP+1 DO 90 I = 1,NP1 DO 70 J = 1,NOP 70 P(J) = G(I,J) CALL FUNCTN(NOP, P, H(I)) NEVAL = NEVAL+1 IF(IPRINT.LE.0) GO TO 90 WRITE(LOUT,1010) NEVAL,H(I),(P(J),J=1,NOP) 1010 FORMAT(/1X, I4, 2X, G12.5, 2X, 5G11.4, 3(/21X, 5G11.4)) 90 CONTINUE C C START OF MAIN CYCLE. C C FIND MAX. & MIN. VALUES FOR CURRENT SIMPLEX (HMAX & HMIN). C 100 LOOP = LOOP+1 IMAX = 1 IMIN = 1 HMAX = H(1) HMIN = H(1) DO 120 I = 2,NP1 IF(H(I).LE.HMAX) GO TO 110 IMAX = I HMAX = H(I) GO TO 120 110 IF(H(I).GE.HMIN) GO TO 120 IMIN = I HMIN = H(I) 120 CONTINUE C C FIND THE CENTROID OF THE VERTICES OTHER THAN P(IMAX) C DO 130 I = 1,NOP 130 PBAR(I) = 0.D0 DO 150 I = 1,NP1 IF(I.EQ.IMAX) GO TO 150 DO 140 J = 1,NOP 140 PBAR(J) = PBAR(J)+G(I,J) 150 CONTINUE DO 160 J = 1,NOP 160 PBAR(J) = PBAR(J)/FLOAT(NAP) C C REFLECT MAXIMUM THROUGH PBAR TO PSTAR, C HSTAR = FUNCTION VALUE AT PSTAR. C DO 170 I = 1,NOP 170 PSTAR(I) = A*(PBAR(I)-G(IMAX,I))+PBAR(I) CALL FUNCTN(NOP, PSTAR, HSTAR) NEVAL = NEVAL+1 IF(IPRINT.LE.0) GO TO 180 IF(MOD(NEVAL,IPRINT).EQ.0) WRITE(LOUT,1010) NEVAL,HSTAR, 1 (PSTAR(J),J=1,NOP) C C IF HSTAR < HMIN, REFLECT PBAR THROUGH PSTAR, C HSTST = FUNCTION VALUE AT PSTST. C 180 IF(HSTAR.GE.HMIN) GO TO 220 DO 190 I = 1,NOP 190 PSTST(I) = C*(PSTAR(I)-PBAR(I))+PBAR(I) CALL FUNCTN(NOP, PSTST, HSTST) NEVAL = NEVAL+1 IF(IPRINT.LE.0) GO TO 200 IF(MOD(NEVAL,IPRINT).EQ.0) WRITE(LOUT,1010) NEVAL,HSTST, 1 (PSTST(J),J=1,NOP) C C IF HSTST < HMIN REPLACE CURRENT MAXIMUM POINT BY PSTST AND C HMAX BY HSTST, THEN TEST FOR CONVERGENCE. C 200 IF(HSTST.GE.HMIN) GO TO 320 DO 210 I = 1,NOP IF(STEP(I).NE.0.D0) G(IMAX,I) = PSTST(I) 210 CONTINUE H(IMAX) = HSTST GO TO 340 C C HSTAR IS NOT < HMIN. C TEST WHETHER IT IS < FUNCTION VALUE AT SOME POINT OTHER THAN C P(IMAX). IF IT IS REPLACE P(IMAX) BY PSTAR & HMAX BY HSTAR. C 220 DO 230 I = 1,NP1 IF(I.EQ.IMAX) GO TO 230 IF(HSTAR.LT.H(I)) GO TO 320 230 CONTINUE C C HSTAR > ALL FUNCTION VALUES EXCEPT POSSIBLY HMAX. C IF HSTAR <= HMAX, REPLACE P(IMAX) BY PSTAR & HMAX BY HSTAR. C IF(HSTAR.GT.HMAX) GO TO 260 DO 250 I = 1,NOP IF(STEP(I).NE.0.D0) G(IMAX,I) = PSTAR(I) 250 CONTINUE HMAX = HSTAR H(IMAX) = HSTAR C C CONTRACTED STEP TO THE POINT PSTST, C HSTST = FUNCTION VALUE AT PSTST. C 260 DO 270 I = 1,NOP 270 PSTST(I) = B*G(IMAX,I)+(1.0-B)*PBAR(I) CALL FUNCTN(NOP, PSTST, HSTST) NEVAL = NEVAL+1 IF(IPRINT.LE.0) GO TO 280 IF(MOD(NEVAL,IPRINT).EQ.0) WRITE(LOUT,1010) NEVAL,HSTST, 1 (PSTST(J),J=1,NOP) C C IF HSTST < HMAX REPLACE P(IMAX) BY PSTST & HMAX BY HSTST. C 280 IF(HSTST.GT.HMAX) GO TO 300 DO 290 I = 1,NOP IF(STEP(I).NE.0.D0) G(IMAX,I) = PSTST(I) 290 CONTINUE H(IMAX) = HSTST GO TO 340 C C HSTST > HMAX. C SHRINK THE SIMPLEX BY REPLACING EACH POINT, OTHER THAN THE CURRENT C MINIMUM, BY A POINT MID-WAY BETWEEN ITS CURRENT POSITION AND THE C MINIMUM. C 300 DO 315 I = 1,NP1 IF(I.EQ.IMIN) GO TO 315 DO 310 J = 1,NOP IF(STEP(J).NE.0.D0) G(I,J) = (G(I,J)+G(IMIN,J))*0.5 310 P(J) = G(I,J) CALL FUNCTN(NOP, P, H(I)) NEVAL = NEVAL+1 IF(IPRINT.LE.0) GO TO 315 IF(MOD(NEVAL,IPRINT).EQ.0) WRITE(LOUT,1010) NEVAL,H(I), 1 (P(J),J=1,NOP) 315 CONTINUE GO TO 340 C C REPLACE MAXIMUM POINT BY PSTAR & H(IMAX) BY HSTAR. C 320 DO 330 I = 1,NOP IF(STEP(I).NE.0.D0) G(IMAX,I) = PSTAR(I) 330 CONTINUE H(IMAX) = HSTAR C C IF LOOP = NLOOP TEST FOR CONVERGENCE, OTHERWISE REPEAT MAIN CYCLE. C 340 IF(LOOP.LT.NLOOP) GO TO 100 C C CALCULATE MEAN & STANDARD DEVIATION OF FUNCTION VALUES FOR THE C CURRENT SIMPLEX. C HSTD = 0.D0 HMEAN = 0.D0 DO 350 I = 1,NP1 350 HMEAN = HMEAN+H(I) HMEAN = HMEAN/FLOAT(NP1) DO 360 I = 1,NP1 360 HSTD = HSTD+(H(I)-HMEAN)**2 HSTD = SQRT(HSTD/FLOAT(NP1)) C C IF THE RMS > STOPCR, SET IFLAG & LOOP TO ZERO AND GO TO THE C START OF THE MAIN CYCLE AGAIN. C IF(HSTD.LE.STOPCR.OR.NEVAL.GT.MAX) GO TO 410 IFLAG = 0 LOOP = 0 GO TO 100 C C FIND THE CENTROID OF THE CURRENT SIMPLEX AND THE FUNCTION VALUE THERE. C 410 DO 380 I = 1,NOP IF(STEP(I).EQ.0.D0) GO TO 380 P(I) = 0.D0 DO 370 J = 1,NP1 370 P(I) = P(I)+G(J,I) P(I) = P(I)/FLOAT(NP1) 380 CONTINUE CALL FUNCTN(NOP, P, FUNC) NEVAL = NEVAL+1 IF(IPRINT.LE.0) GO TO 390 IF(MOD(NEVAL,IPRINT).EQ.0) WRITE(LOUT,1010) NEVAL,FUNC, 1 (P(J),J=1,NOP) C C TEST WHETHER THE NO. OF FUNCTION VALUES ALLOWED, MAX, HAS BEEN C OVERRUN; IF SO, EXIT WITH IFAULT = 1. C 390 IF(NEVAL.LE.MAX) GO TO 420 IFAULT = 1 IF(IPRINT.LT.0) RETURN WRITE(LOUT,1020) MAX 1020 FORMAT(' No. of function evaluations > ',I5) WRITE(LOUT,1030) HSTD 1030 FORMAT(' RMS of function values of last simplex =',G14.6) WRITE(LOUT,1040)(P(I),I=1,NOP) 1040 FORMAT(' Centroid of last simplex =',4(/1X,6G13.5)) WRITE(LOUT,1050) FUNC 1050 FORMAT(' Function value at centroid =',G14.6) RETURN C C CONVERGENCE CRITERION SATISFIED. C IF IFLAG = 0, SET IFLAG & SAVE HMEAN. C IF IFLAG = 1 & CHANGE IN HMEAN <= STOPCR THEN SEARCH IS COMPLETE. C 420 IF(IPRINT.LT.0) GO TO 430 WRITE(LOUT,1060) 1060 FORMAT(/' EVIDENCE OF CONVERGENCE') WRITE(LOUT,1040)(P(I),I=1,NOP) WRITE(LOUT,1050) FUNC 430 IF(IFLAG.GT.0) GO TO 450 IFLAG = 1 440 SAVEMN = HMEAN LOOP = 0 GO TO 100 450 IF(ABS(SAVEMN-HMEAN) .GE. STOPCR) GO TO 440 IF(IPRINT.LT.0) GO TO 460 WRITE(LOUT,1070) NEVAL 1070 FORMAT(//' Minimum found after',I5,' function evaluations') WRITE(LOUT,1080)(P(I),I=1,NOP) 1080 FORMAT(' Minimum at',4(/1X,6G13.6)) WRITE(LOUT,1090) FUNC 1090 FORMAT(' Function value at minimum =',G14.6) 460 IF(IQUAD.LE.0) RETURN C C------------------------------------------------------------------ C C QUADRATIC SURFACE FITTING C IF(IPRINT.GE.0) WRITE(LOUT,1110) 1110 FORMAT(/' Fitting quadratic surface about supposed minimum'/) C C EXPAND THE FINAL SIMPLEX, IF NECESSARY, TO OVERCOME ROUNDING C ERRORS. C HMIN = FUNC NMORE = 0 DO 490 I = 1,NP1 470 TEST = ABS(H(I)-FUNC) IF(TEST.GE.SIMP) GO TO 490 DO 480 J = 1,NOP IF(STEP(J).NE.0.D0) G(I,J) = (G(I,J)-P(J))+G(I,J) 480 PSTST(J) = G(I,J) CALL FUNCTN(NOP, PSTST, H(I)) NMORE = NMORE + 1 NEVAL = NEVAL+1 IF(H(I) .GE. HMIN) GO TO 470 HMIN = H(I) IF(IPRINT.GE.0) WRITE(LOUT, 1010) NEVAL, HMIN, + (PSTST(J),J=1,NOP) GO TO 470 490 CONTINUE C C FUNCTION VALUES ARE CALCULATED AT AN ADDITIONAL NAP POINTS. C DO 510 I = 1,NAP I1 = I+1 DO 500 J = 1,NOP 500 PSTAR(J) = (G(1,J)+G(I1,J))*0.5 CALL FUNCTN(NOP, PSTAR, AVAL(I)) NMORE = NMORE + 1 NEVAL = NEVAL+1 510 CONTINUE C C THE MATRIX OF ESTIMATED SECOND DERIVATIVES IS CALCULATED AND ITS C LOWER TRIANGLE STORED IN BMAT. C A0 = H(1) DO 540 I = 1,NAP I1 = I-1 I2 = I+1 IF(I1.LT.1) GO TO 540 DO 530 J = 1,I1 J1 = J+1 DO 520 K = 1,NOP 520 PSTST(K) = (G(I2,K)+G(J1,K))*0.5 CALL FUNCTN(NOP, PSTST, HSTST) NMORE = NMORE + 1 NEVAL = NEVAL+1 L = I*(I-1)/2+J BMAT(L) = 2.0*(HSTST+A0-AVAL(I)-AVAL(J)) 530 CONTINUE 540 CONTINUE L = 0 DO 550 I = 1,NAP I1 = I+1 L = L+I BMAT(L) = 2.0*(H(I1)+A0-2.0*AVAL(I)) 550 CONTINUE C C THE VECTOR OF ESTIMATED FIRST DERIVATIVES IS CALCULATED AND C STORED IN AVAL. C DO 560 I = 1,NAP I1 = I+1 560 AVAL(I) = 2.0*AVAL(I)-(H(I1)+3.0*A0)*0.5 C C THE MATRIX Q OF NELDER & MEAD IS CALCULATED AND STORED IN G. C DO 570 I = 1,NOP 570 PMIN(I) = G(1,I) DO 580 I = 1,NAP I1 = I+1 DO 580 J = 1,NOP G(I1,J) = G(I1,J)-G(1,J) 580 CONTINUE DO 590 I = 1,NAP I1 = I+1 DO 590 J = 1,NOP G(I,J) = G(I1,J) 590 CONTINUE C C INVERT BMAT C CALL SYMINV(BMAT, NAP, BMAT, TEMP, NULLTY, IFAULT, RMAX) IF(IFAULT.NE.0) GO TO 600 IRANK = NAP-NULLTY GO TO 610 600 IF(IPRINT.GE.0) WRITE(LOUT,1120) 1120 FORMAT(/' MATRIX OF ESTIMATED SECOND DERIVATIVES NOT +VE DEFN.'/ 1 ' MINIMUM PROBABLY NOT FOUND'/) IFAULT = 2 IF(NEVAL .GT. MAX) RETURN WRITE(LOUT, 1130) 1130 FORMAT(/10X, 'Search restarting'/) DO 605 I = 1,NOP 605 STEP(I) = 0.5 * STEP(I) GO TO 30 C C BMAT*A/2 IS CALCULATED AND STORED IN H. C 610 DO 650 I = 1,NAP H(I) = 0.D0 DO 640 J = 1,NAP IF(J.GT.I) GO TO 620 L = I*(I-1)/2+J GO TO 630 620 L = J*(J-1)/2+I 630 H(I) = H(I)+BMAT(L)*AVAL(J) 640 CONTINUE 650 CONTINUE C C FIND THE POSITION, PMIN, & VALUE, YMIN, OF THE MINIMUM OF THE C QUADRATIC. C YMIN = 0.D0 DO 660 I = 1,NAP 660 YMIN = YMIN+H(I)*AVAL(I) YMIN = A0-YMIN DO 670 I = 1,NOP PSTST(I) = 0.D0 DO 670 J = 1,NAP 670 PSTST(I) = PSTST(I)+H(J)*G(J,I) DO 680 I = 1,NOP 680 PMIN(I) = PMIN(I)-PSTST(I) IF(IPRINT.LT.0) GO TO 690 WRITE(LOUT,1140) YMIN,(PMIN(I),I=1,NOP) 1140 FORMAT(' Minimum of quadratic surface =',G14.6,' at', 1 4(/1X,6G13.5)) WRITE(LOUT,1150) 1150 FORMAT(' IF THIS DIFFERS BY MUCH FROM THE MINIMUM ESTIMATED', 1 1X,'FROM THE MINIMIZATION,'/ 2 ' THE MINIMUM MAY BE FALSE &/OR THE INFORMATION MATRIX MAY BE', 3 1X,'INACCURATE'/) C C Q*BMAT*Q'/2 IS CALCULATED & ITS LOWER TRIANGLE STORED IN VC C 690 DO 760 I = 1,NOP DO 730 J = 1,NAP H(J) = 0.D0 DO 720 K = 1,NAP IF(K.GT.J) GO TO 700 L = J*(J-1)/2+K GO TO 710 700 L = K*(K-1)/2+J 710 H(J) = H(J)+BMAT(L)*G(K,I)*0.5 720 CONTINUE 730 CONTINUE DO 750 J = I,NOP L = J*(J-1)/2+I VC(L) = 0.D0 DO 740 K = 1,NAP 740 VC(L) = VC(L)+H(K)*G(K,J) 750 CONTINUE 760 CONTINUE C C THE DIAGONAL ELEMENTS OF VC ARE COPIED INTO VAR. C J = 0 DO 770 I = 1,NOP J = J+I 770 VAR(I) = VC(J) IF(IPRINT.LT.0) RETURN WRITE(LOUT,1160) IRANK 1160 FORMAT(' Rank of information matrix =',I3/ 1 ' Inverse of information matrix:-') IJK = 1 GO TO 880 C 790 CONTINUE WRITE(LOUT,1170) 1170 FORMAT(/' If the function minimized was -LOG(LIKELIHOOD),'/ 1 ' this is the covariance matrix of the parameters.'/ 2 ' If the function was a sum of squares of residuals,'/ 3 ' this matrix must be multiplied by twice the estimated', 4 1X, 'residual variance'/' to obtain the covariance matrix.'/) CALL SYMINV(VC, NAP, BMAT, TEMP, NULLTY, IFAULT, RMAX) C C BMAT NOW CONTAINS THE INFORMATION MATRIX C WRITE(LOUT,1190) 1190 FORMAT(' INFORMATION MATRIX:-'/) IJK = 3 GO TO 880 C 800 IJK = 2 II = 0 IJ = 0 DO 840 I = 1,NOP II = II+I IF(VC(II).GT.0.D0) THEN VC(II) = 1.0/SQRT(VC(II)) ELSE VC(II) = 0.D0 END IF JJ = 0 DO 830 J = 1,I-1 JJ = JJ+J IJ = IJ+1 VC(IJ) = VC(IJ)*VC(II)*VC(JJ) 830 CONTINUE IJ = IJ+1 840 CONTINUE C WRITE(LOUT,1200) 1200 FORMAT(//' CORRELATION MATRIX:-') II = 0 DO 850 I = 1,NOP II = II+I IF(VC(II).NE.0.D0) VC(II) = 1.D0 850 CONTINUE GO TO 880 C C Exit, on successful termination. C 860 WRITE(LOUT,1210) NMORE 1210 FORMAT(/' A further',I4,' function evaluations have been used'/) RETURN C 880 L = 1 890 IF(L.GT.NOP) GO TO (790, 860, 800),IJK II = L*(L-1)/2 DO 910 I = L,NOP I1 = II+L II = II+I I2 = MIN0(II,I1+5) IF(IJK.EQ.3) GO TO 900 WRITE(LOUT,1230)(VC(J),J=I1,I2) GO TO 910 900 WRITE(LOUT,1230)(BMAT(J),J=I1,I2) 910 CONTINUE 1230 FORMAT(1X,6G13.5) WRITE(LOUT,1240) 1240 FORMAT(/) L = L+6 GO TO 890 END SUBROUTINE SYMINV(A, N, C, W, NULLTY, IFAULT, RMAX) C C ALGORITHM AS7, APPLIED STATISTICS, VOL.17, 1968. C C ARGUMENTS:- C A() = INPUT, THE SYMMETRIC MATRIX TO BE INVERTED, STORED IN C LOWER TRIANGULAR FORM C N = INPUT, ORDER OF THE MATRIX C C() = OUTPUT, THE INVERSE OF A (A GENERALIZED INVERSE IF C IS C SINGULAR), ALSO STORED IN LOWER TRIANGULAR. C C AND A MAY OCCUPY THE SAME LOCATIONS. C W() = WORKSPACE, DIMENSION AT LEAST N. C NULLTY = OUTPUT, THE RANK DEFICIENCY OF A. C IFAULT = OUTPUT, ERROR INDICATOR C = 1 IF N < 1 C = 2 IF A IS NOT +VE SEMI-DEFINITE C = 0 OTHERWISE C RMAX = OUTPUT, APPROXIMATE BOUND ON THE ACCURACY OF THE DIAGONAL C ELEMENTS OF C. E.G. IF RMAX = 1.E-04 THEN THE DIAGONAL C ELEMENTS OF C WILL BE ACCURATE TO ABOUT 4 DEC. DIGITS. C C LATEST REVISION - 1 April 1985 C C*************************************************************************** C IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION A(1),C(1),W(N) NROW = N IFAULT = 1 IF(NROW.LE.0) GO TO 100 IFAULT = 0 C C CHOLESKY FACTORIZATION OF A, RESULT IN C C CALL CHOLA(A,NROW,C,NULLTY,IFAULT,RMAX,W) IF(IFAULT.NE.0) GO TO 100 C C INVERT C & FORM THE PRODUCT (CINV)'*CINV, WHERE CINV IS THE INVERSE C OF C, ROW BY ROW STARTING WITH THE LAST ROW. C IROW = THE ROW NUMBER, NDIAG = LOCATION OF LAST ELEMENT IN THE ROW. C NN = NROW*(NROW+1)/2 IROW = NROW NDIAG = NN 10 IF(C(NDIAG).EQ.0.D0) GO TO 60 L = NDIAG DO 20 I = IROW,NROW W(I) = C(L) L = L+I 20 CONTINUE ICOL = NROW JCOL = NN MDIAG = NN 30 L = JCOL X = 0.D0 IF(ICOL.EQ.IROW) X = 1.D0/W(IROW) K = NROW 40 IF(K.EQ.IROW) GO TO 50 X = X-W(K)*C(L) K = K-1 L = L-1 IF(L.GT.MDIAG) L = L-K+1 GO TO 40 50 C(L) = X/W(IROW) IF(ICOL.EQ.IROW) GO TO 80 MDIAG = MDIAG-ICOL ICOL = ICOL-1 JCOL = JCOL-1 GO TO 30 60 L = NDIAG DO 70 J = IROW,NROW C(L) = 0.D0 L = L+J 70 CONTINUE 80 NDIAG = NDIAG-IROW IROW = IROW-1 IF(IROW.NE.0) GO TO 10 100 RETURN END SUBROUTINE CHOLA(A, N, U, NULLTY, IFAULT, RMAX, R) C C ALGORITHM AS6, APPLIED STATISTICS, VOL.17, 1968, WITH C MODIFICATIONS BY A.J.MILLER C C ARGUMENTS:- C A() = INPUT, A +VE DEFINITE MATRIX STORED IN LOWER-TRIANGULAR C FORM. C N = INPUT, THE ORDER OF A C U() = OUTPUT, A LOWER TRIANGULAR MATRIX SUCH THAT U*U' = A. C A & U MAY OCCUPY THE SAME LOCATIONS. C NULLTY = OUTPUT, THE RANK DEFICIENCY OF A. C IFAULT = OUTPUT, ERROR INDICATOR C = 1 IF N < 1 C = 2 IF A IS NOT +VE SEMI-DEFINITE C = 0 OTHERWISE C RMAX = OUTPUT, AN ESTIMATE OF THE RELATIVE ACCURACY OF THE C DIAGONAL ELEMENTS OF U. C R() = OUTPUT, ARRAY CONTAINING BOUNDS ON THE RELATIVE ACCURACY C OF EACH DIAGONAL ELEMENT OF U. C C LATEST REVISION - 1 April 1985 C C*************************************************************************** C IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION A(1),U(1),R(N) C ETA=.5D0 110 ETA=ETA/2.D0 IF (1.D0 + ETA .GT. 1.D0) GOTO 110 ETA=ETA*2.D0 IFAULT = 1 IF(N.LE.0) GO TO 100 IFAULT = 2 NULLTY = 0 RMAX = ETA R(1) = ETA J = 1 K = 0 C C FACTORIZE COLUMN BY COLUMN, ICOL = COLUMN NO. C DO 80 ICOL = 1,N L = 0 C C IROW = ROW NUMBER WITHIN COLUMN ICOL C DO 40 IROW = 1,ICOL K = K+1 W = A(K) IF(IROW.EQ.ICOL) RSQ = (W*ETA)**2 M = J DO 10 I = 1,IROW L = L+1 IF(I.EQ.IROW) GO TO 20 W = W-U(L)*U(M) IF(IROW.EQ.ICOL) RSQ = RSQ+(U(L)**2*R(I))**2 M = M+1 10 CONTINUE 20 IF(IROW.EQ.ICOL) GO TO 50 IF(U(L).EQ.0.D0) GO TO 30 U(K) = W/U(L) GO TO 40 30 U(K) = 0.D0 IF(ABS(W).GT.ABS(RMAX*A(K))) GO TO 100 40 CONTINUE C C END OF ROW, ESTIMATE RELATIVE ACCURACY OF DIAGONAL ELEMENT. C 50 RSQ = SQRT(RSQ) IF(ABS(W).LE.5.*RSQ) GO TO 60 IF(W.LT.0.D0) GO TO 100 U(K) = SQRT(W) R(I) = RSQ/W IF(R(I).GT.RMAX) RMAX = R(I) GO TO 70 60 U(K) = 0.D0 NULLTY = NULLTY+1 70 J = J+ICOL 80 CONTINUE IFAULT = 0.D0 100 RETURN END