SUBROUTINE MINIMIZE_2D(X0,FCN,X,F,Info) IMPLICIT NONE ! Arguments REAL*8, INTENT(INOUT) :: X0(2),X(2),F INTEGER, INTENT(INOUT) :: Info EXTERNAL FCN ! Local Variables INTEGER Nop, MaxF, IPrint, NLoop,IQuad REAL*8 Step(2),StopCryterion,Simplex,Var(2) ! Unpack and call the minimizer X(1) = X0(1) X(2) = X0(2) Step(1) = 1d-2 Step(2) = 1d-2 Nop = 2 MaxF = 400 IPrint = -1 StopCryterion = 1.d-10 NLoop = 4 IQuad = 1 Simplex = 1.d-4 CALL minim(X, Step, Nop, F, MaxF, IPrint, StopCryterion, NLoop, & IQuad, Simplex, Var, FCN, Info) END SUBROUTINE MINIMIZE_2D SUBROUTINE MINIMIZE_1D(X0,FCN,X,F,Info) IMPLICIT NONE ! Arguments REAL*8, INTENT(INOUT) :: X0,X,F INTEGER, INTENT(INOUT) :: Info EXTERNAL FCN ! Local Variables INTEGER Nop, MaxF, IPrint, NLoop,IQuad REAL*8 Step,StopCryterion,Simplex,Var ! Unpack and call the minimizer X = X0 Step = 0.4D0 Nop = 1 MaxF = 400 IPrint = -1 StopCryterion = 1.d-10 NLoop = 4 IQuad = 0 Simplex = 1.d-5 CALL minim(X, Step, Nop, F, MaxF, IPrint, StopCryterion, NLoop, & IQuad, Simplex, Var, FCN, Info) END SUBROUTINE MINIMIZE_1D SUBROUTINE minim(p, step, nop, func, maxfn, iprint, stopcr, nloop, iquad, & simp, var, functn, ifault) ! A PROGRAM FOR FUNCTION MINIMIZATION USING THE SIMPLEX METHOD. ! FOR DETAILS, SEE NELDER & MEAD, THE COMPUTER JOURNAL, JANUARY 1965 ! PROGRAMMED BY D.E.SHAW, ! CSIRO, DIVISION OF MATHEMATICS & STATISTICS ! P.O. BOX 218, LINDFIELD, N.S.W. 2070 ! WITH AMENDMENTS BY R.W.M.WEDDERBURN ! ROTHAMSTED EXPERIMENTAL STATION ! HARPENDEN, HERTFORDSHIRE, ENGLAND ! Further amended by Alan Miller ! CSIRO Division of Mathematics & Statistics ! Private Bag 10, CLAYTON, VIC. 3169 ! Fortran 90 conversion by Alan Miller, June 1995 ! Latest revision - 14 September 1995 ! ARGUMENTS:- ! P() = INPUT, STARTING VALUES OF PARAMETERS ! OUTPUT, FINAL VALUES OF PARAMETERS ! STEP() = INPUT, INITIAL STEP SIZES ! NOP = INPUT, NO. OF PARAMETERS, INCL. ANY TO BE HELD FIXED ! FUNC = OUTPUT, THE FUNCTION VALUE CORRESPONDING TO THE FINAL ! PARAMETER VALUES. ! maxfn = INPUT, THE MAXIMUM NO. OF FUNCTION EVALUATIONS ALLOWED. ! Say, 20 times the number of parameters, NOP. ! IPRINT = INPUT, PRINT CONTROL PARAMETER ! < 0 NO PRINTING ! = 0 PRINTING OF PARAMETER VALUES AND THE FUNCTION ! VALUE AFTER INITIAL EVIDENCE OF CONVERGENCE. ! > 0 AS FOR IPRINT = 0 PLUS PROGRESS REPORTS AFTER ! EVERY IPRINT EVALUATIONS, PLUS PRINTING FOR THE ! INITIAL SIMPLEX. ! STOPCR = INPUT, STOPPING CRITERION. ! The criterion is applied to the standard deviation of ! the values of FUNC at the points of the simplex. ! NLOOP = INPUT, THE STOPPING RULE IS APPLIED AFTER EVERY NLOOP ! FUNCTION EVALUATIONS. Normally NLOOP should be slightly ! greater than NOP, say NLOOP = 2*NOP. ! IQUAD = INPUT, = 1 IF FITTING OF A QUADRATIC SURFACE IS REQUIRED ! = 0 IF NOT ! N.B. The fitting of a quadratic surface is strongly ! recommended, provided that the fitted function is ! continuous in the vicinity of the minimum. It is often ! a good indicator of whether a premature termination of ! the search has occurred. ! SIMP = INPUT, CRITERION FOR EXPANDING THE SIMPLEX TO OVERCOME ! ROUNDING ERRORS BEFORE FITTING THE QUADRATIC SURFACE. ! The simplex is expanded so that the function values at ! the points of the simplex exceed those at the supposed ! minimum by at least an amount SIMP. ! VAR() = OUTPUT, CONTAINS THE DIAGONAL ELEMENTS OF THE INVERSE OF ! THE INFORMATION MATRIX. ! FUNCTN = INPUT, NAME OF THE USER'S SUBROUTINE - ARGUMENTS (P,FUNC) ! WHICH RETURNS THE FUNCTION VALUE FOR A GIVEN SET OF ! PARAMETER VALUES IN ARRAY P. !**** FUNCTN MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM. ! IFAULT = OUTPUT, = 0 FOR SUCCESSFUL TERMINATION ! = 1 IF MAXIMUM NO. OF FUNCTION EVALUATIONS EXCEEDED ! = 2 IF INFORMATION MATRIX IS NOT +VE SEMI-DEFINITE ! = 3 IF NOP < 1 ! = 4 IF NLOOP < 1 ! N.B. P, STEP AND VAR (IF IQUAD = 1) MUST HAVE DIMENSION AT LEAST NOP ! IN THE CALLING PROGRAM. !***************************************************************************** IMPLICIT DOUBLE PRECISION (a-h,o-z) INTEGER, INTENT(IN) :: nop, maxfn, iprint, nloop, iquad INTEGER, INTENT(OUT) :: ifault DOUBLE PRECISION, INTENT(IN) :: stopcr, simp DOUBLE PRECISION, INTENT(INOUT) :: p(nop), step(nop) DOUBLE PRECISION, INTENT(OUT) :: var(nop), func EXTERNAL functn ! Local variables DOUBLE PRECISION :: g(nop+1,nop), h(nop+1), pbar(nop), pstar(nop), & pstst(nop), aval(nop), pmin(nop), temp(nop), & bmat(nop*(nop+1)/2), vc(nop*(nop+1)/2), & zero = 0.d0, half = 0.5d0, one = 1.d0, two = 2.d0 ! A = REFLECTION COEFFICIENT, B = CONTRACTION COEFFICIENT, AND ! C = EXPANSION COEFFICIENT. DOUBLE PRECISION, PARAMETER :: a = 1.d0, b = 0.5d0, c = 2.d0 ! SET LOUT = LOGICAL UNIT NO. FOR OUTPUT INTEGER, PARAMETER :: lout = 6 ! IF PROGRESS REPORTS HAVE BEEN REQUESTED, PRINT HEADING IF (iprint.GT.0) WRITE (lout,5000) iprint ! CHECK INPUT ARGUMENTS ifault = 0 IF (nop.LE.0) ifault = 3 IF (nloop.LE.0) ifault = 4 IF (ifault.NE.0) RETURN ! SET NAP = NO. OF PARAMETERS TO BE VARIED, I.E. WITH STEP.NE.0 nap = COUNT(step.NE.zero) neval = 0 loop = 0 iflag = 0 ! IF NAP = 0 EVALUATE FUNCTION AT THE STARTING POINT AND RETURN IF (nap.LE.0) THEN CALL functn(p,func) RETURN END IF ! SET UP THE INITIAL SIMPLEX 20 g(1,:) = p irow = 2 DO i = 1, nop IF (step(i).NE.zero) THEN g(irow,:) = p g(irow,i) = p(i) + step(i) irow = irow + 1 END IF END DO np1 = nap + 1 DO i = 1, np1 p = g(i,:) CALL functn(p,h(i)) neval = neval + 1 IF (iprint.GT.0) THEN WRITE (lout,5100) neval, h(i), p END IF END DO ! START OF MAIN CYCLE. ! FIND MAX. & MIN. VALUES FOR CURRENT SIMPLEX (HMAX & HMIN). Main_loop: DO loop = loop + 1 imax = 1 imin = 1 hmax = h(1) hmin = h(1) DO i = 2, np1 IF (h(i).GT.hmax) THEN imax = i hmax = h(i) ELSE IF (h(i).LT.hmin) THEN imin = i hmin = h(i) END IF END IF END DO ! FIND THE CENTROID OF THE VERTICES OTHER THAN P(IMAX) pbar = zero DO i = 1, np1 IF (i.NE.imax) THEN pbar = pbar + g(i,:) END IF END DO pbar = pbar / FLOAT(nap) ! REFLECT MAXIMUM THROUGH PBAR TO PSTAR, ! HSTAR = FUNCTION VALUE AT PSTAR. pstar = a * (pbar - g(imax,:)) + pbar CALL functn(pstar,hstar) neval = neval + 1 IF (iprint.GT.0) THEN IF (MOD(neval,iprint).EQ.0) WRITE (lout,5100) neval, hstar, pstar END IF ! IF HSTAR < HMIN, REFLECT PBAR THROUGH PSTAR, ! HSTST = FUNCTION VALUE AT PSTST. IF (hstar.LT.hmin) THEN pstst = c * (pstar - pbar) + pbar CALL functn(pstst,hstst) neval = neval + 1 IF (iprint.GT.0) THEN IF (MOD(neval,iprint).EQ.0) WRITE (lout,5100) neval, hstst, pstst END IF ! IF HSTST < HMIN REPLACE CURRENT MAXIMUM POINT BY PSTST AND ! HMAX BY HSTST, THEN TEST FOR CONVERGENCE. IF (hstst.GE.hmin) THEN ! REPLACE MAXIMUM POINT BY PSTAR & H(IMAX) BY HSTAR. g(imax,:) = pstar h(imax) = hstar ELSE g(imax,:) = pstst h(imax) = hstst END IF GO TO 250 END IF ! HSTAR IS NOT < HMIN. ! TEST WHETHER IT IS < FUNCTION VALUE AT SOME POINT OTHER THAN ! P(IMAX). IF IT IS REPLACE P(IMAX) BY PSTAR & HMAX BY HSTAR. DO i = 1, np1 IF (i.NE.imax) THEN IF (hstar.LT.h(i)) THEN ! REPLACE MAXIMUM POINT BY PSTAR & H(IMAX) BY HSTAR. g(imax,:) = pstar h(imax) = hstar GO TO 250 END IF END IF END DO ! HSTAR > ALL FUNCTION VALUES EXCEPT POSSIBLY HMAX. ! IF HSTAR <= HMAX, REPLACE P(IMAX) BY PSTAR & HMAX BY HSTAR. IF (hstar.LE.hmax) THEN g(imax,:) = pstar hmax = hstar h(imax) = hstar END IF ! CONTRACTED STEP TO THE POINT PSTST, ! HSTST = FUNCTION VALUE AT PSTST. pstst = b * g(imax,:) + (one-b) * pbar CALL functn(pstst,hstst) neval = neval + 1 IF (iprint.GT.0) THEN IF (MOD(neval,iprint).EQ.0) WRITE (lout,5100) neval, hstst, pstst END IF ! IF HSTST < HMAX REPLACE P(IMAX) BY PSTST & HMAX BY HSTST. IF (hstst.LE.hmax) THEN g(imax,:) = pstst h(imax) = hstst GO TO 250 END IF ! HSTST > HMAX. ! SHRINK THE SIMPLEX BY REPLACING EACH POINT, OTHER THAN THE CURRENT ! MINIMUM, BY A POINT MID-WAY BETWEEN ITS CURRENT POSITION AND THE ! MINIMUM. DO i = 1, np1 IF (i.NE.imin) THEN DO j = 1, nop IF (step(j).NE.zero) g(i,j) = (g(i,j) + g(imin,j)) * half p(j) = g(i,j) END DO CALL functn(p,h(i)) neval = neval + 1 IF (iprint.GT.0) THEN IF (MOD(neval,iprint).EQ.0) WRITE (lout,5100) neval, h(i), p END IF END IF END DO ! IF LOOP = NLOOP TEST FOR CONVERGENCE, OTHERWISE REPEAT MAIN CYCLE. 250 IF (loop.LT.nloop) CYCLE Main_loop ! CALCULATE MEAN & STANDARD DEVIATION OF FUNCTION VALUES FOR THE ! CURRENT SIMPLEX. hmean = SUM( h(1:np1) ) / FLOAT(np1) hstd = SUM( (h(1:np1) - hmean) ** 2 ) hstd = SQRT(hstd / FLOAT(np1)) ! IF THE RMS > STOPCR, SET IFLAG & LOOP TO ZERO AND GO TO THE ! START OF THE MAIN CYCLE AGAIN. IF (hstd.GT.stopcr .AND. neval.LE.maxfn) THEN iflag = 0 loop = 0 CYCLE Main_loop END IF ! FIND THE CENTROID OF THE CURRENT SIMPLEX AND THE FUNCTION VALUE THERE. DO i = 1, nop IF (step(i).NE.zero) THEN p(i) = SUM( g(1:np1,i) ) / FLOAT(np1) END IF END DO CALL functn(p,func) neval = neval + 1 IF (iprint.GT.0) THEN IF (MOD(neval,iprint).EQ.0) WRITE (lout,5100) neval, func, p END IF ! TEST WHETHER THE NO. OF FUNCTION VALUES ALLOWED, maxfn, HAS BEEN ! OVERRUN; IF SO, EXIT WITH IFAULT = 1. IF (neval.GT.maxfn) THEN ifault = 1 IF (iprint.LT.0) RETURN WRITE (lout,5200) maxfn WRITE (lout,5300) hstd WRITE (lout,5400) p WRITE (lout,5500) func RETURN END IF ! CONVERGENCE CRITERION SATISFIED. ! IF IFLAG = 0, SET IFLAG & SAVE HMEAN. ! IF IFLAG = 1 & CHANGE IN HMEAN <= STOPCR THEN SEARCH IS COMPLETE. IF (iprint.GE.0) THEN WRITE (lout,5600) WRITE (lout,5400) p WRITE (lout,5500) func END IF IF (iflag.EQ.0 .OR. ABS(savemn-hmean).GE.stopcr) THEN iflag = 1 savemn = hmean loop = 0 ELSE EXIT Main_loop END IF END DO Main_loop IF (iprint.GE.0) THEN WRITE (lout,5700) neval WRITE (lout,5800) p WRITE (lout,5900) func END IF IF (iquad.LE.0) RETURN !------------------------------------------------------------------ ! QUADRATIC SURFACE FITTING IF (iprint.GE.0) WRITE (lout,6000) ! EXPAND THE FINAL SIMPLEX, IF NECESSARY, TO OVERCOME ROUNDING ! ERRORS. hmin = func nmore = 0 DO i = 1, np1 DO test = ABS(h(i)-func) IF (test.LT.simp) THEN DO j = 1, nop IF (step(j).NE.zero) g(i,j) = (g(i,j)-p(j)) + g(i,j) pstst(j) = g(i,j) END DO CALL functn(pstst,h(i)) nmore = nmore + 1 neval = neval + 1 IF (h(i).GE.hmin) CYCLE hmin = h(i) IF (iprint.GE.0) WRITE (lout,5100) neval, hmin, pstst ELSE EXIT END IF END DO END DO ! FUNCTION VALUES ARE CALCULATED AT AN ADDITIONAL NAP POINTS. DO i = 1, nap i1 = i + 1 pstar = (g(1,:) + g(i1,:)) * half CALL functn(pstar,aval(i)) nmore = nmore + 1 neval = neval + 1 END DO ! THE MATRIX OF ESTIMATED SECOND DERIVATIVES IS CALCULATED AND ITS ! LOWER TRIANGLE STORED IN BMAT. a0 = h(1) DO i = 1, nap i1 = i - 1 i2 = i + 1 DO j = 1, i1 j1 = j + 1 pstst = (g(i2,:) + g(j1,:)) * half CALL functn(pstst,hstst) nmore = nmore + 1 neval = neval + 1 l = i * (i-1) / 2 + j bmat(l) = two * (hstst + a0 - aval(i) - aval(j)) END DO END DO l = 0 DO i = 1, nap i1 = i + 1 l = l + i bmat(l) = two * (h(i1) + a0 - two*aval(i)) END DO ! THE VECTOR OF ESTIMATED FIRST DERIVATIVES IS CALCULATED AND ! STORED IN AVAL. DO i = 1, nap i1 = i + 1 aval(i) = two * aval(i) - (h(i1) + 3.d0*a0) * half END DO ! THE MATRIX Q OF NELDER & MEAD IS CALCULATED AND STORED IN G. pmin = g(1,:) DO i = 1, nap i1 = i + 1 g(i1,:) = g(i1,:) - g(1,:) END DO DO i = 1, nap i1 = i + 1 g(i,:) = g(i1,:) END DO ! INVERT BMAT CALL syminv(bmat, nap, bmat, temp, nullty, ifault, rmax) IF (ifault.EQ.0) THEN irank = nap - nullty ELSE ! BMAT not +ve definite ! Resume search for the minimum IF (iprint.GE.0) WRITE (lout,6100) ifault = 2 IF (neval.GT.maxfn) RETURN IF (iprint.GE.0) WRITE (lout,6200) step = half * step GO TO 20 END IF ! BMAT*A/2 IS CALCULATED AND STORED IN H. DO i = 1, nap h(i) = zero DO j = 1, nap IF (j.LE.i) THEN l = i * (i-1) / 2 + j ELSE l = j * (j-1) / 2 + i END IF h(i) = h(i) + bmat(l) * aval(j) END DO END DO ! FIND THE POSITION, PMIN, & VALUE, YMIN, OF THE MINIMUM OF THE ! QUADRATIC. ymin = DOT_PRODUCT( h(1:nap), aval(1:nap) ) ymin = a0 - ymin DO i = 1, nop pstst(i) = DOT_PRODUCT( h(1:nap), g(1:nap,i) ) END DO pmin = pmin - pstst IF (iprint.GE.0) THEN WRITE (lout,6300) ymin, pmin WRITE (lout,6400) END IF ! Q*BMAT*Q'/2 IS CALCULATED & ITS LOWER TRIANGLE STORED IN VC DO i = 1, nop DO j = 1, nap h(j) = zero DO k = 1, nap IF (k.LE.j) THEN l = j * (j-1) / 2 + k ELSE l = k * (k-1) / 2 + j END IF h(j) = h(j) + bmat(l) * g(k,i) * half END DO END DO DO j = i, nop l = j * (j-1) / 2 + i vc(l) = DOT_PRODUCT( h(1:nap), g(1:nap,j) ) END DO END DO ! THE DIAGONAL ELEMENTS OF VC ARE COPIED INTO VAR. j = 0 DO i = 1, nop j = j + i var(i) = vc(j) END DO IF (iprint.LT.0) RETURN WRITE (lout,6500) irank CALL print_tri_matrix(vc, nop, lout) WRITE (lout,6600) CALL syminv(vc, nap, bmat, temp, nullty, ifault, rmax) ! BMAT NOW CONTAINS THE INFORMATION MATRIX WRITE (lout,6700) CALL print_tri_matrix(bmat, nop, lout) ii = 0 ij = 0 DO i = 1, nop ii = ii + i IF (vc(ii).GT.zero) THEN vc(ii) = one / SQRT(vc(ii)) ELSE vc(ii) = zero END IF jj = 0 DO j = 1, i - 1 jj = jj + j ij = ij + 1 vc(ij) = vc(ij) * vc(ii) * vc(jj) END DO ij = ij + 1 END DO WRITE (lout,6800) ii = 0 DO i = 1, nop ii = ii + i IF (vc(ii).NE.zero) vc(ii) = one END DO CALL print_tri_matrix(vc, nop, lout) ! Exit, on successful termination. 650 WRITE (lout,6900) nmore RETURN 5000 FORMAT (' Progress Report every',i4,' function evaluations'/, & ' EVAL. FUNC.VALUE.',10X,'PARAMETER VALUES') 5100 FORMAT (/1X,i4,2X,g12.5,2X,5G11.4,3(/21X,5G11.4)) 5200 FORMAT (' No. of function evaluations > ',i5) 5300 FORMAT (' RMS of function values of last simplex =',g14.6) 5400 FORMAT (' Centroid of last simplex =',4(/1X,6G13.5)) 5500 FORMAT (' Function value at centroid =',g14.6) 5600 FORMAT (/' EVIDENCE OF CONVERGENCE') 5700 FORMAT (/' Minimum found after',i5,' function evaluations') 5800 FORMAT (' Minimum at',4(/1X,6G13.6)) 5900 FORMAT (' Function value at minimum =',g14.6) 6000 FORMAT (/' Fitting quadratic surface about supposed minimum'/) 6100 FORMAT (/' MATRIX OF ESTIMATED SECOND DERIVATIVES NOT +VE DEFN.'/ & ' MINIMUM PROBABLY NOT FOUND'/) 6200 FORMAT (/10X,'Search restarting'/) 6300 FORMAT (' Minimum of quadratic surface =',g14.6,' at',4(/1X,6G13.5)) 6400 FORMAT (' IF THIS DIFFERS BY MUCH FROM THE MINIMUM ESTIMATED',1X, & 'FROM THE MINIMIZATION,'/' THE MINIMUM MAY BE FALSE &/OR THE ' & 'INFORMATION MATRIX MAY BE',1X,'INACCURATE'/) 6500 FORMAT (' Rank of information matrix =',i3/ & ' Inverse of information matrix:-') 6600 FORMAT (/' If the function minimized was -LOG(LIKELIHOOD),'/ & ' this is the covariance matrix of the parameters.'/ & ' If the function was a sum of squares of residuals,'/ & ' this matrix must be multiplied by twice the estimated ', & 'residual variance'/' to obtain the covariance matrix.'/) 6700 FORMAT (' INFORMATION MATRIX:-'/) 6800 FORMAT (/' CORRELATION MATRIX:-') 6900 FORMAT (/' A further',i4,' function evaluations have been used'/) END SUBROUTINE minim SUBROUTINE syminv(a, n, c, w, nullty, ifault, rmax) ! ALGORITHM AS7, APPLIED STATISTICS, VOL.17, 1968. ! ARGUMENTS:- ! A() = INPUT, THE SYMMETRIC MATRIX TO BE INVERTED, STORED IN ! LOWER TRIANGULAR FORM ! N = INPUT, ORDER OF THE MATRIX ! C() = OUTPUT, THE INVERSE OF A (A GENERALIZED INVERSE IF C IS ! SINGULAR), ALSO STORED IN LOWER TRIANGULAR. ! C AND A MAY OCCUPY THE SAME LOCATIONS. ! W() = WORKSPACE, DIMENSION AT LEAST N. ! NULLTY = OUTPUT, THE RANK DEFICIENCY OF A. ! IFAULT = OUTPUT, ERROR INDICATOR ! = 1 IF N < 1 ! = 2 IF A IS NOT +VE SEMI-DEFINITE ! = 0 OTHERWISE ! RMAX = OUTPUT, APPROXIMATE BOUND ON THE ACCURACY OF THE DIAGONAL ! ELEMENTS OF C. E.G. IF RMAX = 1.E-04 THEN THE DIAGONAL ! ELEMENTS OF C WILL BE ACCURATE TO ABOUT 4 DEC. DIGITS. ! LATEST REVISION - 1 April 1985 !*************************************************************************** IMPLICIT DOUBLE PRECISION (a-h,o-z) DIMENSION a(*), c(*), w(n) DOUBLE PRECISION :: zero = 0.d0, one = 1.d0 nrow = n ifault = 1 IF (nrow.GT.0) THEN ifault = 0 ! CHOLESKY FACTORIZATION OF A, RESULT IN C CALL chola(a, nrow, c, nullty, ifault, rmax, w) IF (ifault.EQ.0) THEN ! INVERT C & FORM THE PRODUCT (CINV)'*CINV, WHERE CINV IS THE INVERSE ! OF C, ROW BY ROW STARTING WITH THE LAST ROW. ! IROW = THE ROW NUMBER, NDIAG = LOCATION OF LAST ELEMENT IN THE ROW. nn = nrow * (nrow+1) / 2 irow = nrow ndiag = nn 10 IF (c(ndiag).NE.zero) THEN l = ndiag DO i = irow, nrow w(i) = c(l) l = l + i END DO icol = nrow jcol = nn mdiag = nn 30 l = jcol x = zero IF (icol.EQ.irow) x = one / w(irow) k = nrow 40 IF (k.NE.irow) THEN x = x - w(k) * c(l) k = k - 1 l = l - 1 IF (l.GT.mdiag) l = l - k + 1 GO TO 40 END IF c(l) = x / w(irow) IF (icol.EQ.irow) GO TO 60 mdiag = mdiag - icol icol = icol - 1 jcol = jcol - 1 GO TO 30 END IF ! (c(ndiag).NE.zero) l = ndiag DO j = irow, nrow c(l) = zero l = l + j END DO 60 ndiag = ndiag - irow irow = irow - 1 IF (irow.NE.0) GO TO 10 END IF END IF RETURN END SUBROUTINE syminv SUBROUTINE chola(a, n, u, nullty, ifault, rmax, r) ! ALGORITHM AS6, APPLIED STATISTICS, VOL.17, 1968, WITH ! MODIFICATIONS BY A.J.MILLER ! ARGUMENTS:- ! A() = INPUT, A +VE DEFINITE MATRIX STORED IN LOWER-TRIANGULAR ! FORM. ! N = INPUT, THE ORDER OF A ! U() = OUTPUT, A LOWER TRIANGULAR MATRIX SUCH THAT U*U' = A. ! A & U MAY OCCUPY THE SAME LOCATIONS. ! NULLTY = OUTPUT, THE RANK DEFICIENCY OF A. ! IFAULT = OUTPUT, ERROR INDICATOR ! = 1 IF N < 1 ! = 2 IF A IS NOT +VE SEMI-DEFINITE ! = 0 OTHERWISE ! RMAX = OUTPUT, AN ESTIMATE OF THE RELATIVE ACCURACY OF THE ! DIAGONAL ELEMENTS OF U. ! R() = OUTPUT, ARRAY CONTAINING BOUNDS ON THE RELATIVE ACCURACY ! OF EACH DIAGONAL ELEMENT OF U. ! LATEST REVISION - 1 April 1985 !*************************************************************************** IMPLICIT DOUBLE PRECISION (a-h,o-z) DIMENSION a(*), u(*), r(n) ! ETA SHOULD BE SET EQUAL TO THE SMALLEST +VE VALUE SUCH THAT ! 1.D0 + ETA IS CALCULATED AS BEING GREATER THAN 1.D0 IN THE ACCURACY ! BEING USED. DOUBLE PRECISION :: eta = 1.d-16, zero = 0.d0 ifault = 1 IF (n.GT.0) THEN ifault = 2 nullty = 0 rmax = eta r(1) = eta j = 1 k = 0 ! FACTORIZE COLUMN BY COLUMN, ICOL = COLUMN NO. DO 50 icol = 1, n l = 0 ! IROW = ROW NUMBER WITHIN COLUMN ICOL DO 30 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) EXIT w = w - u(l) * u(m) IF (irow.EQ.icol) rsq = rsq + (u(l)**2*r(i)) ** 2 m = m + 1 10 CONTINUE IF (irow.EQ.icol) GO TO 40 IF (u(l).NE.zero) THEN u(k) = w / u(l) ELSE u(k) = zero IF (ABS(w).GT.ABS(rmax*a(k))) GO TO 60 END IF 30 CONTINUE ! END OF ROW, ESTIMATE RELATIVE ACCURACY OF DIAGONAL ELEMENT. 40 rsq = SQRT(rsq) IF (ABS(w).GT.5.*rsq) THEN IF (w.LT.zero) RETURN u(k) = SQRT(w) r(i) = rsq / w IF (r(i).GT.rmax) rmax = r(i) ELSE u(k) = zero nullty = nullty + 1 END IF j = j + icol 50 CONTINUE ifault = zero END IF 60 RETURN END SUBROUTINE chola SUBROUTINE print_tri_matrix(a, n, lout) IMPLICIT NONE INTEGER, INTENT(IN) :: n, lout DOUBLE PRECISION, INTENT(IN) :: a(*) ! Local variables INTEGER :: i, ii, i1, i2, l l = 1 DO l = 1, n, 6 ii = l * (l-1) / 2 DO i = l, n i1 = ii + l ii = ii + i i2 = MIN(ii,i1+5) WRITE (lout,'(1X,6G13.5)') a(i1:i2) END DO END DO RETURN END SUBROUTINE print_tri_matrix