1: <?php
  2: /**
  3:  *  @package JAMA
  4:  *
  5:  *  Cholesky decomposition class
  6:  *
  7:  *  For a symmetric, positive definite matrix A, the Cholesky decomposition
  8:  *  is an lower triangular matrix L so that A = L*L'.
  9:  *
 10:  *  If the matrix is not symmetric or positive definite, the constructor
 11:  *  returns a partial decomposition and sets an internal flag that may
 12:  *  be queried by the isSPD() method.
 13:  *
 14:  *  @author Paul Meagher
 15:  *  @author Michael Bommarito
 16:  *  @version 1.2
 17:  */
 18: class CholeskyDecomposition {
 19: 
 20:     /**
 21:      *  Decomposition storage
 22:      *  @var array
 23:      *  @access private
 24:      */
 25:     private $L = array();
 26: 
 27:     /**
 28:      *  Matrix row and column dimension
 29:      *  @var int
 30:      *  @access private
 31:      */
 32:     private $m;
 33: 
 34:     /**
 35:      *  Symmetric positive definite flag
 36:      *  @var boolean
 37:      *  @access private
 38:      */
 39:     private $isspd = true;
 40: 
 41: 
 42:     /**
 43:      *  CholeskyDecomposition
 44:      *
 45:      *  Class constructor - decomposes symmetric positive definite matrix
 46:      *  @param mixed Matrix square symmetric positive definite matrix
 47:      */
 48:     public function __construct($A = null) {
 49:         if ($A instanceof Matrix) {
 50:             $this->L = $A->getArray();
 51:             $this->m = $A->getRowDimension();
 52: 
 53:             for($i = 0; $i < $this->m; ++$i) {
 54:                 for($j = $i; $j < $this->m; ++$j) {
 55:                     for($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) {
 56:                         $sum -= $this->L[$i][$k] * $this->L[$j][$k];
 57:                     }
 58:                     if ($i == $j) {
 59:                         if ($sum >= 0) {
 60:                             $this->L[$i][$i] = sqrt($sum);
 61:                         } else {
 62:                             $this->isspd = false;
 63:                         }
 64:                     } else {
 65:                         if ($this->L[$i][$i] != 0) {
 66:                             $this->L[$j][$i] = $sum / $this->L[$i][$i];
 67:                         }
 68:                     }
 69:                 }
 70: 
 71:                 for ($k = $i+1; $k < $this->m; ++$k) {
 72:                     $this->L[$i][$k] = 0.0;
 73:                 }
 74:             }
 75:         } else {
 76:             throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException));
 77:         }
 78:     }   //  function __construct()
 79: 
 80: 
 81:     /**
 82:      *  Is the matrix symmetric and positive definite?
 83:      *
 84:      *  @return boolean
 85:      */
 86:     public function isSPD() {
 87:         return $this->isspd;
 88:     }   //  function isSPD()
 89: 
 90: 
 91:     /**
 92:      *  getL
 93:      *
 94:      *  Return triangular factor.
 95:      *  @return Matrix Lower triangular matrix
 96:      */
 97:     public function getL() {
 98:         return new Matrix($this->L);
 99:     }   //  function getL()
100: 
101: 
102:     /**
103:      *  Solve A*X = B
104:      *
105:      *  @param $B Row-equal matrix
106:      *  @return Matrix L * L' * X = B
107:      */
108:     public function solve($B = null) {
109:         if ($B instanceof Matrix) {
110:             if ($B->getRowDimension() == $this->m) {
111:                 if ($this->isspd) {
112:                     $X  = $B->getArrayCopy();
113:                     $nx = $B->getColumnDimension();
114: 
115:                     for ($k = 0; $k < $this->m; ++$k) {
116:                         for ($i = $k + 1; $i < $this->m; ++$i) {
117:                             for ($j = 0; $j < $nx; ++$j) {
118:                                 $X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k];
119:                             }
120:                         }
121:                         for ($j = 0; $j < $nx; ++$j) {
122:                             $X[$k][$j] /= $this->L[$k][$k];
123:                         }
124:                     }
125: 
126:                     for ($k = $this->m - 1; $k >= 0; --$k) {
127:                         for ($j = 0; $j < $nx; ++$j) {
128:                             $X[$k][$j] /= $this->L[$k][$k];
129:                         }
130:                         for ($i = 0; $i < $k; ++$i) {
131:                             for ($j = 0; $j < $nx; ++$j) {
132:                                 $X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i];
133:                             }
134:                         }
135:                     }
136: 
137:                     return new Matrix($X, $this->m, $nx);
138:                 } else {
139:                     throw new PHPExcel_Calculation_Exception(JAMAError(MatrixSPDException));
140:                 }
141:             } else {
142:                 throw new PHPExcel_Calculation_Exception(JAMAError(MatrixDimensionException));
143:             }
144:         } else {
145:             throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException));
146:         }
147:     }   //  function solve()
148: 
149: }   //  class CholeskyDecomposition
150: