C MRN.for:
C
C Contains routines that implement the modified Newton-Raphson method with line search
C
C The module MNR stores default values, the subroutine FixvMNR does the
C iterations and the subroutine MNRStep does the line search.
C
C Output options use the commands provided by Compaq Fortran.
C
C Fortran 95 Language options are used in the program code
C
C If you use a different compiler, out-comment the GetTextposition, SetTextposition
C command and write or print to the standard output device.


      Module MNR

      Integer(4), Parameter :: MNR_Maxit=500 ! maximum number of iterations
	Real(8), Parameter :: MNR_ParTol=1.0e-11  ! convergence criterium 2
	Real(8), Parameter :: MNR_GradTol=1.0-7   ! convergence criterium 3
	Real(8), Parameter :: MNR_FITol=1.0e-8  ! convergence criterium 1
	
	Integer(4) Jac
	Character(50) MNR_Message(0:6)

      Data Jac/3/ ! Central difference approximation of Jacobian
	Data MNR_Message(0)/'Normal Termination'/
	Data MNR_Message(1)/'Not able to evaluate F'/
	Data MNR_Message(2)/'Not able to evaluate Jacobian'/
	Data MNR_Message(3)/'Not able to solve for dx'/
      Data MNR_Message(4)/'Minimal step size in line search'/
	Data MNR_Message(5)/'Convergence to minimum'/
	Data MNR_Message(6)/'Maximum number of iterations exceeded'/

	Logical MNR_Print

	Data MNR_Print /.true./

	End Module MNR


C*********************************************************************
C
C FixvMN1
C 
C Purpose: Find the solution to f(x)=0
C
C Usage: Call FixvMN1(N,x1,F,x2,crit)
C
C Input:  N, Integer(4), number of elements in x
C        x1, Real(8) x1(N) column vector
C         F, pointer the function that computes f(x),
C         the call must be CALL F(N,N,x,fx), where fx is N x 1
C
C Output: x2, Real(8), N x 1 vector with solution.
C         rc, Integer(4), return code (see below)
C        crit, Real(8), 6 times 1 vector, where:
C              crit(1)=exit condition:
C                                      0: normal termination
C                                      1: not able to evaluate F
C                                      2: not able to evaluate Jacobian
C                                      3: matrix singular
C                                      4: minimal step size reached
C                                      5: parameter convergence without F(x)=0
C                                      6: maximal number of iterations reached
C
C              crit(2)=maximum absolute value of F at x2
C              crit(3)=maximum absolute value of relative change of x
C              crit(4)=maximum absolute value of relative gradient at x2
C              crit(5)=0.5*fx'*fx
C              crit(6)=number of iterations
C
C Remarks: If an error occurs during execution, the program
C          terminates.
C
C          The code that evaluated F must include the statement
C          use Errors and must set FERROR to a non-zero value, if
C          it is not possible to evaluate F at x1.
C
C          At minimum, the module Errors must be defined as follows
C
C          Module Errors
C             Integer(4) FError
C          End Module Errors
C
C The program uses the DLSARG routine from the IMSL library. If you cannot
C link to this library, use the routine DGESV from the LAPACK,
C which is freely available from www.netlib.org.


      Subroutine FixvMN1(N,x0,F,x1,crit)
 
      use MNR
      use IMSLF90 
	use DFLIB
      use Errors

      Integer(4) N
	Real(8) x0(N), x1(N), crit(1:6)

	Integer(4) fh, iercd
	Real(8) xscale(N), epsfcn, df(N,N), fx(N), dx(N)
      Real(8) step, mstep, dg(N), g1
      Real(8) MinStep2, GradTest, ParTest

      Parameter (IPATH=1)

      TYPE (rccoord) curpos	
		
	External F

C	Initialize 
      xscale=1.0     ! no scaling
	epsfcn=0.0
	x1=x0      
	FError=0 ! Program stops if a call to F produces an error code	
	
	fh=GetActiveQQ()

