// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  TU Dresden                                                            ==
// ==                                                                        ==
// ==  Institut f�r Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
// ==                                                                        ==
// ============================================================================

/** \file Cholesky.h */

#ifndef CHOLESKY_H
#define CHOLESKY_H

#include "FixVec.h"

using namespace std;
using namespace AMDiS;

/**
 * \ingroup Solvers
 *
 * \brief
 * Solves a symmetric positive definite system using Cholesky decomposition.
 * Returns true if the matrix in (numerically) positive definite.
 */

class Cholesky
{
 public:
  /** \brief
   * Computes the Cholesky factor of a positive definite matrix A.
   * On input, only the upper triangle of A need be given; it is not modified.
   * The Cholesky factor is returned in the lower triangle of A, except for its
   * diagonal elements which are returned in P.
   */
  static bool factorization(Matrix<double> *A, Vector<double> *p);

  /** \brief
   * Solves system A*X=B, where A is a positive definite matrix.
   * If P=NULL; A is assumed to be positive definite, and a Cholesky
   * decomposition is computed using the previous routine.
   * If P is given, A and P are assumed to be already given as the output of
   * the previous routine.
   */
  static bool solve(Matrix<double> *A, Vector<double> *b, Vector<double> *x,
		    Vector<double> *p = NULL);
  static bool solve(Matrix<double> *A, Vector<WorldVector<double> > *b,
		    Vector<WorldVector<double> > *x,
		    Vector<double> *p = NULL);
};

#endif