program newtonma

  USE HSL_ZD11_double
  use hsl_ma48_double
  use bmproblem, only: pdata_type,iniprob,endprob,evalr,evalj
  use iso_c_binding, only: c_ptr,c_loc

  implicit none

  ! LOCAL SCALARS
  integer :: allocstat,INFO,ierr,istop,iter,jnnz,maxit,n
  real(kind=8) :: eps,finish,start
  type(pdata_type), target :: pdata
  
  TYPE(ZD11_TYPE) :: MATRIX
  TYPE(MA48_FACTORS) FACTORS
  TYPE(MA48_CONTROL) CONTROL

  ! LOCAL ARRAYS
  real(kind=8), allocatable :: r(:),x(:)

  call iniprob(n,jnnz,c_loc(pdata))

  allocate(x(n),r(n),stat=allocstat)
  if ( allocstat .ne. 0 ) then
     write(*,*) 'Error in allocate.'
     stop
  end if
  
  x(1:n) = pdata%x0(1:n)
  
  matrix%n = n
  matrix%m = n
  matrix%ne = jnnz
  
  allocate(matrix%val(jnnz),matrix%row(jnnz),matrix%col(jnnz),stat=allocstat)
  if ( allocstat .ne. 0 ) then
     write(*,*) 'Error in allocate.'
     stop
  end if
  
  CALL MA48_INITIALIZE(FACTORS,CONTROL)
  
  eps = 1.0d-06 * sqrt( real( n ) )
  maxit = 10

  call cpu_time(start)

  call newton(n,x,r,eps,maxit,jnnz,iter,ierr,istop,c_loc(pdata))

  call cpu_time(finish)

  write(*,*) 'If this value is not zero, something went wrong = ',ierr
  write(*,*) 'Satisfied stopping criterion = ',istop
  write(*,*) 'Residual Euclidian norm = ',norm2( r(1:n) )
  write(*,*) 'Number of iterations = ',iter
  write(*,*) 'Number of functional evaluations = ',pdata%counters(1)
  write(*,*) 'Number of Jacobian evaluations = ',pdata%counters(2)
  write(*,*) 'CPU time in seconds = ',finish - start

  open(10,file='tabline.txt')
  write(10,9000) pdata%np,n,ierr,istop,norm2(r(1:n)),iter,pdata%counters(1),pdata%counters(2),finish-start
  close(10)
    
  deallocate(x,r,matrix%val,matrix%row,matrix%col,stat=allocstat)
  if ( allocstat .ne. 0 ) then
     write(*,*) 'Error in deallocate.'
     stop
  end if
  
  CALL MA48_FINALIZE(FACTORS,CONTROL,INFO)
  
  call endprob(c_loc(pdata))

  stop
  
9000 format(1X,I3,1X,I8,1X,I2,1X,I2,1X,1P,D7.1,3(1X,I6),0P,F12.6)

contains
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine newton(n,x,r,eps,maxit,jnnz,iter,ierr,istop,pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: jnnz,maxit,n
    integer, intent(out) :: ierr,istop,iter
    real(kind=8), intent(in) :: eps
    type(c_ptr), optional, intent(in) :: pdataptr

    ! ARRAY ARGUMENTS
    real(kind=8), intent(inout) :: x(n)
    real(kind=8), intent(out) :: r(n)

    ! LOCAL SCALARS
    logical :: firstcallma48
    integer :: FAST,i,iini,iend
    real(kind=8) :: rnorm
    
    TYPE(MA48_AINFO) AINFO
    TYPE(MA48_FINFO) FINFO
    TYPE(MA48_SINFO) SINFO
   
    ! LOCAL ARRAYS
    integer :: jsta(n),jlen(n)
    real(kind=8) :: rhs(n),s(n)
    
    ierr = 0
    iter = 0
    
    firstcallma48 = .true.
    
100 continue

    call evalr(n,x,r,ierr,pdataptr)
    if ( ierr .ne. 0 ) return
    
    rnorm = dot_product( r(1:n), r(1:n) )

    write(*,*) 'iter = ',iter,' rnorm = ',rnorm
    
    if ( rnorm .le. eps ** 2 ) then
       istop = 0
       return
    end if

    if ( iter .ge. maxit ) then
       istop = 1
       return
    end if

    iter = iter + 1
    
    call evalj(n,x,jsta,jlen,jnnz,matrix%col,matrix%val,ierr,pdataptr)
    if ( ierr .ne. 0 ) return

    do i = 1,n
       iini = jsta(i)
       iend = jsta(i) + jlen(i) - 1
       matrix%row(iini:iend) = i
    end do

    if ( firstcallma48 ) then
       CALL MA48_ANALYSE(MATRIX,FACTORS,CONTROL,AINFO,FINFO)
       IF(AINFO%FLAG .ne. 0) THEN
          WRITE(6,'(A,I2)') &
               ' Failure of MA48_ANALYSE with AINFO%FLAG=', AINFO%FLAG
          STOP
       END IF
    else
       CALL MA48_FACTORIZE(MATRIX,FACTORS,CONTROL,FINFO,FAST)
       IF(FINFO%FLAG .ne. 0) THEN
          WRITE(6,'(A,I2)') &
               ' Failure of MA48_FACTORIZE with FINFO%FLAG=', FINFO%FLAG
          STOP
       END IF
    end if

    rhs(1:n) = - r(1:n)
    
    CALL MA48_SOLVE(MATRIX,FACTORS,rhs,s,CONTROL,SINFO)
    IF(SINFO%FLAG .ne. 0) THEN
       WRITE(6,'(A,I2)') &
            ' Failure of MA48_SOLVE with SINFO%FLAG=', SINFO%FLAG
       STOP
    END IF

    x(1:n) = x(1:n) + s(1:n)

    go to 100
    
  end subroutine newton

end program newtonma
