← Back to branch summary

~ubuntu-branches/ubuntu/trusty/moodle/trusty-proposed

« back to all changes in this revision

Viewing changes to lib/phpexcel/PHPExcel/Shared/JAMA/LUDecomposition.php

  • Committer: Package Import Robot
  • Author(s): Thijs Kinkhorst
  • Date: 2013-07-19 08:52:46 UTC
  • mfrom: (1.1.10)
  • Revision ID: package-import@ubuntu.com-20130719085246-yebwditc2exoap2r
Tags: 2.5.1-1
http://bugs.debian.org/429339
http://bugs.debian.org/716972
http://bugs.debian.org/716986
http://bugs.debian.org/717080
http://bugs.debian.org/717108
http://bugs.debian.org/717278
* New upstream version: 2.5.1.
  - Fixes security issues:
    CVE-2013-2242 CVE-2013-2243 CVE-2013-2244 CVE-2013-2245
    CVE-2013-2246
* Depend on apache2 instead of obsolete apache2-mpm-prefork.
* Use packaged libphp-phpmailer (closes: #429339), adodb,
  HTMLPurifier, PclZip.
* Update debconf translations, thanks Salvatore Merone, Pietro Tollot,
  Joe Hansen, Yuri Kozlov, Holger Wansing, Américo Monteiro,
  Adriano Rafael Gomes, victory, Michał Kułach.
  (closes: #716972, #716986, #717080, #717108, #717278)

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/phpexcel/PHPExcel/Shared/JAMA/LUDecomposition.php
 
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 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 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 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
         *      @exception  IllegalArgumentException Matrix row dimensions must agree.
 
222
         *      @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 Exception(self::MatrixSingularException);
 
252
                        }
 
253
                } else {
 
254
                        throw new Exception(self::MatrixSquareException);
 
255
                }
 
256
        }       //      function solve()
 
257
 
 
258
}       //      class PHPExcel_Shared_JAMA_LUDecomposition