module bmsesem

  use bmqr, only: qr,extractQ,extractR,qrupdate,qrsolve
  use bmproblem, only: pdata_type,evalr
  use iso_c_binding, only: c_ptr

  implicit none
  
  public :: sesem

contains

  ! ******************************************************************
  ! ******************************************************************

  subroutine sesem(n,x,l,u,m,res,ftarget,maxit,maxnprogr,nred,redscheme, &
       seed,ndiis,nhlim,iter,fcnt,fbest,iterbest,fcntbest,accacc,bobfai, &
       istop,pdataptr)

    implicit none
    
    ! SCALAR ARGUMENTS
    integer, intent(in) :: m,maxit,maxnprogr,n,nred,ndiis,nhlim,redscheme
    integer, intent(out) :: accacc,bobfai,fcnt,fcntbest,istop,iter,iterbest
    real(kind=8), intent(in) :: ftarget
    real(kind=8), intent(inout) :: seed
    real(kind=8), intent(out) :: fbest
    type(c_ptr), optional, intent(in) :: pdataptr

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: l(n),u(n)
    real(kind=8), intent(inout) :: x(n)
    real(kind=8), intent(out) :: res(m)
    
    ! PARAMETERS
    logical, parameter :: tabline = .false.
    integer, parameter :: iprint = 2
    real(kind=8), parameter :: macheps = epsilon( 1.0d0 ), hsmall = 1.0d-04, hlarge = 0.1d0

    ! LOCAL SCALARS
    logical :: diis,extra
    integer :: jtbr,maxrank,ncol,nh,nprogr,rank,nexti
    real(kind=8) :: alpha,f,f0,facc,ftrial
    
    ! LOCAL ARRAYS
    integer :: perm(nhlim-1)
    real(kind=8) :: col(m),resacc(m),restmp(m),restrial(m),xacc(max(n,nhlim-1)),xtmp(n), &
         xtrial(n),xhist(n,nhlim),rhist(m,nhlim),q(m,nhlim-1),r(nhlim-1,nhlim-1),gam(nhlim-1)

    ! EXTERNAL FUNCTIONS
    real(kind=8), external :: drand

    integer :: bobyqaiprint,i,j,k,maxfun,mred,npt
    real(kind=8) :: fbobyqa,rhobeg,rhoend,smllstre
    integer :: rind(m)
    real(kind=8) :: d(n),lred(nred),ured(nred),xred(nred),xt(n,nred),w(10000000)

    if ( any( l(1:n) .gt. - 1.0d0+20 .or. u(1:n) .lt. 1.0d+20 ) ) then
       write(*,*) 'This version of sesem is not able to solve problems with box constraints!'
       stop
    end if

    diis = .true.
    if ( ndiis .eq. 0 ) then
       diis = .false.
    end if
    
    iter = 0
    fcnt = 0
    nexti = 1
    nprogr = 0
    accacc = 0
    bobfai = 0
    
    call evalr(n,x,m,res,pdataptr)
    fcnt = fcnt + 1

    f = dot_product( res(1:m), res(1:m) )

    fbest = f
    fcntbest = fcnt
    iterbest = iter

    f0 = f
    
    nh = 1
    xhist(1:n,1) = x(1:n)
    rhist(1:m,1) = res(1:m)

    maxrank = 0

    smllstre = 0.1d0
    
