C - stack smashing detected

Question

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;
}

Answer 1

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.