1<?php
2
3namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;
4
5use 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 */
24class 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     * @throws CalculationException illegalArgumentException Matrix row dimensions must agree
246     * @throws CalculationException runtimeException  Matrix is singular
247     *
248     * @return Matrix X so that L*U*X = B(piv,:)
249     */
250    public function solve($B)
251    {
252        if ($B->getRowDimension() == $this->m) {
253            if ($this->isNonsingular()) {
254                // Copy right hand side with pivoting
255                $nx = $B->getColumnDimension();
256                $X = $B->getMatrix($this->piv, 0, $nx - 1);
257                // Solve L*Y = B(piv,:)
258                for ($k = 0; $k < $this->n; ++$k) {
259                    for ($i = $k + 1; $i < $this->n; ++$i) {
260                        for ($j = 0; $j < $nx; ++$j) {
261                            $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
262                        }
263                    }
264                }
265                // Solve U*X = Y;
266                for ($k = $this->n - 1; $k >= 0; --$k) {
267                    for ($j = 0; $j < $nx; ++$j) {
268                        $X->A[$k][$j] /= $this->LU[$k][$k];
269                    }
270                    for ($i = 0; $i < $k; ++$i) {
271                        for ($j = 0; $j < $nx; ++$j) {
272                            $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
273                        }
274                    }
275                }
276
277                return $X;
278            }
279
280            throw new CalculationException(self::MATRIX_SINGULAR_EXCEPTION);
281        }
282
283        throw new CalculationException(self::MATRIX_SQUARE_EXCEPTION);
284    }
285}
286