1

I have written the following code to calculate the determinant of a N*N matrix. it works perfectly for the matrix 4*4 and 5*5. but it could not find the determinant of the 40*40 matrix named Z_z. The elements of matrix Z_z are presented herein.

#include <iostream>
int const N_M=40;
void getCofactor(double A[N_M][N_M], double A_rc[N_M][N_M], int r,int c, int n) 
{ int i = 0, j = 0; 
// Looping for each element of the matrix 
for (int row = 0; row < n; row++) 
{ 
    for (int col = 0; col < n; col++) 
    { 
        //  Copying into temporary matrix only those element 
        //  which are not in given row and column 
        if (row != r && col != c) 
        { 
            A_rc[i][j] = A[row][col]; 
            j=j+1;
            // Row is filled, so increase row index and 
            // reset col index 
            if (j == n - 1) 
            { 
                j = 0; 
                i=i+1; 
            } 
        } 
    } 
} 
}

double determinant(double A[N_M][N_M], int n) 
{ double D = 0.0; // Initialize result

//  Base case : if matrix contains single element
if (n==1) return A[0][0];
else if (n == 2) return (A[0][0]*A[1][1])-(A[0][1]*A[1][0]); 
else {
double sub_Matrix_A_0c[N_M][N_M]; // To store cofactors 

 // Iterate for each element of first row 
for (int c = 0; c < n; c++) 
{ 
    // Getting Cofactor of A[0][f] 
    getCofactor(A, sub_Matrix_A_0c, 0, c, n); 
    D += pow(-1.0,c) * A[0][c] * determinant(sub_Matrix_A_0c, n - 1); 


} 

return D;} 
}

int main () {
double Z_z[N_M][N_M]=

{{-64,16,-4,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {-43.7213019529827,12.4106746539480,-3.52287874528034,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,-43.7213019529827,12.4106746539480,-3.52287874528034,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,-27,9,-3,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,-27,9,-3,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,-16.0579142269798,6.36491716338729,-2.52287874528034,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,-16.0579142269798,6.36491716338729,-2.52287874528034,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,-8,4,-2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-8,4,-2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3.53179897265895,2.31915967282662,-1.52287874528034,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3.53179897265895,2.31915967282662,-1.52287874528034,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,1,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,1,-1,1,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.142956190020121,0.273402182265940,-0.522878745280338,1,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.142956190020121,0.273402182265940,-0.522878745280338,1,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2.20222621658426,1.69267904961742,1.30102999566398,1}, {37.2320239618439,-7.04575749056068,1,0,-37.2320239618439,7.04575749056068,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,27,-6,1,0,-27,6,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,19.0947514901619,-5.04575749056068,1,0,-19.0947514901619,5.04575749056068,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,12,-4,1,0,-12,4,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6.95747901847985,-3.04575749056068,1,0,-6.95747901847985,3.04575749056068,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,-2,1,0,-3,2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.820206546797821,-1.04575749056068,1,0,-0.820206546797821,1.04575749056068,-1,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,-1,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,2,1,0,-3,-2,-1,0}, {-21.1372724716820,2,0,0,21.1372724716820,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,-18,2,0,0,18,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,-15.1372724716820,2,0,0,15.1372724716820,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,-12,2,0,0,12,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-9.13727247168203,2,0,0,9.13727247168203,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,2,0,0,6,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3.13727247168203,2,0,0,3.13727247168203,-2,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,2,0,0,-6,-2,0,0}, {24,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-7.80617997398389,-2,0,0}};

double det=determinant(Z_z, 40); 
cout<<det;

system ("pause");
return 0;}
  • 1
    What does it mean not to find the determinant? Also, be wary of accumulation of numerical errors. – Javier Silva Ortíz Dec 01 '19 at 06:24
  • it means, it gives me the determinant of 4*4 matrix, but for 40*40, it gives me nothing without error. nothing is printed by the COUT command. – Seyyed Kazem Razavi Dec 01 '19 at 06:39
  • If `A[0,c]` is equal to 0, you don't need to calculate the corresponding submatrix nor the corresponding determinant.This method cannot work for such large matrices, except if you have many zeros. So use this characteristic – Damien Dec 01 '19 at 06:48
  • @Damien I added your suggestion to the code by an If. However, it still does not calculate the determinant. – Seyyed Kazem Razavi Dec 01 '19 at 06:52
  • Your submatrix does not need to be so large. Did you check the global amount of memory that your programme is using ? – Damien Dec 01 '19 at 06:59
  • In the past, I used this method for very large sparse matrices. The issue is that you lost the sparsity characteristic with classical methods. I don't know if your matrix is sparse enough, but this method is not *a priori* stupid for this kind of matrices – Damien Dec 01 '19 at 08:00

1 Answers1

2

You recursively call determinant() functuon n times at the first stage, then n - 1 times for each of n calls, etc. So total number of call would be closed to n! (factorial).

When n = 4 or n = 5 the number of calls is still acceptable, but at n = 40 you try to make 40! = 815915283247897734345611269596115894272000000000 virtual calls to say nothing about so many operations of any kind. I don't think you can find a machine to handle that.

Anatoliy R
  • 1,749
  • 2
  • 14
  • 20
  • No, I just use python to calculate the det, which is -4.637771011488998 – Bill Chen Dec 01 '19 at 06:32
  • @Anatoliy so this explains why it doesn't work, but do you have a suggestion to fix it? – harold Dec 01 '19 at 06:37
  • @BillChen Well, not the worst choice if you want just to find a value for given matrix, but if you want to solve, you probably have to look at the matrix. Of course there are faster methods for large matrices, but in general it will be still slow. But this matrix also has many zeros. – Anatoliy R Dec 01 '19 at 06:39
  • @AnatoliyR, I agree, it is too slow to use a loop. – Bill Chen Dec 01 '19 at 06:40
  • @harold using this approach - no. maybe with another approach when you convert matrix to a different form with the same determinant, also the given matrix is "good" because of so many zeros. – Anatoliy R Dec 01 '19 at 06:42
  • @BillChen could you please suggest me the other methods, I really need to find the detriment of such large matrices in a fast way. Also, I didn't understand the exact error here. – Seyyed Kazem Razavi Dec 01 '19 at 06:43
  • 1
    @SeyyedKazemRazavi, [LU_decomposition](https://en.wikipedia.org/wiki/LU_decomposition#Computing_the_determinant). – Evg Dec 01 '19 at 06:48
  • @SeyyedKazemRazavi there is no error in your code, it's just not possible to run it because modern computers cannot handle that with memory and speed resources they have. Consider to find some library in C++ or other language (BTW python is a good choice) or look at the method Evg proposed. Also look [here](https://en.wikipedia.org/wiki/Gaussian_elimination) – Anatoliy R Dec 01 '19 at 06:56
  • @SeyyedKazemRazavi - even with techniques like LU decomposition, bear in mind that the computation time increases with order of the matrix - complexity of the algorithm is at least `O(n*n)`. – Peter Dec 01 '19 at 06:57
  • You can use [Gaussian elimination](https://en.m.wikipedia.org/wiki/Gaussian_elimination) and then get the determinant by multiplying the element of the main diagonal. – Costantino Grana Dec 01 '19 at 07:27
  • Try this one: https://stackoverflow.com/questions/1886280/how-to-find-determinant-of-large-matrix. I will try to implement myself, if I finished I will post the code in answer. – Bill Chen Dec 01 '19 at 18:47