!     Last change:  J    17 Aug 2005    8:09 pm
! *** filename nllsqu1.f95              *** John F. Monahan **
!                                       **********************
!
!  SOME COMPILERS REQUIRE MODULE TO PREVIOUSLY EXIST OR BE
!   AT THE BEGINNING
!
      module passinfo
!                          pass info through module
      integer, parameter :: n = 12
      integer, parameter :: p = 3
      real(KIND=8), dimension(n) :: x,y
      end module passinfo
!
      program nllsqu1
!  test problem 1 of nonlinear least squares algorithms
!  problem is from Fuller, p 228
      use passinfo
      implicit none
      external res,jac
      real(KIND=8), dimension(n) :: r
      real(KIND=8), dimension(p) :: th
      real(KIND=8), dimension( (p*(p+1))/2 ) :: h
      real(KIND=8) sse,s
      integer pp,termcd,i,j
!
  20  format(f6.1,f4.1)
  21  format(' Termination Code',i3,5x,' SSE=',f16.9,' Var est ',     &
     &        f16.9)
  22  format(' Parameter Estimates',3f12.7)
  23  format(' Covariance Matrix of Estimates '/6x,f12.7              &
     &   /6x,2f12.7/6x,3f12.7)
  29  format(' Test Problem 1 -- Fuller example, n=',i3,'  p=',i2/)
!
      open( unit=5, file='fulnls.dat')
      open( unit=6, file='nllsqu1.out')
!
      pp = (p*(p+1)) / 2
      write(6,29) n,p
!                  starting values
      th(1) = 140.d0
      th(2) = -85.d0
      th(3) = 1.50d0
!                  read the data
      do i = 1,n
      read(5,20) y(i),x(i)
      write(6,20) y(i),x(i)
      end do  ! loop on i
!                  call nllsq1
      call nllsq(res,jac,n,p,th,r,10,sse,h,6,termcd)
!                  write out results
      write(6,22) th
      s = sse/real(n-p,8)     !KIND
      write(6,21) termcd,sse,s
!                  compute covariance matrix of estimates
      call chlzoi(h,p)
      do j = 1,pp
      h(j) = h(j)*s
      end do  ! loop on j
      write(6,23) h
      stop
      end program nllsqu1
      subroutine res(th,r)
!  compute residuals for nonlinear least squares example
      use passinfo
      implicit none
      real(KIND=8), dimension(p) :: th
      real(KIND=8), dimension(n) :: r
      real(KIND=8) s
      integer i
  21  format(' in res, sse=',f16.8)
      s = 0.d0
      do i = 1,n
      r(i) = y(i) - th(1) - th(2)*exp( - th(3)*x(i) )
      s = s + r(i)*r(i)
      end do  ! loop on i
      write(6,21) s  
      return
      end subroutine res
      subroutine jac(th,ws)
!  compute Jacobian for nonlinear least squares example
      use passinfo
      implicit none
      real(KIND=8), dimension(p) :: th
      real(KIND=8), dimension(n,p) :: ws
      real(KIND=8) s
      integer i
      do i = 1,n
      ws(i,1) = -1.d0                   ! notice negative signs
      s = exp( - th(3)*x(i) )
      ws(i,2) = -s                      ! we need Jacobian of residuals
      ws(i,3) = s*th(2)*x(i)            !  wrto parameters
      end do  ! loop on i
      return
      end subroutine jac
      subroutine nllsq(res,jac,n,p,x,r,mit,sse,h,iout,termcd)
