source: sourcecode/application/libraries/PHPExcel/Shared/JAMA/LUDecomposition.php @ 1

Last change on this file since 1 was 1, checked in by dungnv, 11 years ago
File size: 6.1 KB
RevLine 
[1]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 */
21class PHPExcel_Shared_JAMA_LUDecomposition {
22
23        const MatrixSingularException   = "Can only perform operation on singular matrix.";
24        const MatrixSquareException             = "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        /**
58         *      LU Decomposition constructor.
59         *
60         *      @param $A Rectangular matrix
61         *      @return Structure to access L, U and piv.
62         */
63        public function __construct($A) {
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::ArgumentTypeException);
119                }
120        }       //      function __construct()
121
122
123        /**
124         *      Get lower triangular factor.
125         *
126         *      @return array Lower triangular factor
127         */
128        public function getL() {
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        /**
145         *      Get upper triangular factor.
146         *
147         *      @return array Upper triangular factor
148         */
149        public function getU() {
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        /**
164         *      Return pivot permutation vector.
165         *
166         *      @return array Pivot vector
167         */
168        public function getPivot() {
169                return $this->piv;
170        }       //      function getPivot()
171
172
173        /**
174         *      Alias for getPivot
175         *
176         *      @see getPivot
177         */
178        public function getDoublePivot() {
179                return $this->getPivot();
180        }       //      function getDoublePivot()
181
182
183        /**
184         *      Is the matrix nonsingular?
185         *
186         *      @return true if U, and hence A, is nonsingular.
187         */
188        public function isNonsingular() {
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        /**
199         *      Count determinants
200         *
201         *      @return array d matrix deterninat
202         */
203        public function det() {
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::MatrixDimensionException);
212                }
213        }       //      function det()
214
215
216        /**
217         *      Solve A*X = B
218         *
219         *      @param  $B  A Matrix with as many rows as A and any number of columns.
220         *      @return  X so that L*U*X = B(piv,:)
221         *      @PHPExcel_Calculation_Exception  IllegalArgumentException Matrix row dimensions must agree.
222         *      @PHPExcel_Calculation_Exception  RuntimeException  Matrix is singular.
223         */
224        public function solve($B) {
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::MatrixSingularException);
252                        }
253                } else {
254                        throw new PHPExcel_Calculation_Exception(self::MatrixSquareException);
255                }
256        }       //      function solve()
257
258}       //      class PHPExcel_Shared_JAMA_LUDecomposition
Note: See TracBrowser for help on using the repository browser.