gmres_rdp Subroutine

    subroutine gmres_rdp(A, b, x, info, rtol, atol, preconditioner, options, transpose)
        class(abstract_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 (optional).
        class(abstract_opts), optional, intent(in) :: options
        !! GMRES options.   
        logical, optional, intent(in) :: transpose
        !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used.

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

        ! Options.
        integer :: kdim, maxiter
        real(dp) :: tol, rtol_, atol_
        logical :: trans
        type(gmres_dp_opts) :: opts

        ! Krylov subspace
        class(abstract_vector_rdp), allocatable :: V(:)
        ! Hessenberg matrix.
        real(dp), allocatable :: H(:, :)
        ! Least-squares variables.
        real(dp), allocatable :: y(:), e(:)
        real(dp) :: beta

        ! Preconditioner
        logical :: has_precond
        class(abstract_precond_rdp), allocatable :: precond

        ! Miscellaneous.
        integer :: i, k, niter
        real(dp), allocatable :: alpha(:)
        class(abstract_vector_rdp), allocatable :: dx, wrk
        character(len=256) :: msg

        ! Deals with the optional args.
        rtol_ = optval(rtol, rtol_dp)
        atol_ = optval(atol, atol_dp)
        if (present(options)) then
            select type (options)
            type is (gmres_dp_opts)
                opts = options
            end select
        else
            opts = gmres_dp_opts()
        endif

        kdim = opts%kdim ; maxiter = opts%maxiter
        tol = atol_ + rtol_ * b%norm()
        trans = optval(transpose, .false.)

        ! Deals with the preconditioner.
        if (present(preconditioner)) then
            has_precond = .true.
            allocate(precond, source=preconditioner)
        else
            has_precond = .false.
        endif

        ! Initialize working variables.
        allocate(wrk, source=b) ; call wrk%zero()
        allocate(V(kdim+1), source=b) ; call zero_basis(V)
        allocate(H(kdim+1, kdim)) ; H = 0.0_dp
        allocate(y(kdim)) ; y = 0.0_dp
        allocate(alpha(kdim)) ; alpha = 0.0_dp
        allocate(e(kdim+1)) ; e = 0.0_dp

        info = 0 ; niter = 0

        ! Initial Krylov vector.
        if (x%norm() > 0) then
            if (trans) then
                call A%rmatvec(x, V(1))
            else
                call A%matvec(x, V(1))
            endif
        endif

        call V(1)%sub(b) ; call V(1)%chsgn()
        beta = V(1)%norm() ; call V(1)%scal(one_rdp/beta)

        ! Iterative solver.
        gmres_iter : do i = 1, maxiter
            ! Zero-out variables.
            H = 0.0_dp ; y = 0.0_dp ; e = 0.0_dp ; e(1) = beta
            call zero_basis(V(2:))

            ! Arnoldi factorization.
            arnoldi_fact: do k = 1, kdim
                ! Preconditioner.
                wrk = V(k) ; if (has_precond) call precond%apply(wrk)

                ! Matrix-vector product.
                if (trans) then
                    call A%rmatvec(wrk, V(k+1))
                else
                    call A%matvec(wrk, V(k+1))
                endif

                ! Double Gram-Schmid orthogonalization
                call double_gram_schmidt_step(V(k+1), V(:k), info, if_chk_orthonormal=.false., beta=H(:k, k))
                call check_info(info, 'double_gram_schmidt_step', module=this_module, procedure='gmres_rdp')

                ! Update Hessenberg matrix and normalize residual Krylov vector.
                H(k+1, k) = V(k+1)%norm()
                if (abs(H(k+1, k)) > tol) then
                    call V(k+1)%scal(one_rdp / H(k+1, k))
                endif

                ! Least-squares problem.
                y(:k) = lstsq(H(:k+1, :k), e(:k+1))

                ! Compute residual.
                beta = norm2(abs(e(:k+1) - matmul(H(:k+1, :k), y(:k))))

                ! Current number of iterations performed.
                niter = niter + 1

                ! Check convergence.
                write(msg,'(A,I3,2(A,E9.2))') 'GMRES(k)   inner step ', k, ': |res|= ', &
                            & abs(beta), ', tol= ', tol
                call logger%log_information(msg, module=this_module, procedure='gmres_rdp')
                if (abs(beta) <= tol) then
                    exit arnoldi_fact
                endif
            enddo arnoldi_fact

            ! Update solution.
            k = min(k, kdim) ; call linear_combination(dx, V(:k), y(:k))
            if (has_precond) call precond%apply(dx) ; call x%add(dx)

            ! Recompute residual for sanity check.
            if (trans) then
                call A%rmatvec(x, v(1))
            else
                call A%matvec(x, v(1))
            endif
            call v(1)%sub(b) ; call v(1)%chsgn()

            ! Initialize new starting Krylov vector if needed.
            beta = v(1)%norm() ; call v(1)%scal(one_rdp / beta)

            write(msg,'(A,I3,2(A,E9.2))') 'GMRES(k) outer step   ', i, ': |res|= ', &
                            & abs(beta), ', tol= ', tol
            call logger%log_information(msg, module=this_module, procedure='gmres_rdp')

            ! Exit gmres if desired accuracy is reached.
            if (abs(beta) <= tol) exit gmres_iter

        enddo gmres_iter

        ! Returns the number of iterations.
        info = niter

        return
    end subroutine gmres_rdp