// @(#)root/smatrix:$Id$
// Authors: T. Glebe, L. Moneta    2005  

#ifndef ROOT_Math_Dfactir
#define ROOT_Math_Dfactir
// ********************************************************************
//
// source:
//
// type:      source code
//
// created:   02. Apr 2001
//
// author:    Thorsten Glebe
//            HERA-B Collaboration
//            Max-Planck-Institut fuer Kernphysik
//            Saupfercheckweg 1
//            69117 Heidelberg
//            Germany
//            E-mail: T.Glebe@mpi-hd.mpg.de
//
// Description: Determinant of a square matrix, needed for Dfinv()
//              Code was taken from CERNLIB::kernlib dfact function, translated
//              from FORTRAN to C++ and optimized.
//
// changes:
// 02 Apr 2001 (TG) creation
//
// ********************************************************************

#include <cmath>

namespace ROOT { 

  namespace Math { 


/** Dfactir.
    Function to compute the determinant from a square matrix, Det(A) of
    dimension idim and order n. A working area ir is returned which is
    needed by the Dfinv routine.

    @author T. Glebe
*/
template <class Matrix, unsigned int n, unsigned int idim>
bool Dfactir(Matrix& rhs, typename Matrix::value_type& det, unsigned int* ir)
  // int *n, float *a, int *idim, int *ir, int *ifail, float *det, int *jfail)
{

#ifdef XXX
  if (idim < n || n <= 0) {
    return false;
  }
#endif


  /* Initialized data */
  typename Matrix::value_type* a = rhs.Array();

  /* Local variables */
  static unsigned int nxch, i, j, k, l;
  static typename Matrix::value_type p, q, tf;
  
  /* Parameter adjustments */
  a -= idim + 1;
  --ir;

  /* Function Body */
  
  // fact.inc
  nxch = 0;
  det = 1.;
  for (j = 1; j <= n; ++j) {
    const unsigned int ji = j * idim;
    const unsigned int jj = j + ji;

    k = j;
    p = std::abs(a[jj]);

    if (j != n) {
      for (i = j + 1; i <= n; ++i) {
	q = std::abs(a[i + ji]);
	if (q > p) {
	  k = i;
	  p = q;
	}
      } // for i

      if (k != j) {
	for (l = 1; l <= n; ++l) {
	  const unsigned int li = l*idim;
	  const unsigned int jli = j + li;
	  const unsigned int kli = k + li;
	  tf = a[jli];
	  a[jli] = a[kli];
	  a[kli] = tf;
	} // for l
	++nxch;
	ir[nxch] = (j << 12) + k;
      } // if k != j
    } // if j!=n

    if (p <= 0.) {
      det = 0;
      return false;
    }

    det *= a[jj];
#ifdef XXX
    t = std::abs(det);
    if (t < 1e-19 || t > 1e19) {
      det = 0;
      return false;
    }
#endif

    a[jj] = 1. / a[jj];
    if (j == n) {
      continue;
    }

    const unsigned int jm1 = j - 1;
    const unsigned int jpi = (j + 1) * idim;
    const unsigned int jjpi = j + jpi;

    for (k = j + 1; k <= n; ++k) {
      const unsigned int ki  = k * idim;
      const unsigned int jki = j + ki;
      const unsigned int kji = k + jpi;
      if (j != 1) {
	for (i = 1; i <= jm1; ++i) {
	  const unsigned int ii = i * idim;
	  a[jki] -= a[i + ki] * a[j + ii];
	  a[kji] -= a[i + jpi] * a[k + ii];
	} // for i
      }
      a[jki] *= a[jj];
      a[kji] -= a[jjpi] * a[k + ji];
    } // for k
  } // for j

  if (nxch % 2 != 0) {
    det = -(det);
  }
  ir[n] = nxch;
  return true;
} // end of Dfact


  }  // namespace Math

}  // namespace ROOT
          


#endif /* ROOT_Math_Dfactir */