...
    r = b-matmul(A,x)
    p = r
    k = 1
    do
       Ap = matmul(A,p)
       rtr = dot_product(r,r)
       alpha = rtr / dot_product(p,Ap)
       x = x + alpha * p
       r = r - alpha * Ap
       if (sqrt(sum(r**2)) <= tolerance .or. k == max_iteration) exit
       beta = dot_product(r,r) / rtr
       p = r + beta * p
       k = k + 1
    end do
  ...