/*************************************************************************
Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).

>>> SOURCE LICENSE >>>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation (www.fsf.org); either version 2 of the 
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

A copy of the GNU General Public License is available at
http://www.fsf.org/licensing/licenses

>>> END OF LICENSE >>>
*************************************************************************/


namespace alglib {
  public class sdet {
    /*************************************************************************
    Determinant calculation of the matrix given by LDLT decomposition.

    Input parameters:
        A       -   LDLT-decomposition of the matrix,
                    output of subroutine SMatrixLDLT.
        Pivots  -   table of permutations which were made during
                    LDLT decomposition, output of subroutine SMatrixLDLT.
        N       -   size of matrix A.
        IsUpper -   matrix storage format. The value is equal to the input
                    parameter of subroutine SMatrixLDLT.

    Result:
        matrix determinant.

      -- ALGLIB --
         Copyright 2005-2008 by Bochkanov Sergey
    *************************************************************************/
    public static double smatrixldltdet(ref double[,] a,
        ref int[] pivots,
        int n,
        bool isupper) {
      double result = 0;
      int k = 0;

      result = 1;
      if (isupper) {
        k = 0;
        while (k < n) {
          if (pivots[k] >= 0) {
            result = result * a[k, k];
            k = k + 1;
          } else {
            result = result * (a[k, k] * a[k + 1, k + 1] - a[k, k + 1] * a[k, k + 1]);
            k = k + 2;
          }
        }
      } else {
        k = n - 1;
        while (k >= 0) {
          if (pivots[k] >= 0) {
            result = result * a[k, k];
            k = k - 1;
          } else {
            result = result * (a[k - 1, k - 1] * a[k, k] - a[k, k - 1] * a[k, k - 1]);
            k = k - 2;
          }
        }
      }
      return result;
    }


    /*************************************************************************
    Determinant calculation of the symmetric matrix

    Input parameters:
        A       -   matrix. Array with elements [0..N-1, 0..N-1].
        N       -   size of matrix A.
        IsUpper -   if IsUpper = True, then symmetric matrix A is given by its
                    upper triangle, and the lower triangle isn’t used by
                    subroutine. Similarly, if IsUpper = False, then A is given
                    by its lower triangle.

    Result:
        determinant of matrix A.

      -- ALGLIB --
         Copyright 2005-2008 by Bochkanov Sergey
    *************************************************************************/
    public static double smatrixdet(double[,] a,
        int n,
        bool isupper) {
      double result = 0;
      int[] pivots = new int[0];

      a = (double[,])a.Clone();

      ldlt.smatrixldlt(ref a, n, isupper, ref pivots);
      result = smatrixldltdet(ref a, ref pivots, n, isupper);
      return result;
    }


    public static double determinantldlt(ref double[,] a,
        ref int[] pivots,
        int n,
        bool isupper) {
      double result = 0;
      int k = 0;

      result = 1;
      if (isupper) {
        k = 1;
        while (k <= n) {
          if (pivots[k] > 0) {
            result = result * a[k, k];
            k = k + 1;
          } else {
            result = result * (a[k, k] * a[k + 1, k + 1] - a[k, k + 1] * a[k, k + 1]);
            k = k + 2;
          }
        }
      } else {
        k = n;
        while (k >= 1) {
          if (pivots[k] > 0) {
            result = result * a[k, k];
            k = k - 1;
          } else {
            result = result * (a[k - 1, k - 1] * a[k, k] - a[k, k - 1] * a[k, k - 1]);
            k = k - 2;
          }
        }
      }
      return result;
    }


    public static double determinantsymmetric(double[,] a,
        int n,
        bool isupper) {
      double result = 0;
      int[] pivots = new int[0];

      a = (double[,])a.Clone();

      ldlt.ldltdecomposition(ref a, n, isupper, ref pivots);
      result = determinantldlt(ref a, ref pivots, n, isupper);
      return result;
    }
  }
}