!
!  minimize nonlinear sum of squares using Levenberg-Marquardt
!   algorithm with analytic Jacobian matrix
!
!  method follows algorithms by
!  Jorge J. More (1977) The Levenberg-Marquardt Algorithm:
!   Implementation and Theory, in 'Numerical Analysis'
!   Lecture Notes in Mathematics 630, Springer-Verlag, pp. 105-116
!
! Argument List
!
!  res             subroutine to compute residuals
!  jac             subroutine to compute Jacobian matrix
!  n     integer   number of observations
!  p     integer   number of parameters
!  x     real*8    (in) start point,    (out) nlls estimates
!  r     real*8    (out) residual vector at end
!  mit   integer   max number of iterations
!  sse   real*8    (out) error sum of squares
!  h     real*8    (out) estimate of (Cholesky factor of) Hessian at end
! iout   integer   output (print) unit number (like 6)
! termcd integer   termination code (1,2)=success, (4,5,6)=failure
!              termcd=1   stopped when residuals orthog to fitted vals
!              termcd=2   stopped when step size too small
!              termcd=4   stopped after too many (mit) iterations
!              termcd=5   stopped when step too small/not local quadrat
!              termcd=6   stopped when SSE goes up, no local improvement
!
!              required routines
!  user supplied
!   res, jac   compute residuals, Jacobian matrix, taking the form
!     subroutine res(x,r)    in-vector x, out-vector r residuals
!     subroutine jac(x,gbig) in-vector x, out-array Jacobian matrix
!                                          gbig(i,j) holds Jac(i,j)
!
!  auxiliary
!   rot734     computes Givens transformation
!   chlzhi     solves L x = y for x, given rhs y and Cholesky factor L
!   chlzih     solves L'x = y for x, given rhs y and Cholesky factor L
!
!  J F Monahan (Dec, 1993) Dept of Statistics, NCSU, Raleigh, NC 27695
!     Recoded November 1999, April 2000 for Fortran 95
!     Two corrections, August 2005
!
!              argument declarations
!
      implicit none
      integer, INTENT(IN) :: n,p,mit,iout
      integer, INTENT(OUT) :: termcd
      real(KIND=8), dimension(p), INTENT(IN OUT) :: x
      real(KIND=8), dimension(n), INTENT(OUT) :: r
      real(KIND=8), dimension( (p*(p+1))/2 ), INTENT(OUT) :: h
      real(KIND=8), INTENT(OUT) :: sse
!
!              local declarations
!
!                          stopping constants
      real(KIND=8), parameter :: epsc = 1.d-6
      real(KIND=8), parameter :: epsx = 1.d-6
!
      real(KIND=8), dimension(n,p) :: gbig
      real(KIND=8), dimension(p) :: t,tcopy,xscale,xold,xtrow
      real(KIND=8), dimension( (p*(p+1))/2 ) :: hcopy
      real(KIND=8) pxr,s,sset,sr,lambda,costh,delta,rho,v,  &
     &      gam,mu,srlam,a,b,c,d,z,lwrb,uprb,phi,phip
!
      integer, parameter :: kupmax = 9  ! tries to get good lambda
      integer i,j,k,kit,kup,kcollin,jm1,jp1,ip1
