#include "thinfilm_BLS.h"

void tf_det(double **M,int dimension,Huge_real *det)
{
  int i, j, k, h, pj, pk;
  int c[2*tf_MAXL +1], r[2*tf_MAXL +1];       /* index for column and row */
  double pivot, pt, Sign;
  double determinant,scale;
  double **m;

  m = matrix(1, dimension, 1, dimension);
  cp_matrix(m, M, dimension);

  pj=pk=0;

//  printf("\ndimension=%d\n",dimension);  
   
  for (i = 1; i <=dimension; i++)
    c[i] = r[i] = i;
  Sign = 1;       

    /* Apply gauss algorithm to m */


  for (i = 1; i <=(dimension - 1); i++)
    {
      /* Search pivot */
      pivot = 0;
      for (j = i; j <=dimension; j++)
        {
          for (k = i; k <=dimension; k++)
            {
              pt = fabs (m[r[j]][c[k]]);
              if (pt > pivot)
                {
                  pivot = pt;
                  pj = j;
                  pk = k;
                }
            }
        }
      if (pivot == 0)
        {
          fprintf (stderr,
                   "\nwarning: matrix is degenerate!\n\
this shouldnt happen often -- otherwise check your settings!\n");
          det->log = -HUGE_VAL;
          det->sign = 0;
          return ;   /* Matrix is degenerate */
        }
     /* Pivot to left upper corner (i,i) of m */
      if (pj != i)
        {
          Sign *= -1;
          h = r[i];
          r[i] = r[pj];
          r[pj] = h;
        }
      if (pk != i)
        {
          Sign *= -1;
          h = c[i];
          c[i] = c[pk];
          c[pk] = h;
        }

      /* Add row i to each lower */
      for (j = i + 1; j <=dimension; j++)
        {
          scale = m[r[j]][c[i]] / m[r[i]][c[i]];
          if (!(scale == 0))
            {
              for (k = i + 1; k <=dimension; k++)
                {
                  m[r[j]][c[k]] = m[r[j]][c[k]] - m[r[i]][c[k]] * scale;
                }
            }
        }
    }

  /* Calculate determinant */

  determinant = 0;
//  printf("\ncalculating det\n");
  for (i = dimension; i >= 1; i--)
    {
        determinant += log(fabs (m[r[i]][c[i]]));
        Sign *= Sign(m[r[i]][c[i]]);         
    }
    
  det->log = determinant;
  det->sign = Sign;

//  printf("\t %e",determinant\n); 

 free_matrix(m,1,dimension,1,dimension);

  return ;
}

void cp_matrix(double **dest, double **source,int N)
{
int i,j;
for(i = 1;i<=N;i++)
	for(j = 1;j<=N;j++)  
	dest[i][j] = source[i][j];

}

