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.colt;
025    
026    import org.ujmp.core.Matrix;
027    import org.ujmp.core.doublematrix.stub.AbstractDenseDoubleMatrix2D;
028    import org.ujmp.core.exceptions.MatrixException;
029    import org.ujmp.core.interfaces.Wrapper;
030    
031    import cern.colt.matrix.DoubleFactory2D;
032    import cern.colt.matrix.DoubleMatrix2D;
033    import cern.colt.matrix.impl.DenseDoubleMatrix2D;
034    import cern.colt.matrix.linalg.Algebra;
035    import cern.colt.matrix.linalg.CholeskyDecomposition;
036    import cern.colt.matrix.linalg.EigenvalueDecomposition;
037    import cern.colt.matrix.linalg.LUDecomposition;
038    import cern.colt.matrix.linalg.QRDecomposition;
039    import cern.colt.matrix.linalg.SingularValueDecomposition;
040    import cern.colt.matrix.linalg.SmpBlas;
041    import cern.jet.math.Functions;
042    
043    public class ColtDenseDoubleMatrix2D extends AbstractDenseDoubleMatrix2D
044                    implements Wrapper<DenseDoubleMatrix2D> {
045            private static final long serialVersionUID = -3223474248020842822L;
046    
047            private static ColtDenseDoubleMatrix2DFactory factory = new ColtDenseDoubleMatrix2DFactory();
048    
049            private DenseDoubleMatrix2D matrix = null;
050    
051            public ColtDenseDoubleMatrix2D(long... size) {
052                    this.matrix = new DenseDoubleMatrix2D((int) size[ROW],
053                                    (int) size[COLUMN]);
054            }
055    
056            public ColtDenseDoubleMatrix2D(DoubleMatrix2D m) {
057                    if (m instanceof DenseDoubleMatrix2D) {
058                            this.matrix = (DenseDoubleMatrix2D) m;
059                    } else {
060                            this.matrix = new DenseDoubleMatrix2D(m.toArray());
061                    }
062            }
063    
064            public ColtDenseDoubleMatrix2D(DenseDoubleMatrix2D m) {
065                    this.matrix = m;
066            }
067    
068            public ColtDenseDoubleMatrix2D(Matrix source) throws MatrixException {
069                    this(source.getSize());
070                    for (long[] c : source.availableCoordinates()) {
071                            setDouble(source.getAsDouble(c), c);
072                    }
073            }
074    
075            public double getDouble(long row, long column) {
076                    return matrix.getQuick((int) row, (int) column);
077            }
078    
079            public double getDouble(int row, int column) {
080                    return matrix.getQuick(row, column);
081            }
082    
083            public long[] getSize() {
084                    return new long[] { matrix.rows(), matrix.columns() };
085            }
086    
087            public void setDouble(double value, long row, long column) {
088                    matrix.setQuick((int) row, (int) column, value);
089            }
090    
091            public void setDouble(double value, int row, int column) {
092                    matrix.setQuick(row, column, value);
093            }
094    
095            public DenseDoubleMatrix2D getWrappedObject() {
096                    return matrix;
097            }
098    
099            public void setWrappedObject(DenseDoubleMatrix2D object) {
100                    this.matrix = object;
101            }
102    
103            public Matrix transpose() {
104                    return new ColtDenseDoubleMatrix2D((DenseDoubleMatrix2D) matrix
105                                    .viewDice().copy());
106            }
107    
108            public Matrix inv() {
109                    return new ColtDenseDoubleMatrix2D(
110                                    (DenseDoubleMatrix2D) Algebra.DEFAULT.inverse(matrix));
111            }
112    
113            public Matrix solve(Matrix b) {
114                    if (b instanceof ColtDenseDoubleMatrix2D) {
115                            ColtDenseDoubleMatrix2D b2 = (ColtDenseDoubleMatrix2D) b;
116                            if (isSquare()) {
117                                    DoubleMatrix2D ret = new LUDecomposition(matrix)
118                                                    .solve(b2.matrix);
119                                    return new ColtDenseDoubleMatrix2D(ret);
120                            } else {
121                                    DoubleMatrix2D ret = new QRDecomposition(matrix)
122                                                    .solve(b2.matrix);
123                                    return new ColtDenseDoubleMatrix2D(ret);
124                            }
125                    } else {
126                            return super.solve(b);
127                    }
128            }
129    
130            public Matrix solveSPD(Matrix b) {
131                    if (b instanceof ColtDenseDoubleMatrix2D) {
132                            ColtDenseDoubleMatrix2D b2 = (ColtDenseDoubleMatrix2D) b;
133                            DoubleMatrix2D ret = new CholeskyDecomposition(matrix)
134                                            .solve(b2.matrix);
135                            return new ColtDenseDoubleMatrix2D(ret);
136                    } else {
137                            return super.solve(b);
138                    }
139            }
140    
141            public Matrix invSPD() {
142                    DoubleMatrix2D ret = new CholeskyDecomposition(matrix)
143                                    .solve(DoubleFactory2D.dense.identity(matrix.rows()));
144                    return new ColtDenseDoubleMatrix2D(ret);
145            }
146    
147            public Matrix plus(double value) {
148                    return new ColtDenseDoubleMatrix2D((DenseDoubleMatrix2D) matrix.copy()
149                                    .assign(Functions.plus(value)));
150            }
151    
152            public Matrix minus(double value) {
153                    return new ColtDenseDoubleMatrix2D((DenseDoubleMatrix2D) matrix.copy()
154                                    .assign(Functions.minus(value)));
155            }
156    
157            public Matrix times(double value) {
158                    return new ColtDenseDoubleMatrix2D((DenseDoubleMatrix2D) matrix.copy()
159                                    .assign(Functions.mult(value)));
160            }
161    
162            public Matrix divide(double value) {
163                    return new ColtDenseDoubleMatrix2D((DenseDoubleMatrix2D) matrix.copy()
164                                    .assign(Functions.div(value)));
165            }
166    
167            public Matrix mtimes(Matrix m) {
168                    if (m instanceof ColtDenseDoubleMatrix2D) {
169                            DenseDoubleMatrix2D ret = new DenseDoubleMatrix2D(
170                                            (int) getRowCount(), (int) m.getColumnCount());
171                            matrix.zMult(((ColtDenseDoubleMatrix2D) m).matrix, ret);
172                            return new ColtDenseDoubleMatrix2D(ret);
173                    } else {
174                            return super.mtimes(m);
175                    }
176            }
177    
178            public Matrix plus(Matrix m) {
179                    if (m instanceof ColtDenseDoubleMatrix2D) {
180                            DenseDoubleMatrix2D ret = new DenseDoubleMatrix2D(
181                                            (int) getRowCount(), (int) m.getColumnCount());
182                            ret.assign(matrix);
183                            SmpBlas.smpBlas.daxpy(1, ((ColtDenseDoubleMatrix2D) m)
184                                            .getWrappedObject(), ret);
185                            return new ColtDenseDoubleMatrix2D(ret);
186                    } else {
187                            return super.plus(m);
188                    }
189            }
190    
191            public Matrix minus(Matrix m) {
192                    if (m instanceof ColtDenseDoubleMatrix2D) {
193                            DenseDoubleMatrix2D ret = new DenseDoubleMatrix2D(
194                                            (int) getRowCount(), (int) m.getColumnCount());
195                            ret.assign(matrix);
196                            SmpBlas.smpBlas.daxpy(-1, ((ColtDenseDoubleMatrix2D) m)
197                                            .getWrappedObject(), ret);
198                            return new ColtDenseDoubleMatrix2D(ret);
199                    } else {
200                            return super.minus(m);
201                    }
202            }
203    
204            public Matrix[] svd() {
205                    if (getColumnCount() > getRowCount()) {
206                            SingularValueDecomposition svd = new SingularValueDecomposition(
207                                            matrix.viewDice());
208                            Matrix u = new ColtDenseDoubleMatrix2D(svd.getV());
209                            Matrix s = new ColtDenseDoubleMatrix2D(svd.getS());
210                            Matrix v = new ColtDenseDoubleMatrix2D(svd.getU());
211                            return new Matrix[] { u, s, v };
212                    } else {
213                            SingularValueDecomposition svd = new SingularValueDecomposition(
214                                            matrix);
215                            Matrix u = new ColtDenseDoubleMatrix2D(svd.getU());
216                            Matrix s = new ColtDenseDoubleMatrix2D(svd.getS());
217                            Matrix v = new ColtDenseDoubleMatrix2D(svd.getV());
218                            return new Matrix[] { u, s, v };
219                    }
220            }
221    
222            public Matrix[] qr() {
223                    if (getColumnCount() > getRowCount()) {
224                            throw new MatrixException("matrix size must be m>=n");
225                    }
226                    QRDecomposition qr = new QRDecomposition(matrix);
227                    Matrix q = new ColtDenseDoubleMatrix2D(qr.getQ());
228                    Matrix r = new ColtDenseDoubleMatrix2D(qr.getR());
229                    return new Matrix[] { q, r };
230            }
231    
232            public Matrix[] eig() {
233                    EigenvalueDecomposition eig = new EigenvalueDecomposition(matrix);
234                    Matrix v = new ColtDenseDoubleMatrix2D(eig.getV());
235                    Matrix d = new ColtDenseDoubleMatrix2D(eig.getD());
236                    return new Matrix[] { v, d };
237            }
238    
239            public Matrix chol() {
240                    CholeskyDecomposition chol = new CholeskyDecomposition(matrix);
241                    Matrix r = new ColtDenseDoubleMatrix2D(chol.getL());
242                    return r;
243            }
244    
245            public Matrix[] lu() {
246                    if (getColumnCount() > getRowCount()) {
247                            throw new MatrixException("only supported for m>=n");
248                    }
249                    LUDecomposition lu = new LUDecomposition(matrix);
250                    Matrix l = new ColtDenseDoubleMatrix2D(lu.getL());
251                    Matrix u = new ColtDenseDoubleMatrix2D(lu.getU().viewPart(0, 0,
252                                    (int) getColumnCount(), (int) getColumnCount()));
253                    int[] piv = lu.getPivot();
254                    int m = (int) getRowCount();
255                    Matrix p = new ColtDenseDoubleMatrix2D(m, m);
256                    for (int i = 0; i < m; i++) {
257                            p.setAsDouble(1, i, piv[i]);
258                    }
259                    return new Matrix[] { l, u, p };
260            }
261    
262            public Matrix copy() {
263                    Matrix m = new ColtDenseDoubleMatrix2D((DenseDoubleMatrix2D) matrix
264                                    .copy());
265                    if (getAnnotation() != null) {
266                            m.setAnnotation(getAnnotation().clone());
267                    }
268                    return m;
269            }
270    
271            public ColtDenseDoubleMatrix2DFactory getFactory() {
272                    return factory;
273            }
274    
275    }