[Summary view] [Print] [Text view]

   1  <?php
   2  /**
   3   *    @package JAMA
   4   *
   5   *    For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
   6   *    unit lower triangular matrix L, an n-by-n upper triangular matrix U,
   7   *    and a permutation vector piv of length m so that A(piv,:) = L*U.
   8   *    If m < n, then L is m-by-m and U is m-by-n.
   9   *
  10   *    The LU decompostion with pivoting always exists, even if the matrix is
  11   *    singular, so the constructor will never fail. The primary use of the
  12   *    LU decomposition is in the solution of square systems of simultaneous
  13   *    linear equations. This will fail if isNonsingular() returns false.
  14   *
  15   *    @author Paul Meagher
  16   *    @author Bartosz Matosiuk
  17   *    @author Michael Bommarito
  18   *    @version 1.1
  19   *    @license PHP v3.0
  20   */
  21  class PHPExcel_Shared_JAMA_LUDecomposition
  22  {
  23      const MATRIX_SINGULAR_EXCEPTION    = "Can only perform operation on singular matrix.";
  24      const MATRIX_SQUARE_EXCEPTION      = "Mismatched Row dimension";
  25  
  26      /**
  27       *    Decomposition storage
  28       *    @var array
  29       */
  30      private $LU = array();
  31  
  32      /**
  33       *    Row dimension.
  34       *    @var int
  35       */
  36      private $m;
  37  
  38      /**
  39       *    Column dimension.
  40       *    @var int
  41       */
  42      private $n;
  43  
  44      /**
  45       *    Pivot sign.
  46       *    @var int
  47       */
  48      private $pivsign;
  49  
  50      /**
  51       *    Internal storage of pivot vector.
  52       *    @var array
  53       */
  54      private $piv = array();
  55  
  56      /**
  57       *    LU Decomposition constructor.
  58       *
  59       *    @param $A Rectangular matrix
  60       *    @return Structure to access L, U and piv.
  61       */
  62      public function __construct($A)
  63      {
  64          if ($A instanceof PHPExcel_Shared_JAMA_Matrix) {
  65              // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
  66              $this->LU = $A->getArray();
  67              $this->m  = $A->getRowDimension();
  68              $this->n  = $A->getColumnDimension();
  69              for ($i = 0; $i < $this->m; ++$i) {
  70                  $this->piv[$i] = $i;
  71              }
  72              $this->pivsign = 1;
  73              $LUrowi = $LUcolj = array();
  74  
  75              // Outer loop.
  76              for ($j = 0; $j < $this->n; ++$j) {
  77                  // Make a copy of the j-th column to localize references.
  78                  for ($i = 0; $i < $this->m; ++$i) {
  79                      $LUcolj[$i] = &$this->LU[$i][$j];
  80                  }
  81                  // Apply previous transformations.
  82                  for ($i = 0; $i < $this->m; ++$i) {
  83                      $LUrowi = $this->LU[$i];
  84                      // Most of the time is spent in the following dot product.
  85                      $kmax = min($i, $j);
  86                      $s = 0.0;
  87                      for ($k = 0; $k < $kmax; ++$k) {
  88                          $s += $LUrowi[$k] * $LUcolj[$k];
  89                      }
  90                      $LUrowi[$j] = $LUcolj[$i] -= $s;
  91                  }
  92                  // Find pivot and exchange if necessary.
  93                  $p = $j;
  94                  for ($i = $j+1; $i < $this->m; ++$i) {
  95                      if (abs($LUcolj[$i]) > abs($LUcolj[$p])) {
  96                          $p = $i;
  97                      }
  98                  }
  99                  if ($p != $j) {
 100                      for ($k = 0; $k < $this->n; ++$k) {
 101                          $t = $this->LU[$p][$k];
 102                          $this->LU[$p][$k] = $this->LU[$j][$k];
 103                          $this->LU[$j][$k] = $t;
 104                      }
 105                      $k = $this->piv[$p];
 106                      $this->piv[$p] = $this->piv[$j];
 107                      $this->piv[$j] = $k;
 108                      $this->pivsign = $this->pivsign * -1;
 109                  }
 110                  // Compute multipliers.
 111                  if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
 112                      for ($i = $j+1; $i < $this->m; ++$i) {
 113                          $this->LU[$i][$j] /= $this->LU[$j][$j];
 114                      }
 115                  }
 116              }
 117          } else {
 118              throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ARGUMENT_TYPE_EXCEPTION);
 119          }
 120      }    //    function __construct()
 121  
 122      /**
 123       *    Get lower triangular factor.
 124       *
 125       *    @return array Lower triangular factor
 126       */
 127      public function getL()
 128      {
 129          for ($i = 0; $i < $this->m; ++$i) {
 130              for ($j = 0; $j < $this->n; ++$j) {
 131                  if ($i > $j) {
 132                      $L[$i][$j] = $this->LU[$i][$j];
 133                  } elseif ($i == $j) {
 134                      $L[$i][$j] = 1.0;
 135                  } else {
 136                      $L[$i][$j] = 0.0;
 137                  }
 138              }
 139          }
 140          return new PHPExcel_Shared_JAMA_Matrix($L);
 141      }    //    function getL()
 142  
 143      /**
 144       *    Get upper triangular factor.
 145       *
 146       *    @return array Upper triangular factor
 147       */
 148      public function getU()
 149      {
 150          for ($i = 0; $i < $this->n; ++$i) {
 151              for ($j = 0; $j < $this->n; ++$j) {
 152                  if ($i <= $j) {
 153                      $U[$i][$j] = $this->LU[$i][$j];
 154                  } else {
 155                      $U[$i][$j] = 0.0;
 156                  }
 157              }
 158          }
 159          return new PHPExcel_Shared_JAMA_Matrix($U);
 160      }    //    function getU()
 161  
 162      /**
 163       *    Return pivot permutation vector.
 164       *
 165       *    @return array Pivot vector
 166       */
 167      public function getPivot()
 168      {
 169          return $this->piv;
 170      }    //    function getPivot()
 171  
 172      /**
 173       *    Alias for getPivot
 174       *
 175       *    @see getPivot
 176       */
 177      public function getDoublePivot()
 178      {
 179          return $this->getPivot();
 180      }    //    function getDoublePivot()
 181  
 182      /**
 183       *    Is the matrix nonsingular?
 184       *
 185       *    @return true if U, and hence A, is nonsingular.
 186       */
 187      public function isNonsingular()
 188      {
 189          for ($j = 0; $j < $this->n; ++$j) {
 190              if ($this->LU[$j][$j] == 0) {
 191                  return false;
 192              }
 193          }
 194          return true;
 195      }    //    function isNonsingular()
 196  
 197      /**
 198       *    Count determinants
 199       *
 200       *    @return array d matrix deterninat
 201       */
 202      public function det()
 203      {
 204          if ($this->m == $this->n) {
 205              $d = $this->pivsign;
 206              for ($j = 0; $j < $this->n; ++$j) {
 207                  $d *= $this->LU[$j][$j];
 208              }
 209              return $d;
 210          } else {
 211              throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MATRIX_DIMENSION_EXCEPTION);
 212          }
 213      }    //    function det()
 214  
 215      /**
 216       *    Solve A*X = B
 217       *
 218       *    @param  $B  A Matrix with as many rows as A and any number of columns.
 219       *    @return  X so that L*U*X = B(piv,:)
 220       *    @PHPExcel_Calculation_Exception  IllegalArgumentException Matrix row dimensions must agree.
 221       *    @PHPExcel_Calculation_Exception  RuntimeException  Matrix is singular.
 222       */
 223      public function solve($B)
 224      {
 225          if ($B->getRowDimension() == $this->m) {
 226              if ($this->isNonsingular()) {
 227                  // Copy right hand side with pivoting
 228                  $nx = $B->getColumnDimension();
 229                  $X  = $B->getMatrix($this->piv, 0, $nx-1);
 230                  // Solve L*Y = B(piv,:)
 231                  for ($k = 0; $k < $this->n; ++$k) {
 232                      for ($i = $k+1; $i < $this->n; ++$i) {
 233                          for ($j = 0; $j < $nx; ++$j) {
 234                              $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
 235                          }
 236                      }
 237                  }
 238                  // Solve U*X = Y;
 239                  for ($k = $this->n-1; $k >= 0; --$k) {
 240                      for ($j = 0; $j < $nx; ++$j) {
 241                          $X->A[$k][$j] /= $this->LU[$k][$k];
 242                      }
 243                      for ($i = 0; $i < $k; ++$i) {
 244                          for ($j = 0; $j < $nx; ++$j) {
 245                              $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
 246                          }
 247                      }
 248                  }
 249                  return $X;
 250              } else {
 251                  throw new PHPExcel_Calculation_Exception(self::MATRIX_SINGULAR_EXCEPTION);
 252              }
 253          } else {
 254              throw new PHPExcel_Calculation_Exception(self::MATRIX_SQUARE_EXCEPTION);
 255          }
 256      }
 257  }