/** @(#) factorizecholesky.cpp */
//#include "factorizecholesky.h"

/** @(#) factorizecholesky.h */
#ifndef FACTORIZECHOLESKY_H
#define FACTORIZECHOLESKY_H
//----------------------------------------------------------------------------
#include "../matrixdouble/squaredmatrix.h"
#include <iostream>
#include <stdexcept>	// runtime_exception
using namespace std;

/**
* Obtiene la factorizacion de Cholesky.
* Logica uno.
* @param A Matriz cuadrada simetrica a factorizar.
* No se requiere U, porque U = Lt
* Matriz positiva definida. det(A) = 324
*	| 1	2	4  |
* A =	| 3	13	23 |
*	| 4	23	77 |
*
*	| 1	0	0 |
* L =	| 2	3	0 |
*	| 4	5	6 |
*
*       | 60 30 20 |
*  A =  | 30 20 15 |   is positive definite, non symmetric
*       | 20 15 12 |
*
*        | -60  30  20 |
*  A' =  |  30 -20  15 |   is not positive definite
*        |  20  15 -12 |
*
* @param L return Factor en forma de matriz triangular inferior.
* @author Daniel Alba Cuellar, Omar Posada Villarreal
* @return 1=true si la factorización tiene éxito, 0=false en caso contrario
* @version 1.0, 02/03/2002
*/
bool factorizeCholesky(int n, SquaredMatrix &A, SquaredMatrix &L)
{
   // { posada
   A.setLogic(1);
   L.setLogic(1);
   // } posada

   int k, s, i;
   double sum;

   for (k = 1; k <= n; k++)
   {
      sum = 0.0;
      for (s = 1; s <= k-1; s++)
         sum += L[k][s] * L[k][s];

      if (A[k][k]-sum < 0.0) {
         cout << "factorizeCholesky failed. Matrix is not positive definite.";
         return false;
      }

      L[k][k] = sqrt( A[k][k] - sum );

      for (i = k+1; i <= n; i++)
      {
         sum = 0.0;
         for (s = 1; s <= k-1; s++)
            sum += L[i][s] * L[k][s];

         L[i][k] = (A[i][k] - sum)/L[k][k];
      }
   }
   return true;
}

/** Checa parametros: dimensiones, rangos, divisiones entre cero */
bool checkFactorizeCholesky(int n, SquaredMatrix &A, SquaredMatrix &L)
                throw(invalid_argument) {
	// Condiciones invalidas
        if ( (n < 1)
        	|| (A.getSize() != n)
        	|| (L.getSize() != n)
        	) {
        	throw invalid_argument("factorizeCholesky: arg invalidos.");
        }
	return factorizeCholesky(n, A, L);
}

//----------------------------------------------------------------------------
#endif
// Fin------------------------------------------------------------------------