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.mtj;
025    
026    import java.io.IOException;
027    import java.io.ObjectInputStream;
028    import java.io.ObjectOutputStream;
029    
030    import no.uib.cipr.matrix.DenseCholesky;
031    import no.uib.cipr.matrix.DenseLU;
032    import no.uib.cipr.matrix.DenseMatrix;
033    import no.uib.cipr.matrix.EVD;
034    import no.uib.cipr.matrix.Matrices;
035    import no.uib.cipr.matrix.QR;
036    
037    import org.ujmp.core.Matrix;
038    import org.ujmp.core.calculation.Calculation.Ret;
039    import org.ujmp.core.doublematrix.stub.AbstractDenseDoubleMatrix2D;
040    import org.ujmp.core.exceptions.MatrixException;
041    import org.ujmp.core.interfaces.Wrapper;
042    import org.ujmp.mtj.calculation.Inv;
043    import org.ujmp.mtj.calculation.SVD;
044    
045    public class MTJDenseDoubleMatrix2D extends AbstractDenseDoubleMatrix2D
046                    implements Wrapper<DenseMatrix> {
047            private static final long serialVersionUID = -2386081646062313108L;
048    
049            private transient DenseMatrix matrix = null;
050    
051            public MTJDenseDoubleMatrix2D(DenseMatrix m) {
052                    this.matrix = m;
053            }
054    
055            public MTJDenseDoubleMatrix2D(no.uib.cipr.matrix.Matrix m) {
056                    this.matrix = new DenseMatrix(m);
057            }
058    
059            public MTJDenseDoubleMatrix2D(Matrix m) throws MatrixException {
060                    if (m instanceof MTJDenseDoubleMatrix2D) {
061                            this.matrix = ((MTJDenseDoubleMatrix2D) m).matrix.copy();
062                    } else {
063                            this.matrix = new DenseMatrix(m.toDoubleArray());
064                    }
065            }
066    
067            public MTJDenseDoubleMatrix2D(long... size) {
068                    this.matrix = new DenseMatrix((int) size[ROW], (int) size[COLUMN]);
069            }
070    
071            public Matrix[] svd() throws MatrixException {
072                    return SVD.INSTANCE.calc(this);
073            }
074    
075            public Matrix[] qr() throws MatrixException {
076                    if (getRowCount() >= getColumnCount()) {
077                            try {
078                                    QR qr = QR.factorize(getWrappedObject());
079                                    Matrix q = new MTJDenseDoubleMatrix2D(qr.getQ());
080                                    Matrix r = new MTJDenseDoubleMatrix2D(qr.getR());
081                                    return new Matrix[] { q, r };
082                            } catch (Exception e) {
083                                    throw new MatrixException(e);
084                            }
085                    } else {
086                            throw new MatrixException("only allowed for matrices m>=n");
087                    }
088            }
089    
090            public Matrix[] lu() throws MatrixException {
091                    try {
092                            DenseLU lu = DenseLU.factorize(getWrappedObject());
093                            Matrix l = new MTJDenseDoubleMatrix2D(lu.getL());
094                            Matrix u = new MTJDenseDoubleMatrix2D(lu.getU());
095                            int m = (int) getRowCount();
096                            int[] piv = lu.getPivots();
097                            Matrix p = new MTJDenseDoubleMatrix2D(m, m);
098                            // pivots seem to be broken
099                            // http://code.google.com/p/matrix-toolkits-java/issues/detail?id=1
100                            // for (int i = 0; i < m; i++) {
101                            // p.setAsDouble(1, i, piv[i]);
102                            // }
103                            p.eye(Ret.ORIG);
104                            return new Matrix[] { l, u, p };
105                    } catch (Exception e) {
106                            throw new MatrixException(e);
107                    }
108            }
109    
110            public Matrix chol() throws MatrixException {
111                    try {
112                            DenseCholesky chol = DenseCholesky.factorize(getWrappedObject());
113                            Matrix l = new MTJDenseDoubleMatrix2D(chol.getL());
114                            return l;
115                    } catch (Exception e) {
116                            throw new MatrixException(e);
117                    }
118            }
119    
120            public Matrix[] eig() throws MatrixException {
121                    try {
122                            EVD evd = EVD.factorize(getWrappedObject());
123                            Matrix v = new MTJDenseDoubleMatrix2D(evd.getRightEigenvectors());
124                            int m = (int) getRowCount();
125                            double[] evds = evd.getRealEigenvalues();
126                            Matrix d = new MTJDenseDoubleMatrix2D(m, m);
127                            for (int i = 0; i < m; i++) {
128                                    d.setAsDouble(evds[i], i, i);
129                            }
130                            return new Matrix[] { v, d };
131                    } catch (Exception e) {
132                            throw new MatrixException(e);
133                    }
134            }
135    
136            public double getDouble(long row, long column) {
137                    return matrix.getData()[(int) (row + column * matrix.numRows())];
138            }
139    
140            public double getDouble(int row, int column) {
141                    return matrix.getData()[(row + column * matrix.numRows())];
142            }
143    
144            public long[] getSize() {
145                    return new long[] { matrix.numRows(), matrix.numColumns() };
146            }
147    
148            public void setDouble(double value, long row, long column) {
149                    matrix.getData()[(int) (row + column * matrix.numRows())] = value;
150            }
151    
152            public void setDouble(double value, int row, int column) {
153                    matrix.getData()[(row + column * matrix.numRows())] = value;
154            }
155    
156            public Matrix transpose() {
157                    DenseMatrix ret = new DenseMatrix((int) getColumnCount(),
158                                    (int) getRowCount());
159                    return new MTJDenseDoubleMatrix2D((DenseMatrix) matrix.transpose(ret));
160            }
161    
162            public Matrix inv() {
163                    return new Inv(this).calcNew();
164            }
165    
166            public DenseMatrix getWrappedObject() {
167                    return matrix;
168            }
169    
170            public void setWrappedObject(DenseMatrix object) {
171                    this.matrix = object;
172            }
173    
174            private void readObject(ObjectInputStream s) throws IOException,
175                            ClassNotFoundException {
176                    s.defaultReadObject();
177                    double[][] values = (double[][]) s.readObject();
178                    matrix = new DenseMatrix(values);
179            }
180    
181            private void writeObject(ObjectOutputStream s) throws IOException,
182                            MatrixException {
183                    s.defaultWriteObject();
184                    s.writeObject(this.toDoubleArray());
185            }
186    
187            public Matrix mtimes(Matrix m2) throws MatrixException {
188                    if (m2 instanceof MTJDenseDoubleMatrix2D) {
189                            DenseMatrix a = matrix;
190                            DenseMatrix b = ((MTJDenseDoubleMatrix2D) m2).getWrappedObject();
191                            DenseMatrix c = new DenseMatrix(a.numRows(), b.numColumns());
192                            try {
193                                    a.mult(b, c);
194                                    return new MTJDenseDoubleMatrix2D(c);
195                            } catch (Exception e) {
196                                    // sometimes BLAS cannot be found. Don't know why. Use direct
197                                    // method instead.
198                                    double[] Bd = ((DenseMatrix) b).getData(), Cd = ((DenseMatrix) c)
199                                                    .getData();
200                                    org.netlib.blas.Dgemm.dgemm("N", "N", c.numRows(), c
201                                                    .numColumns(), a.numColumns(), 1, a.getData(), 0, Math
202                                                    .max(1, a.numRows()), Bd, 0, Math.max(1, b.numRows()),
203                                                    1, Cd, 0, Math.max(1, c.numRows()));
204                                    return new MTJDenseDoubleMatrix2D(c);
205                            }
206                    }
207                    return super.mtimes(m2);
208            }
209    
210            public Matrix plus(Matrix m2) throws MatrixException {
211                    if (m2 instanceof MTJDenseDoubleMatrix2D) {
212                            DenseMatrix ret = matrix.copy();
213                            ret.add(((MTJDenseDoubleMatrix2D) m2).getWrappedObject());
214                            return new MTJDenseDoubleMatrix2D(ret);
215                    } else {
216                            return super.plus(m2);
217                    }
218            }
219    
220            public Matrix times(double f) throws MatrixException {
221                    DenseMatrix ret = matrix.copy();
222                    ret.scale(f);
223                    return new MTJDenseDoubleMatrix2D(ret);
224            }
225    
226            public Matrix divide(double f) throws MatrixException {
227                    DenseMatrix ret = matrix.copy();
228                    ret.scale(1.0 / f);
229                    return new MTJDenseDoubleMatrix2D(ret);
230            }
231    
232            public Matrix copy() {
233                    Matrix m = new MTJDenseDoubleMatrix2D(matrix.copy());
234                    if (getAnnotation() != null) {
235                            m.setAnnotation(getAnnotation().clone());
236                    }
237                    return m;
238            }
239    
240            public Matrix solve(Matrix b) {
241                    if (b instanceof MTJDenseDoubleMatrix2D) {
242                            MTJDenseDoubleMatrix2D b2 = (MTJDenseDoubleMatrix2D) b;
243                            DenseMatrix x = new DenseMatrix((int) getColumnCount(), (int) b2
244                                            .getColumnCount());
245                            matrix.solve(b2.matrix, x);
246                            return new MTJDenseDoubleMatrix2D(x);
247                    } else {
248                            return super.solve(b);
249                    }
250            }
251    
252            public Matrix invSPD() {
253                    DenseCholesky chol = DenseCholesky.factorize(getWrappedObject());
254                    return new MTJDenseDoubleMatrix2D(chol.solve(Matrices.identity(matrix
255                                    .numRows())));
256            }
257    }