C Start Iterations
      crit=1.0
	
	Call Erset(0,0,0)	
      FError=0
      CALL F(N,N,x1,fx)
	if (FError.ne.0) then
	   crit(1)=1
	   return
	Endif
	crit(2)=MaxVal(dabs(fx))
	g1=0.5*Dot_Product(fx,fx)	
	crit(5)=g1
      if (crit(2).lt.MNR_FITol) then
	   crit(1)=0.0
         crit(3)=0.0
	   crit(6)=0.0
	   return
      endif
	
	Call ClearScreen($GWINDOW)
	CALL SETTEXTPOSITION(1,1,curpos)
	WRITE(fh,*) 'Algorithm: modified Newton-Rhapson with line search'
	
	do while (crit(6) .lt. dble(MNR_MaxIt))

		Call SetTextPosition(3,2,curpos)
		WRITE(fh,"('Iteration No.: ',I4,' FError=',I2)") IDINT(crit(6)), FError
		if (MNR_Print) then
			WRITE(11,"('MNR-Iteration #: ',I4,' FError=',I2)") IDINT(crit(6)), FError
		Endif
		Call SetTextPosition(4,2,curpos)
		WRITE(fh,"('max|fx|= ',E15.8,', g(x)=',E15.8,' max|dx/x|=',E15.8)")  crit(2), crit(5), crit(3)
		if (MNR_Print) then
	     	   WRITE(11,"('max|fx|= ',E15.8,', g(x)=',E15.8,' max|dx/x|=',E15.8)") crit(2), crit(5), crit(3)
		Endif
		if (FError .ne. 0) then
			crit(1)=1.
			exit
		endif
		Select Case (Jac)
			case(1)
				Call DFDJAC(F,N,N,x1,xscale,fx,epsfcn,df,N)	
			case(2)
				Call MyJac(n,n,x1,f,df)
			case(3)
 				Call CDJac(n,n,x1,f,df)
		End Select	   
		if (MNR_Print) then
			WRITE(11,"('Error Code from F after call to CDFJAC ',I2)") FError	   
		Endif
		if (FError .ne. 0) then
			crit(1)=2.
			exit
		endif
		Call DLSARG(N,df,N,-fx,IPATH,dx)     	   
		if (iercd() .eq. 2) then
			crit(1)=3.
			exit
		endif
		dg=MatMul(transpose(df),fx)	   
		crit(4)=GradTest(n,dg,x1,crit(5))
		mstep=MinStep2(n,x1,dx)
		if (MNR_Print) then
			WRITE(11,"('GradTest=',F15.10)") crit(4)
			WRITE(11,"('MinStep =',E15.4)") mstep
		Endif
C Compute step size
		Call MNRStep(n,x1,dx,dg,g1,mstep,F,step,fx)
		Call SetTextPosition(5,2,curpos)
		WRITE(fh,"('Step= ',E15.8)") step
		if (MNR_Print) then
			WRITE(11,"('Step= ',E15.8)") step	     
		Endif
		if ((mstep .le. 1.0) .and. (step .lt. mstep)) then
			if (crit(2) .lt. MNR_FITol) then
				crit(1)=0.0
			elseif (crit(4) .lt. MNR_GradTol) then
				crit(1)=5.0
			else
				crit(1)=4.0
			endif
			exit
		endif			  	
		crit(3)=ParTest(n,x1,x1+step*dx) 
		x1=x1+step*dx         
		Call F(N,N,x1,fx)
		crit(2)=MaxVal(dabs(fx))         
		crit(5)=0.5*Dot_Product(fx,fx)
		g1=crit(5)
		crit(6)=crit(6)+1.0
		if (crit(2) .lt. MNR_FITol) then
			crit(1)=0.0
			exit
		Endif
      End Do

      if (crit(6) .ge. dble(MNR_maxit)) then
		crit(1)=6.
      endif

      Call Erset(0,2,2) !restore default error behavior of dlsarg
	return
      End Subroutine


