cg_rdp Subroutine

    subroutine cg_rdp(A, b, x, info, rtol, atol, preconditioner, options)
        class(abstract_sym_linop_rdp), intent(in) :: A
        !! Linear operator to be inverted.
        class(abstract_vector_rdp), intent(in) :: b
        !! Right-hand side vector.
        class(abstract_vector_rdp), intent(inout) :: x
        !! Solution vector.
        integer, intent(out) :: info
        !! Information flag.
        real(dp), optional, intent(in) :: rtol
        !! Relative solver tolerance
        real(dp), optional, intent(in) :: atol
        !! Absolute solver tolerance
        class(abstract_precond_rdp), optional, intent(in) :: preconditioner
        !! Preconditioner (not yet supported).
        type(cg_dp_opts), optional, intent(in) :: options
        !! Options for the conjugate gradient solver.

        !----------------------------------------
        !-----     Internal variables      ------
        !----------------------------------------

        ! Options.
        integer :: maxiter
        real(dp) :: tol, rtol_, atol_
        type(cg_dp_opts) :: opts

        ! Working variables.
        class(abstract_vector_rdp), allocatable :: r, p, Ap
        real(dp) :: alpha, beta, r_dot_r_old, r_dot_r_new
        real(dp) :: residual

        ! Miscellaneous.
        integer :: i
        character(len=256) :: msg

        ! Deals with the optional args.
        rtol_ = optval(rtol, rtol_dp)
        atol_ = optval(atol, atol_dp)
        if (present(options)) then
            opts = options
        else
            opts = cg_dp_opts()
        endif
        tol = atol_ + rtol_ * b%norm() ; maxiter = opts%maxiter

        ! Initialize vectors.
        allocate(r, source=b)  ; call r%zero()
        allocate(p, source=b)  ; call p%zero()
        allocate(Ap, source=b) ; call Ap%zero()

        info = 0

        ! Compute initial residual r = b - Ax.
        if (x%norm() > 0) call A%matvec(x, r)
        call r%sub(b) ; call r%chsgn()

        ! Initialize direction vector.
        p = r

        ! Initialize dot product of residual r_dot_r_old = r' * r
        r_dot_r_old = r%dot(r)

        ! Conjugate gradient iteration.
        cg_loop: do i = 1, maxiter
            ! Compute A @ p
            call A%matvec(p, Ap)
            ! Compute step size.
            alpha = r_dot_r_old / p%dot(Ap)
            ! Update solution x = x + alpha*p
            call x%axpby(one_rdp, p, alpha)
            ! Update residual r = r - alpha*Ap
            call r%axpby(one_rdp, Ap, -alpha)
            ! Compute new dot product of residual r_dot_r_new = r' * r.
            r_dot_r_new = r%dot(r)
            ! Check for convergence.
            residual = sqrt(r_dot_r_new)
            ! Current number of iterations performed.
            info = info + 1

            if (residual < tol) exit cg_loop

            ! Compute new direction beta = r_dot_r_new / r_dot_r_old.
            beta = r_dot_r_new / r_dot_r_old
            ! Update direction p = beta*p + r
            call p%axpby(beta, r, one_rdp)
            ! Update r_dot_r_old for next iteration.
            r_dot_r_old = r_dot_r_new

            write(msg,'(A,I3,2(A,E9.2))') 'CG step ', i, ': res= ', residual, ', tol= ', tol
            call logger%log_information(msg, module=this_module, procedure='cg_rdp')
        enddo cg_loop
        
        return
    end subroutine cg_rdp