001    /*
002     * Copyright (C) 2008-2010 by Holger Arndt
003     *
004     * This file is part of the Universal Java Matrix Package (UJMP).
005     * See the NOTICE file distributed with this work for additional
006     * information regarding copyright ownership and licensing.
007     *
008     * UJMP is free software; you can redistribute it and/or modify
009     * it under the terms of the GNU Lesser General Public License as
010     * published by the Free Software Foundation; either version 2
011     * of the License, or (at your option) any later version.
012     *
013     * UJMP is distributed in the hope that it will be useful,
014     * but WITHOUT ANY WARRANTY; without even the implied warranty of
015     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
016     * GNU Lesser General Public License for more details.
017     *
018     * You should have received a copy of the GNU Lesser General Public
019     * License along with UJMP; if not, write to the
020     * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
021     * Boston, MA  02110-1301  USA
022     */
023    
024    package org.ujmp.core.doublematrix.calculation.general.decomposition;
025    
026    import org.ujmp.core.Matrix;
027    import org.ujmp.core.MatrixFactory;
028    import org.ujmp.core.util.DecompositionOps;
029    import org.ujmp.core.util.UJMPSettings;
030    
031    /**
032     * Cholesky Decomposition.
033     * <P>
034     * For a symmetric, positive definite matrix A, the Cholesky decomposition is an
035     * lower triangular matrix L so that A = L*L'.
036     * <P>
037     * If the matrix is not symmetric or positive definite, the constructor returns
038     * a partial decomposition and sets an internal flag that may be queried by the
039     * isSPD() method.
040     */
041    
042    public interface Chol<T> {
043    
044            public T calc(T source);
045    
046            public T solve(T source, T b);
047    
048            public static int THRESHOLD = 100;
049    
050            public static final Chol<Matrix> MATRIX = new Chol<Matrix>() {
051    
052                    public final Matrix calc(Matrix source) {
053                            if (UJMPSettings.getNumberOfThreads() == 1) {
054                                    if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) {
055                                            return MATRIXLARGESINGLETHREADED.calc(source);
056                                    } else {
057                                            return MATRIXSMALLSINGLETHREADED.calc(source);
058                                    }
059                            } else {
060                                    if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) {
061                                            return MATRIXLARGEMULTITHREADED.calc(source);
062                                    } else {
063                                            return MATRIXSMALLMULTITHREADED.calc(source);
064                                    }
065                            }
066                    }
067    
068                    public final Matrix solve(Matrix source, Matrix b) {
069                            if (UJMPSettings.getNumberOfThreads() == 1) {
070                                    if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) {
071                                            return MATRIXLARGESINGLETHREADED.solve(source, b);
072                                    } else {
073                                            return MATRIXSMALLSINGLETHREADED.solve(source, b);
074                                    }
075                            } else {
076                                    if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) {
077                                            return MATRIXLARGEMULTITHREADED.solve(source, b);
078                                    } else {
079                                            return MATRIXSMALLMULTITHREADED.solve(source, b);
080                                    }
081                            }
082                    }
083            };
084    
085            public static final Chol<Matrix> INSTANCE = MATRIX;
086    
087            public static final Chol<Matrix> UJMP = new Chol<Matrix>() {
088    
089                    public final Matrix calc(Matrix source) {
090                            CholMatrix chol = new CholMatrix(source);
091                            return chol.getL();
092                    }
093    
094                    public final Matrix solve(Matrix source, Matrix b) {
095                            CholMatrix chol = new CholMatrix(source);
096                            return chol.solve(b);
097                    }
098            };
099    
100            public static final Chol<Matrix> MATRIXSMALLMULTITHREADED = UJMP;
101    
102            public static final Chol<Matrix> MATRIXSMALLSINGLETHREADED = UJMP;
103    
104            public static final Chol<Matrix> MATRIXLARGESINGLETHREADED = new Chol<Matrix>() {
105                    public Matrix calc(Matrix source) {
106                            Chol<Matrix> chol = null;
107                            if (UJMPSettings.isUseJBlas()) {
108                                    chol = DecompositionOps.CHOL_JBLAS;
109                            }
110                            if (chol == null) {
111                                    chol = UJMP;
112                            }
113                            return chol.calc(source);
114                    }
115    
116                    public final Matrix solve(Matrix source, Matrix b) {
117                            Chol<Matrix> chol = null;
118                            if (UJMPSettings.isUseJBlas()) {
119                                    chol = DecompositionOps.CHOL_JBLAS;
120                            }
121                            if (chol == null) {
122                                    chol = UJMP;
123                            }
124                            return chol.solve(source, b);
125                    }
126            };
127    
128            public static final Chol<Matrix> MATRIXLARGEMULTITHREADED = new Chol<Matrix>() {
129                    public Matrix calc(Matrix source) {
130                            Chol<Matrix> chol = null;
131                            if (UJMPSettings.isUseJBlas()) {
132                                    chol = DecompositionOps.CHOL_JBLAS;
133                            }
134                            if (chol == null && UJMPSettings.isUseOjalgo()) {
135                                    chol = DecompositionOps.CHOL_OJALGO;
136                            }
137                            if (chol == null) {
138                                    chol = UJMP;
139                            }
140                            return chol.calc(source);
141                    }
142    
143                    public final Matrix solve(Matrix source, Matrix b) {
144                            Chol<Matrix> chol = null;
145                            if (UJMPSettings.isUseJBlas()) {
146                                    chol = DecompositionOps.CHOL_JBLAS;
147                            }
148                            if (chol == null && UJMPSettings.isUseOjalgo()) {
149                                    chol = DecompositionOps.CHOL_OJALGO;
150                            }
151                            if (chol == null) {
152                                    chol = UJMP;
153                            }
154                            return chol.solve(source, b);
155                    }
156            };
157    
158            final class CholMatrix {
159                    private static final long serialVersionUID = 400514872358216115L;
160    
161                    /**
162                     * Array for internal storage of decomposition.
163                     * 
164                     * @serial internal array storage.
165                     */
166                    private final double[][] L;
167    
168                    /**
169                     * Row and column dimension (square matrix).
170                     * 
171                     * @serial matrix dimension.
172                     */
173                    private final int n;
174    
175                    /**
176                     * Symmetric and positive definite flag.
177                     * 
178                     * @serial is symmetric and positive definite flag.
179                     */
180                    private boolean isspd;
181    
182                    /*
183                     * ------------------------ Constructor ------------------------
184                     */
185    
186                    /**
187                     * Cholesky algorithm for symmetric and positive definite matrix.
188                     * 
189                     * @param A
190                     *            Square, symmetric matrix.
191                     * @return Structure to access L and isspd flag.
192                     */
193    
194                    public CholMatrix(Matrix Arg) {
195                            final double[][] A = Arg.toDoubleArray();
196                            n = (int) Arg.getRowCount();
197                            L = new double[n][n];
198                            isspd = (Arg.getColumnCount() == n);
199                            // Main loop.
200                            double[] Lrowj = null;
201                            double[] Lrowk = null;
202                            double[] Aj = null;
203                            for (int j = 0; j < n; j++) {
204                                    Lrowj = L[j];
205                                    Aj = A[j];
206                                    double d = 0.0;
207                                    for (int k = 0; k < j; k++) {
208                                            Lrowk = L[k];
209                                            double s = 0.0;
210                                            for (int i = 0; i < k; i++) {
211                                                    s += Lrowk[i] * Lrowj[i];
212                                            }
213                                            Lrowj[k] = s = (Aj[k] - s) / Lrowk[k];
214                                            d = d + s * s;
215                                            isspd = isspd & (A[k][j] == Aj[k]);
216                                    }
217                                    d = Aj[j] - d;
218                                    isspd = isspd & (d > 0.0);
219                                    Lrowj[j] = Math.sqrt(Math.max(d, 0.0));
220                                    for (int k = j + 1; k < n; k++) {
221                                            Lrowj[k] = 0.0;
222                                    }
223                            }
224                    }
225    
226                    /*
227                     * ------------------------ Temporary, experimental code.
228                     * ------------------------ *\
229                     * 
230                     * \** Right Triangular Cholesky Decomposition. <P> For a symmetric,
231                     * positive definite matrix A, the Right Cholesky decomposition is an
232                     * upper triangular matrix R so that A = R'*R. This constructor computes
233                     * R with the Fortran inspired column oriented algorithm used in LINPACK
234                     * and MATLAB. In Java, we suspect a row oriented, lower triangular
235                     * decomposition is faster. We have temporarily included this
236                     * constructor here until timing experiments confirm this suspicion.\
237                     * 
238                     * \** Array for internal storage of right triangular decomposition. **\
239                     * private transient double[][] R;
240                     * 
241                     * \** Cholesky algorithm for symmetric and positive definite matrix.
242                     * 
243                     * @param A Square, symmetric matrix.
244                     * 
245                     * @param rightflag Actual value ignored.
246                     * 
247                     * @return Structure to access R and isspd flag.\
248                     * 
249                     * public CholeskyDecomposition (Matrix Arg, int rightflag) { //
250                     * Initialize. double[][] A = Arg.getArray(); n =
251                     * Arg.getColumnDimension(); R = new double[n][n]; isspd =
252                     * (Arg.getColumnDimension() == n); // Main loop. for (int j = 0; j < n;
253                     * j++) { double d = 0.0; for (int k = 0; k < j; k++) { double s =
254                     * A[k][j]; for (int i = 0; i < k; i++) { s = s - R[i][k]*R[i][j]; }
255                     * R[k][j] = s = s/R[k][k]; d = d + s*s; isspd = isspd & (A[k][j] ==
256                     * A[j][k]); } d = A[j][j] - d; isspd = isspd & (d > 0.0); R[j][j] =
257                     * Math.sqrt(Math.max(d,0.0)); for (int k = j+1; k < n; k++) { R[k][j] =
258                     * 0.0; } } }
259                     * 
260                     * \** Return upper triangular factor.
261                     * 
262                     * @return R\
263                     * 
264                     * public Matrix getR () { return new Matrix(R,n,n); }
265                     * 
266                     * \* ------------------------ End of temporary code.
267                     * ------------------------
268                     */
269    
270                    /*
271                     * ------------------------ Public Methods ------------------------
272                     */
273    
274                    /**
275                     * Is the matrix symmetric and positive definite?
276                     * 
277                     * @return true if A is symmetric and positive definite.
278                     */
279    
280                    public final boolean isSPD() {
281                            return isspd;
282                    }
283    
284                    /**
285                     * Return triangular factor.
286                     * 
287                     * @return L
288                     */
289    
290                    public final Matrix getL() {
291                            return MatrixFactory.linkToArray(L);
292                    }
293    
294                    /**
295                     * Solve A*X = B
296                     * 
297                     * @param B
298                     *            A Matrix with as many rows as A and any number of columns.
299                     * @return X so that L*L'*X = B
300                     * @exception IllegalArgumentException
301                     *                Matrix row dimensions must agree.
302                     * @exception RuntimeException
303                     *                Matrix is not symmetric positive definite.
304                     */
305    
306                    public final Matrix solve(Matrix B) {
307                            if (B.getRowCount() != n) {
308                                    throw new IllegalArgumentException("Matrix row dimensions must agree.");
309                            }
310                            if (!isspd) {
311                                    throw new RuntimeException("Matrix is not symmetric positive definite.");
312                            }
313    
314                            // Copy right hand side.
315                            final double[][] X = B.toDoubleArray();
316                            final int nx = (int) B.getColumnCount();
317    
318                            // Solve L*Y = B;
319                            for (int k = 0; k < n; k++) {
320                                    for (int j = 0; j < nx; j++) {
321                                            for (int i = 0; i < k; i++) {
322                                                    X[k][j] -= X[i][j] * L[k][i];
323                                            }
324                                            X[k][j] /= L[k][k];
325                                    }
326                            }
327    
328                            // Solve L'*X = Y;
329                            for (int k = n - 1; k >= 0; k--) {
330                                    for (int j = 0; j < nx; j++) {
331                                            for (int i = k + 1; i < n; i++) {
332                                                    X[k][j] -= X[i][j] * L[i][k];
333                                            }
334                                            X[k][j] /= L[k][k];
335                                    }
336                            }
337    
338                            return MatrixFactory.linkToArray(X);
339                    }
340    
341            }
342    }