module brown_m
!! Contains routines for evaluating **B**._d_**W** in the SDE for Brownian Dynamics
!! simulation.

use iso_c_binding, only: c_loc, c_ptr, c_f_pointer
use mkl_blas
use mkl_lapack
use constants_m
use random_m
use strings_m
use logger_m

implicit none

private
public :: brn_init, brn_finish, brn_calc_dw, brn_calc_bdw

character(4) :: meth
    !!Method to use for generating Brownian terms. Possible values are
    !!'CHOL' or 'LANC'.
integer :: sizm          
    !!Size of the mobility matrix
integer :: s             
    !!Number of Brownian vectors to be generated; equivalent to block size in the
    !!block Lanczos algorithm.
integer :: dimk_min = 2
    !!Minimum dimension of KSP
integer :: dimk_max      
    !!Maximum dimension of KSP
real(rp) :: ethres = 1.0e-3_rp  
    !!Error threshold for KSP-based method.
real(rp), dimension(:,:), allocatable :: bdw_old
    !!(`sizm`, `s`) array. Stores B.dw at previous iteration.
real(rp), dimension(:,:), allocatable :: lncv
    !!(`sizm`, `s*dimk_max`) array. Stores the Lanczos vectors.
real(rp), dimension(:), allocatable :: eigvlh
    !!(`s*dimk_max`,) array. Eigenvalues of the tridiagonal/block triadiagonal
    !! matrix H generated by the Lanczos algorithm.
real(rp), dimension(:,:), allocatable :: eigvch
    !!(`s*dimk_max`, `s*dimk_max`) array. Eigenvectors of the matrix H.
real(rp), dimension(:,:), allocatable :: w
    !!(`sizm`, `s`) array. Used for generating the Lanczos vectors.
real(rp), dimension(:), allocatable :: h_d
    !!(`dimk_max`,) array. Diagonal of the symmetric tridiagonal matrix H
    !! generated in case of `s = 1`.
real(rp), dimension(:), allocatable :: h_sd
    !!(`dimk_max`,) array. Subdiagonal of the symmetric tridiagonal matrix H
    !! generated in case of `s = 1`.
real(rp), dimension(:), allocatable :: tau
    !! (`s`,) array. Used in the QR decomposition of `w`.
real(rp), dimension(:,:), allocatable :: matr
    !! (`s`,`s`) array. Used in the QR decomposition of `w`.
real(rp), dimension(:,:), allocatable :: h_bd
    !!(`s`, `s*dimk_max`) array. Diagonal blocks of the matrix H.
    !!Used for KSP with `s > 1`.
real(rp), dimension(:,:), allocatable :: h_bsd
    !!(`s`, `s*dimk_max`) array. Subdiagonal blocks of the matrix H.
    !!Used for KSP with `s > 1`.
real(rp), dimension(:,:), allocatable :: h_bnd
    !!(`2*s`, `s*dimk_max`) array. Matrix H in banded storage.
    !!Used for KSP with `s > 1`.
real(rp), dimension(:,:), allocatable :: tmp
    !!(`s*dimk_max`, `2*s`) array used as temporary storage space.
real(rp), dimension(:), allocatable :: work
    !!Workspace array for LAPACK routines.
integer, dimension(:), allocatable :: iwork
    !!Workspace array for LAPACK routines.

contains

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

subroutine brn_init(n, nblks, mth, f, e)
    !! Performs initial setup.

    integer, intent(in) :: n
        !! Size of the mobility matrix.
    integer, intent(in) :: nblks
        !! Number of **B**._d_**W** vectors to generate for the same mobility matrix.
        !! `1 <= nblks < n`.
    character(4), intent(in) :: mth
        !! Method for calculating the Brownian terms. `mth = 'CHOL'` for
        !! Cholesky decomposition; `mth = 'LANC'` for KSP-based method.
    integer, intent(in), optional :: f
        !! For KSP-based method, maximum number of iterations.
        !! Must be present if `mth = 'LANC'`. `f` must be less than `n`.
    real(rp), intent(in), optional :: e
        !! For KSP-based method, error threshold for convergence. 
        !! Must be present if `mth = 'LANC'`.

    sizm = n; s = nblks; meth = mth

    if (meth == 'CHOL') then
        !Nothing to do here.
    else if (meth == 'LANC') then
        dimk_max = max(dimk_min, f) + 1
        !Adding 1 to ensure dimk_max > dimk_min
        ethres = e

        allocate( bdw_old(sizm, s) )
        allocate( lncv(sizm,s*dimk_max) )
        allocate( eigvlh(s*dimk_max) )
        allocate( eigvch(s*dimk_max,s*dimk_max) )
        allocate( w(sizm,s))
        allocate( tmp(s*dimk_max,2*s))
        if (s == 1) then
            allocate( h_d(dimk_max) )
            allocate( h_sd(dimk_max) )
        else
            allocate( tau(s))
            allocate( matr(s,s))
            allocate( h_bd(s,s*dimk_max) )
            allocate( h_bsd(s,s*dimk_max) )
            allocate( h_bnd(2*s,s*dimk_max))
        end if
    end if

    end subroutine

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

