...
    r = b - A*x
    p = r
    k = 1
    do
       rtr = dot(r, r)
       Ap = A*p
       alpha = rtr / dot(p, Ap)
       x = x + alpha*p
       r = r - alpha*Ap
       if (vec_norm(r) <= tolerance .or. k == max_iters) exit
       beta = dot(r, r) / rtr
       p = r + beta*p
       k = k + 1
    end do
...