001    package org.ujmp.core.doublematrix.calculation.general.decomposition;
002    
003    import org.ujmp.core.Matrix;
004    import org.ujmp.core.MatrixFactory;
005    import org.ujmp.core.doublematrix.DenseDoubleMatrix2D;
006    import org.ujmp.core.util.DecompositionOps;
007    import org.ujmp.core.util.UJMPSettings;
008    
009    /**
010     * LU Decomposition.
011     * <P>
012     * For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n unit
013     * lower triangular matrix L, an n-by-n upper triangular matrix U, and a
014     * permutation vector piv of length m so that A(piv,:) = L*U. If m < n, then L
015     * is m-by-m and U is m-by-n.
016     * <P>
017     * The LU decompostion with pivoting always exists, even if the matrix is
018     * singular, so the constructor will never fail. The primary use of the LU
019     * decomposition is in the solution of square systems of simultaneous linear
020     * equations. This will fail if isNonsingular() returns false.
021     */
022    
023    public interface LU<T> {
024    
025            public T[] calc(T source);
026    
027            public T solve(T source, T b);
028    
029            public static int THRESHOLD = 100;
030    
031            public static final LU<Matrix> MATRIX = new LU<Matrix>() {
032    
033                    public final Matrix[] calc(Matrix source) {
034                            if (UJMPSettings.getNumberOfThreads() == 1) {
035                                    if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) {
036                                            return MATRIXLARGESINGLETHREADED.calc(source);
037                                    } else {
038                                            return MATRIXSMALLSINGLETHREADED.calc(source);
039                                    }
040                            } else {
041                                    if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) {
042                                            return MATRIXLARGEMULTITHREADED.calc(source);
043                                    } else {
044                                            return MATRIXSMALLMULTITHREADED.calc(source);
045                                    }
046                            }
047                    }
048    
049                    public final Matrix solve(Matrix source, Matrix b) {
050                            if (UJMPSettings.getNumberOfThreads() == 1) {
051                                    if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) {
052                                            return MATRIXLARGESINGLETHREADED.solve(source, b);
053                                    } else {
054                                            return MATRIXSMALLSINGLETHREADED.solve(source, b);
055                                    }
056                            } else {
057                                    if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) {
058                                            return MATRIXLARGEMULTITHREADED.solve(source, b);
059                                    } else {
060                                            return MATRIXSMALLMULTITHREADED.solve(source, b);
061                                    }
062                            }
063                    }
064            };
065    
066            public static final LU<Matrix> INSTANCE = MATRIX;
067    
068            public static final LU<Matrix> UJMP = new LU<Matrix>() {
069    
070                    public final Matrix[] calc(Matrix source) {
071                            LUMatrix lu = new LUMatrix(source);
072                            return new Matrix[] { lu.getL(), lu.getU(), lu.getP() };
073                    }
074    
075                    public final Matrix solve(Matrix source, Matrix b) {
076                            LUMatrix lu = new LUMatrix(source);
077                            return lu.solve(b);
078                    }
079            };
080    
081            public static final LU<Matrix> MATRIXSMALLMULTITHREADED = UJMP;
082    
083            public static final LU<Matrix> MATRIXSMALLSINGLETHREADED = UJMP;
084    
085            public static final LU<Matrix> MATRIXLARGESINGLETHREADED = new LU<Matrix>() {
086                    public final Matrix[] calc(Matrix source) {
087                            LU<Matrix> lu = null;
088                            if (UJMPSettings.isUseJBlas()) {
089                                    lu = DecompositionOps.LU_JBLAS;
090                            }
091                            if (lu == null) {
092                                    lu = UJMP;
093                            }
094                            return lu.calc(source);
095                    }
096    
097                    public final Matrix solve(Matrix source, Matrix b) {
098                            LU<Matrix> lu = null;
099                            if (UJMPSettings.isUseJBlas()) {
100                                    lu = DecompositionOps.LU_JBLAS;
101                            }
102                            if (lu == null && UJMPSettings.isUseOjalgo()) {
103                                    lu = DecompositionOps.LU_OJALGO;
104                            }
105                            if (lu == null) {
106                                    lu = UJMP;
107                            }
108                            return lu.solve(source, b);
109                    }
110            };
111    
112            public static final LU<Matrix> MATRIXLARGEMULTITHREADED = new LU<Matrix>() {
113                    public final Matrix[] calc(Matrix source) {
114                            LU<Matrix> lu = null;
115                            if (UJMPSettings.isUseJBlas()) {
116                                    lu = DecompositionOps.LU_JBLAS;
117                            }
118                            if (lu == null && UJMPSettings.isUseOjalgo()) {
119                                    lu = DecompositionOps.LU_OJALGO;
120                            }
121                            if (lu == null) {
122                                    lu = UJMP;
123                            }
124                            return lu.calc(source);
125                    }
126    
127                    public final Matrix solve(Matrix source, Matrix b) {
128                            LU<Matrix> lu = null;
129                            if (UJMPSettings.isUseJBlas()) {
130                                    lu = DecompositionOps.LU_JBLAS;
131                            }
132                            if (lu == null && UJMPSettings.isUseOjalgo()) {
133                                    lu = DecompositionOps.LU_OJALGO;
134                            }
135                            if (lu == null) {
136                                    lu = UJMP;
137                            }
138                            return lu.solve(source, b);
139                    }
140            };
141    
142            final class LUMatrix {
143                    /**
144                     * Array for internal storage of decomposition.
145                     * 
146                     * @serial internal array storage.
147                     */
148                    private final double[][] LU;
149    
150                    /**
151                     * Row and column dimensions, and pivot sign.
152                     * 
153                     * @serial column dimension.
154                     * @serial row dimension.
155                     * @serial pivot sign.
156                     */
157                    private final int m, n;
158                    private int pivsign;
159    
160                    /**
161                     * Internal storage of pivot vector.
162                     * 
163                     * @serial pivot vector.
164                     */
165                    private final int[] piv;
166    
167                    /*
168                     * ------------------------ Constructor ------------------------
169                     */
170    
171                    /**
172                     * LU Decomposition
173                     * 
174                     * @param A
175                     *            Rectangular matrix
176                     * @return Structure to access L, U and piv.
177                     */
178                    public LUMatrix(Matrix A) {
179    
180                            // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
181    
182                            LU = A.toDoubleArray();
183                            m = (int) A.getRowCount();
184                            n = (int) A.getColumnCount();
185                            piv = new int[m];
186                            for (int i = 0; i < m; i++) {
187                                    piv[i] = i;
188                            }
189                            pivsign = 1;
190    
191                            final double[] LUcolj = new double[m];
192                            double[] LUrowi = null;
193    
194                            // Outer loop.
195    
196                            for (int j = 0; j < n; j++) {
197    
198                                    // Make a copy of the j-th column to localize references.
199                                    for (int i = 0; i < m; i++) {
200                                            LUcolj[i] = LU[i][j];
201                                    }
202    
203                                    // Apply previous transformations.
204    
205                                    for (int i = 0; i < m; i++) {
206                                            LUrowi = LU[i];
207    
208                                            // Most of the time is spent in the following dot product.
209    
210                                            final int kmax = Math.min(i, j);
211                                            double s = 0.0;
212                                            for (int k = 0; k < kmax; k++) {
213                                                    s += LUrowi[k] * LUcolj[k];
214                                            }
215    
216                                            LUrowi[j] = LUcolj[i] -= s;
217                                    }
218    
219                                    int p = j;
220                                    for (int i = j + 1; i < m; i++) {
221                                            if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) {
222                                                    p = i;
223                                            }
224                                    }
225                                    if (p != j) {
226                                            for (int k = 0; k < n; k++) {
227                                                    double t = LU[p][k];
228                                                    LU[p][k] = LU[j][k];
229                                                    LU[j][k] = t;
230                                            }
231                                            int k = piv[p];
232                                            piv[p] = piv[j];
233                                            piv[j] = k;
234                                            pivsign = -pivsign;
235                                    }
236    
237                                    // Compute multipliers.
238    
239                                    // http://cio.nist.gov/esd/emaildir/lists/jama/msg01498.html
240                                    if (j < m && LU[j][j] != 0.0) {
241                                            for (int i = j + 1; i < m; i++) {
242                                                    LU[i][j] /= LU[j][j];
243                                            }
244                                    }
245                            }
246                    }
247    
248                    /*
249                     * ------------------------ Temporary, experimental code.
250                     * ------------------------ *\
251                     * 
252                     * \** LU Decomposition, computed by Gaussian elimination. <P> This
253                     * constructor computes L and U with the "daxpy"-based elimination
254                     * algorithm used in LINPACK and MATLAB. In Java, we suspect the
255                     * dot-product, Crout algorithm will be faster. We have temporarily
256                     * included this constructor until timing experiments confirm this
257                     * suspicion. <P>
258                     * 
259                     * @param A Rectangular matrix
260                     * 
261                     * @param linpackflag Use Gaussian elimination. Actual value ignored.
262                     * 
263                     * @return Structure to access L, U and piv.\
264                     * 
265                     * public LUDecomposition (Matrix A, int linpackflag) { // Initialize.
266                     * LU = A.getArrayCopy(); m = A.getRowDimension(); n =
267                     * A.getColumnDimension(); piv = new int[m]; for (int i = 0; i < m; i++)
268                     * { piv[i] = i; } pivsign = 1; // Main loop. for (int k = 0; k < n;
269                     * k++) { // Find pivot. int p = k; for (int i = k+1; i < m; i++) { if
270                     * (Math.abs(LU[i][k]) > Math.abs(LU[p][k])) { p = i; } } // Exchange if
271                     * necessary. if (p != k) { for (int j = 0; j < n; j++) { double t =
272                     * LU[p][j]; LU[p][j] = LU[k][j]; LU[k][j] = t; } int t = piv[p]; piv[p]
273                     * = piv[k]; piv[k] = t; pivsign = -pivsign; } // Compute multipliers
274                     * and eliminate k-th column. if (LU[k][k] != 0.0) { for (int i = k+1; i
275                     * < m; i++) { LU[i][k] /= LU[k][k]; for (int j = k+1; j < n; j++) {
276                     * LU[i][j] -= LU[i][k]*LU[k][j]; } } } } }
277                     * 
278                     * \* ------------------------ End of temporary code.
279                     * ------------------------
280                     */
281    
282                    /*
283                     * ------------------------ Public Methods ------------------------
284                     */
285    
286                    /**
287                     * Is the matrix nonsingular?
288                     * 
289                     * @return true if U, and hence A, is nonsingular.
290                     */
291    
292                    public final boolean isNonsingular() {
293                            for (int j = 0; j < n; j++) {
294                                    if (LU[j][j] == 0)
295                                            return false;
296                            }
297                            return true;
298                    }
299    
300                    /**
301                     * Return lower triangular factor
302                     * 
303                     * @return L
304                     */
305    
306                    public final Matrix getL() {
307                            final int min = Math.min(m, n);
308                            final double[][] L = new double[m][min];
309                            for (int i = 0; i < m; i++) {
310                                    for (int j = 0; j < min; j++) {
311                                            if (i > j) {
312                                                    L[i][j] = LU[i][j];
313                                            } else if (i == j) {
314                                                    L[i][j] = 1.0;
315                                            }
316                                    }
317                            }
318                            return MatrixFactory.linkToArray(L);
319                    }
320    
321                    /**
322                     * Return upper triangular factor
323                     * 
324                     * @return U
325                     */
326    
327                    public final Matrix getU() {
328                            final int min = Math.min(m, n);
329                            final double[][] U = new double[min][n];
330                            for (int i = 0; i < min; i++) {
331                                    for (int j = 0; j < n; j++) {
332                                            if (i <= j) {
333                                                    U[i][j] = LU[i][j];
334                                            }
335                                    }
336                            }
337                            return MatrixFactory.linkToArray(U);
338                    }
339    
340                    /**
341                     * Return pivot permutation vector
342                     * 
343                     * @return piv
344                     */
345    
346                    public final int[] getPivot() {
347                            final int[] p = new int[m];
348                            for (int i = 0; i < m; i++) {
349                                    p[i] = piv[i];
350                            }
351                            return p;
352                    }
353    
354                    public final Matrix getP() {
355                            final DenseDoubleMatrix2D p = DenseDoubleMatrix2D.factory.zeros(m, m);
356                            for (int i = 0; i < m; i++) {
357                                    p.setDouble(1, i, piv[i]);
358                            }
359                            return p;
360                    }
361    
362                    /**
363                     * Return pivot permutation vector as a one-dimensional double array
364                     * 
365                     * @return (double) piv
366                     */
367    
368                    public final double[] getDoublePivot() {
369                            final double[] vals = new double[m];
370                            for (int i = 0; i < m; i++) {
371                                    vals[i] = (double) piv[i];
372                            }
373                            return vals;
374                    }
375    
376                    /**
377                     * Determinant
378                     * 
379                     * @return det(A)
380                     * @exception IllegalArgumentException
381                     *                Matrix must be square
382                     */
383    
384                    public final double det() {
385                            if (m != n) {
386                                    throw new IllegalArgumentException("Matrix must be square.");
387                            }
388                            double d = (double) pivsign;
389                            for (int j = 0; j < n; j++) {
390                                    d *= LU[j][j];
391                            }
392                            return d;
393                    }
394    
395                    /**
396                     * Solve A*X = B
397                     * 
398                     * @param B
399                     *            A Matrix with as many rows as A and any number of columns.
400                     * @return X so that L*U*X = B(piv,:)
401                     * @exception IllegalArgumentException
402                     *                Matrix row dimensions must agree.
403                     * @exception RuntimeException
404                     *                Matrix is singular.
405                     */
406    
407                    public final Matrix solve(Matrix B) {
408                            if (B.getRowCount() != m) {
409                                    throw new IllegalArgumentException("Matrix row dimensions must agree.");
410                            }
411                            if (!this.isNonsingular()) {
412                                    throw new RuntimeException("Matrix is singular.");
413                            }
414    
415                            // Copy right hand side with pivoting
416                            final int nx = (int) B.getColumnCount();
417                            // final Matrix Xmat = B.selectRows(Ret.NEW,
418                            // MathUtil.toLongArray(piv));
419    
420                            final double[][] X = new double[piv.length][(int) B.getColumnCount()];
421                            if (B instanceof DenseDoubleMatrix2D) {
422                                    DenseDoubleMatrix2D m = (DenseDoubleMatrix2D) B;
423                                    for (int c = (int) B.getColumnCount(); --c >= 0;) {
424                                            for (int r = piv.length; --r >= 0;) {
425                                                    X[r][c] = m.getDouble(piv[r], c);
426                                            }
427                                    }
428                            } else {
429                                    for (int c = (int) B.getColumnCount(); --c >= 0;) {
430                                            for (int r = piv.length; --r >= 0;) {
431                                                    X[r][c] = B.getAsDouble(piv[r], c);
432                                            }
433                                    }
434                            }
435    
436                            // Solve L*Y = B(piv,:)
437                            for (int k = 0; k < n; k++) {
438                                    for (int i = k + 1; i < n; i++) {
439                                            for (int j = 0; j < nx; j++) {
440                                                    X[i][j] -= X[k][j] * LU[i][k];
441                                            }
442                                    }
443                            }
444    
445                            for (int k = n - 1; k >= 0; k--) {
446                                    for (int j = 0; j < nx; j++) {
447                                            X[k][j] /= LU[k][k];
448                                    }
449                                    for (int i = 0; i < k; i++) {
450                                            for (int j = 0; j < nx; j++) {
451                                                    X[i][j] -= X[k][j] * LU[i][k];
452                                            }
453                                    }
454                            }
455                            return MatrixFactory.linkToArray(X);
456                    }
457    
458            };
459    }