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.jama; 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 Jama.CholeskyDecomposition; 032 import Jama.EigenvalueDecomposition; 033 import Jama.LUDecomposition; 034 import Jama.QRDecomposition; 035 import Jama.SingularValueDecomposition; 036 037 public class JamaDenseDoubleMatrix2D extends AbstractDenseDoubleMatrix2D 038 implements Wrapper<Jama.Matrix> { 039 private static final long serialVersionUID = -6065454603299978242L; 040 041 private Jama.Matrix matrix = null; 042 043 public JamaDenseDoubleMatrix2D(long... size) { 044 this.matrix = new Jama.Matrix((int) size[ROW], (int) size[COLUMN]); 045 } 046 047 public JamaDenseDoubleMatrix2D(Jama.Matrix matrix) { 048 this.matrix = matrix; 049 } 050 051 public JamaDenseDoubleMatrix2D(Matrix source) throws MatrixException { 052 this(source.getSize()); 053 for (long[] c : source.availableCoordinates()) { 054 setDouble(source.getAsDouble(c), c); 055 } 056 } 057 058 public static Jama.Matrix identity(int m, int n) { 059 Jama.Matrix A = new Jama.Matrix(m, n); 060 double[][] X = A.getArray(); 061 for (int i = 0; i < m; i++) { 062 for (int j = 0; j < n; j++) { 063 X[i][j] = (i == j ? 1.0 : 0.0); 064 } 065 } 066 return A; 067 } 068 069 public Matrix inv() throws MatrixException { 070 return new JamaDenseDoubleMatrix2D(matrix.inverse()); 071 } 072 073 public Matrix invSPD() throws MatrixException { 074 CholeskyDecomposition chol = new CholeskyDecomposition(matrix); 075 return new JamaDenseDoubleMatrix2D(chol.solve(Jama.Matrix.identity( 076 matrix.getRowDimension(), matrix.getRowDimension()))); 077 } 078 079 public Matrix[] svd() throws MatrixException { 080 if (getColumnCount() > getRowCount()) { 081 SingularValueDecomposition svd = new SingularValueDecomposition( 082 matrix.transpose()); 083 Matrix u = new JamaDenseDoubleMatrix2D(svd.getV()); 084 Matrix s = new JamaDenseDoubleMatrix2D(svd.getS().transpose()); 085 Matrix v = new JamaDenseDoubleMatrix2D(svd.getU()); 086 return new Matrix[] { u, s, v }; 087 } else { 088 SingularValueDecomposition svd = new SingularValueDecomposition( 089 matrix); 090 Matrix u = new JamaDenseDoubleMatrix2D(svd.getU()); 091 Matrix s = new JamaDenseDoubleMatrix2D(svd.getS()); 092 Matrix v = new JamaDenseDoubleMatrix2D(svd.getV()); 093 return new Matrix[] { u, s, v }; 094 } 095 } 096 097 public double getDouble(long row, long column) { 098 return matrix.get((int) row, (int) column); 099 } 100 101 public double getDouble(int row, int column) { 102 return matrix.get(row, column); 103 } 104 105 public long[] getSize() { 106 return new long[] { matrix.getRowDimension(), 107 matrix.getColumnDimension() }; 108 } 109 110 public void setDouble(double value, long row, long column) { 111 matrix.set((int) row, (int) column, value); 112 } 113 114 public void setDouble(double value, int row, int column) { 115 matrix.set(row, column, value); 116 } 117 118 public Jama.Matrix getWrappedObject() { 119 return matrix; 120 } 121 122 public void setWrappedObject(Jama.Matrix object) { 123 this.matrix = object; 124 } 125 126 public final Matrix copy() throws MatrixException { 127 Matrix m = new JamaDenseDoubleMatrix2D(matrix.copy()); 128 if (getAnnotation() != null) { 129 m.setAnnotation(getAnnotation().clone()); 130 } 131 return m; 132 } 133 134 public Matrix transpose() { 135 return new JamaDenseDoubleMatrix2D(matrix.transpose()); 136 } 137 138 public Matrix[] qr() { 139 if (getRowCount() >= getColumnCount()) { 140 QRDecomposition qr = new QRDecomposition(matrix); 141 Matrix q = new JamaDenseDoubleMatrix2D(qr.getQ()); 142 Matrix r = new JamaDenseDoubleMatrix2D(qr.getR()); 143 return new Matrix[] { q, r }; 144 } else { 145 throw new MatrixException( 146 "QR decomposition only works for matrices m>=n"); 147 } 148 } 149 150 public Matrix[] lu() { 151 LUDecomposition lu = new LUDecomposition(matrix); 152 Matrix l = new JamaDenseDoubleMatrix2D(lu.getL()); 153 Matrix u = new JamaDenseDoubleMatrix2D(lu.getU()); 154 int m = (int) getRowCount(); 155 int[] piv = lu.getPivot(); 156 Matrix p = new JamaDenseDoubleMatrix2D(m, m); 157 for (int i = 0; i < m; i++) { 158 p.setAsDouble(1, i, piv[i]); 159 } 160 return new Matrix[] { l, u, p }; 161 } 162 163 public Matrix[] eig() { 164 EigenvalueDecomposition eig = new EigenvalueDecomposition(matrix); 165 Matrix v = new JamaDenseDoubleMatrix2D(eig.getV()); 166 Matrix d = new JamaDenseDoubleMatrix2D(eig.getD()); 167 return new Matrix[] { v, d }; 168 } 169 170 public Matrix chol() { 171 CholeskyDecomposition chol = new CholeskyDecomposition(matrix); 172 Matrix r = new JamaDenseDoubleMatrix2D(chol.getL()); 173 return r; 174 } 175 176 public Matrix mtimes(Matrix m) { 177 if (m instanceof JamaDenseDoubleMatrix2D) { 178 return new JamaDenseDoubleMatrix2D(matrix 179 .times(((JamaDenseDoubleMatrix2D) m).matrix)); 180 } else { 181 return super.mtimes(m); 182 } 183 } 184 185 public Matrix times(double value) { 186 return new JamaDenseDoubleMatrix2D(matrix.times(value)); 187 } 188 189 public Matrix divide(double value) { 190 return new JamaDenseDoubleMatrix2D(matrix.times(1.0 / value)); 191 } 192 193 public double det() { 194 return matrix.det(); 195 } 196 197 public Matrix plus(Matrix m) { 198 if (m instanceof JamaDenseDoubleMatrix2D) { 199 return new JamaDenseDoubleMatrix2D(matrix 200 .plus(((JamaDenseDoubleMatrix2D) m).matrix)); 201 } else { 202 return super.plus(m); 203 } 204 } 205 206 public Matrix minus(Matrix m) { 207 if (m instanceof JamaDenseDoubleMatrix2D) { 208 return new JamaDenseDoubleMatrix2D(matrix 209 .minus(((JamaDenseDoubleMatrix2D) m).matrix)); 210 } else { 211 return super.minus(m); 212 } 213 } 214 215 public Matrix solve(Matrix b) { 216 if (b instanceof JamaDenseDoubleMatrix2D) { 217 JamaDenseDoubleMatrix2D b2 = (JamaDenseDoubleMatrix2D) b; 218 Jama.Matrix x = matrix.solve(b2.matrix); 219 return new JamaDenseDoubleMatrix2D(x); 220 } else { 221 return super.solve(b); 222 } 223 } 224 225 public Matrix solveSPD(Matrix b) { 226 if (b instanceof JamaDenseDoubleMatrix2D) { 227 JamaDenseDoubleMatrix2D b2 = (JamaDenseDoubleMatrix2D) b; 228 CholeskyDecomposition chol = new CholeskyDecomposition(matrix); 229 return new JamaDenseDoubleMatrix2D(chol.solve(b2.matrix)); 230 } else { 231 return super.solve(b); 232 } 233 } 234 235 }