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 }