You are here

CholeskyDecomposition.php in Loft Data Grids 6.2

File

vendor/phpoffice/phpexcel/Classes/PHPExcel/Shared/JAMA/CholeskyDecomposition.php
View source
<?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(ArgumentTypeException));
    }
  }

  //	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(MatrixDimensionException));
      }
    }
    else {
      throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException));
    }
  }

}

//	class CholeskyDecomposition

Classes

Namesort descending Description
CholeskyDecomposition @package JAMA