Reputation: 2977
I need to implement a pretty easy in-place LU-decomposition of matrix A. I'm using Gaussian elimination and I want to test it with a 3x3 matrix. The problem is, I keep getting stack smashing
error and I don't have any idea why. I don't see any problems in my code, which could do this. Do you have any idea?
The problem is probably in the Factorization block.
###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: %d\n", 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;
}
Upvotes: 0
Views: 502
Reputation: 141554
Your original Matlab code was (Matlab has 1-based indexing):
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
What rows = k + 1:n
this does is that it sets rows
to represent a range. The expression A(rows, k)
is actually a reference to a vector-shaped slice of the matrix, and Matlab can divide a vector by a scalar.
On the last line, A(rows, rows)
is a matrix-shaped slice , and A(rows, k) * A(k, rows)
is a matrix multiplication, e.g. multiplying matrices of dimension (1,3) and (3,1) to get one of (3,3).
In C you can't do that using the builtin =
and /
operators.
The C equivalent is:
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];
}
(disclaimer: untested!)
The first part is straightforward: every value in a vector is being divided by a scalar.
However, the second line is more complicated. The Matlab code includes a matrix multiplication and a matrix subtraction ; and also the operation of extracting a sub-matrix from a matrix. If we tried to write a direct translation of that to C, it is very complicated.
We need to use two nested loops to iterate over the rows and columns to perform this operation on the square matrix.
Upvotes: 2
Reputation: 25855
Your error is most likely here:
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];
This means that rows
goes through the values 1, 3, 5; and is then used to access an array with only three elements. That would, indeed, overflow, as the only valid offset among those is 1.
EDIT: Looking at your Matlab code, it is doing something completely different, as rows = k + 1:n
sets rows
to a small vector, which it then uses the splice the matrix, something C does not support as a primitive. You would need to reimplement both that and the matrix multiplication A(rows, k) * A(k, rows)
using explicit loops.
Upvotes: 6