subroutine brn_finish()
    !! Performs cleanup.

    if (allocated(bdw_old)) deallocate(bdw_old)
    if (allocated(lncv))    deallocate(lncv)
    if (allocated(eigvlh))  deallocate(eigvlh)
    if (allocated(eigvch))  deallocate(eigvch)
    if (allocated(w))       deallocate(w)
    if (allocated(tmp))     deallocate(tmp)

    if (allocated(h_d))     deallocate(h_d)
    if (allocated(h_sd))    deallocate(h_sd)
    if (allocated(tau))     deallocate(tau)
    if (allocated(matr))    deallocate(matr)
    if (allocated(h_bd))    deallocate(h_bd)
    if (allocated(h_bsd))   deallocate(h_bsd)
    if (allocated(h_bnd))   deallocate(h_bnd)

    if (allocated(work))    deallocate(work)
    if (allocated(iwork))   deallocate(iwork)

    end subroutine

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

subroutine brn_calc_bdw(mob, bdw, ierr, lconv, f)
    !! Main driver routine for calculating **B**._d_**W**.

    real(rp), dimension(:,:), intent(in out) :: mob
        !! (`sizm`,`sizm`) symmetric positive definite matrix.
        !! On _entry_, contains the mobility matrix.  On _return_, if `meth = 'CHOL'`,
        !! the upper triangular part is overwritten with the result of Cholesky
        !! decomposition.  The strictly lower triangular part is left intact.
    real(rp), dimension(:,:), intent(in out) :: bdw
        !! (`sizm`,`s`) matrix.
        !! On _entry_, contains `s` column vectors drawn from standard normal 
        !! distribution. On _return_, the `s` columns are overwritten with the 
        !! result of **B**.d**W**.
    integer, intent(out) :: ierr
        !! Error flag.
    logical, intent(out), optional :: lconv
        !! If `meth = 'LANC'`, returns *true* if converged, *false* otherwise.
    integer, intent(out), optional :: f
        !! If `meth = 'LANC'`, the number of iterations performed.
    integer :: f_
    logical :: lconv_

    if (meth == 'CHOL') then
        call calc_bdw_cholesky(mob, bdw, ierr)
    else if (meth == 'LANC') then
        if (s == 1) then
            call calc_bdw_lanc(mob, bdw, lconv_, f_, ierr)
        else
            call calc_bdw_blanc(mob, bdw, lconv_, f_, ierr)
        end if
    end if
    if (present(lconv)) lconv = lconv_
    if (present(f)) f = f_

    end subroutine

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

