[Summary view]

Differences Between: [Versions 310 and 401] [Versions 311 and 401] [Versions 39 and 401]

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