Determinant of a matrix by Gaussian elimination C++

Question

I was trying to find the code for finding the determinant of a square matrix , and I came across this code.

    int det(vector<vector<int> > mat) {
        int n = mat.size();


        for(int col = 0; col < n; ++col) {
            bool found = false;
            for(int row = col; row < n; ++row) {
                if(mat[row][col]) {
                    mat[row].swap(mat[col]);
                    found = true;
                    break;
                }
            }
            if(!found) {
                return 0;
            }
            for(int row = col + 1; row < n; ++row) {
                while(true) {
                    int del = mat[row][col] / mat[col][col];
                    for (int j = col; j < n; ++j) {
                        mat[row][j] -= del * mat[col][j];
                    }
                    if (mat[row][col] == 0)
                        break;
                    else
                        mat[row].swap(mat[col]);
                }
            }
        }

        li res = 1;

        for(int i = 0; i < n; ++i) {
            res *= mat[i][i];
        }
        return abs(res);
    }

But I am having trouble understanding line 20-29 i.e where the subtraction of row from multiple of another row is performed. I mean why the while loop is required here ? as I am subtracting the quotient*dividend , It should always be 0 , right ? So I think it should be just one iteration. So , why we need to perform this mat[row].swap(mat[col]); operation ? Thanks in advance.

Answer 1

There is some strange logic in your code to account for the fact that you are performing your calculations using integer arithmetic.

Say you have a 3x3 matrix in which the first two rows are:

4 6 5
1 2 3

When you compute del for col=0 and row=1, you will get:

del = 1/4 = 0

With that, when you compute:

mat[row][j] -= del * mat[col][j];

mat[row][j] doesn't change at all.

To account for that, you swap the rows. Now the first two rows are:

1 2 3
4 6 5

With the rows swapped like that, the value of del is 4/1 = 4. Now the line:

mat[row][j] -= del * mat[col][j];

does make a difference. The value of mat[1][0] ends up being zero, which is what you need. So you break out of the while loop.

Here's an instrumented version of your function that produces a lot of debug output, with a helper function to print the matrix and the main function to test the code.

#include <iostream>
#include <vector>
#include <stdlib.h>

using namespace std;

void printMatrix(vector<vector<int> > const& mat)
{
   int n = mat.size();
   for(int row = 0; row < n; ++row) {
      for(int col = 0; col < n; ++col) {
         cout << mat[row][col] << " ";
      }
      cout << "\n";
   }
   cout << "\n";
}

int det(vector<vector<int> > mat) {
   int n = mat.size();

   for(int col = 0; col < n; ++col) {
      cout << "Column: " << col << "\n";
      printMatrix(mat);
      bool found = false;
      for(int row = col; row < n; ++row) {
         if(mat[row][col]) {
            cout << "Got non-zero value for row " << row << " and col " << col << "\n";
            if ( row != col )
            {
               cout << "(1) Swapping rows " << col << " and " << row << "\n";
               mat[row].swap(mat[col]);
               printMatrix(mat);
            }
            else
            {
               cout << "Not swapping rows\n";
            }
            found = true;
            break;
         }
      }

      if(!found) {
         cout << "Did not find a non-zero row. Column: " << col << "\n";
         return 0;
      }

      for(int row = col + 1; row < n; ++row) {
         while(true) {
            int del = mat[row][col] / mat[col][col];
            cout << "del: " << del << "\n";
            for (int j = col; j < n; ++j) {
               mat[row][j] -= del * mat[col][j];
            }
            if (mat[row][col] == 0)
            {
               break;
            }
            else
            {
               cout << "(2) Swapping rows " << col << " and " << row << "\n";
               mat[row].swap(mat[col]);
               printMatrix(mat);
            }
         }
      }
   }

   printMatrix(mat);
   long res = 1;

   for(int i = 0; i < n; ++i) {
      res *= mat[i][i];
   }
   return abs(res);
}

int main()
{
   vector<vector<int> > mat = { {4, 6, 5}, {1, 2, 3}, {8, 10, 9} };
   int r = det(mat);
   cout << "Determinant: " << r << endl;
   return 0;
}