MODULE util_minim
        USE nrtype
        USE error
        IMPLICIT NONE

        INTERFACE assert_eq
                MODULE PROCEDURE assert_eq2,assert_eq3,assert_eqn
        END INTERFACE
        INTERFACE diagmult
                MODULE PROCEDURE diagmult_rv,diagmult_r
        END INTERFACE
!        INTERFACE outerprod
!                MODULE PROCEDURE outerprod_r,outerprod_d
!        END INTERFACE
        INTERFACE imaxloc
                MODULE PROCEDURE imaxloc_r,imaxloc_i
        END INTERFACE
        INTERFACE swap
                MODULE PROCEDURE swap_r,swap_rv
        END INTERFACE        
        CONTAINS

!BL
!        SUBROUTINE assert_v(n,string)
!        CHARACTER(LEN=*), INTENT(IN) :: string
!        LOGICAL, DIMENSION(:), INTENT(IN) :: n
!        if (.not. all(n)) then
!                write (*,*) 'nrerror: an assertion failed with this tag:', &
!                        string
!                STOP 'program terminated by assert_v'
!        end if
!        END SUBROUTINE assert_v
!BL
        FUNCTION assert_eq2(n1,n2,string)
        CHARACTER(LEN=*), INTENT(IN) :: string
        INTEGER, INTENT(IN) :: n1,n2
        INTEGER :: assert_eq2
        if (n1 == n2) then
                assert_eq2=n1
        else
                write (*,*) 'nrerror: an assert_eq failed with this tag:', &
                        string
!                STOP 'program terminated by assert_eq2'
                CALL error_msg('program terminated by assert_eq2')
        end if
        END FUNCTION assert_eq2
!BL
        FUNCTION assert_eq3(n1,n2,n3,string)
        CHARACTER(LEN=*), INTENT(IN) :: string
        INTEGER, INTENT(IN) :: n1,n2,n3
        INTEGER :: assert_eq3
        if (n1 == n2 .and. n2 == n3) then
                assert_eq3=n1
        else
                write (*,*) 'nrerror: an assert_eq failed with this tag:', &
                        string
!                STOP 'program terminated by assert_eq3'
                CALL error_msg('program terminated by assert_eq3')
        end if
        END FUNCTION assert_eq3
!BL
        FUNCTION assert_eqn(nn,string)
        CHARACTER(LEN=*), INTENT(IN) :: string
        INTEGER, DIMENSION(:), INTENT(IN) :: nn
        INTEGER :: assert_eqn
        if (all(nn(2:) == nn(1))) then
                assert_eqn=nn(1)
        else
                write (*,*) 'nrerror: an assert_eq failed with this tag:', &
                        string
!                STOP 'program terminated by assert_eqn'
                CALL error_msg('program terminated by assert_eqn')
        end if
        END FUNCTION assert_eqn
!BL
!        SUBROUTINE nrerror(string)
!        CHARACTER(LEN=*), INTENT(IN) :: string
!        write (*,*) 'nrerror: ',string
!        STOP 'program terminated by nrerror'
!        END SUBROUTINE nrerror
!BL
        SUBROUTINE swap_r(a,b)
        REAL(DP), INTENT(INOUT) :: a,b
        REAL(DP) :: dum
        dum=a
        a=b
        b=dum
        END SUBROUTINE swap_r
!BL
        SUBROUTINE swap_rv(a,b)
        REAL(DP), DIMENSION(:), INTENT(INOUT) :: a,b
        REAL(DP), DIMENSION(SIZE(a)) :: dum
        dum=a
        a=b
        b=dum
        END SUBROUTINE swap_rv
!BL
        SUBROUTINE diagmult_rv(mat,diag)
        REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: mat
        REAL(DP), DIMENSION(:), INTENT(IN) :: diag
        INTEGER(I4B) :: j,n
        n = assert_eq2(size(diag),min(size(mat,1),size(mat,2)),'diagmult_rv')
        do j=1,n
                mat(j,j)=mat(j,j)*diag(j)
        end do
        END SUBROUTINE diagmult_rv
!BL
        SUBROUTINE diagmult_r(mat,diag)
        REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: mat
        REAL(DP), INTENT(IN) :: diag
        INTEGER(I4B) :: j,n
        n = min(size(mat,1),size(mat,2))
        do j=1,n
                mat(j,j)=mat(j,j)*diag
        end do
        END SUBROUTINE diagmult_r
!BL
        FUNCTION outerand(a,b)
        LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: a,b
        LOGICAL(LGT), DIMENSION(size(a),size(b)) :: outerand
        outerand = spread(a,dim=2,ncopies=size(b)) .and. &
                spread(b,dim=1,ncopies=size(a))
        END FUNCTION outerand
!BL
        FUNCTION outerprod_r(a,b)
        REAL(DP), DIMENSION(:), INTENT(IN) :: a,b
        REAL(DP), DIMENSION(size(a),size(b)) :: outerprod_r
        outerprod_r = spread(a,dim=2,ncopies=size(b)) * &
                spread(b,dim=1,ncopies=size(a))
        END FUNCTION outerprod_r
