<?php
/**
 *    @package JAMA
 *
 *    Cholesky decomposition class
 *
 *    For a symmetric, positive definite matrix A, the Cholesky decomposition
 *    is an lower triangular matrix L so that A = L*L'.
 *
 *    If the matrix is not symmetric or positive definite, the constructor
 *    returns a partial decomposition and sets an internal flag that may
 *    be queried by the isSPD() method.
 *
 *    @author Paul Meagher
 *    @author Michael Bommarito
 *    @version 1.2
 */
class CholeskyDecomposition
{
  /**
   * Decomposition storage
   * 
   * @var array
   * @access private
   */
  private $L = array ();
  
  /**
   * Matrix row and column dimension
   * 
   * @var int
   * @access private
   */
  private $m;
  
  /**
   * Symmetric positive definite flag
   * 
   * @var boolean
   * @access private
   */
  private $isspd = true;
  
  /**
   * CholeskyDecomposition
   *
   * Class constructor - decomposes symmetric positive definite matrix
   * 
   * @param
   *          mixed Matrix square symmetric positive definite matrix
   */
  public function __construct($A = null)
  {
    if ($A instanceof Matrix)
    {
      $this->L = $A->getArray ();
      $this->m = $A->getRowDimension ();
      
      for($i = 0; $i < $this->m; ++ $i)
      {
        for($j = $i; $j < $this->m; ++ $j)
        {
          for($sum = $this->L [$i] [$j], $k = $i - 1; $k >= 0; -- $k)
          {
            $sum -= $this->L [$i] [$k] * $this->L [$j] [$k];
          }
          if ($i == $j)
          {
            if ($sum >= 0)
            {
              $this->L [$i] [$i] = sqrt ( $sum );
            }
            else
            {
              $this->isspd = false;
            }
          }
          else
          {
            if ($this->L [$i] [$i] != 0)
            {
              $this->L [$j] [$i] = $sum / $this->L [$i] [$i];
            }
          }
        }
        
        for($k = $i + 1; $k < $this->m; ++ $k)
        {
          $this->L [$i] [$k] = 0.0;
        }
      }
    }
    else
    {
      throw new PHPExcel_Calculation_Exception ( JAMAError ( ARGUMENT_TYPE_EXCEPTION ) );
    }
  } // function __construct()
  
  /**
   * Is the matrix symmetric and positive definite?
   *
   * @return boolean
   */
  public function isSPD()
  {
    return $this->isspd;
  } // function isSPD()
  
  /**
   * getL
   *
   * Return triangular factor.
   * 
   * @return Matrix Lower triangular matrix
   */
  public function getL()
  {
    return new Matrix ( $this->L );
  } // function getL()
  
  /**
   * Solve A*X = B
   *
   * @param $B Row-equal
   *          matrix
   * @return Matrix L * L' * X = B
   */
  public function solve($B = null)
  {
    if ($B instanceof Matrix)
    {
      if ($B->getRowDimension () == $this->m)
      {
        if ($this->isspd)
        {
          $X = $B->getArrayCopy ();
          $nx = $B->getColumnDimension ();
          
          for($k = 0; $k < $this->m; ++ $k)
          {
            for($i = $k + 1; $i < $this->m; ++ $i)
            {
              for($j = 0; $j < $nx; ++ $j)
              {
                $X [$i] [$j] -= $X [$k] [$j] * $this->L [$i] [$k];
              }
            }
            for($j = 0; $j < $nx; ++ $j)
            {
              $X [$k] [$j] /= $this->L [$k] [$k];
            }
          }
          
          for($k = $this->m - 1; $k >= 0; -- $k)
          {
            for($j = 0; $j < $nx; ++ $j)
            {
              $X [$k] [$j] /= $this->L [$k] [$k];
            }
            for($i = 0; $i < $k; ++ $i)
            {
              for($j = 0; $j < $nx; ++ $j)
              {
                $X [$i] [$j] -= $X [$k] [$j] * $this->L [$k] [$i];
              }
            }
          }
          
          return new Matrix ( $X, $this->m, $nx );
        }
        else
        {
          throw new PHPExcel_Calculation_Exception ( JAMAError ( MatrixSPDException ) );
        }
      }
      else
      {
        throw new PHPExcel_Calculation_Exception ( JAMAError ( MATRIX_DIMENSION_EXCEPTION ) );
      }
    }
    else
    {
      throw new PHPExcel_Calculation_Exception ( JAMAError ( ARGUMENT_TYPE_EXCEPTION ) );
    }
  } // function solve()
}
