我正在使用自己的例程来计算复杂的正定对称矩阵的 Cholesky 分解。
我的代码如下所示:
void CholeskyDecomposition(Complex **L, Complex **A, const int dim)
{
Complex sum, t;
Complex sum1;
for (int k = 0; k < dim; k++)
{
// First, calculate the diagonal element of a column
sum1 = A[k][k];
for (int j = 0; j < k; j++)
{
Complex product;
Cconjmul(product, L[j][k], L[j][k]);
Csub(sum1, sum1, product);
}
Csqrt(&sum1, sum1);
L[k][k] = sum1;
// Now, calculate the elements right the diagonal element
for (int i = k + 1; i < dim; i++)
{
sum = A[k][i];
for (int j = 0; j < k; j++)
{
Cconjmul(t, L[j][k], L[j][i]);
Csub(sum, sum, t);
}
Cdiv(L[k][i], sum, sum1);
}
}
}
我针对 MATLAB Cholesky-Factorization 测试了这段代码,结果是
A =
13.6393 + 0.0000i 1.8844 + 3.4319i -4.8788 - 4.0195i
1.8844 - 3.4319i 5.7265 + 0.0000i 1.1568 - 0.4937i
-4.8788 + 4.0195i 1.1568 + 0.4937i 6.6512 + 0.0000i
K>> L = chol(A)
L =
3.6931 + 0.0000i 0.5102 + 0.9293i -1.3210 - 1.0884i
0.0000 + 0.0000i 2.1454 + 0.0000i 1.3248 - 0.5435i
0.0000 + 0.0000i 0.0000 + 0.0000i 1.2927 + 0.0000i
我的 C-Implementation 的结果似乎并非完全错误。第一行与 MATLAB 解决方案相同。第二行和第三行在最后一列包含错误:
MATLAB <--> C
1.3248-0.5435i <--> 1.3248+0.08324i
1.2927+0.0000i <--> 1.3998+0.00000i
谁能帮我确定错误来源?