subroutine calc_bdw_lanc(mob, bdw, lconv, f, ierr)
 !! Implements Lanczos algorithm for calculating **B**._d_**W**.

    real(rp), dimension(:,:), intent(in) :: mob
        !! (`sizm`,`sizm`) symmetric positive definite matrix.
    real(rp), dimension(:,:), intent(in out) :: bdw
        !! (`sizm`,1) array.
        !! On _entry_, contains a vector drawn from standard normal 
        !! distribution. On _return_, is overwritten with the result of
        !! **B**._d_**W**.
    logical, intent(out) :: lconv
        !! Returns *true* if converged, *false* otherwise.
    integer, intent(out) :: f
        !! The number of iterations performed.
    integer, intent(out) :: ierr
        !! Error flag.
    integer :: j, info, lwork, liwork
    real(rp), dimension(1) :: qwork  !Size 1 array for workspace query.
    integer, dimension(1)  :: qiwork !Size 1 array for workspace query.
    real(rp) :: normz, err
    character(len=256) :: msg

    lconv = .false.; ierr = 0; f = 0
    lncv = 0.0_rp; h_d = 0.0_rp; h_sd = 0.0_rp
    eigvlh = 0.0_rp; eigvch = 0.0_rp
    w = 0.0_rp; tmp = 0.0_rp; bdw_old = 0.0_rp

    !Workspace query
    if ((.not. allocated(work)) .or. (.not. allocated(iwork))) then
        call dstevd('V', dimk_max, eigvlh, tmp(1:dimk_max,1), eigvch, dimk_max, &
                    qwork, -1, qiwork, -1, info)
        if (info /= 0) then
            write(msg,'(a,i0,a)') '<calc_bdw_lanc> dstevd err = ', info, &
                ' on workspace query'
            call logger%log_msg(msg)
            ierr = 1; return
        end if
        lwork = int(qwork(1)); liwork = qiwork(1)
        allocate( work(lwork) ); allocate( iwork(liwork) )
    end if

    !Lanczos process
    normz = norm2(bdw(:,1))
    lncv(:,1) = bdw(:,1)/normz

    do j = 1, dimk_max
        f = f + 1 !Iteration count
        call dsymv('L', sizm, 1.0_rp, mob, sizm, lncv(:,j), 1, 0.0_rp, w(:,1), 1)
        if (j > 1) then
            w(:,1) = w(:,1) - h_sd(j-1)*lncv(:,j-1)
        end if
        h_d(j) = dot_product(w(:,1), lncv(:,j))

        if (j < dimk_max) then
            w(:,1) = w(:,1) - h_d(j)*lncv(:,j)
            h_sd(j) = norm2(w(:,1))
            lncv(:,j+1) = w(:,1)/h_sd(j)
        end if

        !Calculate B.dW and check for convergence
        if (j >= dimk_min) then
            eigvlh(1:j) = h_d(1:j)
            tmp(1:j,1) = h_sd(1:j)

            !Eigen decomposition of H = P * d * P^T
            !eigvlh[1:j] <- d; eigvch[1:j,1:j] <- P
            call dstevd('V', j, eigvlh(1:j), tmp(1:j,1), eigvch, dimk_max, &
                    work, size(work), iwork, size(iwork), info)
            if (info /= 0) then
                write(msg,'(a,i0,a,i0)') '<calc_bdw_lanc> dstevd err = ', info, &
                    ' on iteration ', j
                call logger%log_msg(msg)
                ierr = 1; return
            end if

            !eigvlh(1:j) <- d^(1/2)
            eigvlh(1:j) = sqrt(eigvlh(1:j))

            !tmp[1:j,1] <- d^(1/2) * P^T * e_1, where e_1 if the first column of
            !an identity matrix of size j.
            tmp(1:j,1) = eigvch(1,1:j)
            tmp(1:j,1) = eigvlh(1:j)*tmp(1:j,1)

            !tmp[1:j,2] <- P * tmp[1:j,1]
            call dgemv('N', j, j, 1.0_rp, eigvch, dimk_max, tmp(1:j,1), 1, &
                        0.0_rp, tmp(1:j,2), 1)

            !bdw <- normz * lncv[:,1:j] * tmp[1:j,2]
            call dgemv('N', sizm, j, normz, lncv, sizm, tmp(1:j,2), 1, &
                        0.0_rp, bdw, 1)

            !Error calculation for j > dimk_min. For j = dimk_min just save
            !bdw for error calculation in the next iteration.
            if (j > dimk_min) then
                err = norm2(bdw(:,1) - bdw_old(:,1))/norm2(bdw_old(:,1))
                if (err <= ethres) then
                    lconv = .true.; exit
                end if
            end if
            bdw_old = bdw
        end if
    end do

    end subroutine

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

