我试着用C++写一个简单的共轭梯度法函数--很简单吧？但它确实不能解任何方程组。它收敛到完全错误的值上。
我已经把函数写在这里了，并提供了主函数，其中我使用对称矩阵M生成正定矩阵A。我已经提供了所有的代码，看看是否有人可以发现任何错误。请不要因为我这样做而嘲笑我--这是一个很好理解的算法，不提供所有的代码是没有意义的，对吗？请假设ApplyLHS和innerProduct函数工作正常。

template <typename TData>
void DoConjGrad(
    std::vector<TData> &in,
    std::vector<TData> &out,
    std::vector<TData> &A
)
{
    TData tol = 10e-8;

    size_t N = in.size();

    std::vector<TData> r(N);
    std::vector<TData> p(N);
    std::vector<TData> w(N);
    size_t k = 0;

    TData rTr;
    TData rTr_prev;
    TData pTw;

    TData alpha;
    TData beta;

    // initial guess for x0
    std::fill(out.begin(), out.end(), 0);

    // r0 = b - Ax0
    ApplyLHS(out, r, A);
    std::transform(in.cbegin(), in.cend(), r.cbegin(), r.begin(), 
        [](const TData &bElem, const TData &rElem){ return bElem - rElem; } );

    // p0 := r0
    std::copy(r.cbegin(), r.cend(), p.begin());

    // calculate (rTr)0
    rTr = innerProduct(r, r);

    for (int i = 0; i < 100; ++i)
    {
        ApplyLHS(p, w, A);
        pTw = innerProduct(p, w);

        alpha = rTr / pTw;

        std::transform(out.cbegin(), out.cend(), p.cbegin(), out.begin(),
            [alpha](const TData &xElem, const TData &pElem)
            {
                return xElem + alpha*pElem;
            });

        std::transform(r.cbegin(), r.cend(), w.cbegin(), r.begin(),
            [alpha](const TData &rElem, const TData &wElem)
            {
                return rElem - alpha*wElem;
            });

        if (rTr < tol)
            break;

        rTr_prev = rTr;
        rTr = innerProduct(r, r);

        beta = rTr / rTr_prev;

        std::transform(r.cbegin(), r.cend(), p.cbegin(), p.begin(),
            [beta](const TData &rElem, const TData &pElem)
            {
                return rElem + beta*pElem;
            });
    }
}

主要：

int main()
{
    srand(100);
    size_t N = 2;

    std::vector<double> A(N * N);
    std::vector<double> M(N * N);
    std::vector<double> x_true(N);
    std::vector<double> x_calc(N);
    std::vector<double> b(N);
    
    for (int i = 0; i < N; ++i)
        for (int j = 0; j <= i; ++j)
        {
            double r = 2*((double)rand())/((double)RAND_MAX) - 1;
            M[i*N + j] = r;
            M[j*N + i] = r;
        }
    
    // get positive semi-definite matrix = M.T @ M
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
            for (int k = 0; k < N; ++k)
                A[i*N + j] = M[k*N + i] * M[k*N + j];

    // generate random x
    for (int i = 0; i < N; ++i)
        x_true[i] = 2*((double)rand())/((double)RAND_MAX) - 1;

    // calculate b
    ApplyLHS(x_true, b, A);

    // solve for x calc
    DoConjGrad<double>(b, x_calc, A);
}

辅助函数：

template <typename TData>
void printMat(
    size_t nr,
    size_t nc,
    std::vector<TData> &mat
)
{
    for (int i = 0; i < nr; ++i)
    {
        for (int j = 0; j < nc; ++j)
        {
            std::cout.precision(5);
            std::cout.width(10);
            std::cout << mat[i*nc + j] << "\t";
        }
        std::cout << "\n";
    }
}

template <typename TData>
void ApplyLHS(
    std::vector<TData> &in,
    std::vector<TData> &out,
    std::vector<TData> &A
)
{
    size_t N = in.size();
    std::fill(out.begin(), out.end(), 0);

    for (int i = 0; i < N; ++i)
    for (int j = 0; j < N; ++j)
        out[i] += in[j] * A[i*N + j];
}

template <typename TData>
TData innerProduct(
    std::vector<TData> &a,
    std::vector<TData> &b
)
{
    TData result = 0;
    for (int i = 0; i < a.size(); ++i)
    {
        result += a[i] * b[i];
    }
    return result;
}