!BL
!        FUNCTION outerprod_d(a,b)
!        REAL(DP), DIMENSION(:), INTENT(IN) :: a,b
!        REAL(DP), DIMENSION(size(a),size(b)) :: outerprod_d
!        outerprod_d = spread(a,dim=2,ncopies=size(b)) * &
!                spread(b,dim=1,ncopies=size(a))
!        END FUNCTION outerprod_d
!BL
        FUNCTION imaxloc_r(arr)
        REAL(DP), DIMENSION(:), INTENT(IN) :: arr
        INTEGER(I4B) :: imaxloc_r
        INTEGER(I4B), DIMENSION(1) :: imax
        imax=maxloc(arr(:))
        imaxloc_r=imax(1)
        END FUNCTION imaxloc_r
!BL
        FUNCTION imaxloc_i(iarr)
        INTEGER(I4B), DIMENSION(:), INTENT(IN) :: iarr
        INTEGER(I4B), DIMENSION(1) :: imax
        INTEGER(I4B) :: imaxloc_i
        imax=maxloc(iarr(:))
        imaxloc_i=imax(1)
        END FUNCTION imaxloc_i
!BL
        FUNCTION iminloc(arr)
        REAL(DP), DIMENSION(:), INTENT(IN) :: arr
        INTEGER(I4B), DIMENSION(1) :: imin
        INTEGER(I4B) :: iminloc
        imin=minloc(arr(:))
        iminloc=imin(1)
        END FUNCTION iminloc
!BL
        SUBROUTINE covsrt(covar,maska)
        USE nrtype
!        USE util_minim, ONLY : assert_eq,swap_rv
        IMPLICIT NONE
        REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: covar
        LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: maska
        INTEGER(I4B) :: ma,mfit,j,k
        ma=assert_eq(size(covar,1),size(covar,2),size(maska),'covsrt')
        mfit=count(maska)
        covar(mfit+1:ma,1:ma)=0.0
        covar(1:ma,mfit+1:ma)=0.0
        k=mfit
        do j=ma,1,-1
                if (maska(j)) then
                        call swap_rv(covar(1:ma,k),covar(1:ma,j))
                        call swap_rv(covar(k,1:ma),covar(j,1:ma))
                        k=k-1
                end if
        end do
        END SUBROUTINE covsrt
!BL
        SUBROUTINE gaussj(a,b)
        USE nrtype 
!        USE util_minim, ONLY : assert_eq,nrerror,outerand,outerprod,swap_rv
        USE error
        IMPLICIT NONE
        REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: a,b
        INTEGER(I4B), DIMENSION(size(a,1)) :: ipiv,indxr,indxc
        LOGICAL(LGT), DIMENSION(size(a,1)) :: lpiv
        REAL(DP) :: pivinv
        REAL(DP), DIMENSION(size(a,1)) :: dumc
        INTEGER(I4B), TARGET :: irc(2)
        INTEGER(I4B) :: i,l,n
        INTEGER(I4B), POINTER :: irow,icol
        n=assert_eq(size(a,1),size(a,2),size(b,1),'gaussj')
        irow => irc(1)
        icol => irc(2)
        ipiv=0
        do i=1,n
                lpiv = (ipiv == 0)
                irc=maxloc(abs(a),outerand(lpiv,lpiv))
                ipiv(icol)=ipiv(icol)+1
!                if (ipiv(icol) > 1) call nrerror('gaussj: singular matrix (1)')
                if (ipiv(icol) > 1) CALL error_msg('gaussj: singular matrix (1)')
                if (irow /= icol) then
                        call swap_rv(a(irow,:),a(icol,:))
                        call swap_rv(b(irow,:),b(icol,:))
                end if
                indxr(i)=irow
                indxc(i)=icol
                !added by corrado to avoid singular matrix when the
                !derivatives are zero
                if (a(icol,icol) == 0.0) a(icol,icol)=1.0
!               !#######################################
!                if (a(icol,icol) == 0.0) &
!                        call nrerror('gaussj: singular matrix (2)')
!                        CALL error_msg('gaussj: singular matrix (2)')
                pivinv=1.0_dp/a(icol,icol)
                a(icol,icol)=1.0
                a(icol,:)=a(icol,:)*pivinv
                b(icol,:)=b(icol,:)*pivinv
                dumc=a(:,icol)
                a(:,icol)=0.0
                a(icol,icol)=pivinv
                a(1:icol-1,:)=a(1:icol-1,:)-outerprod_r(dumc(1:icol-1),a(icol,:))
                b(1:icol-1,:)=b(1:icol-1,:)-outerprod_r(dumc(1:icol-1),b(icol,:))
                a(icol+1:,:)=a(icol+1:,:)-outerprod_r(dumc(icol+1:),a(icol,:))
                b(icol+1:,:)=b(icol+1:,:)-outerprod_r(dumc(icol+1:),b(icol,:))
        end do
        do l=n,1,-1
                call swap_rv(a(:,indxr(l)),a(:,indxc(l)))
        end do
        END SUBROUTINE gaussj

END MODULE util_minim