subroutine calc_bdw_blanc(mob, bdw, lconv, f, ierr)
    !! Implements block Lanczos algorithm for calculating **B**._d_**W**.

    real(rp), dimension(:,:), intent(in) :: mob
        !! (`sizm`,`sizm`) symmetric positive definite matrix.
    real(rp), dimension(:,:), intent(in out) :: bdw
        !! (`sizm`,`s`) matrix.
        !! On _entry_, contains `s` column vectors drawn from standard normal 
        !! distribution. On _return_, the `s` columns are overwritten with the 
        !! result of **B**._d_**W**.
    logical, intent(out) :: lconv
        !! Returns *true* if converged, *false* otherwise.
    integer, intent(out) :: f
        !! The number of iterations performed.
    integer, intent(out) :: ierr
        !! Error flag.
    real(rp) :: err
    integer :: j, jb, irow, jcol
    integer :: jgbeg, jgend, igbeg, ig, jg
    integer :: info, lwork, liwork 
    real(rp), dimension(1) :: qwork
    integer, dimension(1)  :: qiwork
    character(len=256) :: msg

    lconv = .false.; ierr = 0; f = 0
    lncv = 0.0_rp;  h_bd = 0.0_rp; h_bsd = 0.0_rp
    eigvlh = 0.0_rp; eigvch = 0.0_rp
    w = 0.0_rp; tmp = 0.0_rp; bdw_old = 0.0_rp
    tau = 0.0_rp; matr = 0.0_rp; h_bnd = 0.0_rp

    !Workspace query
    if (.not. allocated(work)) then
        call dgeqp3(sizm, s, bdw, sizm, qiwork, tau, qwork, -1, info)
        if (info /= 0) then
            write(msg,'(a,i0,a)') '<calc_bdw_blanc> dgeqp3 err = ', info, &
                ' on workspace query'
            call logger%log_msg(msg)
            ierr = 1; return
        end if
        lwork = int( qwork(1) )

        call dorgqr(sizm, s, s, bdw, sizm, tau, qwork, -1, info)
        if (info /= 0) then
            write(msg,'(a,i0,a)') '<calc_bdw_blanc> dorgqr err = ', info, &
                ' on workspace query'
            call logger%log_msg(msg)
            ierr = 1; return
        end if
        lwork = max( lwork, int(qwork(1)) )

        call dsbevd('V', 'L', dimk_max*s, 2*s-1, h_bnd, 2*s, &
            eigvlh, eigvch, s*dimk_max, qwork, -1, qiwork, -1, info)
        if (info /= 0) then
            write(msg,'(a,i0,a)') '<calc_bdw_blanc> dsbevd err = ', info, &
                ' on workspace query'
            call logger%log_msg(msg)
            ierr = 1; return
        end if
        lwork = max( lwork, int(qwork(1)) )
        liwork = qiwork(1)

        allocate(work(lwork)); allocate(iwork(liwork))
    end if

    !Lanczos process
    !Reduced QR factorization of bdw with pivoting. The results overwrite bdw.
    iwork = 0; qwork = 0.0_rp
    call dgeqp3(sizm, s, bdw, sizm, iwork, tau, work, size(work), info)
    if (info /= 0) then
       write(msg,'(a,i0,a)') '<calc_bdw_blanc> dgeqp3 err = ', info, &
           ' on iteration 0'
       call logger%log_msg(msg)
       ierr = 1; return
    end if

    !Extract the (s,s) permuted upper triangular matrix R and store it in matr.
    do jcol = 1, s
        matr(1:jcol, iwork(jcol)) = bdw(1:jcol, jcol)
    end do

    !Extract the (sizm,s) matrix Q and store it in lncv[:,1:s].
    call dorgqr(sizm, s, s, bdw, sizm, tau, work, size(work), info)
    if (info /= 0) then
        write(msg,'(a,i0,a)') '<calc_bdw_blanc> dorgqr err = ', info, &
            ' on iteration 0'
        call logger%log_msg(msg)
        ierr = 1; return
    end if
    lncv(:,1:s) = bdw

    do j = 1, dimk_max
        f = f + 1 !Iteration count
        !W <- mob * V_j
        jgbeg = (j-1)*s+1; jgend = j*s
        call dsymm('L', 'L', sizm, s, 1.0_rp, mob, sizm, &
                    lncv(:,jgbeg:jgend), sizm, 0.0_rp, w, sizm )

        if (j > 1) then
            !W -> W - V_(j-1) * H_(j-1,j)
            !   = W - V_(j-1) * H_(j,j-1)^T, using the subdiagonal block
            jgbeg = (j-2)*s+1; jgend = (j-1)*s
            call dgemm('N', 'T', sizm, s, s, -1.0_rp, lncv(:,jgbeg:jgend), sizm, &
                        h_bsd(:,jgbeg:jgend), s, 1.0_rp, w, sizm)
        end if

        !H_(j,j) <- V_j^T * W 
        jgbeg = (j-1)*s+1; jgend = j*s
        call dgemm('T', 'N', s, s, sizm, 1.0_rp, lncv(:,jgbeg:jgend), sizm, &
                    w, sizm, 0.0_rp, h_bd(:,jgbeg:jgend), s)

        if (j < dimk_max) then
            !W -> W - V_(j) * H_(j,j)
            !Note: H_(j,j) is symmetric.
            jgbeg = (j-1)*s+1; jgend = j*s
            call dsymm('R', 'L', sizm, s, -1.0_rp, h_bd(:,jgbeg:jgend), &
                        s, lncv(:,jgbeg:jgend), sizm, 1.0_rp, w, sizm)

            !Reduced QR factorization of w.
            iwork = 0
            call dgeqp3(sizm, s, w, sizm, iwork, tau, work, size(work), info)
            if (info /= 0) then
                write(msg,'(a,i0,a,i0)') '<calc_bdw_blanc> dgeqp3 err = ', info, &
                    ' on iteration ', j
                call logger%log_msg(msg)
                ierr = 1; return
            end if

            !Extract the (s,s) permuted upper triangular matrix R and store
            !it in H_(j+1,j).
            do jcol = 1, s
                h_bsd(1:jcol, (j-1)*s+iwork(jcol)) = w(1:jcol, jcol)
            end do

            !Extract the (sizm,s) matrix Q and store it in V_(j+1).
            call dorgqr(sizm, s, s, w, sizm, tau, work, size(work), info)
            if (info /= 0) then
                write(msg,'(a,i0,a,i0)') '<calc_bdw_blanc> dorgqr err = ', info, &
                    ' on iteration ', j
                call logger%log_msg(msg)
                ierr = 1; return
            end if
            jgbeg = j*s+1; jgend = (j+1)*s
            lncv(:,jgbeg:jgend) = w
        end if

        !Calculate B.dW and check for convergence
        if (j >= dimk_min) then
            jgbeg = (j-1)*s+1; jgend = j*s
            !Create matrix H in band storage format, in preparation for eigen
            !decomposition with symmetric band matrix routine.
            h_bnd = 0.0_rp
            do jb = 1, j
                jgbeg = (jb-1)*s; igbeg = jgbeg
                do jcol = 1,s
                    jg = jgbeg + jcol
                    do irow = jcol,s
                        ig = igbeg + irow
                        h_bnd(ig-jg+1,jg) = h_bd(irow, jg)
                    end do
                end do
            end do

            do jb = 1, (j-1)
                jgbeg = (jb-1)*s; igbeg = jb*s
                do jcol = 1,s
                    jg = jgbeg + jcol
                    do irow = 1,s
                        ig = igbeg + irow
                        h_bnd(ig-jg+1,jg) = h_bsd(irow, jg)
                    end do
                end do
            end do

            !Eigen decomposition of H = P * d * P^T
            !eigvlh <- d; eigvch <- P
            call dsbevd('V', 'L', j*s, 2*s-1, h_bnd(:,1:j*s), 2*s, eigvlh, &
                        eigvch, s*dimk_max, work, size(work), iwork, &
                        size(iwork), info)
            if (info /= 0) then
                write(msg,'(a,i0,a,i0)') '<calc_bdw_blanc> dsbevd err = ', info, &
                    ' on iteration ', j
                call logger%log_msg(msg)
                ierr = 1; return
            end if

            !Square root of the eigen values
            eigvlh(1:j*s) = sqrt(eigvlh(1:j*s))

            !tmp[1:j*s,1:s] <- eigvch[1:s,1:j*s]^T * matr
            call dgemm('T', 'N', j*s, s, s, 1.0_rp, eigvch, s*dimk_max, &
                        matr, s, 0.0_rp, tmp(:,1:s), s*dimk_max)

            !tmp[1:j*s,1:s] <- eigvlh[1:j*s] * tmp[1:j*s,1:s]
            do jcol = 1, s
                tmp(1:j*s,jcol) = eigvlh(1:j*s)*tmp(1:j*s,jcol) 
            end do

            !tmp[1:j*s,s+1:2*s] <- eigvch[1:j*s,1:j*s] * tmp[1:j*s,1:s]
            call dgemm('N', 'N', j*s, s, j*s, 1.0_rp, eigvch, s*dimk_max, &
                        tmp(:,1:s), s*dimk_max, 0.0_rp, tmp(:,s+1:2*s), s*dimk_max)

            !bdw <- lncv[:,1:j*s] * tmp[1:j*s,s+1:2*s]
            call dgemm('N', 'N', sizm, s, j*s, 1.0_rp, lncv, sizm, &
                    tmp(:,s+1:2*s), s*dimk_max, 0.0_rp, bdw, sizm)

            !Error calculation for j > dimk_min. For j == dimk_min just save
            !bdw for error calculation in the next iteration. Using 1-norm as
            !the matrix norm here.
            if (j > dimk_min) then
                err = maxval(sum(abs(bdw-bdw_old), 1)) &
                    / maxval(sum(abs(bdw_old), 1))
                if (err <= ethres) then
                    lconv = .true.; f = j
                    exit
                end if
            end if
            bdw_old = bdw
        end if
    end do

    end subroutine

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

