! Last change: DOS 29 Jul 2000 4:57 pm ! *** copyright 2000 *** ! *** filename plum2t.f95 *** John F. Monahan ** ! ********************** program pplum2t ! test problem for plum2t -- second version of plumit ! general minimizer of functions -- requires gradient & hessian ! implicit none real(kind=8), external :: f1,f2,f3,f4 external gh1,gh2,gh3,gh4 real(kind=8), dimension(10) :: x,typx,grad real(kind=8), dimension(55) :: hess real(kind=8) fnew integer i,j,n,termcd ! data x/ -1.2d0, 1.d0, -1.2d0, 1.d0, 6*0.d0 / data typx/ 10*1.d0 / ! 21 format(//' Termination code',i4,3x,'Function at min',f12.6) 22 format(' Solution Gradient',10x,'Hessian') 23 format(f11.6,f10.6,2x,5f11.4/25x,5f11.4) 25 format(//'First Test of plum2t -- Rosenbrock function, n=4') 26 format(//'Second Test of plum2t -- Rosenbrock function, n=10') 27 format(//'Third Test of plum2t -- Cox likelihood, n=2') 28 format(//'Fourth Test of plum2t -- D&S Example 9.4.1') ! OPEN( UNIT=6, FILE='plum2t.out' ) ! first problem in 4 dimensions write(6,25) n = 4 call plum2t(f1,gh1,x,n,40,typx,fnew,grad,hess,6,termcd) write(6,25) write(6,21) termcd,fnew write(6,22) do i = 1,n write(6,23) x(i),grad(i),(hess((i*(i-1))/2+j),j=1,i) end do ! loop on i ! second problem in 10 dimensions write(6,26) n =10 do i = 1,5 x(2*i-1) = -1.2d0 x(2*i) = 1.d0 end do ! loop on i ! call plum2t(f2,gh2,x,n,40,typx,fnew,grad,hess,6,termcd) write(6,26) write(6,21) termcd,fnew write(6,22) do i = 1,n write(6,23) x(i),grad(i),(hess((i*(i-1))/2+j),j=1,i) end do ! loop on i ! third problem is Cox loglikelihood write(6,27) n = 2 x(1) = 0.d0 x(2) = 0.d0 call plum2t(f3,gh3,x,n,40,typx,fnew,grad,hess,6,termcd) write(6,27) write(6,21) termcd,fnew write(6,22) do i = 1,n write(6,23) x(i),grad(i),(hess((i*(i-1))/2+j),j=1,i) end do ! loop on i ! fourth problem write(6,28) n = 2 x(1) = 0.d0 x(2) = 0.d0 call plum2t(f4,gh4,x,n,40,typx,fnew,grad,hess,6,termcd) write(6,28) write(6,21) termcd,fnew write(6,22) do i = 1,n write(6,23) x(i),grad(i),(hess((i*(i-1))/2+j),j=1,i) end do ! loop on i stop end program pplum2t real(kind=8) function f1(x) implicit none real(kind=8), dimension(4), intent(in) :: x ! extended rosenbrock function in 4 dimensions f1 = 100.d0*(x(2)-x(1)*x(1))*(x(2)-x(1)*x(1)) & & + (1.d0-x(1))*(1.d0-x(1)) & & + 100.d0*(x(4)-x(3)*x(3))*(x(4)-x(3)*x(3)) & & + (1.d0-x(3))*(1.d0-x(3)) return end function f1 subroutine gh1(x,g,h) implicit none real(kind=8), dimension(4), intent(in) :: x real(kind=8), dimension(4), intent(out) :: g real(kind=8), dimension(10), intent(out) :: h ! gradient and hessian for rosenbrock function ! first gradient g(1) = -400.*x(1)*(x(2)-x(1)*x(1)) - 2.*(1.-x(1)) g(2) = 200.*(x(2)-x(1)*x(1)) g(3) = -400.*x(3)*(x(4)-x(3)*x(3)) - 2.*(1.-x(3)) g(4) = 200.*(x(4)-x(3)*x(3)) ! hessian -- symmetric storage h(1) = 1200.*x(1)*x(1) - 400.*x(2) + 2. h(2) = -400.*x(1) h(3) = 200. h(4) = 0.d0 h(5) = 0.d0 h(6) = 1200.*x(3)*x(3) - 400.*x(4) + 2. h(7) = 0. h(8) = 0. h(9) = -400.*x(3) h(10) = 200. return end subroutine gh1 real(kind=8) function f2(x) implicit none real(kind=8), dimension(10), intent(in) :: x integer i,i2,i2m ! extended rosenbrock function in 10 dimensions f2 = 0.d0 do i = 1,5 i2 = i + i i2m = i2 - 1 f2 = f2 + 100.d0*(x(i2)-x(i2m)*x(i2m))*(x(i2)-x(i2m)*x(i2m)) & & + (1.d0-x(i2m))*(1.d0-x(i2m)) end do ! loop on i return end function f2 subroutine gh2(x,g,h) implicit none real(kind=8), dimension(10), intent(in) :: x real(kind=8), dimension(10), intent(out) :: g real(kind=8), dimension(55), intent(out) :: h integer i,i2,i2m ! gradient and hessian for rosenbrock function ! fill hessian with zeros first do i = 1,55 h(i) = 0.d0 end do ! loop on i ! do i = 1,5 i2 = i + i i2m = i2 - 1 ! gradient g(i2m) = - 400.*x(i2m)*(x(i2)-x(i2m)*x(i2m)) - 2.*(1.-x(i2m)) g(i2) = 200.*(x(i2)-x(i2m)*x(i2m)) ! hessian -- in symmetric storage mode ! it is block diagonal h( (i2m*i2)/2 ) = 1200.*x(i2m)*x(i2m) -400.*x(i2)+2. h( (i2*(i2+1))/2 ) = 200. h( (i2*(i2+1))/2 -1 ) = -400.*x(i2m) end do ! loop on i return end subroutine gh2 real(kind=8) function f3(x) implicit none real(kind=8), dimension(2), intent(in) :: x real(kind=8), dimension(4) :: t real(kind=8) v integer, dimension(4) :: n,s integer i ! Cox logistic regression example data n/ 55, 157, 159, 16 / data s/ 0, 2, 7, 3 /, t/ 7.d0, 14.d0, 27.d0, 51.d0 / f3 = 0.d0 do i = 1,4 v = x(1) + x(2)*t(i) f3 = f3 - v*s(i) + n(i)*log(1.d0+exp(v) ) end do return end function f3 subroutine gh3(x,g,h) implicit none real(kind=8), dimension(2), intent(in) :: x real(kind=8), dimension(2), intent(out) :: g real(kind=8), dimension(3), intent(out) :: h real(kind=8), dimension(4) :: t real(kind=8) v,p integer i,n(4),s(4) ! Cox logistic regression example data n/ 55, 157, 159, 16 / data s/ 0, 2, 7, 3 /, t/ 7.d0, 14.d0, 27.d0, 51.d0 / g(1) = 0.d0 g(2) = 0.d0 h(1) = 0.d0 h(2) = 0.d0 h(3) = 0.d0 do i = 1,4 v = x(1) + x(2)*t(i) p = 1.d0/(1.d0+exp(-v)) ! gradient pieces g(1) = g(1) -s(i) + n(i)*p g(2) = g(2) -s(i)*t(i) + n(i)*t(i)*p ! hessian pieces h(1) = h(1) + n(i)*p*(1.d0-p) h(2) = h(2) + n(i)*t(i)*p*(1.d0-p) h(3) = h(3) + n(i)*t(i)*t(i)*p*(1.d0-p) end do ! loop on i return end subroutine gh3 real(kind=8) function f4(x) implicit none real(kind=8), dimension(2), intent(in) :: x ! D&S Example 9.4.1 f4 = (x(1)-2.)**4 + (x(1)-2.d0)*(x(1)-2.d0)*x(2)*x(2) & & + (x(2)+1.d0)*(x(2)+1.d0) return end function f4 subroutine gh4(x,g,h) implicit none real(kind=8), dimension(2), intent(in) :: x real(kind=8), dimension(2), intent(out) :: g real(kind=8), dimension(3), intent(out) :: h ! gradient for D&S example g(1) = 4.d0*(x(1)-2.)*(x(1)-2.)*(x(1)-2.) & & + 2.d0*(x(1)-2.d0)*x(2)*x(2) g(2) = 2.d0*(x(1)-2.d0)*(x(1)-2.d0)*x(2) + 2.d0*(x(2)+1.d0) ! Hessian h(1) = 12.d0*(x(1)-2.)*(x(1)-2.) & & + 2.d0*x(2)*x(2) h(2) = 4.d0*(x(1)-2.d0)*x(2) h(3) = 12.d0*(x(1)-2.)*(x(1)-2.) + 2.d0 return end subroutine gh4 subroutine plum2t(f,gh,x,n,mit,typx,fnew,grad,hess, & & outcod,termcd) ! subroutine to minimize f ! requiring both function f(x) and routine to compute gradient ! and Hessian gh(x,grad,hess) ! designed following recommendations of Dennis and Schnabel (1983) ! Numerical Method for Unconstrained Optimization and Nonlinear ! Equations, Prentice-Hall ! ! f function to be minimized (function subprogram) ! gh subroutine to compute gradient and Hessian ! x real vector in: starting point for search ! out: minimum ! n integer dimension of x ! mit integer maximum number of iterations ! typx real vector typical values for |x|, used for scaling ! fnew real out: value of function at minimum ! grad real vector out: gradient of f at minimum ! hess real vector out: Hessian matrix at minimum in symmetric ! storage mode ** declare n*(n+1)/2 length ** ! outcod integer output code for writing messages ! termcd integer termination code ! termcd=1 means stopped on small gradient (success) ! termcd=2 means stop for small step size ! termcd=3 means stop for failing backtracks ! termcd=4 means stop -- too many iterations ! termcd=5 means too many big steps, maybe unbounded ! termcd=6 means Hessian not positive definite ! termcd=7 means descent direction goes up (not posdf) ! 1 is success, 2 or 3 maybe done, perhaps limited by accuracy of f ! 4, 5, 6, 7 are failures ! ! required external routines: ! del12f computes numerical gradient and Hessian ! chlzky computes Cholesky decomposition of pos def matrix ! chlzhi solves system of equations using Cholesky ! chlzih factor as matrix ! adjust called by chlzky to control size of determinant ! ! j f monahan (july, 1993) ! recoded October 1999, April 2000 for Fortran 95 ! IMPLICIT NONE REAL(KIND=8) f INTEGER, INTENT(IN) :: n,mit REAL(KIND=8), DIMENSION( n ), INTENT(IN OUT) :: x,typx,grad REAL(KIND=8), DIMENSION( (n*(n+1))/2 ), INTENT(OUT) :: hess REAL(KIND=8), INTENT(IN OUT) :: fnew INTEGER, INTENT(IN) :: outcod INTEGER, INTENT(OUT) :: termcd REAL(KIND=8), DIMENSION( n ) :: step,bstep REAL(KIND=8), DIMENSION( (n*(n+1))/2 ) :: thess REAL(KIND=8) stepmx,stsize,fold,flast,gradmx,det,c,slope,rlleng,& & lam2sm,sl,sltry,slast,v1,v2,a,b,disc ! integer i,nn,kit,idet,obig,mbig ! ! *** these constants should change with computer *** ! designed for ieee binary double precision ! machine epsilon REAL(KIND=8), PARAMETER :: eta = 4.44d-16 ! if the evaluations of function f are only accurate to d digits, ! then a better value would be eta = 10.**(-d) ! ! other constants are functions of this machine epsilon ! ! gradtl used for stopping on small gradient REAL(KIND=8), PARAMETER :: gradtl = 7.63d-6 !**d gradtl = cbrt(eta) ! steptl used for stopping on small steps REAL(KIND=8), PARAMETER :: steptl = 5.82d-11 !**d steptl = cbrt(eta)*cbrt(eta) ! ddiff used for finite differences REAL(KIND=8), PARAMETER :: ddiff = 7.63d-6 !**d ddiff = gradtl ! alphmn is required improvement in step REAL(KIND=8), PARAMETER :: alphmn = 1.d-3 ! *** done with constants *** ! 21 format(' begin iteration ',i4,' out of max ',i4,' current x'& & /2x,6f12.6/4x,6f12.6/4x,6f12.6/6x,6f12.6) 22 format(' stepsize',f10.4,' stepmx ',f10.4,' lambda ',f5.3, & & ' fnew ',f10.4) 23 format(' quadratic backtrack') 24 format(' cubic backtrack ') 43 format(' stopped -- too many iterations') 44 format(' stopped -- backtracking failed, but might be done') 45 format(' stopped -- too many big steps, probably unbounded') 46 format(' stopped -- taking steps too small') 47 format(' stopped -- Hessian no longer positive definite') 48 format(' stopped -- gradient small -- success') 49 format(' stopped -- descent direction going up ') ! ! need this later nn = (n*(n+1))/2 ! initialize iteration count kit = 0 ! find maximum stepsize stepmx = 0. stsize = 0. do i = 1,n typx(i) = abs(typx(i)) stepmx = stepmx + ( x(i)/typx(i) ) * ( x(i)/typx(i) ) stsize = stsize + 1./( typx(i) * typx(i) ) end do ! loop on i stepmx = 400.*sqrt( max( stepmx, stsize ) ) ! stepmx limits size of any step ! and keeps from running far away ! ! first evaluate function here fold = f(x) ! get gradients (and Hessian) ! now call gh call gh(x,grad,hess) do i = 1,nn thess(i) = hess(i) ! save Hessian end do ! loop on i ! have gradient and Hessian ! test if gradient too small to continue gradmx = 0. do i = 1,n ! also copy gradient for next step step(i) = - grad(i) c = max( abs(x(i)), typx(i) )/abs( fold ) gradmx = max( gradmx, abs(grad(i))*c ) end do ! loop on i ! are we done already? if( gradmx .lt. gradtl/1000. ) then ! successful completion ! stop on gradient too small write(outcod,48) termcd = 1 return end if ! (gradmx.lt.gradtl/1000) ! divide by 1000 -- stricter at start than finish ! ! gradient not zero -- where to move? ! compute tentative step using Newton ! ! first Cholesky factorization of hess call chlzky(thess,n,det,idet) ! if not pos def then die quickly if( idet .lt. -2000000000 ) then ! Hessian not positive definite! write(outcod,47) termcd = 6 return end if ! ( Cholesky failed ) ! compute step call chlzhi(thess,n,step) call chlzih(thess,n,step) ! ! now the big section of line search code ! ! start main loop do kit = 1,mit ! write(outcod,21) kit,mit,(x(i),i=1,n) ! first step -- see how well Newton step does ! ! how big is this step stsize = 0. do i = 1,n stsize = stsize + ( step(i)/typx(i) ) * ( step(i)/typx(i) ) end do ! loop on i stsize = sqrt( stsize ) ! if step is too big, shorten it if( stsize .gt. stepmx ) then c = stepmx/stsize do i = 1,n step(i) = c * step(i) end do ! loop on i stsize = stepmx ! end if ! ( stsize .gt. stepmx ) ! compute slope of f in step direction slope = 0. rlleng = 0. do i = 1,n slope = slope + step(i)*grad(i) rlleng = max( rlleng, abs(step(i))/max(abs(x(i)),typx(i)) ) end do ! loop on i if( slope .gt. 0. ) then ! steepest descent seems to go up ! quadratic form in inverse of Hessian < 0 ! really step*grad > 0 write(outcod,49) termcd = 7 return end if ! ( slope .gt. 0. ) ! temporary branch if something is goofy ! lam2sm is smallest relative stepsize lam2sm = steptl / rlleng ! ! go to new point and test the water ! ! sl = 1 means full step sl = 1. ! make new point do ! *** backtrack repeat do *** do i = 1,n bstep(i) = x(i) + sl*step(i) end do ! loop on i ! evaluate function fnew = f(bstep) write(outcod,22) stsize,stepmx,sl,fnew ! how much improvement?? if( fnew .lt. fold + alphmn*slope*sl ) exit ! from backtracking ! if not an improvement, then backtrack ! ! have backtracked too much already? if( sl .lt. lam2sm ) then ! stop because of failing backtracks write(outcod,44) termcd = 3 return end if ! ( sl .lt. lam2sm ) ! backtracking didn't work -- quit ! ! backtrack code ! if( sl .eq. 1. ) then ! first fit quadratic write(outcod,23) sltry = - slope/( 2.*(fnew - fold - slope) ) flast = fnew slast = sl sl = max( sl/10., sltry ) ! go try this one else write(outcod,24) ! next try with cubic if failed before ! this code is messy v1 = fnew - fold - sl*slope v2 = flast - fold - slast*slope a = ( v1/(sl*sl) - v2/(slast*slast) )/(sl-slast) b = ( v2*sl/(slast*slast) - v1*slast/(sl*sl) )/(sl-slast) disc = b*b - 3.*a*slope if ( a .eq. 0. ) then ! cubic is a quadratic sltry = - slope/(2.*b) else ! solve cubic sltry = (sqrt(disc)-b)/(3.*a) end if ! ( a .eq. 0 ) ! shorten at least by half sltry = min( sltry, sl/2. ) flast = fnew slast = sl sl = max( sl/10., sltry ) ! go and try this one end if ! ( sl .eq. 1. ) end do ! *** backtrack loop **** ! ! last step was fruitful, where next? ! was the last step as big as could? obig = 0 if( (sl .eq. 1.) .and. (stsize .gt. .99*stepmx ) ) obig = 1 ! several things at once ! copy new point, new differences, stepsize stsize = 0. do i = 1,n step(i) = x(i) - bstep(i) x(i) = bstep(i) c = max( abs(x(i)), typx(i) ) stsize = max( stsize, abs(step(i))/c ) end do ! loop on i ! ! get gradient and Hessian at new point call gh(x,grad,hess) do i = 1,nn thess(i) = hess(i) ! save Hessian end do ! loop on i ! do I stop now? ! first see if gradient too small to continue gradmx = 0. do i = 1,n ! also copy gradient for next step step(i) = - grad(i) c = max( abs(x(i)), typx(i) )/abs(fold) gradmx = max( gradmx, abs(grad(i)*c) ) end do ! loop on i fold = fnew if ( gradmx .lt. gradtl ) then ! successful completion ! stop on gradient too small write(outcod,48) termcd = 1 return end if ! ( gradmx .lt. gradtl ) ! gradient not too small, what about step? if ( stsize .lt. steptl ) then ! stop because step size is too small write(outcod,46) termcd = 2 return end if ! ( stsize .lt. steptl ) ! too many maxsteps? are we running away? if ( obig .gt. 0 ) then ! things are ok -- work on next step mbig = mbig + 1 ! we've done mbig maxsteps in a row if ( mbig .eq. 5 ) then ! stop because of too many maxsteps ! probably unbounded minimum write(outcod,45) termcd = 5 return end if ! ( mbig .eq. 5 ) ! work on next step end if ! ( obig .gt. 0 ) if( obig .eq. 0 ) mbig = 0 ! reset count if last one was ok ! ! compute tenative step using Newton call chlzky(thess,n,det,idet) ! if not pos def die quickly if( idet .lt. -2000000000 ) then ! Hessian not positive definite! write(outcod,47) termcd = 6 return end if ! ( Cholesky failed ) ! compute step call chlzhi(thess,n,step) call chlzih(thess,n,step) ! now go back to the line search code ! too many iterations? end do ! loop on kit ! ! stop because of too many iterations write(outcod,43) termcd = 4 return end subroutine plum2t SUBROUTINE ADJUST(D,I) ! NORMALIZES D WHILE KEEPING CONSTANT VALUE OF D * (2**I) INTEGER, PARAMETER :: IBIG = -2147483644 REAL(KIND=8), INTENT (IN OUT) :: D INTEGER, INTENT(IN OUT) :: I IF ( D .GT. 0.0 ) THEN DO WHILE ( D .LT. 1.0 ) D = D*16. I = I-4 ENDDO DO WHILE ( D .GT. 16.0 ) D = D/16. I = I+4 ENDDO ELSE I = IBIG ! IF D < 0 THEN I = - BIG ENDIF RETURN END SUBROUTINE ADJUST SUBROUTINE CHLZKY(A,N,DET,IDET) ! CHLSKY COMPUTES THE CHOLESKY (SQUARE-ROOT) FACTORIZATION ! A = L * ( L - TRANSPOSE ) WHERE L IS LOWER TRIANGULAR ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! A IS OVERWRITTEN BY L IN ITS LOWER TRIANGULAR PART ! SUBROUTINE ADJUST KEEPS DETERMINANT FROM EXPLODING USING ! 1 LE DET LE 16 AND DETERMINANT OF A IS DET*2**IDET ! ! J F MONAHAN (DEC,1983) DEPT OF STAT, N C S U, RALEIGH, N C 27650 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL(KIND=8), DIMENSION( (N*(N+1))/2 ), INTENT(IN OUT) :: A REAL(KIND=8), INTENT(OUT) :: DET INTEGER, INTENT(OUT) :: IDET REAL(KIND=8) T,S INTEGER I,J,K,IM1,JM1 ! INITIALIZE FOR DETERMINANT DET = 1. IDET = 0 ! DO I-TH ROW DO I = 1,N T = A( (I*(I+1))/2 ) ! FIRST ROW IS TRIVIAL IF(I .GT. 1) THEN IM1 = I-1 DO J=1,IM1 ! WORK ON (I,J)-TH ELEMENT S = A( (I*(I-1))/2 + J ) IF(J .GT. 1) THEN JM1 = J-1 DO K = 1,JM1 S = S - A( (J*(J-1))/2 + K ) * A( (I*(I-1))/2 + K ) ENDDO ! LOOP ON K ENDIF A( (I*(I-1))/2 + J )= S / A( (J*(J+1))/2 ) ! WORK ON DIAGONAL ELEMENT T = T - A( (I*(I-1))/2 + J ) * A( (I*(I-1))/2 + J ) ENDDO ! LOOP ON J ! FINISHED WITH LOOP ON J ! UPDATE DETERMINANT WITH DIAGONAL ENDIF DET = DET * T ! KEEP DET FROM EXPLODING CALL ADJUST(DET,IDET) ! ABORT IF DET IS NEGATIVE OR ZERO IF(IDET.LT.-2000000000) RETURN ! FINISH A POSITIVE DIAGONAL A( (I*(I+1))/2 ) = SQRT(T) ENDDO ! LOOP ON I RETURN END SUBROUTINE CHLZKY SUBROUTINE CHLZHI(A,N,Y) ! ! OVERWRITES Y WITH INVERSE( L ) * Y ! WHERE L IS A LOWER TRIANGULAR MATRIX (N BY N) STORED IN A ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! TO COMPUTE QUADRATIC FORM IN THE INVERSE OF A POSITIVE DEFINITE ! MATRIX A, THAT IS Y' INVERSE( A ) Y THEN ! FIRST: FACTOR A BY CHLSKY(A,N,DET,IDET) ! THIS COMPUTES LOWER TRIANGULAR L SO THAT A = L * L' ! SECOND: COMPUTE INVERSE(L) * Y BY CHLSHI(A,N,Y) ! THIRD: COMPUTE THE SUM OF SQUARES OF THE ELEMENTS OF Y ! ! J F MONAHAN (JULY,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL(KIND=8), DIMENSION( (N*(N+1))/2 ), INTENT(IN) :: A REAL(KIND=8), DIMENSION( N ), INTENT(IN OUT) :: Y INTEGER I,J,IM1 ! ! SINCE THE MATRIX IS LOWER TRIANGULAR, SOLVE FROM THE TOP DOWN ! DO I = 1,N ! BRANCH AROUND ON THE FIRST ONE IF( I .GT. 1 ) THEN IM1 = I - 1 DO J = 1,IM1 Y(I) = Y(I) - A( (I*(I-1))/2 + J )*Y(J) ENDDO ENDIF Y(I) = Y(I) / A( (I*(I+1))/2 ) ENDDO ! LOOP ON I RETURN END SUBROUTINE CHLZHI SUBROUTINE CHLZIH(A,N,Y) ! *** NOTICE TRANSPOSE *** ! OVERWRITES Y WITH INVERSE( L' ) * Y ! WHERE L IS A LOWER TRIANGULAR MATRIX (N BY N) STORED IN A ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! TO COMPUTE THE PRODUCT OF THE INVERSE OF A POSITIVE DEFINITE MATRIX ! AND A VECTOR Y, THAT IS INVERSE( A ) * Y THEN ! FIRST: FACTOR A BY CHLSKY(A,N,DET,IDET) ! THIS COMPUTES LOWER TRIANGULAR L SO THAT A = L * L' ! SECOND: COMPUTE INVERSE(L) * Y BY CHLSHI(A,N,Y) ! THIRD: COMPUTE INVERSE(L') * ( INVERSE(L)*Y ) BY CHLSIH(A,N,Y) ! ! J F MONAHAN (JULY,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL(KIND=8), DIMENSION( (N*(N+1))/2 ), INTENT(IN) :: A REAL(KIND=8), DIMENSION( N ), INTENT(IN OUT) :: Y INTEGER I,II,IP1,J ! ! SINCE THE MATRIX IS UPPER TRIANGULAR, SOLVE FROM THE BOTTOM UP ! ! DO THE LAST ONE SEPARATELY Y(N) = Y(N) / A( (N*(N+1))/2 ) ! IF N = 1 THEN QUIT ELSE START LOOPING IF( N .EQ. 1 ) RETURN DO II = 2,N ! I BELOW NOW GOES FROM N-1 DOWN TO 1 I = N + 1 - II IP1 = I + 1 DO J = IP1,N Y(I) = Y(I) - A( (J*(J-1))/2 + I )*Y(J) ENDDO ! LOOP ON J Y(I) = Y(I) / A( (I*(I+1))/2 ) ENDDO ! LOOP ON II,I RETURN END SUBROUTINE CHLZIH ! *** end of filename plum2t.f95 *********************