LUDecomposition::getDoublePivot() - Code Metrics - Inspection of "Support COUNTIFS multiple arguments" - PHPOffice/PhpSpreadsheet - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Push — develop ( 50ed76...6a48b5 )

by Adrien

created 2019-01-01 05:39 UTC

LUDecomposition::getDoublePivot() A

↳ Parent: LUDecomposition

Complexity

Conditions	1
Paths	1

Size

Total Lines	3
Code Lines	1

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	1
eloc	1
nc	1
nop	0
dl	0
loc	3
rs	10
c	0
b	0
f	0

<?php

namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;

use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;

/**
 *    For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
 *    unit lower triangular matrix L, an n-by-n upper triangular matrix U,
 *    and a permutation vector piv of length m so that A(piv,:) = L*U.
 *    If m < n, then L is m-by-m and U is m-by-n.
 *
 *    The LU decompostion with pivoting always exists, even if the matrix is
 *    singular, so the constructor will never fail. The primary use of the
 *    LU decomposition is in the solution of square systems of simultaneous
 *    linear equations. This will fail if isNonsingular() returns false.
 *
 *    @author Paul Meagher
 *    @author Bartosz Matosiuk
 *    @author Michael Bommarito
 *
 *    @version 1.1
 */
class LUDecomposition
{
    const MATRIX_SINGULAR_EXCEPTION = 'Can only perform operation on singular matrix.';
    const MATRIX_SQUARE_EXCEPTION = 'Mismatched Row dimension';

    /**
     * Decomposition storage.
     *
     * @var array
     */
    private $LU = [];

    /**
     * Row dimension.
     *
     * @var int
     */
    private $m;

    /**
     * Column dimension.
     *
     * @var int
     */
    private $n;

    /**
     * Pivot sign.
     *
     * @var int
     */
    private $pivsign;

    /**
     * Internal storage of pivot vector.
     *
     * @var array
     */
    private $piv = [];