100 continue

    if ( iprint .ge. 0 ) then
     ! write(*,*) 'iter = ',iter,' (fcnt = ',fcnt,') f (sum of squares) = ',f,' x = ',x(1:n),' res = ',res(1:m)
     ! write(*,*) 'iter = ',iter,' (fcnt = ',fcnt,') f (sum of squares) = ',f,' x = ',x(1:n)
       write(*,*) 'iter = ',iter,' (fcnt = ',fcnt,') f (sum of squares) = ',f
    end if

    if ( tabline ) then
       open(10,file='tabline-interrupted.txt')
       write(10,9000) n,m,9,f,0.0d0,iter,fcnt,fbest,iterbest,fcntbest
       close(10)
    end if

    if ( f .le. ftarget ) then
       if ( iprint .ge. 0 ) then
          write(*,*) 'success!'
       end if
       istop = 0
       return
    end if

    if ( iter .ge. maxit ) then
       if ( iprint .ge. 0 ) then
          write(*,*) 'maximum number of iterations reached!'
       end if
       istop = 1
       return
    end if

    if ( nprogr .gt. maxnprogr ) then
       if ( iprint .ge. 0 ) then
          write(*,*) 'lack of progress!'
       end if
       istop = 2
       return
    end if
       
    iter = iter + 1

    mred = m
    rind(1:mred) = (/ (j,j=1,mred) /)

    npt = 2 * nred + 1
    maxfun = 1000 * ( nred + 1 )

    rhobeg = 0.1d0
    rhoend = smllstre
    
    bobyqaiprint = 0
    
    if ( redscheme .eq. 1 ) then
       ! Computes the displacement in an affine subspace
       do j = 1,nred
          do i = 1,n
             xt(i,j) = 2.0d0 * drand(seed) - 1.0d0
          end do
       end do
       
       xred(1:nred) =   0.0d0
       lred(1:nred) = - 1.0d+20
       ured(1:nred) =   1.0d+20

    else if ( redscheme .eq. 2 ) then
       ! Computes the displacement using linear-splines with variable nodes
       if ( mod( nred, 2 ) .ne. 0 .or. nred .lt. 4 ) then
          write(*,*) 'In sesem, nred must be even and at least 4 when redscheme=2.'
          write(*,*) 'Value of nred = ',nred
          stop
       end if
       
       k = ( nred - 2 ) / 2

       do i = 1,k
          xred(i) = drand(seed)
       end do

       lred(1:k) = 0.0d0
       ured(1:k) = 1.0d0
       
       xred(k+1:2*k+2) =   0.0d0
       lred(k+1:2*k+2) = - 1.0d+20
       ured(k+1:2*k+2) =   1.0d+20
       
    else
       ! Computes the displacement using linear-splines with fixed nodes
       if ( nred .lt. 2 ) then
          write(*,*) 'In sesem, nred must be at least 2 when redscheme=3.'
          write(*,*) 'Value of nred = ',nred
          stop
       end if
       
       do i = 1,nred
          xred(i) = drand(seed)
       end do

       lred(1:nred) = - 1.0d+20
       ured(1:nred) =   1.0d+20
    end if
    
200 continue

    call bobyqa(nred,npt,xred,lred,ured,rhobeg,rhoend,bobyqaiprint,maxfun,w,fbobyqa,ftarget,calfun)
    smllstre = min( smllstre, rhoend )
        
    if ( redscheme .eq. 1 ) then
       d(1:n) = matmul( xt(1:n,1:nred), xred(1:nred) )
    else if ( redscheme .eq. 2 ) then
       call frompvtox2(nred,xred,n,d)
    else
       call frompvtox3(nred,xred,n,d)
    end if
    
    xtrial(1:n) = x(1:n) + d(1:n)

    call evalr(n,xtrial,m,restrial,pdataptr)
    fcnt = fcnt + 1

    ftrial = dot_product( restrial(1:m), restrial(1:m) )

    if (abs( ftrial - fbobyqa ) .gt. 1.0d-12 ) then
       write(*,*) 'There is something wrong: ftrial should coincide with fbobyqa!'
       stop
    end if
    
    if ( iprint .ge. 1 ) then
     ! write(*,*) '    rhobeg = ',rhobeg,' rhoend = ',rhoend,' ftrial+ = ',ftrial,' xtrial = ',xtrial(1:n),' restrial = ',restrial(1:m)
     ! write(*,*) '    rhobeg = ',rhobeg,' rhoend = ',rhoend,' ftrial- = ',ftrial,' xtrial = ',xtrial(1:n)
       write(*,*) '    rhobeg = ',rhobeg,' rhoend = ',rhoend,' ftrial+ = ',ftrial
    end if

    if ( .not. ftrial .le. f - sqrt( macheps ) * abs( f ) ) then
       if ( rhobeg .ge. 10.0d0 * sqrt( macheps ) .and. rhoend .ge. 10.0d0 * macheps ) then
          rhobeg = 0.1d0 * rhoend
          rhoend = rhobeg
       
          if ( iprint .ge. 1 ) then
             write(*,*) 'A new run of bobyqa with a smaller rhoend will be tried.'
          end if
          go to 200
       end if
    end if

    eta = max( 1.0d0, ( f0 - ftarget ) ) / ( 2 ** ( iter - 1 ) )
    
    if ( ftrial .le. max( ftarget, f + eta - sqrt( macheps ) * ( f - ftarget ) ) ) then
       if ( iprint .ge. 1 ) then
          write(*,*) 'The reduction algorithm computed xtrial with sufficient descent.'
       end if
       go to 300
    end if

    if ( iprint .ge. 1 ) then
       write(*,*) 'The reduction algorithm has failed.'
    end if

    bobfai = bobfai + 1
    
    do i = 1,n
       d(i) = 2.0d0 * drand(seed) - 1.0d0
    end do

    d(1:n) = 10.0d0 * d(1:n) / norm2( d(1:n) )

    alpha = 1.0d0

