1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math.linear;
19  
20  import org.apache.commons.math.linear.InvalidMatrixException;
21  import org.apache.commons.math.linear.LUDecomposition;
22  import org.apache.commons.math.linear.LUDecompositionImpl;
23  import org.apache.commons.math.linear.MatrixUtils;
24  import org.apache.commons.math.linear.RealMatrix;
25  
26  import junit.framework.Test;
27  import junit.framework.TestCase;
28  import junit.framework.TestSuite;
29  
30  public class LUDecompositionImplTest extends TestCase {
31      private double[][] testData = {
32              { 1.0, 2.0, 3.0},
33              { 2.0, 5.0, 3.0},
34              { 1.0, 0.0, 8.0}
35      };
36      private double[][] testDataMinus = {
37              { -1.0, -2.0, -3.0},
38              { -2.0, -5.0, -3.0},
39              { -1.0,  0.0, -8.0}
40      };
41      private double[][] luData = {
42              { 2.0, 3.0, 3.0 },
43              { 0.0, 5.0, 7.0 },
44              { 6.0, 9.0, 8.0 }
45      };
46      
47      // singular matrices
48      private double[][] singular = {
49              { 2.0, 3.0 },
50              { 2.0, 3.0 }
51      };
52      private double[][] bigSingular = {
53              { 1.0, 2.0,   3.0,    4.0 },
54              { 2.0, 5.0,   3.0,    4.0 },
55              { 7.0, 3.0, 256.0, 1930.0 },
56              { 3.0, 7.0,   6.0,    8.0 }
57      }; // 4th row = 1st + 2nd
58  
59      private static final double entryTolerance = 10e-16;
60  
61      private static final double normTolerance = 10e-14;
62  
63      public LUDecompositionImplTest(String name) {
64          super(name);
65      }
66  
67      public static Test suite() {
68          TestSuite suite = new TestSuite(LUDecompositionImplTest.class);
69          suite.setName("LUDecompositionImpl Tests");
70          return suite;
71      }
72  
73      /** test dimensions */
74      public void testDimensions() {
75          RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
76          LUDecomposition LU = new LUDecompositionImpl(matrix);
77          assertEquals(testData.length, LU.getL().getRowDimension());
78          assertEquals(testData.length, LU.getL().getColumnDimension());
79          assertEquals(testData.length, LU.getU().getRowDimension());
80          assertEquals(testData.length, LU.getU().getColumnDimension());
81          assertEquals(testData.length, LU.getP().getRowDimension());
82          assertEquals(testData.length, LU.getP().getColumnDimension());
83  
84      }
85  
86      /** test non-square matrix */
87      public void testNonSquare() {
88          try {
89              new LUDecompositionImpl(MatrixUtils.createRealMatrix(new double[3][2]));
90          } catch (InvalidMatrixException ime) {
91              // expected behavior
92          } catch (Exception e) {
93              fail("wrong exception caught");
94          }
95      }
96  
97      /** test PA = LU */
98      public void testPAEqualLU() {
99          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 }