C MNRStep
C
C Computes the step size so that the Newton-Raphson algorithm
C always moves in the direction of a (local) minimum of (1/2)(f(x)'f(x))
C
C usage: Call MNRStep(N,x,dx,dg,g0,sl,F,s,fx)
C
C Input: N, Integer(4) the number of elements in x,
C        x, Real(8) N times 1 vector, the current parameter estimate
C       dx, Real(8) N times 1 vector, the current change of the parameter estimate
C       dg, Real(8) N times 1 vector, the gradient of (1/2)(f(x)'f(x)) at x
C       g0, Real(8) (1/2)*f(x)'f(x)
C       sl, Real(8), the minimal step size
C        F, pointer the function that returns the right hand side
C           of the system of functions whose zero is to be located,
C           the call must be F(N,N,x,fx), where fx is N times 1
C           Real(8) and holds the right hand side of the system of equations
C
C Output: s, Real(8) the step size found by the algorithm
C        fx, Real(8) N times 1 vector, f(x+s*dx) at return

                 
      Subroutine MNRStep(N,x,dx,dg,g0,sl,F,s,fx)

      use DFLIB
      use MNR
	use Errors

	Integer(4) N
	Real(8) x(N), dx(N), dg(N), sl, s, g0, fx(n)

      Integer(4) fh
      Real(8) smult, smin, smax, disc, s1, s2, amat(2,2), 
     +        bvec(1:2,1), ab(1:2,1), x1(N), g1, g2, dgdx

      External F

      Type (rccoord) curpos		

	fh=GetActiveQQ()
C  Initialize 
      smult=1.0E-4
	smin=0.1E0
	smax=0.5E0
     
	dgdx=Dot_Product(dg,dx)

      s1=1.0
	s=1.0

C Reduce step size until no function error occurs
      FError=0
      Call F(N,N,x+dx,fx)
	Do while ((FError .ne. 0) .and. (s1 .gt. sl))
	   Call SetTextPosition(3,29,curpos)
	   Write(fh,"(I2)") FError
         s1=0.75*s1
         FError=0
	   Call F(N,N,x+s1*dx,fx)
	End Do
	If (s1 .lt. sl) return
      dx=s1*dx
	s1=1.0

      g1=0.5*Dot_Product(fx,fx)

C Try the full Newton step s=1 
      if (g1 .le. (g0 + smult*dgdx)) then
        s=s1
	  return
      else
        s=-dgdx/(2.0*(g1-g0-dgdx))
        if (s .lt. smin) s=smin
        if (s  .gt. smax) s=smax
        x1=x+s*dx

        Call F(N,N,x1,fx)
	  g2=0.5*Dot_Product(fx,fx)

      endif

      s2=s
	
C  Reduce s2 further unless g2 < g0 + s2*smult*dgdx 
      do while (g2 .gt. (g0 + smult*s2*dgdx) )

        amat(1,1)=1.0/(s2**2)
	  amat(1,2)=-1.0/(s1**2)
	  amat(2,1)=-s1/(s2**2)
	  amat(2,2)=s2/(s1**2)
        bvec(1,1)=g2-s2*dgdx-g0
	  bvec(2,1)=g1-s1*dgdx-g0
	  ab=MatMul(amat,bvec)
	  ab=ab/(s2-s1)
	  if (ab(1,1) .eq. 0.0) then
            s=-dgdx/(2*ab(2,1))
	  else
	      disc=(ab(2,1)**2)-3*ab(1,1)*dgdx	
	      if (disc .lt. 0.0) then
	         s=s2*smax
	       elseif (ab(2,1) .le. 0.0) then                 
	         s=(-ab(2,1)+dsqrt(disc))/(3*ab(1,1))
	       else
	         s=-dgdx/(ab(2,1)+dsqrt(disc))
	       endif
	  endif	          
      
        if (s .lt. s2*smin) s=s2*smin
        if (s .gt. s2*smax) s=s2*smax

        if (s .lt. sl) return
	
        s1=s2
        s2=s
        g1=g2
        x1=x+s2*dx
C	  write(11,*) x1
        Call F(N,N,x1,fx)
        g2=0.5*Dot_Product(fx,fx)

      End Do

      Return

      End Subroutine
      
C ***********************************************************************************
C
C MyJac: Compute Jacobian of f using forward differences
C
C Usage: Call MyJac(m,n,x,f,df)
C
C Input: m, integer(4), the number of elements in x
C        n, integer(4), the number of elemtend returned by f(m,n,x,fx) in fx
C        x, real(8), m times 1 vector
C        f, pointer to the function f(m,n,x,fx)
C
C Output: df, real(8) n times n matrix
C
C Remarks: Code adapted from Press et al. (2001),p. 381

      Subroutine MyJac(m,n,x,f,df)

	Integer(4) n, m
	Real(8) x(m), df(n,m)

      Integer(4) i, j
	Real(8) eps, eps1, h, temp, fx1(n), fx2(n), x1(m)
	External f

      x1=x
      Call f(n,n,x,fx1)

      eps1=Epsilon(eps1)
	eps=Epsilon(eps)**(0.5)
      
      do j=1,m
	   temp=x1(j)
	   h=eps*abs(temp)
C	   if (h .eq. 0.0) h=eps
         if (h.le.eps1) h=eps
	   x1(j)=temp+h
	   h=x1(j)-temp
	   Call F(m,n,x1,fx2)
         x1(j)=temp
	   do i=1,n
	      df(i,j)=(fx2(i)-fx1(i))/h
	   end do
	end do

	return

	End Subroutine

C ***********************************************************************************
C
C CDJac: Compute Jacobian of f using central differences
C
C Usage: Call CDJac(m,n,x,f,df)
C
C Input: m, integer(4), the number of elements in x
C        n, integer(4), the number of elements returned by f(m,n,x,fx) in fx
C        x, real(8), m times 1 vector
C        f, pointer to the function f(m,n,x,fx)
C
C Output: df, real(8) n times n matrix
C
C Remarks: This is algorithm A.5.6.4. from Dennis and Schnabel (1983), p. 323

      Subroutine CDJac(m,n,x,f,df)

      use Errors
	Integer(4) n, m
	Real(8) x(m), df(n,m)

      Integer(4) i, j
	Real(8) eps, eps1, h, temp, fx1(n), fx2(n), x1(m), x2(m)
	External f

      x1=x      
      x2=x

      eps1=Epsilon(eps1)
	eps=Epsilon(eps)**(1.0/3.0)
      
      do j=1,m
	   temp=x1(j)
	   h=eps*dabs(temp)
	   if (h .le. eps1) h=eps
	   x1(j)=temp+h
	   x2(j)=temp-h
	   h=x1(j)-temp
	   Call F(m,n,x1,fx1)
	   if (FError.ne.0) return
	   Call F(m,n,x2,fx2)     
	   if (FError.ne.0) return    
	   do i=1,n
	      df(i,j)=(fx1(i)-fx2(i))/(2.*h)
	   end do
	   x1(j)=x(j)
	   x2(j)=x(j)
	end do

	return

	End Subroutine


C***********************************************************************
C
C MinStep2: Returns the minimal step size sl so that x1=x0+s*dx=0
C           (in terms of the parameter tolerance criterium)
C          
C Usage: sl=MinStep2(nx,x0,dx)
C
C Input: nx, Integer(4), the number of elements in x0
C        x0, Real(8), nx times 1 vector
C        dx, Real(8), nx times 1 vector, change of x
C
C Output: sl, Real(8)

	Real(8) Function MinStep2(nx,x0,dx) 

      use MNR
	Integer(4) nx
	Real(8) x0(nx), dx(nx)

	Integer(4) i
	Real(8) converge, temp

      converge=0.0

	Do i=1,nx
	   temp=dabs(dx(i))/dmax1(dabs(x0(i)),1.0)
	   if (temp .gt. converge) converge=temp
      end do

	MinStep2=MNR_ParTol/converge

	Return

	End Function
                    
C ParTest: implements the parameter convergence criterium
C          from Dennis and Schnabel (1983), p. 160
C
C usage: converge=ParTest(nx,x0,x1)
C
C Input: nx, Integer(4), the number of elements in x0, x1, and typx
C        x0, Real(8) vector with nx elements, the first value of x
C        x1, Real(8) vector with nx elements, the second value of x
C
C Output: Real(8), the minimal value of the convergence criterium

      Real(8) Function ParTest(nx,x0,x1) 

	Integer(4) nx
	Real(8) x0(nx), x1(nx)

	Integer(4) i
	Real(8) converge, temp

      converge=0.0

	Do i=1,nx
	   temp=dabs(x1(i)-x0(i))/dmax1(dabs(x1(i)),1.0)
	   if (temp .gt. converge) converge=temp
      end do

	ParTest=converge

	Return

	End Function


C*********************************************************************
C
C FixvMN2
C 
C Purpose: Find the solution to f(x)=0
C
C Usage: Call FixvMN2(N,x1,bounds,F,x2)
C
C Input:  N, Integer(4), number of elements in x
C        x1, Real(8) x1(N) column vector
C    bounds, Real(8), N x 2 matrix with lower and upper bounds on x
C         F, pointer the function that computes f(x),
C         the call must be CALL F(N,x,fx), where fx is N x 1
C
C Output: x2, Real(8), N x 1 vector with solution.

      Subroutine FixvMN2(N,x0,bounds,F,x2)

      use IMSLF90
	use Errors
	use MNR
	 
      Integer(4) N
	Real(8) x0(N), x1(N), bounds(1:N,1:2), x2(N)

	Integer(4) itn
	Real(8) crit, xscale(N), epsfcn, df(N,N), fx(N), dx(N),
     +        step, CheckBounds

	External F

C	Initialize             
      xscale=1.0     ! no scaling
	epsfcn=0.0
	x1=x0
      itn=1           
C Start Iterations
      crit=1.0
      do while ((crit .gt. MNR_FITol) .and. (itn .lt. MNR_Maxit))
	   CALL F(N,N,x1,fx)
	   if (FError.ne.0) return
         Call DFDJAC(F,N,N,x1,xscale,fx,epsfcn,df,N)		  
	   if (FError.ne.0) return
         Call DLSARG(N,df,N,-fx,1,dx)      
	   x2=x1+dx
         step=1
	   step=CheckBounds(n,x1,dx,bounds)
         x2=x1+step*dx
         CALL F(N,N,x2,fx)
	   if (FError.ne.0) return	   
	   crit=MaxVal(abs(fx))
         x1=x2;
        itn=itn+1
      End Do
      if (itn .gt. MNR_Maxit) then
	   print *, 'Maximum number of iterations exceeded'
	   print *, 'Returned solution may be missleaden'
      endif

      End Subroutine

C********************************************************************
C CheckBounds
C
C Purpose: Reduce the Newton-Raphson step so that
C          the point x_0 + s*dx remains in the domain
C          of the vector valued function f
C
C Usage:   s=CheckBoundsn,x,dx,bounds)
C
C Input:    n := integer(4), the number of variables
C           x := n times 1 real(8) vector, the inital point
C          dx := n times 1 real(8) vector, the Newton-Raphson step
C      bounds := n times 2 real(8) matrix , whose first (second) column
C                holds the lower (upper) bounds for x
C Output    s := Real(8) scalar
C

      Real(8) Function CheckBounds(n,x,dx,bounds)

      Integer(4) n, i
	Real(8) x(1:n), dx(1:n), bounds(1:n,1:2), l1(1:n), l2(1:n), s1, s2

	Do i=1,n
	   l1(i)=(bounds(i,1)-x(i))/dx(i)
	   l2(i)=(bounds(i,2)-x(i))/dx(i)
	End Do
      s1=MinVal(l1,mask=l1 .gt. 0.0)
	s2=MinVal(l2,mask=l2 .gt. 0.0)

      s1=Min(s1,s2)
	if (s1 .lt. 1.0) then
        CheckBounds=0.95*s1
	else
	  CheckBounds=1.0
	endif

	Return

	End Function	


C***************************************************************************
C
C GradTest: implements the gradient criterium
C           from Dennis and Schnabel (1983), p. 160
C
C usage: converge=GradTest(nx,df,x,fx)
C
C Input:   nx, Integer(4), the number of elements in x0, x1
C          df, Real(8) vector with nx elements, the gradient of
C              f at x
C           x,  Real(8) vector with nx elements, the first value of x
C          fx, the value of f at x
C
C Output: Real(8), the maximal value of the convergence criterium

      Real(8) Function GradTest(nx,df,x,fx)                        

	Integer(4) nx
	Real(8) df(nx), x(nx), fx

	Integer(4) i
	Real(8) converge, temp

      converge=0.0
	Do i=1,nx
	   temp=dabs(df(i))*dmax1(dabs(x(i)),1.0)
	   temp=temp/dmax1(dabs(fx),1.0)
         if (temp .gt. converge) converge=temp
      End Do

	GradTest=converge

	Return

	End Function