250 continue
    
    xtrial(1:n) = x(1:n) + alpha * d(1:n)
    
    call evalr(n,xtrial,m,restrial,pdataptr)
    fcnt = fcnt + 1

    ftrial = dot_product( restrial(1:m), restrial(1:m) )

    if ( iprint .ge. 1 ) then
       write(*,*) ' alpha = ',alpha,' ftrial = ',ftrial
    end if

    eta = 1.0d0 / ( 2 ** ( iter - 1 ) )
    
    if ( ftrial .le. max( ftarget, f + eta - sqrt( macheps ) * ( f - ftarget ) ) ) then
       if ( iprint .ge. 1 ) then
          write(*,*) 'The linear search was successful.'
       end if
       go to 300
    end if

    alpha = 0.5d0 * alpha
    go to 250
    
300 continue

    if ( ftrial .le. ftarget ) then
       go to 500
    end if

    if ( .not. ( diis .and. mod( iter, ndiis ) .eq. 0 ) ) then
       go to 500
    end if

    if ( iprint .ge. 2 ) then
       write(*,*) 'Including xtrial in history.'
    end if
       
    if ( nh + 1 .le. min( nhlim, min( n + 1, m + 1 ) ) ) then
       if ( iprint .ge. 2 ) then
          write(*,*) 'There is space to include xtrial without removing any other point ', &
               '(nh = ',nh,' min( nhlim, min( n + 1, m + 1 ) ) = ',min( nhlim, min( n + 1, m + 1 ) ),').'
       end if
       
       nh = nh + 1
       jtbr = 0
       
    else
       if ( iprint .ge. 2 ) then
          write(*,*) 'History is full. Oldest point is being removed ', &
               '(nh = ',nh,' min( nhlim, min( n + 1, m + 1 ) ) = ',min( nhlim, min( n + 1, m + 1 ) ),').'
       end if
       
       xhist(1:n,1:nh-1) = xhist(1:n,2:nh)
       rhist(1:m,1:nh-1) = rhist(1:m,2:nh)
       jtbr = 1
    end if

    xhist(1:n,nh) = xtrial(1:n)
    rhist(1:m,nh) = restrial(1:m)
       
    col(1:m) = rhist(1:m,nh) - rhist(1:m,nh-1)

    extra = .false.
       
    if ( nh .eq. 2 ) then
       ncol = 1
       call qr(m,ncol,m,col,gam,perm,rank)
       call extractQ(m,ncol,m,col,gam,rank,m,Q)
       call extractR(ncol,m,col,rank,nhlim-1,R)
       maxrank = max( maxrank, rank )

       if ( iprint .ge. 2 ) then
          write(*,*) 'Number of columns = ',ncol,' rank = ',rank,' (largest achieved rank = ',maxrank,')'
       end if

    else
       call qrupdate(m,ncol,rank,nhlim-1,jtbr,m,Q,nhlim-1,R,perm,col)
       maxrank = max( maxrank, rank )

       if ( iprint .ge. 2 ) then
          write(*,*) 'Number of columns = ',ncol,' rank = ',rank,' (largest achieved rank = ',maxrank,')'
       end if

       if ( ncol .ne. nh - 1 ) then
          write(*,*) 'ERROR in subroutine sesem: ncol is not equal to nh - 1! ncol = ',ncol,' nh = ',nh
          stop
       end if
          
       if ( rank .lt. maxrank ) then
          if ( iprint .ge. 2 ) then
             write(*,*) 'Current rank is smaller that largest achieved rank. Thus, an extra point will be considerd.'
          end if
             
          extra = .true.

          xtmp(1:n) = x(1:n)
          xtmp(nexti) = x(nexti) + hsmall * max( 1.0d0, abs( x(nexti) ) )
          
          nexti = mod( nexti, n ) + 1

          call evalr(n,xtmp,m,restmp,pdataptr)
             
          if ( iprint .ge. 2 ) then
             write(*,*) 'Including extra point in history. (Considered ith canonical array with i = ',nexti,')'
          end if
             
          if ( nh + 1 .le. min( nhlim, min( n + 1, m + 1 ) ) ) then
             if ( iprint .ge. 2 ) then
                write(*,*) 'There is space to include the extra point without removing any other point ', &
                     '(nh = ',nh,' min( nhlim, min( n + 1, m + 1 ) ) = ',min( nhlim, min( n + 1, m + 1 ) ),').'
             end if
             
             nh = nh + 1
             jtbr = 0
          else
             if ( iprint .ge. 2 ) then
                write(*,*) 'History is full. Oldest point is being removed ', &
                     '(nh = ',nh,' min( nhlim, min( n + 1, m + 1 ) ) = ',min( nhlim, min( n + 1, m + 1 ) ),').'
             end if
                
             xhist(1:n,1:nh-1) = xhist(1:n,2:nh)
             rhist(1:m,1:nh-1) = rhist(1:m,2:nh)
             jtbr = 1
          end if

          xhist(1:n,nh) = xtmp(1:n)
          rhist(1:m,nh) = restmp(1:m)

          col(1:m) = rhist(1:m,nh) - rhist(1:m,nh-1)

          call qrupdate(m,ncol,rank,nhlim-1,jtbr,m,Q,nhlim-1,R,perm,col)
          maxrank = max( maxrank, rank )

          if ( iprint .ge. 2 ) then
             write(*,*) 'Number of columns = ',ncol,' rank = ',rank,' (largest achieved rank = ',maxrank,')'
          end if
          
          if ( ncol .ne. nh - 1 ) then
             write(*,*) 'ERROR in subroutine sesem: ncol is not equal to nh - 1! ncol = ',ncol,' nh = ',nh
             stop
          end if
          
       end if
    end if

    if ( rank .gt. 0 ) then
       call qrsolve(m,ncol,rank,m,Q,nhlim-1,R,perm,restrial,xacc)
       
       xacc(1:n) = xtrial(1:n) - matmul( xhist(1:n,2:nh) - xhist(1:n,1:nh-1), xacc(1:nh-1) )
       
       if ( extra ) then
          if ( iprint .ge. 2 ) then
             write(*,*) 'The extra point is being removed.'
          end if
          
          nh = nh - 1
          
          jtbr = ncol
          
          call qrupdate(m,ncol,rank,nhlim-1,jtbr,m,Q,nhlim-1,R,perm)
          maxrank = max( maxrank, rank )
          
          if ( iprint .ge. 2 ) then
             write(*,*) 'Number of columns = ',ncol,' rank = ',rank,' (largest achieved rank = ',maxrank,')'
          end if
          
          if ( ncol .ne. nh - 1 ) then
             write(*,*) 'ERROR in subroutine sesem: ncol is not equal to nh - 1! ncol = ',ncol,' nh = ',nh
             stop
          end if
       end if

       if ( norm2( xacc(1:n) ) .le. 10.0d0 * max( 1.0d0, norm2( x(1:n) ) ) ) then
          call evalr(n,xacc,m,resacc,pdataptr)
          fcnt = fcnt + 1
          
          facc = dot_product( resacc(1:m), resacc(1:m) )
          
          if ( iprint .ge. 1 ) then
           ! write(*,*) '    facc = ',facc,' xacc = ',xacc(1:n),' resacc = ',resacc(1:m)
           ! write(*,*) '    facc = ',facc,' xacc = ',xacc(1:n)
             write(*,*) '    facc = ',facc
          end if

          if ( facc .lt. ftrial .and. norm2( xacc(1:n) - x(1:n) ) .gt. macheps * max( 1.0d0, norm2( x(1:n) ) ) ) then
             if ( iprint .ge. 2 ) then
                write(*,*) 'Accelerated point is being accepted; and history will be updated substituting xtrial by xacc.'
             end if

             accacc = accacc + 1
             
             xhist(1:n,nh) = xacc(1:n)
             rhist(1:m,nh) = resacc(1:m)
             
             jtbr = ncol
             
             col(1:m) = rhist(1:m,nh) - rhist(1:m,nh-1)
             
             call qrupdate(m,ncol,rank,nhlim-1,jtbr,m,Q,nhlim-1,R,perm,col)
             maxrank = max( maxrank, rank )
             
             if ( iprint .ge. 2 ) then
                write(*,*) 'Number of columns = ',ncol,' rank = ',rank,' (largest achieved rank = ',maxrank,')'
             end if
             
             if ( ncol .ne. nh - 1 ) then
                write(*,*) 'ERROR in subroutine sesem: ncol is not equal to nh - 1! ncol = ',ncol,' nh = ',nh
                stop
             end if
             
             xtrial(1:n) = xacc(1:n)
             restrial(1:m) = resacc(1:m)
             ftrial = facc
          end if
       end if
    end if
    
    if ( rank .eq. 0 ) then
       if ( iprint .ge. 0 ) then
          write(*,*) 'Acceleration matrix has null rank; thus, it will be rebuild from scratch.'
          write(*,*) 'It will be populated with neighbours of the current point taking a small step into coordinate directions.'
       end if

       nh = 0
          