!              
!              formats
!
  21  format(/' Iteration',i4,'  SSE',d18.8)
  22  format(' Current Estimates ',5f12.7/2x,5f12.7)
  23  format(' Step size, Delta',2f16.8)
  24  format(' Step size too small (failure) -- terminated '//)  
  25  format('  Warning -- have collinearity problem',i3,2f12.7)
  26  format(' SSE going up too much -- terminated'//)
  27  format(' Cosine of angle between fitted and resid',f8.4)
  28  format(' Too many (',i6,') iterations -- terminated'//)
  29  format(' Residuals orthogonal to fitted surface -- terminated'//)
  30  format(' Step size too small (success) -- terminated'//)
!
!                          initial range of trust region
      delta = 100.d0
      h = 0.d0          !  initialize Hessian *** first fix, Aug 2005
!                          initialize relative change for use as flag
      rho = 1.d0
!                          start here -- get residuals
      call res(x,r)
!                          compute sse here
      sse = 0.d0
      do i = 1,n
      sse = sse + r(i)*r(i)
      end do  ! loop on i
!                           initialize scale to zero
      xscale = 0.d0
!                           start main loop *************************
      do kit = 1,mit
      write(iout,21) kit,sse
      write(iout,22) (x(j),j=1,p)
      if( rho .gt. 1.d-4 ) then         ! update jacobian
!                           get jacobian matrix
      call jac(x,gbig)
!                           initialize L-M parameter lambda
      lambda = 1.d-6
!                           initialize count of L-M iterations
      kup = 0
!                           set instability flag to none
      kcollin = 0
!                           do Modified Gram-Schmidt 
!
      do j = 1,p
!                      normalize
      s = 0.d0
      do i = 1,n
      s = s + gbig(i,j)*gbig(i,j)
      end do  ! loop on i
!                      compute other ss from this col
      sr = 0.d0
      if( j .gt. 1 ) then
      jm1 = j - 1
      do k = 1,jm1
      sr = sr + h( (j*(j-1))/2 + k )*h( (j*(j-1))/2 + k )
      end do  ! loop on k
                     end if  ! ( j .gt. 1 )
!                      if not much variation left in 
!                       this col (high R**2) then collinear
      if( sr .gt. 10000.d0*s ) write(iout,25) j,sr,s
      if( sr .gt. 10000.d0*s ) kcollin = kcollin + 1
!                      if R**2 > .9999 write warning
      sr = sr + s
!                      use this for scaling
      xscale( j ) = max( sqrt(sr), xscale(j) )
      s = sqrt(s)
      h( (j*(j+1))/2 ) = s 
!                      have norm, normalize
      do i = 1,n
      gbig(i,j) = gbig(i,j)/s
      end do  ! loop on i
!                      MGS -- work on other cols
      if( j .eq. p ) exit
      jp1 = j + 1
      do k = jp1,p
!                      work on column k -- first inner prod
      s = 0.d0
      do i = 1,n
      s = s + gbig(i,j)*gbig(i,k)
      end do  ! loop on i
!                      store in h
      h( (k*(k-1))/2 + j ) = s
!                      now take out of col k
      do i = 1,n
      gbig(i,k) = gbig(i,k) - s*gbig(i,j)
      end do  ! loop on i
!                      done with column k
      end do  ! loop on k
!                      done with MGS -- have R in h, Q in ws
      end do  ! loop on j
!                      next get Q*resids
      pxr = 0.d0
      do j = 1,p
      s = 0.d0
      do i = 1,n
      s = s + gbig(i,j)*r(i)
      end do  ! loop on i
      t(j) = s
!                      save an extra copy
      tcopy(j) = s
!                      used for termination criterion
      pxr = pxr + s*s
!
      end do  ! loop on j
!                      cosine of angle between resids and fitted
      costh = sqrt(pxr/sse)
      write(iout,27) costh
!                      if near zero -- stop -- done
      if( costh .lt. epsc ) then
!                      stop -- resids orthogonal to surface
      termcd = 1
      write(iout,29)
      return
                            end if  ! ( costh .lt. epsc )
!                      if no collinearity problem, set lambda to zero
      if( kcollin .eq. 0 ) lambda = 0.d0
!                      copy h for safe keeping
      hcopy = h
!                      end of update of jacobian
                           end if       ! ( rho .le. 1.d-4 )
!                      loop for iterating for lambda
      do kup = 1,kupmax
!
!                      solve ridge-regression-like normal eqns
!
      srlam = sqrt(lambda)
!                      this will be numerator for More's gamma 
      gam = 0.d0
      do j = 1,p
      do i = j,p
!                      add another observation row -- from lambda*D
      xtrow(i) = 0.d0
      end do  ! loop on i
!                      put scale part here
      xtrow(j) = srlam*xscale(j)
!                      need new value for rhs
      z = 0.d0
!                      rotate out new obs start with diagonal j
      do i = j,p
!                      get rotation to eliminate in col i
      a = h( (i*(i+1))/2 )
      b = xtrow(i)
      call rot734(a,b,c,s,d)
      h( (i*(i+1))/2 ) = d
      xtrow( i ) = 0.d0
!                      apply rotation to remainder of row i
!                       and rest of new observation's vector
      if( i .lt. p ) then
      ip1 = i + 1
      do k = ip1,p
      a = c*h( (k*(k-1))/2 + i ) + s*xtrow(k)
      b = c*xtrow(k) - s*h( (k*(k-1))/2 + i )
      h( (k*(k-1))/2 + i ) = a
      xtrow(k) = b
      end do  ! loop on k
                     end if ! ( i .lt. p )
!                      apply to rhs
      a = c*t(i) + s*z
      b = c*z - s*t(i)
      t(i) = a
      z = b
!                      done with col j, go to next col in row j
      end do  ! loop on i
!                      done with new observation row j
      gam = gam + t(j)*t(j)
      end do  ! loop on j
!                      now solve for next step t
      call chlzih(h,p,t)
!                      have step, now what?
!                        copy and move
!                        old x to end of ws, new x is x + step
      s = 0.d0
      do j = 1,p
!                        copy old x
      xold(j) = x(j)
!                        new x
      x(j) = x(j) - t(j)
!                        compute length
      s = s + (t(j)*xscale(j))*(t(j)*xscale(j))
      t(j) = - t(j)
      end do  ! loop on j
      write(iout,22) (x(j),j=1,p)
      sr = sqrt(s)
      write(iout,23) sr,delta
!
!     *** big branch ***
!                       is step length the right size?
!     *** big branch ***
!
      if( (sr .le. 1.5*delta) .and. (lambda .le. 1.1d-4) ) exit
      if( (sr .le. 1.5*delta) .and. (sr .ge. .75*delta) ) exit
!
!                       *** find new value for lambda ***
!
!                       compute phi and phip
      phi = sr - delta
      do j = 1,p
      t(j) = xscale(j)*xscale(j)*t(j)
      end do  ! loop on j
      call chlzhi(h,p,t)
      phip = 0.d0
      do j = 1,p
      phip = phip + t(j)*t(j)
      end do  ! loop on j
      phip = phip / sr
!                       first time through or not
      if( kup .eq. 1 ) then
!                       initialize upper and lower bounds
      lwrb = lambda
!                       l0 easy, u0 = norm Dinv*J'r = norm Dinv*R'Q'r
      uprb = 0.d0
      do j = 1,p
      v = 0.d0
      do k = 1,j
      v = v + hcopy(((j-1)*j)/2 + k)*tcopy(k)
      end do  ! loop on k
      v = v / xscale(j)
      uprb = uprb + v*v
      end do  ! loop on j
      uprb = sqrt(uprb) / delta
!                       done with initial values, skip to lambda
                       else ! ( kup .eq. 1 )
!                       update upper and lower bounds
      lwrb = max( lwrb, lambda+(phi/phip) )
      if( phi .lt. 0.d0 ) uprb = lambda
                       end if ! ( kup .eq. 1 )
!                       get new lambda
      lambda = lambda + (sr/delta)*(phi/phip)
      if( (lambda .gt. uprb) .or. (lambda .lt. lwrb) ) lambda =       &
     &                       max( 1.d-3*uprb, sqrt( lwrb*uprb ) )
!                       have new lambda -- go and do again
!                       copy back
      k = 0
      do j = 1,p
      t(j) = tcopy(j)
      x(j) = xold(j)
      do i = 1,j
      k = k + 1
      h(k) = hcopy( k )
      end do  ! loop on i
      end do  ! loop on j
!                      go back and get new step 
      end do  ! loop on kup -- finding new lambda
!                      if too many ups, stop -- flustered
      if( kup .ge. kupmax ) then
!                      stop -- SSE going up instead of down
      termcd = 6
      write(iout,26)
      return
                            end if ! ( kup .ge. kupmax )
!
!                      if step size looks ok, then step ahead
!                      get new residuals
      call res(x,r)
!                      does sse go up or down
      sset = 0.d0
      do i = 1,n
      sset = sset + r(i)*r(i)
      end do  ! loop on i
      write(iout,21) kit,sset
!                      check this out
      rho = 0.d0
      mu = 0.5d0
      if( sset .le. sse ) then
!                      get constant gam (actually More's -gamma)  
      gam = gam/sse
!                      get predicted decrement in ss  
      rho = (1.d0 - (sset/sse)) / ( gam + lambda*(s/sse) )
!                      update for delta
      mu = 0.d0
      if( sset .lt. sse ) mu = .5d0
      if( sset .gt. 10.d0*sse ) mu = .1d0
      if( mu .eq. 0.d0 ) mu = gam/(2.d0*gam - (1.d0 - (sset/sse)) )
                          end if ! ( sset .le. sse )
!                      update delta      
      delta = delta*mu
!                      if things are bad, drop delta and go back
      if( rho .lt. 1.d-4 ) then
      x = xold          !     go back to xold *** second fix, Aug 2005
!                       quit (failure) if delta too small
      if( delta .lt. epsx ) then
!                      stop -- steps too small  -- failure
      termcd = 5
      write(iout,24)
      return
                            end if  ! ( delta .lt. epsx ) failure
                           else     ! ( rho .lt. 1.d-4 )
!                                     normal situation
!                      if things are real good, use length
      sr = sqrt(s)
      if( (rho .gt. 0.25d0) .and. (lambda .eq. 0.d0) ) delta = 2.d0*sr
      if( rho .gt. .75d0 ) delta = 2.d0*sr
      sse = sset
!                      quit if delta too small
      if( delta .lt. epsx ) then
!                      stop -- steps too small -- good
      termcd = 2
      write(iout,30)
      return
                            end if  ! ( delta .lt. epsx )
                           end if     ! ( rho .lt. 1.d-4 )
!
!                      end of main iteration loop *******************
      end do  ! loop on kit
!                      terminated -- too many iterations
      termcd = 4
      write(iout,28) mit
      return
      end subroutine nllsq
      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
      SUBROUTINE CHLZOI(A,N)
!
!   OVERWRITES L WITH INVERSE(L')*INVERSE(L) WHERE L IS A LOWER
!   TRIANGULAR (N BY N) MATRIX STORED IN A
!
!   TO COMPUTE THE INVERSE OF A POSITIVE DEFINITE MATRIX A
!     FIRST:  FACTOR A = L * L' WHERE L IS LOWER TRIANGULAR
!             USING CHLSKY(A,N,D,IDET)
!     SECOND: COMPUTE THE INVERSE(A) = INVERSE(L') * INVERSE(L)
!             USING CHLSOI(A,N)
!
! J F MONAHAN (JULY,1984) DEPT OF STATISTICS, N C S U, RALEIGH, N C 2769
!     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) S
      INTEGER I,J,K,IM1,KP1
!                     FIRST INVERT L IN PLACE
      DO K = 1,N
!                     SOLVE L * X = E(K)
      A( (K*(K+1))/2 ) = 1./A( (K*(K+1))/2 )
!                     BRANCH AROUND IF K EQUALS N
      IF( K .LT. N ) THEN
        KP1 = K + 1
        DO I = KP1,N
        S = 0.
        IM1 = I - 1
          DO J = K,IM1
          S = S - A( (I*(I-1))/2 + J )*A( (J*(J-1))/2 + K )
          ENDDO ! LOOP ON J
        A( (I*(I-1))/2 + K ) = S/A( (I*(I+1))/2 )
        ENDDO ! LOOP ON I
                      ENDIF
      ENDDO         ! LOOP ON K
!                    DONE WITH INVERSION, NOW COMPUTE INNER PRODUCT
      DO J = 1,N
        DO I = J,N
!                    GET (I,J) ELEMENT OF MATRIX IS THE INNER PRODUCT
!                    OF COLUMN I AND COLUMN J
        S = 0.
          DO K = I,N
          S = S + A( (K*(K-1))/2 + I )*A( (K*(K-1))/2 + J )
          ENDDO     ! LOOP ON K
!                    STORE THE INNER PRODUCT
        A( (I*(I-1))/2 + J ) = S
        ENDDO       ! LOOP ON I
      ENDDO         ! LOOP ON J
      RETURN
      END SUBROUTINE CHLZOI
      SUBROUTINE ROT734(ALPHA,BETA,GAMA,SIGMA,RNU)
!  GIVEN ALPHA AND BETA, PRODUCES GIVENS TRANSFORMATION VALUES
!
!      (  GAMA       SIGMA  )  (  ALPHA  )       (  RNU  )
!      (                    )  (         )   =   (       )
!      ( -SIGMA       GAMA  )  (  BETA   )       ( ZERO  )
!
!  ALGORITHM 7.3.4 FROM STEWART
!
      REAL(KIND=8), INTENT(IN) :: ALPHA,BETA
      REAL(KIND=8), INTENT(OUT) :: GAMA,SIGMA,RNU
      REAL(KIND=8) ETA,DELTA
!
      ETA = MAX( ABS(ALPHA), ABS(BETA) )
!                     IF BOTH ALPHA AND BETA ARE ZERO THEN
!                     USE IDENTITY TRANSFORMATION
      IF( ETA .EQ. 0. ) then
!                     SET UP IDENTITY TRANSFORMATION IF NO WORK DONE
      RNU = 0.
      GAMA = 1.                                                         
      SIGMA = 0.                                                        
      RETURN
                        else
      DELTA = (ALPHA/ETA)*(ALPHA/ETA) + (BETA/ETA)*(BETA/ETA)           
      DELTA = SQRT(DELTA)                                               
      RNU = ETA * DELTA                                                 
      GAMA = ALPHA / RNU                                                
      SIGMA = BETA / RNU
                        END if
      RETURN                                                            
      END subroutine rot734
! *** end of filename nllsqu1.f95       *********************