Determinant using Gauss Elimination We had a question in examination to compute determinant of general square matrix.
Now I had made the normal program using Laplace Formula of determinant, ie, recursive solution.

Now I was thinking of reducing it into Upper Triangle Matrix and then the product of diagonal elements gives the determinant. 

The code is following.

#include &lt;iostream&gt;
using namespace std;

const int N = 5;
int main()
{
    double array[N][N], det = 1;

    // initialise
    cout &lt;&lt; "Enter " &lt;&lt; N &lt;&lt; "*" &lt;&lt; N &lt;&lt; " matrix:\n";
    for (int i = 0; i &lt; N; ++i) {
        for (int j = 0; j &lt; N; ++j) {
            cin &gt;&gt; array[i][j];
        }
    }

    for (int i = 0; i &lt; N; ++i) 
    {
        bool flag = false; // flag := are all entries below a[i][i] including it zero?
        if (array[i][i] == 0) // If a[i][i] is zero then check rows below and swap
        {
            flag = true;
            for (int j = i; j &lt; N; ++j)
            {
                if (array[j][i] != 0) {
                    det *= -1;
                    for (int k = 0; k &lt; N; ++k) { 
                        double t = array[i][k]; // swapping
                        array[i][k] = array[j][k];
                        array[j][k] = t;
                        flag = false;
                    }
                }
            }
        }

        if (flag == true) {
            det = 0;
            break;
        }

        else {

            for (int j = i+1; j &lt; N; ++j)
            {
                double store = array[j][i];
                for (int k = i; k &lt; N; ++k) {
                    array[j][k] -= (array[i][k]*store)/array[i][i];

                } 
            }
            det *= array[i][i];
        }
    }

    cout &lt;&lt; "Determinant: " &lt;&lt; det;
    return 0;
}


The problems with this approach is that matrix needs to be of floating point datatypes. Worse, it gives wrong result for least significant place (ones digit) for all integral entries if I use datatype of entries as float. To overcome that, we can use double.

But still I :


need tips on improving performance of this method
need to know whether this method is better than Simple Laplace Formula?
Some general remarks:


Don't use namespace std;, see for example Why is “using namespace std;” considered bad practice?.
Don't put all the code into main(). Computing the determinant in a separate
function increases the overall clarity of the program and makes it easier to add
test cases. In addition, that gives you a function which can be reused in other
programs.
Swapping two values can be done in C++ simply with std::swap.
return 0; at the end of the main program can be omitted.


At some places your code does unnecessary work:


If the diagonal element a[i][i] is zero then all subsequent rows 
with array[j][i] != 0 are swapped with the ith row. It suffices to
swap one row with a non-zero leading element.
When swapping rows i and j, it suffices to swap the elements starting
at column i because the preceding elements are not used anymore.
Similarly, when adding a multiple of row i to row j, it suffices to
update elements starting at column i + 1, because all values of column
i will not be read anymore.


As already said in the comments, the Gaussian elimination is faster than the
Laplace expansion for large matrices (\$ O(N^3) \$ vs \$ O(N!) \$ complexity).
However, the “pivoting” (i.e. which rows to swap if an diagonal element is zero)
can be improved. A common choice is “partial pivoting”:


  In partial pivoting, the algorithm selects the entry with largest absolute value from the column of the matrix that is currently being considered as the pivot element. Partial pivoting is generally sufficient to adequately reduce round-off error. ...


Using a (nested) std::vector for the matrix storage makes it
easier to pass it an argument to another function.
In addition, it allows to swap two rows in constant time with vector::swap.

Putting it all together, the determinant function could look like this:

#include &lt;vector&gt;
#include &lt;cmath&gt;

double determinant(std::vector&lt;std::vector&lt;double&gt;&gt; &amp;matrix) {
    int N = static_cast&lt;int&gt;(matrix.size());
    double det = 1;

    for (int i = 0; i &lt; N; ++i) {

        double pivotElement = matrix[i][i];
        int pivotRow = i;
        for (int row = i + 1; row &lt; N; ++row) {
            if (std::abs(matrix[row][i]) &gt; std::abs(pivotElement)) {
                pivotElement = matrix[row][i];
                pivotRow = row;
            }
        }
        if (pivotElement == 0.0) {
            return 0.0;
        }
        if (pivotRow != i) {
            matrix[i].swap(matrix[pivotRow]);
            det *= -1.0;
        }
        det *= pivotElement;

        for (int row = i + 1; row &lt; N; ++row) {
            for (int col = i + 1; col &lt; N; ++col) {
                matrix[row][col] -= matrix[row][i] * matrix[i][col] / pivotElement;
            }
        }
    }

    return det;
}