[Summary view]

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

   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          for ($i = 0; $i < $this->m; ++$i) {
 139              for ($j = 0; $j < $this->n; ++$j) {
 140                  if ($i > $j) {
 141                      $L[$i][$j] = $this->LU[$i][$j];
 142                  } elseif ($i == $j) {
 143                      $L[$i][$j] = 1.0;
 144                  } else {
 145                      $L[$i][$j] = 0.0;
 146                  }
 147              }
 148          }
 149  
 150          return new Matrix($L);
 151      }
 152  
 153      //    function getL()
 154  
 155      /**
 156       * Get upper triangular factor.
 157       *
 158       * @return Matrix Upper triangular factor
 159       */
 160      public function getU()
 161      {
 162          for ($i = 0; $i < $this->n; ++$i) {
 163              for ($j = 0; $j < $this->n; ++$j) {
 164                  if ($i <= $j) {
 165                      $U[$i][$j] = $this->LU[$i][$j];
 166                  } else {
 167                      $U[$i][$j] = 0.0;
 168                  }
 169              }
 170          }
 171  
 172          return new Matrix($U);
 173      }
 174  
 175      //    function getU()
 176  
 177      /**
 178       * Return pivot permutation vector.
 179       *
 180       * @return array Pivot vector
 181       */
 182      public function getPivot()
 183      {
 184          return $this->piv;
 185      }
 186  
 187      //    function getPivot()
 188  
 189      /**
 190       * Alias for getPivot.
 191       *
 192       *    @see getPivot
 193       */
 194      public function getDoublePivot()
 195      {
 196          return $this->getPivot();
 197      }
 198  
 199      //    function getDoublePivot()
 200  
 201      /**
 202       *    Is the matrix nonsingular?
 203       *
 204       * @return bool true if U, and hence A, is nonsingular
 205       */
 206      public function isNonsingular()
 207      {
 208          for ($j = 0; $j < $this->n; ++$j) {
 209              if ($this->LU[$j][$j] == 0) {
 210                  return false;
 211              }
 212          }
 213  
 214          return true;
 215      }
 216  
 217      //    function isNonsingular()
 218  
 219      /**
 220       * Count determinants.
 221       *
 222       * @return array d matrix deterninat
 223       */
 224      public function det()
 225      {
 226          if ($this->m == $this->n) {
 227              $d = $this->pivsign;
 228              for ($j = 0; $j < $this->n; ++$j) {
 229                  $d *= $this->LU[$j][$j];
 230              }
 231  
 232              return $d;
 233          }
 234  
 235          throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
 236      }
 237  
 238      //    function det()
 239  
 240      /**
 241       * Solve A*X = B.
 242       *
 243       * @param mixed $B a Matrix with as many rows as A and any number of columns
 244       *
 245       * @return Matrix X so that L*U*X = B(piv,:)
 246       */
 247      public function solve($B)
 248      {
 249          if ($B->getRowDimension() == $this->m) {
 250              if ($this->isNonsingular()) {
 251                  // Copy right hand side with pivoting
 252                  $nx = $B->getColumnDimension();
 253                  $X = $B->getMatrix($this->piv, 0, $nx - 1);
 254                  // Solve L*Y = B(piv,:)
 255                  for ($k = 0; $k < $this->n; ++$k) {
 256                      for ($i = $k + 1; $i < $this->n; ++$i) {
 257                          for ($j = 0; $j < $nx; ++$j) {
 258                              $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
 259                          }
 260                      }
 261                  }
 262                  // Solve U*X = Y;
 263                  for ($k = $this->n - 1; $k >= 0; --$k) {
 264                      for ($j = 0; $j < $nx; ++$j) {
 265                          $X->A[$k][$j] /= $this->LU[$k][$k];
 266                      }
 267                      for ($i = 0; $i < $k; ++$i) {
 268                          for ($j = 0; $j < $nx; ++$j) {
 269                              $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
 270                          }
 271                      }
 272                  }
 273  
 274                  return $X;
 275              }
 276  
 277              throw new CalculationException(self::MATRIX_SINGULAR_EXCEPTION);
 278          }
 279  
 280          throw new CalculationException(self::MATRIX_SQUARE_EXCEPTION);
 281      }
 282  }