400    continue
       
       if ( nh + 1 .lt. min( nhlim, min( n + 1, m + 1 ) ) ) then
          xtmp(1:n) = xtrial(1:n)
          xtmp(nexti) = xtrial(nexti) + hlarge * max( 1.0d0, abs( xtrial(nexti) ) )
          
          nexti = mod( nexti, n ) + 1

          call evalr(n,xtmp,m,restmp,pdataptr)

          if ( iprint .ge. 2 ) then
             write(*,*) '(Considered ith canonical array with i = ',nexti,')'
             write(*,*) 'f (sum of squares) at considered neighbour = ',sum( restmp(1:m) ** 2 )
          end if
          
          nh = nh + 1
             
          xhist(1:n,nh) = xtmp(1:n)
          rhist(1:m,nh) = restmp(1:m)

          if ( nh .ge. 2 ) then
             col(1:m) = rhist(1:m,nh) - rhist(1:m,nh-1)

             if ( nh .eq. 2 ) then
                ncol = 1
                call qr(m,ncol,m,col,gam,perm,rank)
                call extractQ(m,ncol,m,col,gam,rank,m,Q)
                call extractR(ncol,m,col,rank,nhlim-1,R)
                maxrank = max( maxrank, rank )
                   
                if ( iprint .ge. 2 ) then
                   write(*,*) 'Number of columns = ',ncol,' rank = ',rank,' (largest achieved rank = ',maxrank,')'
                end if
             else
                jtbr = 0
                call qrupdate(m,ncol,rank,nhlim-1,jtbr,m,Q,nhlim-1,R,perm,col)
                maxrank = max( maxrank, rank )
          
                if ( iprint .ge. 2 ) then
                   write(*,*) 'Number of columns = ',ncol,' rank = ',rank,' (largest achieved rank = ',maxrank,')'
                end if
          
                if ( ncol .ne. nh - 1 ) then
                   write(*,*) 'ERROR in subroutine sesem: ncol is not equal to nh - 1! ncol = ',ncol,' nh = ',nh
                   stop
                end if
             end if
          end if
          
          go to 400
       end if

       if ( iprint .ge. 2 ) then
          write(*,*) 'Include xtrial.'
       end if
       
       if ( iprint .ge. 2 ) then
          write(*,*) 'f (sum of squares) at xtrial = ',sum( restrial(1:m) ** 2 )
       end if
          
       nh = nh + 1
          
       xhist(1:n,nh) = xtrial(1:n)
       rhist(1:m,nh) = restrial(1:m)

       jtbr = 0

       col(1:m) = rhist(1:m,nh) - rhist(1:m,nh-1)

       call qrupdate(m,ncol,rank,nhlim-1,jtbr,m,Q,nhlim-1,R,perm,col)
       maxrank = max( maxrank, rank )

       if ( iprint .ge. 2 ) then
          write(*,*) 'Number of columns = ',ncol,' rank = ',rank,' (largest achieved rank = ',maxrank,')'
       end if
          
       call qrsolve(m,ncol,rank,m,Q,nhlim-1,R,perm,restrial,xacc)
       
       xacc(1:n) = xtrial(1:n) - matmul( xhist(1:n,2:nh) - xhist(1:n,1:nh-1), xacc(1:nh-1) )

       call evalr(n,xacc,m,resacc,pdataptr)
       fcnt = fcnt + 1
       
       facc = dot_product( resacc(1:m), resacc(1:m) )

       if ( iprint .ge. 1 ) then
        ! write(*,*) '    facc = ',facc,' xacc = ',xacc(1:n),' resacc = ',resacc(1:m)
        ! write(*,*) '    facc = ',facc,' xacc = ',xacc(1:n)
          write(*,*) '    facc = ',facc
       end if

       if ( facc .lt. ftrial .and. norm2( xacc(1:n) - x(1:n) ) .gt. macheps * max( 1.0d0, norm2( x(1:n) ) ) ) then

          accacc = accacc + 1

          xhist(1:n,nh) = xacc(1:n)
          rhist(1:m,nh) = resacc(1:m)
       
          jtbr = ncol

          col(1:m) = rhist(1:m,nh) - rhist(1:m,nh-1)

          call qrupdate(m,ncol,rank,nhlim-1,jtbr,m,Q,nhlim-1,R,perm,col)
          maxrank = max( maxrank, rank )
                    
          xtrial(1:n) = xacc(1:n)
          restrial(1:m) = resacc(1:m)
          ftrial = facc
       end if
    end if

