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.jmatrices; 025 026 import org.jmatrices.dbl.MatrixFactory; 027 import org.jmatrices.dbl.decomposition.CholeskyDecomposition; 028 import org.jmatrices.dbl.decomposition.EigenvalueDecomposition; 029 import org.jmatrices.dbl.decomposition.LUDecomposition; 030 import org.jmatrices.dbl.decomposition.QRDecomposition; 031 import org.jmatrices.dbl.decomposition.SingularValueDecomposition; 032 import org.jmatrices.dbl.operator.MatrixOperator; 033 import org.jmatrices.dbl.transformer.MatrixTransformer; 034 import org.ujmp.core.Coordinates; 035 import org.ujmp.core.Matrix; 036 import org.ujmp.core.doublematrix.stub.AbstractDenseDoubleMatrix2D; 037 import org.ujmp.core.exceptions.MatrixException; 038 import org.ujmp.core.interfaces.Wrapper; 039 040 public class JMatricesDenseDoubleMatrix2D extends AbstractDenseDoubleMatrix2D 041 implements Wrapper<org.jmatrices.dbl.Matrix> { 042 private static final long serialVersionUID = 513251881654621L; 043 044 private org.jmatrices.dbl.Matrix matrix = null; 045 046 public JMatricesDenseDoubleMatrix2D(org.jmatrices.dbl.Matrix matrix) { 047 this.matrix = matrix; 048 } 049 050 public JMatricesDenseDoubleMatrix2D(long... size) { 051 if (size[ROW] > 0 && size[COLUMN] > 0) { 052 this.matrix = MatrixFactory.getMatrix((int) size[ROW], 053 (int) size[COLUMN], null); 054 } 055 } 056 057 public JMatricesDenseDoubleMatrix2D(org.ujmp.core.Matrix source) 058 throws MatrixException { 059 this(source.getSize()); 060 for (long[] c : source.availableCoordinates()) { 061 setDouble(source.getAsDouble(c), c); 062 } 063 } 064 065 public double getDouble(long row, long column) { 066 return matrix.get((int) row + 1, (int) column + 1); 067 } 068 069 public double getDouble(int row, int column) { 070 return matrix.get(row + 1, column + 1); 071 } 072 073 public long[] getSize() { 074 return matrix == null ? Coordinates.ZERO2D : new long[] { 075 matrix.rows(), matrix.cols() }; 076 } 077 078 public void setDouble(double value, long row, long column) { 079 matrix.set((int) row + 1, (int) column + 1, value); 080 } 081 082 public void setDouble(double value, int row, int column) { 083 matrix.set(row + 1, column + 1, value); 084 } 085 086 public org.jmatrices.dbl.Matrix getWrappedObject() { 087 return matrix; 088 } 089 090 public void setWrappedObject(org.jmatrices.dbl.Matrix object) { 091 this.matrix = object; 092 } 093 094 public Matrix transpose() { 095 return new JMatricesDenseDoubleMatrix2D(MatrixTransformer 096 .transpose(matrix)); 097 } 098 099 public Matrix inv() { 100 return new JMatricesDenseDoubleMatrix2D(MatrixTransformer 101 .inverse(matrix)); 102 } 103 104 public Matrix[] eig() { 105 EigenvalueDecomposition evd = new EigenvalueDecomposition(matrix); 106 Matrix v = new JMatricesDenseDoubleMatrix2D(evd.getV()); 107 Matrix d = new JMatricesDenseDoubleMatrix2D(evd.getD()); 108 return new Matrix[] { v, d }; 109 } 110 111 public Matrix[] qr() { 112 if (getRowCount() >= getColumnCount()) { 113 QRDecomposition qr = new QRDecomposition(matrix); 114 Matrix q = new JMatricesDenseDoubleMatrix2D(qr.getQ()); 115 Matrix r = new JMatricesDenseDoubleMatrix2D(qr.getR()); 116 return new Matrix[] { q, r }; 117 } else { 118 throw new MatrixException("only allowed for matrices m>=n"); 119 } 120 } 121 122 public Matrix[] svd() { 123 if (isSquare()) { 124 SingularValueDecomposition qr = new SingularValueDecomposition( 125 matrix); 126 Matrix u = new JMatricesDenseDoubleMatrix2D(qr.getU()); 127 Matrix s = new JMatricesDenseDoubleMatrix2D(qr.getS()); 128 Matrix v = new JMatricesDenseDoubleMatrix2D(qr.getV()); 129 return new Matrix[] { u, s, v }; 130 } else { 131 throw new MatrixException("only allowed for square matrices"); 132 } 133 } 134 135 public Matrix chol() { 136 CholeskyDecomposition chol = new CholeskyDecomposition(matrix); 137 Matrix r = new JMatricesDenseDoubleMatrix2D(chol.getL()); 138 return r; 139 } 140 141 // some error 142 // public Matrix invSPD() throws MatrixException { 143 // CholeskyDecomposition chol = new CholeskyDecomposition(matrix); 144 // org.jmatrices.dbl.Matrix eye = MatrixFactory.getIdentityMatrix(matrix 145 // .rows(), null); 146 // return new JMatricesDenseDoubleMatrix2D(chol.solve(eye)); 147 // } 148 149 public Matrix[] lu() { 150 if (getRowCount() >= getColumnCount()) { 151 LUDecomposition lu = new LUDecomposition(matrix); 152 Matrix l = new JMatricesDenseDoubleMatrix2D(lu.getL()); 153 Matrix u = new JMatricesDenseDoubleMatrix2D(lu.getU().getSubMatrix( 154 1, 1, (int) getColumnCount(), (int) getColumnCount())); 155 int m = (int) getRowCount(); 156 int[] piv = lu.getPivot(); 157 Matrix p = new JMatricesDenseDoubleMatrix2D(m, m); 158 for (int i = 0; i < m; i++) { 159 p.setAsDouble(1, i, piv[i] - 1); 160 } 161 return new Matrix[] { l, u, p }; 162 } else { 163 throw new MatrixException("only allowed for matrices m>=n"); 164 } 165 } 166 167 public Matrix mtimes(Matrix m2) { 168 if (m2 instanceof JMatricesDenseDoubleMatrix2D) { 169 return new JMatricesDenseDoubleMatrix2D(MatrixOperator.multiply( 170 matrix, ((JMatricesDenseDoubleMatrix2D) m2).matrix)); 171 } else { 172 return super.mtimes(m2); 173 } 174 } 175 176 public Matrix solve(Matrix b) { 177 if (b instanceof JMatricesDenseDoubleMatrix2D) { 178 JMatricesDenseDoubleMatrix2D b2 = (JMatricesDenseDoubleMatrix2D) b; 179 if (isSquare()) { 180 org.jmatrices.dbl.Matrix x = new LUDecomposition(matrix) 181 .solve(b2.matrix); 182 return new JMatricesDenseDoubleMatrix2D(x); 183 } else { 184 org.jmatrices.dbl.Matrix x = new QRDecomposition(matrix) 185 .solve(b2.matrix); 186 return new JMatricesDenseDoubleMatrix2D(x); 187 } 188 } else { 189 return super.solve(b); 190 } 191 } 192 }