subroutine calc_bdw_cholesky(mob, bdw, ierr)
    !! Calculates **B**._d_**W** using Cholesky decomposition.

    real(rp), dimension(:,:), intent(in out) :: mob
        !! (`sizm`,`sizm`) symmetric positive definite matrix.
        !! On _entry_, contains the mobility matrix.  On _return_,
        !! the upper triangular part is overwritten with the result of Cholesky
        !! decomposition.  The strictly lower triangular part is left intact.
    real(rp), dimension(:,:), intent(in out) :: bdw
        !! (`sizm`,`s`) matrix.
        !! On _entry_, contains `s` column vectors drawn from standard normal 
        !! distribution. On _return_, the `s` columns are overwritten with the 
        !! result of **B**._d_**W**.
    integer, intent(out) :: ierr
        !! Error flag.
    integer :: info
    character(len=256) :: msg

    ierr = 0

    !Cholesky decomposition: mob = U^T * U. The upper triangular part of mob is
    !overwritten with U. The strictly lower triangular part of mob is not
    !accessed.
    call dpotrf('U', sizm, mob, sizm, info)
    if (info /= 0) then
        write(msg,'(a,i0)') '<calc_bdw_cholesky> dpotrf err = ', info
        call logger%log_msg(msg)
        ierr = 1; return
    end if

    !Product of U (stored in the upper triangular part of mob) with each of the
    !random vectors: bdw = U*bdw
    if (s == 1) then
        call dtrmv('U', 'N', 'N', sizm, mob, sizm, bdw(:,1), 1)
    else
        call dtrmm('L', 'U', 'N', 'N', sizm, s, 1.0_rp, mob, sizm, bdw, sizm)
    end if

    end subroutine

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