500 continue
    
    x(1:n) = xtrial(1:n)
    res(1:m) = restrial(1:m)
    f = ftrial

    if ( f .lt. fbest ) then
       nprogr = 0
       fbest = f
       fcntbest = fcnt
       iterbest = iter
    else
       nprogr = nprogr + 1
    end if
    
    go to 100

9000 format(1X,I7,1X,I7,1X,I1,1X,1P,D12.6,1X,1P,D12.6,2(1X,I10),1X,1P,D12.6,2(1X,I10))

  contains

    ! ******************************************************************
    ! ******************************************************************

    subroutine calfun(nred,xred,f)

      implicit none
      
      ! SCALAR ARGUMENTS
      integer, intent(in) :: nred
      real(kind=8), intent(out) :: f
      
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: xred(nred)

      ! LOCAL ARRAYS
      real(kind=8) :: d(n),rfull(m),xfull(n)

      if ( redscheme .eq. 1 ) then
         d(1:n) =  matmul(xt(1:n,1:nred),xred(1:nred))
      else if ( redscheme .eq. 2 ) then
         call frompvtox2(nred,xred,n,d)
      else
         call frompvtox3(nred,xred,n,d)
      end if
      
      if ( any( isnan( d(1:n) ) ) ) then
         write(*,*) 'There is a NaN in d!'
         stop
      end if
      
      xfull(1:n) = x(1:n) + d(1:n)
      xfull(1:n) = max( l(1:n), min( xfull(1:n), u(1:n) ) )
      
      if ( any( isnan( xfull(1:n) ) ) ) then
         write(*,*) 'There is a NaN in xfull!'
         stop
      end if
      
