conjugate_gradient_solver Subroutine

    subroutine conjugate_gradient_solver(n, rhs, x, eel, matvec, precnd, &
                                         arg_tol, arg_n_iter)
        !! Conjugate gradient solver (TODO)
        ! TODO add more printing
    
        use mod_constants, only: eps_rp
        use mod_memory, only: mallocate, mfree

        implicit none

        integer(ip), intent(in) :: n
        !! Size of the matrix
        real(rp), intent(in), optional :: arg_tol
        !! Optional convergence criterion in input, if not present
        !! OMMP_DEFAULT_SOLVER_TOL is used.
        real(rp) :: tol
        !! Convergence criterion, it is required that RMS norm < tol

        integer(ip), intent(in), optional :: arg_n_iter
        !! Optional maximum number of iterations for the solver, if not present
        !! OMMP_DEFAULT_SOLVER_ITER is used.
        integer(ip) :: n_iter
        !! Maximum number of iterations for the solver 

        real(rp), dimension(n), intent(in) :: rhs
        !! Right hand side of the linear system
        real(rp), dimension(n), intent(inout) :: x
        !! In input, initial guess for the solver, in output the solution
        type(ommp_electrostatics_type), intent(in) :: eel
        !! Electrostatics data structure
        external :: matvec
        !! Routine to perform matrix-vector product
        external :: precnd
        !! Preconditioner routine

        integer(ip) :: it
        real(rp) :: rms_norm, alpha, gnew, gold, gama
        real(rp), allocatable :: r(:), p(:), h(:), z(:)
        character(len=OMMP_STR_CHAR_MAX) :: msg

        ! Optional arguments handling
        if(present(arg_tol)) then
            tol = arg_tol
        else
            tol = OMMP_DEFAULT_SOLVER_TOL
        end if

        if(present(arg_n_iter)) then
            n_iter = arg_n_iter
        else
            n_iter = OMMP_DEFAULT_SOLVER_ITER
        end if

        call ommp_message("Solving linear system with CG solver", OMMP_VERBOSE_LOW)
        write(msg, "(A, I4)") "Max iter:", n_iter
        call ommp_message(msg, OMMP_VERBOSE_LOW)
        write(msg, "(A, E8.1)") "Tolerance: ", tol
        call ommp_message(msg, OMMP_VERBOSE_LOW)

        call mallocate('conjugate_gradient_solver [r]', n, r)
        call mallocate('conjugate_gradient_solver [p]', n, p)
        call mallocate('conjugate_gradient_solver [h]', n, h)
        call mallocate('conjugate_gradient_solver [z]', n, z)

        ! compute a guess, if required:
        rms_norm = dot_product(x,x)
        if(rms_norm < eps_rp) then
            call ommp_message("Input guess has zero norm, generating a guess&
                              & from preconditioner.", OMMP_VERBOSE_HIGH)
            call precnd(eel, x, x)
        else
            call ommp_message("Using input guess as a starting point for&
                              & iterative solver.", OMMP_VERBOSE_HIGH)
        end if

        ! compute the residual:
        call matvec(eel, x, z, .true.)
        r = rhs - z
        ! apply the preconditioner and get the first direction:
        call precnd(eel, r, z)
        p = z
        gold = dot_product(r, z)
        gama = 0.0_rp

        do it = 1, n_iter
            ! compute the step:
            call matvec(eel, p, h, .true.)
            gama = dot_product(h, p)

            ! unlikely quick return:
            if(abs(gama) < eps_rp) then
                call ommp_message("Direction vector with zero norm, exiting &
                                  &iterative solver.", OMMP_VERBOSE_HIGH)
                exit
            end if

            alpha = gold / gama
            x = x + alpha * p
            r = r - alpha * h

            ! apply the preconditioner:
            call precnd(eel, r, z)
            gnew = dot_product(r, z)
            rms_norm = sqrt(gnew/dble(n))

            write(msg, "('iter=',i4,' residual rms norm: ', d14.4)") it, rms_norm
            call ommp_message(msg, OMMP_VERBOSE_HIGH)

            ! Check convergence
            if(rms_norm < tol) then
                call ommp_message("Required convergence threshold reached, &
                                  &exiting iterative solver.", OMMP_VERBOSE_HIGH)
                exit
            end if

            ! compute the next direction:
            gama = gnew/gold
            p    = gama*p + z
            gold = gnew
        end do

        call mfree('conjugate_gradient_solver [r]', r)
        call mfree('conjugate_gradient_solver [p]', p)
        call mfree('conjugate_gradient_solver [h]', h)
        call mfree('conjugate_gradient_solver [z]', z)

        if(rms_norm > tol .and. abs(gama) > eps_rp) then
            call fatal_error("Iterative solver did not converged")
        end if

    end subroutine conjugate_gradient_solver