subroutine brn_calc_dw(dw)
    !! Calculates standard normally distributed random vectors.

    real(rp), dimension(:,:), target, intent(out) :: dw
        !! (`n`,`m`) array.
    real(rp), dimension(:), pointer :: vec_dw
    type(c_ptr) :: loc_dw
    integer :: n, m

    n = size(dw,1); m = size(dw,2)
    dw = 0.0_rp
    loc_dw = c_loc(dw)
    call c_f_pointer(loc_dw, vec_dw, [n*m])
    call get_rv_gaussian(0.0_rp, 1.0_rp, vec_dw, 1000000) 

    end subroutine

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

subroutine calc_bdw_lanc_sor(mob, bdw, lconv, f, ierr)

    real(rp), dimension(:,:), intent(in) :: mob
        !! (`sizm`,`sizm`) symmetric positive definite matrix.
    real(rp), dimension(:,:), intent(in out) :: bdw
        !! (`sizm`,1) array.
        !! On _entry_, contains a vector drawn from standard normal 
        !! distribution. On _return_, is overwritten with the result of
        !! **B**.d**W**.
    logical, intent(out) :: lconv
        !! Returns *true* if converged, *false* otherwise.
    integer, intent(out) :: f
        !! The number of iterations performed.
    integer, intent(out) :: ierr
        !! Error flag.
    integer :: j, info, lwork, liwork, k
    real(rp), dimension(1) :: qwork  !Size 1 array for workspace query.
    integer, dimension(1)  :: qiwork !Size 1 array for workspace query.
    real(rp), dimension(dimk_max) :: errbnd
    real(rp), dimension(sizm) :: ritzv
    real(rp) :: normz, err, normm
    character(len=256) :: msg

    lconv = .false.; ierr = 0; f = 0
    tmp = 0.0_rp; eigvlh = 0.0_rp; eigvch = 0.0_rp
    h_d = 0.0_rp; h_sd = 0.0_rp
    bdw_old = 0.0_rp; lncv = 0.0_rp; w = 0.0_rp

    !Workspace query
    if ((.not. allocated(work)) .or. (.not. allocated(iwork))) then
        call dstevd('V', dimk_max, eigvlh, tmp(1:dimk_max,1), eigvch, dimk_max, &
                    qwork, -1, qiwork, -1, info)
        if (info /= 0) then
            write(msg,'(a,i0,a)') '<calc_bdw_lanc_sor> dstevd err = ', info, &
                ' on workspace query'
            call logger%log_msg(msg)
            ierr = 1; return
        end if
        lwork = int(qwork(1)); liwork = qiwork(1)
        allocate( work(lwork) ); allocate( iwork(liwork) )
    end if

    !Lanczos process
    normz = norm2(bdw(:,1))
    normm = maxval(sum(abs(mob),1))
    lncv(:,1) = bdw(:,1)/normz

    do j = 1, dimk_max
        f = f + 1 !Iteration count
        call dsymv('L', sizm, 1.0_rp, mob, sizm, lncv(:,j), 1, 0.0_rp, w(:,1), 1)
        if (j > 1) then
            w(:,1) = w(:,1) - h_sd(j-1)*lncv(:,j-1)
        end if
        h_d(j) = dot_product(w(:,1), lncv(:,j))

        if (j < dimk_max) then
            w(:,1) = w(:,1) - h_d(j)*lncv(:,j)
            do k = 1, j
                if (errbnd(k) < sqrt(epsilon(1.0_rp)*normm)) then
                    call dgemv('N', sizm, j, 1.0_rp, lncv, sizm, eigvch(:,k), 1, &
                                0.0_rp, ritzv, 1)
                    w(:,1) = w(:,1) - dot_product(ritzv,w(:,1))*ritzv
                end if
            end do
            h_sd(j) = norm2(w(:,1))
            lncv(:,j+1) = w(:,1)/h_sd(j)
        end if

        !Compute eigenvalues, eigenvectors and error bounds for H
        if (j > 1) then
            eigvlh(1:j) = h_d(1:j)
            tmp(1:j,1) = h_sd(1:j)

            !Eigen decomposition of H = P * d * P^T
            !eigvlh[1:j] <- d; eigvch[1:j,1:j] <- P
            call dstevd('V', j, eigvlh(1:j), tmp(1:j,1), eigvch, dimk_max, &
                    work, size(work), iwork, size(iwork), info)
            if (info /= 0) then
                write(msg,'(a,i0,a,i0)') '<calc_bdw_lanc_sor> dstevd err = ', info, &
                    ' on iteration ', j
                call logger%log_msg(msg)
                ierr = 1; return
            end if
            errbnd(1:j) = abs(h_sd(j))*eigvch(j,1:j)
        end if

        !Calculate B.dW and check for convergence
        if (j >= dimk_min) then
            !eigvlh(1:j) <- d^(1/2)
            eigvlh(1:j) = sqrt(eigvlh(1:j))

            !tmp[1:j,1] <- d^(1/2) * P^T * e_1, where e_1 if the first column of
            !an identity matrix of size j.
            tmp(1:j,1) = eigvch(1,1:j)
            tmp(1:j,1) = eigvlh(1:j)*tmp(1:j,1)

            !tmp[1:j,2] <- P * tmp[1:j,1]
            call dgemv('N', j, j, 1.0_rp, eigvch, dimk_max, tmp(1:j,1), 1, &
                        0.0_rp, tmp(1:j,2), 1)

            !bdw <- normz * lncv[:,1:j] * tmp[1:j,2]
            call dgemv('N', sizm, j, normz, lncv, sizm, tmp(1:j,2), 1, &
                        0.0_rp, bdw, 1)

            !Error calculation for j > dimk_min. For j = dimk_min just save
            !bdw for error calculation in the next iteration.
            if (j > dimk_min) then
                err = norm2(bdw(:,1) - bdw_old(:,1))/norm2(bdw_old(:,1))
                !write(*,*) j, 'err = ', err
                if (err <= ethres) then
                    lconv = .true.
                    exit
                end if
            end if
            bdw_old = bdw
        end if
    end do

    end subroutine

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

end module brown_m