!!$      write(*,*) 'Computed point: '
!!$      write(*,*) 'x = ',xfull(1:n)
    
      call evalr(n,xfull,m,rfull,pdataptr)
      fcnt = fcnt + 1

      if ( any( isnan( rfull(1:m) ) ) ) then
         write(*,*) 'There is a NaN in rfull!'
         stop
      end if
      
      f = sum( rfull(rind(1:mred)) ** 2 )

!!$      write(*,*) 'Value of f = ',f
    
      if ( isnan( f ) ) then
         write(*,*) 'f is NaN !!!'
         stop
      end if
      
    end subroutine calfun
    
  end subroutine sesem
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine frompvtox2(nred,pv,n,d)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n,nred

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: pv(nred)
    real(kind=8), intent(out) :: d(n)
  
    ! LOCAL SCALARS
    integer :: i,j,k,left,npoints
    real(kind=8) :: ival,p
  
    ! LOCAL ARRAYS
    logical :: remove((nred-2)/2+2)
    real(kind=8) :: snodes((nred-2)/2+2),values((nred-2)/2+2)

    k = ( nred - 2 ) / 2
    
    if ( any( pv(1:k) .lt. 0.0d0 .or. pv(1:k) .gt. 1.0d0 ) ) then
       write(*,*) 'In frompvtox2, there are points out of range.'
       stop
    end if

    ! Extract points and values

    snodes(1:k+2) = (/ 0.0d0, pv(1:k), 1.0d0 /)
    values(1:k+2) = pv(k+1:2*k+2)

