pcgrbpini Subroutine

  subroutine pcgrbpini(r, pk, error, d, xm_in, a_pk)
    use precision, only: dp

    implicit none

    real(kind=dp), intent(out), dimension(:) :: r
    real(kind=dp), intent(out), dimension(:) :: pk
    real(kind=dp), intent(out) :: error
    real(kind=dp), intent(in), dimension(:) :: d
    real(kind=dp), intent(in), dimension(:) :: xm_in
    real(kind=dp), intent(in), dimension(:) :: a_pk

    real(kind=dp) :: beta

  ! R ini and R norm(error)
    r = d-a_pk
    error = dot_product(r, r)

  ! Beta ini
    beta = 1.0_dp/dot_product(r**2, xm_in)

  ! pk ini
    pk = beta*xm_in*r

  end subroutine pcgrbpini