module modminres  

    use iso_c_binding, only: c_ptr

    implicit none

    type :: matrix_type
        integer :: ne,n  
        integer, allocatable, dimension(:) :: row, col
        real(kind=8), allocatable, dimension(:) :: val
    end type

    type lss_type
        integer :: nnzmax
        type(matrix_type) :: matrix
    end type lss_type
  
    real(kind=8), parameter :: macheps = epsilon( 1.0d0 ),macheps12 = sqrt(macheps)

    public minres  

  contains
    subroutine minres (n,nfr,b,eta,y,r,k,kmax,iprint,Dtype,lss)
        implicit none

        !SCALARS ARGUMENTS
        integer, intent(in) :: n,nfr,kmax,iprint
        integer, intent(out):: Dtype,k
        real(kind=8), intent(in) :: eta
    
        !ARRAY ARGUMENTS
        real(kind=8), intent(in)  :: b(n)
        real(kind=8), intent(out) :: y(n),r(n)
        type(lss_type), intent(in) :: lss

        !LOCAL SCALARS
        real(kind=8) :: normb,beta,c,s,phi,tau,delta1,delta2,gamma1,gamma2,alpha,epscurr,epsnext
    
        !LOCAL ARRAYS
        integer :: i,j,z
        real(kind=8) :: vcurr(n),vprev(nfr),dcurr(n),dprev(n),dprev2(n), &
            q(n),p(nfr),yprev(n),rprev(n),vcurrfull(n)

        !INICIALIZATION

        normb = norm2( b(1:nfr) )
        beta  = normb 
        phi   = normb

        r(1:nfr) = b(1:nfr)
        vcurr(1:nfr) = b(1:nfr)/normb 
        vcurr(nfr+1:n) = 0.0d0 

        vprev(1:nfr)   = 0.0d0
        y(1:n)         = 0.0d0
        dprev(1:n)     = 0.0d0 
        dprev2(1:n)    = 0.0d0 
        vcurrfull(1:n) = 0.0d0 

        epscurr = 0.0d0
        c = -1.0d0
        s = 0.0d0

        tau = 0.0d0
        delta1 = 0.0d0

        k = 0
    100 continue 
        if (iprint .GE. 2) then 
            write(*,*)
            write(*,*) "MINRES: iter = ", k, " norm2(r) = ", norm2(r(1:nfr)) 
        end if 

        q(1:n) = 0.0d0
        do z = 1,lss%matrix%ne
            i = lss%matrix%row(z)
            j = lss%matrix%col(z)
            if (j .ne. i) then
                q(i) = q(i) + lss%matrix%val(z) * vcurr(j)
                q(j) = q(j) + lss%matrix%val(z) * vcurr(i)
            else
                q(i) = q(i) + lss%matrix%val(z) * vcurr(j)
            end if 
        end do   

        alpha = dot_product(vcurr(1:nfr),q(1:nfr))

        p(1:nfr) = q(1:nfr) - beta*vprev(1:nfr) - alpha*vcurr(1:nfr)
        beta = norm2(p(1:nfr))

        delta2     =  c*delta1 + s*alpha
        gamma1     =  s*delta1 - c*alpha
        epsnext    =  s*beta
        delta1     = -c*beta

        yprev(1:nfr) = y(1:nfr)
        rprev(1:nfr) = r(1:nfr)
        if (c*gamma1 .GE. 0) then 
            if (iprint .GE. 1) then 
                write(*,*)
                write(*,*) "MINRES: NPC directions returned. Final Iter = ", k
            end if 
            Dtype = 1
            return
        end if 

        if (phi * sqrt(gamma1*gamma1 + delta1*delta1) .LE. eta * sqrt(normb*normb - phi*phi)) then
            if (iprint .GE. 1) then 
                write(*,*)
                write(*,*) "MINRES: Solution found. Final Iter = ", k
            end if 
            Dtype = 0
            return
        end if 

        if (k .GE. kmax) then 
            if (iprint .GE. 1) then 
                write(*,*)
                write(*,*) "MINRES: Maximum iteration reached."
            end if 
            Dtype = 2
            return
        end if 

        gamma2 = sqrt(gamma1*gamma1 + beta*beta)
        if (gamma2 .GT. 0.0d0) then 
            c = gamma1 / gamma2
            s = beta   / gamma2
            tau = c*phi
            phi = s*phi

            dcurr(1:nfr) = (vcurr(1:nfr) - delta2*dprev(1:nfr) - epscurr*dprev2(1:nfr))/gamma2
            y(1:nfr) = yprev(1:nfr) + tau*dcurr(1:nfr)

            if (beta .GT. 0.0d0) then 
                vprev(1:nfr) = vcurr(1:nfr)
                vcurr(1:nfr) = p(1:nfr)/beta

                r(1:nfr) = s*s*r(1:nfr) - phi*c*vcurr(1:nfr)
            end if 
        else 
            c = 0.0d0 
            s = 1.0d0 
            tau = 0.0d0 
            y(1:nfr) = yprev(1:nfr)
        end if 
        dprev2(1:nfr) = dprev(1:nfr)
        dprev(1:nfr)  = dcurr(1:nfr)
        epscurr = epsnext

        k = k + 1 
        go to 100

    end subroutine minres

  end module modminres
  