    /**
     * LU Decomposition constructor.
     *
     * @param Matrix $A Rectangular matrix
     */
    public function __construct($A)
    {
        if ($A instanceof Matrix) {

            // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
            $this->LU = $A->getArray();
            $this->m = $A->getRowDimension();
            $this->n = $A->getColumnDimension();
            for ($i = 0; $i < $this->m; ++$i) {
                $this->piv[$i] = $i;
            }
            $this->pivsign = 1;
            $LUrowi = $LUcolj = [];


            // Outer loop.
            for ($j = 0; $j < $this->n; ++$j) {
                // Make a copy of the j-th column to localize references.
                for ($i = 0; $i < $this->m; ++$i) {
                    $LUcolj[$i] = &$this->LU[$i][$j];
                }
                // Apply previous transformations.
                for ($i = 0; $i < $this->m; ++$i) {
                    $LUrowi = $this->LU[$i];
                    // Most of the time is spent in the following dot product.
                    $kmax = min($i, $j);
                    $s = 0.0;
                    for ($k = 0; $k < $kmax; ++$k) {
                        $s += $LUrowi[$k] * $LUcolj[$k];
                    }
                    $LUrowi[$j] = $LUcolj[$i] -= $s;
                }
                // Find pivot and exchange if necessary.
                $p = $j;
                for ($i = $j + 1; $i < $this->m; ++$i) {
                    if (abs($LUcolj[$i]) > abs($LUcolj[$p])) {
                        $p = $i;
                    }
                }
                if ($p != $j) {
                    for ($k = 0; $k < $this->n; ++$k) {
                        $t = $this->LU[$p][$k];
                        $this->LU[$p][$k] = $this->LU[$j][$k];
                        $this->LU[$j][$k] = $t;
                    }
                    $k = $this->piv[$p];
                    $this->piv[$p] = $this->piv[$j];
                    $this->piv[$j] = $k;
                    $this->pivsign = $this->pivsign * -1;
                }
                // Compute multipliers.
                if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
                    for ($i = $j + 1; $i < $this->m; ++$i) {
                        $this->LU[$i][$j] /= $this->LU[$j][$j];
                    }
                }
            }
        } else {
            throw new CalculationException(Matrix::ARGUMENT_TYPE_EXCEPTION);
        }
    }

    //    function __construct()

    /**
     * Get lower triangular factor.
     *
     * @return Matrix Lower triangular factor
     */
    public function getL()
    {
        for ($i = 0; $i < $this->m; ++$i) {
            for ($j = 0; $j < $this->n; ++$j) {
                if ($i > $j) {
                    $L[$i][$j] = $this->LU[$i][$j];
                } elseif ($i == $j) {
                    $L[$i][$j] = 1.0;
                } else {
                    $L[$i][$j] = 0.0;
                }
            }
        }

        return new Matrix($L);

    }

    //    function getL()

    /**
     * Get upper triangular factor.
     *
     * @return Matrix Upper triangular factor
     */
    public function getU()
    {
        for ($i = 0; $i < $this->n; ++$i) {
            for ($j = 0; $j < $this->n; ++$j) {
                if ($i <= $j) {
                    $U[$i][$j] = $this->LU[$i][$j];
                } else {
                    $U[$i][$j] = 0.0;
                }
            }
        }

        return new Matrix($U);

    }

    //    function getU()

    /**
     * Return pivot permutation vector.
     *
     * @return array Pivot vector
     */
    public function getPivot()
    {
        return $this->piv;
    }

    //    function getPivot()

    /**
     * Alias for getPivot.
     *
     *    @see getPivot
     */
    public function getDoublePivot()
    {
        return $this->getPivot();
    }

    //    function getDoublePivot()

    /**
     *    Is the matrix nonsingular?
     *
     * @return bool true if U, and hence A, is nonsingular
     */
    public function isNonsingular()
    {
        for ($j = 0; $j < $this->n; ++$j) {
            if ($this->LU[$j][$j] == 0) {
                return false;
            }
        }

        return true;
    }

    //    function isNonsingular()

    /**
     * Count determinants.
     *
     * @return array d matrix deterninat
     */
    public function det()
    {
        if ($this->m == $this->n) {
            $d = $this->pivsign;
            for ($j = 0; $j < $this->n; ++$j) {
                $d *= $this->LU[$j][$j];
            }

            return $d;

        }

        throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
    }

    //    function det()

    /**
     * Solve A*X = B.
     *
     * @param mixed $B a Matrix with as many rows as A and any number of columns
     *
     * @throws CalculationException illegalArgumentException Matrix row dimensions must agree
     * @throws CalculationException runtimeException  Matrix is singular
     *
     * @return Matrix X so that L*U*X = B(piv,:)
     */
    public function solve($B)
    {
        if ($B->getRowDimension() == $this->m) {
            if ($this->isNonsingular()) {
                // Copy right hand side with pivoting
                $nx = $B->getColumnDimension();
                $X = $B->getMatrix($this->piv, 0, $nx - 1);
                // Solve L*Y = B(piv,:)
                for ($k = 0; $k < $this->n; ++$k) {
                    for ($i = $k + 1; $i < $this->n; ++$i) {
                        for ($j = 0; $j < $nx; ++$j) {
                            $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
                        }
                    }
                }
                // Solve U*X = Y;
                for ($k = $this->n - 1; $k >= 0; --$k) {
                    for ($j = 0; $j < $nx; ++$j) {
                        $X->A[$k][$j] /= $this->LU[$k][$k];
                    }
                    for ($i = 0; $i < $k; ++$i) {
                        for ($j = 0; $j < $nx; ++$j) {
                            $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
                        }
                    }
                }

                return $X;
            }

            throw new CalculationException(self::MATRIX_SINGULAR_EXCEPTION);
        }

        throw new CalculationException(self::MATRIX_SQUARE_EXCEPTION);
    }
}


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) {
			0 ignored issues – show introduced 2018-10-21 04:55 UTC by Report Bug Copy Issue Report `$A` is always a sub-type of `PhpOffice\PhpSpreadsheet\Shared\JAMA\Matrix`. Loading history...
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 = [];
			0 ignored issues – show Unused Code introduced 2017-12-03 09:09 UTC by Report Bug Copy Issue Report The assignment to `$LUrowi` is dead and can be removed. Loading history...
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);
			0 ignored issues – show Comprehensibility Best Practice introduced 2017-12-03 09:09 UTC by Report Bug Copy Issue Report The variable `$L` does not seem to be defined for all execution paths leading up to this point. Loading history...
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);
			0 ignored issues – show Comprehensibility Best Practice introduced 2017-12-03 09:09 UTC by Report Bug Copy Issue Report The variable `$U` does not seem to be defined for all execution paths leading up to this point. Loading history...
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;
			0 ignored issues – show Bug Best Practice introduced 2017-12-03 09:09 UTC by Report Bug Copy Issue Report The expression `return $d` returns the type `integer` which is incompatible with the documented return type `array`. Loading history...
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 LUX = 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

PHPOffice / PhpSpreadsheet

Push — develop ( 50ed76...6a48b5 )

LUDecomposition::getDoublePivot() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like