/*---------------------------------------------------------------------------*
 *                                   IT++			             *
 *---------------------------------------------------------------------------*
 * Copyright (c) 1995-2004 by Tony Ottosson, Thomas Eriksson, Pål Frenger,   *
 * Tobias Ringström, and Jonas Samuelsson.                                   *
 *                                                                           *
 * Permission to use, copy, modify, and distribute this software and its     *
 * documentation under the terms of the GNU General Public License is hereby *
 * granted. No representations are made about the suitability of this        *
 * software for any purpose. It is provided "as is" without expressed or     *
 * implied warranty. See the GNU General Public License for more details.    *
 *---------------------------------------------------------------------------*/

/*!
  \file
  \brief Implementation of Cholesky factorisation functions
  \author Tony Ottosson

  1.9

  2004/06/16 15:48:33
*/

#include <algorithm>
#include <cassert>
#include "itconfig.h"
#include "base/cholesky.h"
#include "../src/base/lapack.h"

namespace itpp { 

#ifndef NO_LAPACK

  bool chol(const mat &X, mat &F)
  {
    char uplo='U';
    int n, lda, info;
    n = lda = X.rows();

    F = X; // input matrix is overwritten

    dpotrf_(&uplo, &n, F._data(), &lda, &info);

    // Set lower part to zero
    for (int i=0; i<n; i++)
      for(int j=i+1; j<n; j++)
	F(j,i) = 0;

    return (info==0);
  }

  bool chol(const cmat &X, cmat &F)
  {
    char uplo='U';
    int n, lda, info;
    n = lda = X.rows();

    F = X; // input matrix is overwritten

    zpotrf_(&uplo, &n, F._data(), &lda, &info);

    // Set lower part to zero
    for (int i=0; i<n; i++)
      for(int j=i+1; j<n; j++)
	F(j,i) = 0;

    return (info==0);
  }

  cmat chol(const cmat &X)
  {
    cmat F;
    if (!chol(X, F)) {
      it_warning("cholesky factorization didn't succeed");
    }

    return F;
  }


  mat chol(const mat &X)
  {
    mat F;
    if (!chol(X, F)) {
      it_warning("cholesky factorization didn't succeed");
    }

    return F;
  }

#endif //NO_LAPACK

} //namespace itpp