我需要实现矩阵 A 的一个非常简单的就地 LU 分解。我正在使用高斯消除,我想用 3x3 矩阵测试它。问题是,我一直收到stack smashing
错误,我不知道为什么。我在我的代码中没有看到任何问题,它可以做到这一点。你有什么想法吗?
问题可能出在因子分解块中。
###My code:###
#include <stdio.h>
int main() {
int n = 3; // matrix size
int A[3][3] = {
{1, 4, 7},
{2, 5, 8},
{3, 6, 10}
};
printf("Matrix A:n");
for( int i=0; i < n; i++ ) {
for( int j=0; j < n; j++ ) {
printf("%d ", A[i][j]);
if ( j % 2 == 0 && j != 0 ) {
printf("n");
}
}
}
// FACTORIZATION
int k;
int rows;
for( k = 0; k < n; k++ ) {
rows = k + k+1;
A[rows][k] = A[rows][k]/A[k][k];
A[rows][rows] = A[rows][rows] - A[rows][k] * A[k][rows];
printf("k: %dn", k);
}
printf("Matrix after decomp:n");
for( int i=0; i < n; i++ ) {
for( int j=0; j < n; j++ ) {
printf("%d ", A[i][j]);
if ( j % 3 == 0 && j != 0 ) {
printf("n");
}
}
}
return 0;
}
你的错误很可能在这里:
rows = k + k+1;
A[rows][k] = A[rows][k]/A[k][k];
A[rows][rows] = A[rows][rows] - A[rows][k] * A[k][rows];
这意味着rows
遍历值 1、3、5;然后用于访问只有三个元素的数组。这确实会溢出,因为其中唯一有效的偏移量是 1。
编辑:看看你的Matlab代码,它正在做一些完全不同的事情,因为rows = k + 1:n
rows
设置为一个小向量,然后它使用矩阵的拼接,这是C不支持的原语。您需要使用显式循环重新实现该和矩阵乘法A(rows, k) * A(k, rows)
。
你最初的 Matlab 代码是(Matlab 有 1 个基于 1 的索引):
for k = 1:n - 1
rows = k + 1:n
A(rows, k) = A(rows, k) / A(k, k)
A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows)
end
这样做rows = k + 1:n
是将rows
设置为表示范围。表达式A(rows, k)
实际上是对矩阵的矢量形切片的引用,Matlab 可以将向量除以标量。
在最后一行,A(rows, rows)
是一个矩阵形切片,A(rows, k) * A(k, rows)
是矩阵乘法,例如将维度为(1,3)和(3,1)的矩阵相乘得到(3,3)中的一个。
在 C 语言中,您无法使用内置的 =
和 /
运算符来执行此操作。
C 等效项为:
for ( int k = 0; k < n - 1; ++k )
{
// A(rows, k) = A(rows, k) / A(k, k)
for ( int row = k + 1; row < n; ++row )
A[row][k] /= A[k][k];
// A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows)
for ( int row = k + 1; row < n; ++row )
for ( int col = k + 1; col < n; ++col )
A[row][col] -= A[row][k] * A[k][col];
}
(免责声明:未经测试!
第一部分很简单:向量中的每个值都被标量除以。
但是,第二行更复杂。Matlab代码包括矩阵乘法和矩阵减法;以及从矩阵中提取子矩阵的操作。如果我们试图将其直接翻译成 C,那就非常复杂了。
我们需要使用两个嵌套循环来迭代行和列,以在方阵上执行此操作。