!!$    write(*,*) 'In frompvtox2, k = ',k
!!$    write(*,*) 'snodes = ',snodes(1:k+2)
!!$    write(*,*) 'values = ',values(1:k+2)

    ! Sort (both) in non-decreasing order
  
    call dsort(snodes,values,k+2,2)

!!$    write(*,*) 'After sorting.'
!!$    write(*,*) 'snodes = ',snodes(1:k+2)
!!$    write(*,*) 'values = ',values(1:k+2)

    ! Shrink coincident points
  
    i = 1

10  continue
  
    remove(i) = .false.

    ival = values(i)

    j = i

20  continue
  
    if ( j + 1 .le. k + 2 ) then
       if ( snodes(j+1) .eq. snodes(i) ) then
          ival = ival + values(j+1)
          remove(j+1) = .true.
          j = j + 1
          go to 20
       end if
    end if

    values(i) = ival / ( j + 1 - i )
  
    i = j + 1

    if ( i .le. k + 2 ) then
       go to 10
    end if

    npoints = count( .not. remove(1:k+2) )

    if ( npoints .gt. k + 2 ) then
       write(*,*) 'There is something wrong in frompvtox2 because, ', &
            'after removing coincident point, npoints is greater than k+2.'
       stop
    end if
     
    snodes(1:npoints) = pack( snodes(1:k+2), .not. remove(1:k+2) )
    values(1:npoints) = pack( values(1:k+2), .not. remove(1:k+2) )
    
    if ( npoints .lt. 2 ) then
       write(*,*) 'There is something wrong in frompvtox2 because, ', &
            'after removing coincident point, npoints is smaller than 2.'
       stop
    end if
    
!!$    write(*,*) 'After shrinking, npoints = ',npoints
!!$    write(*,*) 'snodes = ',snodes(1:npoints)
!!$    write(*,*) 'values = ',values(1:npoints)
    
    ! Compute d using linear interpolation

    d(1) = values(1)
    d(n) = values(npoints)
  
    left = 1
    do i = 2,n - 1
       p = real(i-1) / (n-1)
30     continue
       if ( left .eq. npoints ) then
          write(*,*) 'In frompvtox2, there is something wrong.'
          stop
       end if
       if ( snodes(left+1) .le. p ) then
          left = left + 1
          go to 30
       end if
       d(i) = values(left) + ( values(left+1) - values(left) ) * ( p - snodes(left) ) / ( snodes(left+1) - snodes(left) )
    end do

!!$    write(*,*) 'Computed direction: '
!!$    write(*,*) 'direction = ',d(1:n)
    
  end subroutine frompvtox2

  ! ******************************************************************
  ! ******************************************************************

  subroutine frompvtox3(nred,values,n,d)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nred,n

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: values(nred)
    real(kind=8), intent(out) :: d(n)
  
    ! LOCAL SCALARS
    integer :: i,left
    real(kind=8) :: p
  
    ! LOCAL ARRAYS
    real(kind=8) :: snodes(nred)

    ! Compute d using linear interpolation

    snodes(1:nred) = (/ (real(i-1)/(nred-1),i=1,nred) /)

    d(1) = values(1)
    d(n) = values(nred)
  
    left = 1
    do i = 2,n - 1
       p = real(i-1) / (n-1)
10     continue
       if ( left .eq. nred ) then
          write(*,*) 'In frompvtox3, there is something wrong.'
          stop
       end if
       if ( snodes(left+1) .le. p ) then
          left = left + 1
          go to 10
       end if
       d(i) = values(left) + ( values(left+1) - values(left) ) * ( p - snodes(left) ) / ( snodes(left+1) - snodes(left) )
    end do

  end subroutine frompvtox3

end module bmsesem