001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.linear;
019    
020    import org.apache.commons.math.linear.InvalidMatrixException;
021    import org.apache.commons.math.linear.LUDecomposition;
022    import org.apache.commons.math.linear.LUDecompositionImpl;
023    import org.apache.commons.math.linear.MatrixUtils;
024    import org.apache.commons.math.linear.RealMatrix;
025    
026    import junit.framework.Test;
027    import junit.framework.TestCase;
028    import junit.framework.TestSuite;
029    
030    public class LUDecompositionImplTest extends TestCase {
031        private double[][] testData = {
032                { 1.0, 2.0, 3.0},
033                { 2.0, 5.0, 3.0},
034                { 1.0, 0.0, 8.0}
035        };
036        private double[][] testDataMinus = {
037                { -1.0, -2.0, -3.0},
038                { -2.0, -5.0, -3.0},
039                { -1.0,  0.0, -8.0}
040        };
041        private double[][] luData = {
042                { 2.0, 3.0, 3.0 },
043                { 0.0, 5.0, 7.0 },
044                { 6.0, 9.0, 8.0 }
045        };
046        
047        // singular matrices
048        private double[][] singular = {
049                { 2.0, 3.0 },
050                { 2.0, 3.0 }
051        };
052        private double[][] bigSingular = {
053                { 1.0, 2.0,   3.0,    4.0 },
054                { 2.0, 5.0,   3.0,    4.0 },
055                { 7.0, 3.0, 256.0, 1930.0 },
056                { 3.0, 7.0,   6.0,    8.0 }
057        }; // 4th row = 1st + 2nd
058    
059        private static final double entryTolerance = 10e-16;
060    
061        private static final double normTolerance = 10e-14;
062    
063        public LUDecompositionImplTest(String name) {
064            super(name);
065        }
066    
067        public static Test suite() {
068            TestSuite suite = new TestSuite(LUDecompositionImplTest.class);
069            suite.setName("LUDecompositionImpl Tests");
070            return suite;
071        }
072    
073        /** test dimensions */
074        public void testDimensions() {
075            RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
076            LUDecomposition LU = new LUDecompositionImpl(matrix);
077            assertEquals(testData.length, LU.getL().getRowDimension());
078            assertEquals(testData.length, LU.getL().getColumnDimension());
079            assertEquals(testData.length, LU.getU().getRowDimension());
080            assertEquals(testData.length, LU.getU().getColumnDimension());
081            assertEquals(testData.length, LU.getP().getRowDimension());
082            assertEquals(testData.length, LU.getP().getColumnDimension());
083    
084        }
085    
086        /** test non-square matrix */
087        public void testNonSquare() {
088            try {
089                new LUDecompositionImpl(MatrixUtils.createRealMatrix(new double[3][2]));
090            } catch (InvalidMatrixException ime) {
091                // expected behavior
092            } catch (Exception e) {
093                fail("wrong exception caught");
094            }
095        }
096    
097        /** test PA = LU */
098        public void testPAEqualLU() {
099            RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
100            LUDecomposition lu = new LUDecompositionImpl(matrix);
101            RealMatrix l = lu.getL();
102            RealMatrix u = lu.getU();
103            RealMatrix p = lu.getP();
104            double norm = l.multiply(u).subtract(p.multiply(matrix)).getNorm();
105            assertEquals(0, norm, normTolerance);
106    
107            matrix = MatrixUtils.createRealMatrix(testDataMinus);
108            lu = new LUDecompositionImpl(matrix);
109            l = lu.getL();
110            u = lu.getU();
111            p = lu.getP();
112            norm = l.multiply(u).subtract(p.multiply(matrix)).getNorm();
113            assertEquals(0, norm, normTolerance);
114    
115            matrix = MatrixUtils.createRealIdentityMatrix(17);
116            lu = new LUDecompositionImpl(matrix);
117            l = lu.getL();
118            u = lu.getU();
119            p = lu.getP();
120            norm = l.multiply(u).subtract(p.multiply(matrix)).getNorm();
121            assertEquals(0, norm, normTolerance);
122    
123            matrix = MatrixUtils.createRealMatrix(singular);
124            lu = new LUDecompositionImpl(matrix);
125            assertFalse(lu.getSolver().isNonSingular());
126            assertNull(lu.getL());
127            assertNull(lu.getU());
128            assertNull(lu.getP());
129    
130            matrix = MatrixUtils.createRealMatrix(bigSingular);
131            lu = new LUDecompositionImpl(matrix);
132            assertFalse(lu.getSolver().isNonSingular());
133            assertNull(lu.getL());
134            assertNull(lu.getU());
135            assertNull(lu.getP());
136    
137        }
138    
139        /** test that L is lower triangular with unit diagonal */
140        public void testLLowerTriangular() {
141            RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
142            RealMatrix l = new LUDecompositionImpl(matrix).getL();
143            for (int i = 0; i < l.getRowDimension(); i++) {
144                assertEquals(l.getEntry(i, i), 1, entryTolerance);
145                for (int j = i + 1; j < l.getColumnDimension(); j++) {
146                    assertEquals(l.getEntry(i, j), 0, entryTolerance);
147                }
148            }
149        }
150    
151        /** test that U is upper triangular */
152        public void testUUpperTriangular() {
153            RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
154            RealMatrix u = new LUDecompositionImpl(matrix).getU();
155            for (int i = 0; i < u.getRowDimension(); i++) {
156                for (int j = 0; j < i; j++) {
157                    assertEquals(u.getEntry(i, j), 0, entryTolerance);
158                }
159            }
160        }
161    
162        /** test that P is a permutation matrix */
163        public void testPPermutation() {
164            RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
165            RealMatrix p   = new LUDecompositionImpl(matrix).getP();
166    
167            RealMatrix ppT = p.multiply(p.transpose());
168            RealMatrix id  = MatrixUtils.createRealIdentityMatrix(p.getRowDimension());
169            assertEquals(0, ppT.subtract(id).getNorm(), normTolerance);
170    
171            for (int i = 0; i < p.getRowDimension(); i++) {
172                int zeroCount  = 0;
173                int oneCount   = 0;
174                int otherCount = 0;
175                for (int j = 0; j < p.getColumnDimension(); j++) {
176                    final double e = p.getEntry(i, j);
177                    if (e == 0) {
178                        ++zeroCount;
179                    } else if (e == 1) {
180                        ++oneCount;
181                    } else {
182                        ++otherCount;
183                    }
184                }
185                assertEquals(p.getColumnDimension() - 1, zeroCount);
186                assertEquals(1, oneCount);
187                assertEquals(0, otherCount);
188            }
189    
190            for (int j = 0; j < p.getColumnDimension(); j++) {
191                int zeroCount  = 0;
192                int oneCount   = 0;
193                int otherCount = 0;
194                for (int i = 0; i < p.getRowDimension(); i++) {
195                    final double e = p.getEntry(i, j);
196                    if (e == 0) {
197                        ++zeroCount;
198                    } else if (e == 1) {
199                        ++oneCount;
200                    } else {
201                        ++otherCount;
202                    }
203                }
204                assertEquals(p.getRowDimension() - 1, zeroCount);
205                assertEquals(1, oneCount);
206                assertEquals(0, otherCount);
207            }
208    
209        }
210    
211    
212        /** test singular */
213        public void testSingular() {
214            LUDecomposition lu =
215                new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData));
216            assertTrue(lu.getSolver().isNonSingular());
217            lu = new LUDecompositionImpl(MatrixUtils.createRealMatrix(singular));
218            assertFalse(lu.getSolver().isNonSingular());
219            lu = new LUDecompositionImpl(MatrixUtils.createRealMatrix(bigSingular));
220            assertFalse(lu.getSolver().isNonSingular());
221        }
222    
223        /** test matrices values */
224        public void testMatricesValues1() {
225           LUDecomposition lu =
226                new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData));
227            RealMatrix lRef = MatrixUtils.createRealMatrix(new double[][] {
228                    { 1.0, 0.0, 0.0 },
229                    { 0.5, 1.0, 0.0 },
230                    { 0.5, 0.2, 1.0 }
231            });
232            RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
233                    { 2.0,  5.0, 3.0 },
234                    { 0.0, -2.5, 6.5 },
235                    { 0.0,  0.0, 0.2 }
236            });
237            RealMatrix pRef = MatrixUtils.createRealMatrix(new double[][] {
238                    { 0.0, 1.0, 0.0 },
239                    { 0.0, 0.0, 1.0 },
240                    { 1.0, 0.0, 0.0 }
241            });
242            int[] pivotRef = { 1, 2, 0 };
243    
244            // check values against known references
245            RealMatrix l = lu.getL();
246            assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
247            RealMatrix u = lu.getU();
248            assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
249            RealMatrix p = lu.getP();
250            assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
251            int[] pivot = lu.getPivot();
252            for (int i = 0; i < pivotRef.length; ++i) {
253                assertEquals(pivotRef[i], pivot[i]);
254            }
255    
256            // check the same cached instance is returned the second time
257            assertTrue(l == lu.getL());
258            assertTrue(u == lu.getU());
259            assertTrue(p == lu.getP());
260            
261        }
262    
263        /** test matrices values */
264        public void testMatricesValues2() {
265           LUDecomposition lu =
266                new LUDecompositionImpl(MatrixUtils.createRealMatrix(luData));
267            RealMatrix lRef = MatrixUtils.createRealMatrix(new double[][] {
268                    {    1.0,    0.0, 0.0 },
269                    {    0.0,    1.0, 0.0 },
270                    { 1.0 / 3.0, 0.0, 1.0 }
271            });
272            RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
273                    { 6.0, 9.0,    8.0    },
274                    { 0.0, 5.0,    7.0    },
275                    { 0.0, 0.0, 1.0 / 3.0 }
276            });
277            RealMatrix pRef = MatrixUtils.createRealMatrix(new double[][] {
278                    { 0.0, 0.0, 1.0 },
279                    { 0.0, 1.0, 0.0 },
280                    { 1.0, 0.0, 0.0 }
281            });
282            int[] pivotRef = { 2, 1, 0 };
283    
284            // check values against known references
285            RealMatrix l = lu.getL();
286            assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
287            RealMatrix u = lu.getU();
288            assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
289            RealMatrix p = lu.getP();
290            assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
291            int[] pivot = lu.getPivot();
292            for (int i = 0; i < pivotRef.length; ++i) {
293                assertEquals(pivotRef[i], pivot[i]);
294            }
295    
296            // check the same cached instance is returned the second time
297            assertTrue(l == lu.getL());
298            assertTrue(u == lu.getU());
299            assertTrue(p == lu.getP());
300            
301        }
302    
303    }