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 }