/branches/v1.6.5/Classes/PHPExcel/Shared/JAMA/CholeskyDecomposition.php

<?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

  */

  var $L = array();

  

  /**

  * Matrix row and column dimension

  * @var int

  * @access private

  */

  var $m;

  

  /**

  * Symmetric positive definite flag

  * @var boolean

  * @access private

  */

  var $isspd = true;

  

  /**

  * CholeskyDecomposition

  * Class constructor - decomposes symmetric positive definite matrix

  * @param mixed Matrix square symmetric positive definite matrix

  */

  function CholeskyDecomposition( $A = null ) {

    if( is_a($A, '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 {

      trigger_error(ArgumentTypeException, ERROR);

    }

  }

  

  /** 

  * Is the matrix symmetric and positive definite?

  * @return boolean

  */

  function isSPD () {

    return $this->isspd;

  }

  

  /** 

  * getL

  * Return triangular factor.

  * @return Matrix Lower triangular matrix

  */

  function getL () {

    return new Matrix($this->L);

  }

  

  /** 

  * Solve A*X = B

  * @param $B Row-equal matrix

  * @return Matrix L * L' * X = B

  */

  function solve ( $B = null ) {

    if( is_a($B, '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 {

          trigger_error(MatrixSPDException, ERROR);

        }

      } else {

        trigger_error(MatrixDimensionException, ERROR);

      }

    } else {

      trigger_error(ArgumentTypeException, ERROR);

    }

  }

}



?>