CG

Matrix<type_numeric,dim,1> r = b - A * x;
Matrix<type_numeric,dim,1> p = r;
type_numeric err_i0 = r | r;

for (uint32_t i=0u;i<max_iters;++i) {
    Matrix<type_numeric,dim,1> Ap = A * p;

    type_numeric alpha = err_i0 / (p | Ap);
    x += alpha * p;
    r -= alpha * Ap;

    type_numeric err_i1 = r | r;
    if (sqrt(err_i1)<tolerance) break;

    p = r + (err_i1/err_i0)*p;

    err_i0 = err_i1;
}