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 <vector>
#include <cmath>
double determinant(std::vector<std::vector<double>> &matrix) {
int N = static_cast<int>(matrix.size());
double det = 1;
for (int i = 0; i < N; ++i) {
double pivotElement = matrix[i][i];
int pivotRow = i;
for (int row = i + 1; row < N; ++row) {
if (std::abs(matrix[row][i]) > std::abs(pivotElement)) {
pivotElement = matrix[row][i];
pivotRow = row;
}
}
if (pivotElement == 0.0) {
return 0.0;
}
if (pivotRow != i) {
for (int col = i; col < N; ++col) {
std::matrix[i].swap(matrix[i][col], matrix[pivotRow][col]matrix[pivotRow]);
}
det *= -1.0;
}
det *= pivotElement;
for (int row = i + 1; row < N; ++row) {
for (int col = i + 1; col < N; ++col) {
matrix[row][col] -= matrix[row][i] * matrix[i][col] / pivotElement;
}
}
}
return det;
}
Here I'm using a (nested) std::vector for the matrix storage because that is
easier to pass as an argument to another function.
Putting it all together, the determinant function could look like this:
#include <vector>
#include <cmath>
double determinant(std::vector<std::vector<double>> &matrix) {
int N = static_cast<int>(matrix.size());
double det = 1;
for (int i = 0; i < N; ++i) {
double pivotElement = matrix[i][i];
int pivotRow = i;
for (int row = i + 1; row < N; ++row) {
if (std::abs(matrix[row][i]) > std::abs(pivotElement)) {
pivotElement = matrix[row][i];
pivotRow = row;
}
}
if (pivotElement == 0.0) {
return 0.0;
}
if (pivotRow != i) {
for (int col = i; col < N; ++col) {
std::swap(matrix[i][col], matrix[pivotRow][col]);
}
det *= -1.0;
}
det *= pivotElement;
for (int row = i + 1; row < N; ++row) {
for (int col = i + 1; col < N; ++col) {
matrix[row][col] -= matrix[row][i] * matrix[i][col] / pivotElement;
}
}
}
return det;
}
Here I'm using a (nested) std::vector for the matrix storage because that is
easier to pass as an argument to another function.
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 <vector>
#include <cmath>
double determinant(std::vector<std::vector<double>> &matrix) {
int N = static_cast<int>(matrix.size());
double det = 1;
for (int i = 0; i < N; ++i) {
double pivotElement = matrix[i][i];
int pivotRow = i;
for (int row = i + 1; row < N; ++row) {
if (std::abs(matrix[row][i]) > 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 < N; ++row) {
for (int col = i + 1; col < N; ++col) {
matrix[row][col] -= matrix[row][i] * matrix[i][col] / pivotElement;
}
}
}
return det;
}
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.
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.
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.
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.
At some places your code does unnecessary work:
If the diagonal element
a[i][i]is zero then all subsequent rows witharray[j][i] != 0are swapped with theith row. It suffices to swap one row with a non-zero leading element.When swapping rows
iandj, it suffices to swap the elements starting at columnibecause the preceding elements are not used anymore.Similarly, when adding a multiple of row
ito rowj, it suffices to update elements starting at columni + 1, because all values of columniwill 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. ...
Putting it all together, the determinant function could look like this:
#include <vector>
#include <cmath>
double determinant(std::vector<std::vector<double>> &matrix) {
int N = static_cast<int>(matrix.size());
double det = 1;
for (int i = 0; i < N; ++i) {
double pivotElement = matrix[i][i];
int pivotRow = i;
for (int row = i + 1; row < N; ++row) {
if (std::abs(matrix[row][i]) > std::abs(pivotElement)) {
pivotElement = matrix[row][i];
pivotRow = row;
}
}
if (pivotElement == 0.0) {
return 0.0;
}
if (pivotRow != i) {
for (int col = i; col < N; ++col) {
std::swap(matrix[i][col], matrix[pivotRow][col]);
}
det *= -1.0;
}
det *= pivotElement;
for (int row = i + 1; row < N; ++row) {
for (int col = i + 1; col < N; ++col) {
matrix[row][col] -= matrix[row][i] * matrix[i][col] / pivotElement;
}
}
}
return det;
}
Here I'm using a (nested) std::vector for the matrix storage because that is
easier to pass as an argument to another function.