//==============================================================================
//
// Calculation of Determinant of square Matrix
//
// Developer: Henry Guennadi Levkin
//
//==============================================================================

#include <stdio.h>

//==============================================================================
void MatrixPrint(int nDim, double* pfMatr)
{
  int i,j;

  printf("------------------------------------------------------------------\n");
  for(i=0; i<nDim; i++)
  {
    for(j=0; j<nDim; j++)
    {
      printf("%9.2lf ", pfMatr[nDim*i + j]);
    }
    printf("\n");
  }
}

//==============================================================================
// function is not optimized for speed
// matrix becomes triangle after function call
//==============================================================================

#include <math.h>

double MatrixDeterminant(int nDim, double* pfMatr)
{
  double fDet = 1.;
  double fMaxElem;
  double fAcc;

  int i , j, k, m;


  for(k=0; k<(nDim-1); k++) // base row of matrix
  {
    // search of line with max element
    fMaxElem = fabs( pfMatr[k*nDim + k] );
    m = k;
    for(i=k+1; i<nDim; i++)
    {
      if(fMaxElem < fabs(pfMatr[i*nDim + k]) )
      {
        fMaxElem = pfMatr[i*nDim + k];
        m = i;
      }
    }
    
    // permutation of base line (index k) and max element line(index m)
    if(m != k)
    {
      for(i=k; i<nDim; i++)
      {
        fAcc               = pfMatr[k*nDim + i];
        pfMatr[k*nDim + i] = pfMatr[m*nDim + i];
        pfMatr[m*nDim + i] = fAcc;
      }
      fDet *= (-1.);
    }

    if( pfMatr[k*nDim + k] == 0.) return 0.0;

    // trianglulation of matrix
    for(j=(k+1); j<nDim; j++) // current row of matrix
    {
      fAcc = - pfMatr[j*nDim + k] / pfMatr[k*nDim + k];
      for(i=k; i<nDim; i++)
      {
        pfMatr[j*nDim + i] = pfMatr[j*nDim + i] + fAcc*pfMatr[k*nDim + i];
      }
    }
  }

  for(i=0; i<nDim; i++)
  {
    fDet *= pfMatr[i*nDim + i]; // diagonal elements multiplication
  }

  return fDet;
}

//==============================================================================
// testing of function MatrixDeterminant()
//==============================================================================

#define MATRIX_DIMENSION 3

int main(int nArgs, char** pArgs)
{
  int nDim = MATRIX_DIMENSION;
  double fMatr[MATRIX_DIMENSION*MATRIX_DIMENSION] =
  {
  1.0, 2.0, 3.0,
  4.0, 5.0, 7.0,
  6.0, 8.0, 9.0,
  };

  double fDeterm;

  MatrixPrint(nDim, fMatr);
  
  fDeterm = MatrixDeterminant(nDim, fMatr);
  
  printf("\nDeterminant = %lf\n", fDeterm);

  return 0;
}