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