001/* 002 * JScience - Java(TM) Tools and Libraries for the Advancement of Sciences. 003 * Copyright (C) 2006 - JScience (http://jscience.org/) 004 * All rights reserved. 005 * 006 * Permission to use, copy, modify, and distribute this software is 007 * freely granted, provided that this notice is preserved. 008 */ 009package org.jscience.mathematics.vector; 010 011import java.util.Iterator; 012import java.util.List; 013 014import javolution.context.ObjectFactory; 015import javolution.lang.MathLib; 016import javolution.util.FastComparator; 017import javolution.util.FastMap; 018import javolution.util.FastTable; 019import javolution.util.Index; 020 021import org.jscience.mathematics.structure.Field; 022 023/** 024 * <p> This class represents a matrix made of {@link SparseVector sparse 025 * vectors} (as rows). To create a sparse matrix made of column vectors the 026 * {@link #transpose} method can be used. 027 * For example:[code] 028 * SparseVector<Rational> column0 = SparseVector.valueOf(...); 029 * SparseVector<Rational> column1 = SparseVector.valueOf(...); 030 * SparseMatrix<Rational> M = SparseMatrix.valueOf(column0, column1).transpose(); 031 * [/code]</p> 032 * <p> As for any concrete {@link org.jscience.mathematics.structure.Structure 033 * structure}, this class is declared <code>final</code> (otherwise most 034 * operations would have to be overridden to return the appropriate type). 035 * Specialized dense matrix should sub-class {@link Matrix} directly. 036 * For example:[code] 037 * // Extension through composition. 038 * final class BandMatrix <F extends Field<F>> extends Matrix<F> { 039 * private SparseMatrix<F> _value; 040 * ... 041 * public BandMatrix opposite() { // Returns the right type. 042 * return BandMatrix.valueOf(_value.opposite()); 043 * } 044 * ... 045 * }[/code] 046 * </p> 047 * @author <a href="mailto:jean-marie@dautelle.com">Jean-Marie Dautelle</a> 048 * @version 3.3, January 2, 2007 049 */ 050public final class SparseMatrix<F extends Field<F>> extends Matrix<F> { 051 052 /** 053 * Holds the number of columns n or the number of rows m if transposed. 054 */ 055 int _n; 056 057 /** 058 * Holds the zero. 059 */ 060 F _zero;; 061 062 /** 063 * Indicates if this matrix is transposed (the rows are then the columns). 064 */ 065 boolean _transposed; 066 067 /** 068 * Holds this matrix rows (or columns when transposed). 069 */ 070 final FastTable<SparseVector<F>> _rows = new FastTable<SparseVector<F>>(); 071 072 /** 073 * Returns the sparse square matrix having the specified diagonal 074 * vector. This method is typically used to create an identity matrix. 075 * For example:[code] 076 * SparseMatrix<Real> IDENTITY = Matrix.valueOf( 077 * DenseVector.valueOf({Real.ONE, Real.ONE, Real.ONE}), Real.ZERO); 078 * [/code] 079 * 080 * @param diagonal the diagonal vector. 081 * @param zero value of non-diagonal elements. 082 * @return a square matrix with <code>diagonal</code> on the diagonal and 083 * <code>zero</code> elsewhere. 084 */ 085 public static <F extends Field<F>> SparseMatrix<F> valueOf( 086 Vector<F> diagonal, F zero) { 087 int n = diagonal.getDimension(); 088 SparseMatrix<F> M = SparseMatrix.newInstance(n, zero, false); 089 for (int i = 0; i < n; i++) { 090 SparseVector<F> row = SparseVector.valueOf(n, zero, i, diagonal 091 .get(i)); 092 M._rows.add(row); 093 } 094 return M; 095 } 096 097 /** 098 * Returns a sparse matrix holding the specified row vectors 099 * (column vectors if {@link #transpose transposed}). 100 * 101 * @param rows the row vectors. 102 * @return the matrix having the specified rows. 103 * @throws DimensionException if the rows do not have the same dimension. 104 */ 105 public static <F extends Field<F>> SparseMatrix<F> valueOf( 106 SparseVector<F>... rows) { 107 final int n = rows[0]._dimension; 108 final F zero = rows[0]._zero; 109 SparseMatrix<F> M = SparseMatrix.newInstance(n, zero, false); 110 for (int i = 0, m = rows.length; i < m; i++) { 111 SparseVector<F> rowi = rows[i]; 112 if (rowi._dimension != n) 113 throw new DimensionException( 114 "All vectors must have the same dimension."); 115 if (!zero.equals(rowi._zero)) 116 throw new DimensionException( 117 "All vectors must have the same zero element."); 118 M._rows.add(rowi); 119 } 120 return M; 121 } 122 123 /** 124 * Returns a sparse matrix holding the row vectors from the specified 125 * collection (column vectors if {@link #transpose transposed}). 126 * 127 * @param rows the list of row vectors. 128 * @return the matrix having the specified rows. 129 * @throws DimensionException if the rows do not have the same dimension. 130 */ 131 public static <F extends Field<F>> SparseMatrix<F> valueOf( 132 List<SparseVector<F>> rows) { 133 final int n = rows.get(0)._dimension; 134 final F zero = rows.get(0)._zero; 135 SparseMatrix<F> M = SparseMatrix.newInstance(n, zero, false); 136 Iterator<SparseVector<F>> iterator = rows.iterator(); 137 for (int i = 0, m = rows.size(); i < m; i++) { 138 SparseVector<F> rowi = iterator.next(); 139 if (rowi.getDimension() != n) 140 throw new DimensionException( 141 "All vectors must have the same dimension."); 142 if (!zero.equals(rowi._zero)) 143 throw new DimensionException( 144 "All vectors must have the same zero element."); 145 M._rows.add(rowi); 146 } 147 return M; 148 } 149 150 /** 151 * Returns a sparse matrix equivalent to the specified matrix but with 152 * the zero elements removed using the default object equality comparator. 153 * 154 * @param that the matrix to convert. 155 * @param zero the zero element for the sparse vector to return. 156 * @return <code>SparseMatrix.valueOf(that, zero, FastComparator.DEFAULT)</code> or a dense matrix holding the same elements 157 */ 158 public static <F extends Field<F>> SparseMatrix<F> valueOf(Matrix<F> that, F zero) { 159 return SparseMatrix.valueOf(that, zero, FastComparator.DEFAULT); 160 } 161 162 /** 163 * Returns a sparse matrix equivalent to the specified matrix but with 164 * the zero elements removed using the specified object equality comparator. 165 * 166 * @param that the matrix to convert. 167 * @param zero the zero element for the sparse vector to return. 168 * @param comparator the comparator used to determinate zero equality. 169 * @return <code>that</code> or a dense matrix holding the same elements 170 * as the specified matrix. 171 */ 172 public static <F extends Field<F>> SparseMatrix<F> valueOf(Matrix<F> that, 173 F zero, FastComparator<? super F> comparator) { 174 if (that instanceof SparseMatrix) 175 return (SparseMatrix<F>) that; 176 int n = that.getNumberOfColumns(); 177 int m = that.getNumberOfRows(); 178 SparseMatrix<F> M = SparseMatrix.newInstance(n, zero, false); 179 for (int i = 0; i < m; i++) { 180 SparseVector<F> rowi = SparseVector.valueOf(that.getRow(i), zero, 181 comparator); 182 M._rows.add(rowi); 183 } 184 return M; 185 } 186 187 /** 188 * Returns the value of the non-set elements for this sparse matrix. 189 * 190 * @return the element corresponding to zero. 191 */ 192 public F getZero() { 193 return _zero; 194 } 195 196 @Override 197 public int getNumberOfRows() { 198 return _transposed ? _n : _rows.size(); 199 } 200 201 @Override 202 public int getNumberOfColumns() { 203 return _transposed ? _rows.size() : _n; 204 } 205 206 @Override 207 public F get(int i, int j) { 208 return _transposed ? _rows.get(j).get(i) : _rows.get(i).get(j); 209 } 210 211 @Override 212 public SparseVector<F> getRow(int i) { 213 if (!_transposed) 214 return _rows.get(i); 215 // Else transposed. 216 int n = _rows.size(); 217 int m = _n; 218 if ((i < 0) || (i >= m)) 219 throw new DimensionException(); 220 SparseVector<F> V = SparseVector.newInstance(n, _zero); 221 for (int j = 0; j < n; j++) { 222 SparseVector<F> row = _rows.get(j); 223 F e = row._elements.get(Index.valueOf(i)); 224 if (e != null) { 225 V._elements.put(Index.valueOf(j), e); 226 } 227 } 228 return V; 229 } 230 231 @Override 232 public SparseVector<F> getColumn(int j) { 233 if (_transposed) 234 return _rows.get(j); 235 int m = _rows.size(); 236 if ((j < 0) || (j >= _n)) 237 throw new DimensionException(); 238 SparseVector<F> V = SparseVector.newInstance(_n, _zero); 239 for (int i = 0; i < m; i++) { 240 SparseVector<F> row = _rows.get(i); 241 F e = row._elements.get(Index.valueOf(j)); 242 if (e != null) { 243 V._elements.put(Index.valueOf(i), e); 244 } 245 } 246 return V; 247 } 248 249 @Override 250 public SparseVector<F> getDiagonal() { 251 int m = this.getNumberOfRows(); 252 int n = this.getNumberOfColumns(); 253 int dimension = MathLib.min(m, n); 254 SparseVector<F> V = SparseVector.newInstance(_n, _zero); 255 for (int i = 0; i < dimension; i++) { 256 SparseVector<F> row = _rows.get(i); 257 F e = row._elements.get(Index.valueOf(i)); 258 if (e != null) { 259 V._elements.put(Index.valueOf(i), e); 260 } 261 } 262 return V; 263 } 264 265 @Override 266 public SparseMatrix<F> opposite() { 267 SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero, _transposed); 268 for (int i = 0, p = _rows.size(); i < p; i++) { 269 M._rows.add(_rows.get(i).opposite()); 270 } 271 return M; 272 } 273 274 @Override 275 public SparseMatrix<F> plus(Matrix<F> that) { 276 if (this.getNumberOfRows() != that.getNumberOfRows()) 277 throw new DimensionException(); 278 SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero, _transposed); 279 for (int i = 0, p = _rows.size(); i < p; i++) { 280 M._rows.add(_rows.get(i).plus( 281 _transposed ? that.getColumn(i) : that.getRow(i))); 282 } 283 return M; 284 } 285 286 @Override 287 public SparseMatrix<F> minus(Matrix<F> that) { // Returns more specialized type. 288 return this.plus(that.opposite()); 289 } 290 291 @Override 292 public SparseMatrix<F> times(F k) { 293 SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero, _transposed); 294 for (int i = 0, p = _rows.size(); i < p; i++) { 295 M._rows.add(_rows.get(i).times(k)); 296 } 297 return M; 298 } 299 300 @Override 301 public SparseVector<F> times(Vector<F> v) { 302 if (v.getDimension() != this.getNumberOfColumns()) 303 throw new DimensionException(); 304 final int m = this.getNumberOfRows(); 305 SparseVector<F> V = SparseVector.newInstance(m, _zero); 306 for (int i = 0; i < m; i++) { 307 F e = this.getRow(i).times(v); 308 if (!_zero.equals(e)) { 309 V._elements.put(Index.valueOf(i), e); 310 } 311 } 312 return V; 313 } 314 315 @Override 316 public SparseMatrix<F> times(Matrix<F> that) { 317 final int m = this.getNumberOfRows(); 318 final int n = this.getNumberOfColumns(); 319 final int p = that.getNumberOfColumns(); 320 if (that.getNumberOfRows() != n) 321 throw new DimensionException(); 322 // Creates a mxp matrix in transposed form (p columns vectors of size m) 323 FastTable<SparseVector<F>> rows = this.getRows(); 324 SparseMatrix<F> M = SparseMatrix.newInstance(m, _zero, true); 325 for (int j = 0; j < p; j++) { 326 Vector<F> thatColj = that.getColumn(j); 327 SparseVector<F> column = SparseVector.newInstance(m, _zero); 328 M._rows.add(column); // M is transposed. 329 for (int i = 0; i < m; i++) { 330 F e = rows.get(i).times(thatColj); 331 if (!_zero.equals(e)) { 332 column._elements.put(Index.valueOf(i), e); 333 } 334 } 335 } 336 return M; 337 } 338 private FastTable<SparseVector<F>> getRows() { 339 if (!_transposed) 340 return _rows; 341 FastTable<SparseVector<F>> rows = FastTable.newInstance(); 342 for (int i = 0; i < _n; i++) { 343 rows.add(this.getRow(i)); 344 } 345 return rows; 346 } 347 348 @Override 349 public SparseMatrix<F> inverse() { 350 if (!isSquare()) 351 throw new DimensionException("Matrix not square"); 352 F detInv = this.determinant().inverse(); 353 SparseMatrix<F> A = this.adjoint(); 354 // Multiply adjoint elements with 1 / determinant. 355 for (int i = 0, m = A._rows.size(); i < m; i++) { 356 SparseVector<F> row = A._rows.get(i); 357 for (FastMap.Entry<Index, F> e = row._elements.head(), end = row._elements 358 .tail(); (e = e.getNext()) != end;) { 359 F element = e.getValue(); 360 e.setValue(detInv.times(element)); 361 } 362 } 363 return A; 364 } 365 366 @Override 367 public F determinant() { 368 if (!isSquare()) 369 throw new DimensionException("Matrix not square"); 370 if (_n == 1) 371 return this.get(0, 0); 372 // Expansion by minors (also known as Laplacian) 373 // This algorithm is division free but too slow for dense matrix. 374 SparseVector<F> row0 = this.getRow(0); 375 F det = null; 376 for (FastMap.Entry<Index, F> e = row0._elements.head(), end = row0._elements 377 .tail(); (e = e.getNext()) != end;) { 378 int i = e.getKey().intValue(); 379 F d = e.getValue().times(cofactor(0, i)); 380 if (i % 2 != 0) { 381 d = d.opposite(); 382 } 383 det = (det == null) ? d : det.plus(d); 384 } 385 return det == null ? _zero : det; 386 } 387 388 @Override 389 public Matrix<F> solve(Matrix<F> y) { 390 return this.inverse().times(y); 391 } 392 393 @Override 394 public SparseMatrix<F> transpose() { 395 SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero, !_transposed); 396 M._rows.addAll(this._rows); 397 return M; 398 } 399 400 @Override 401 public F cofactor(int i, int j) { 402 if (_transposed) { 403 int k = i; 404 i = j; 405 j = k; // Swaps i,j 406 } 407 int m = _rows.size(); 408 SparseMatrix<F> M = SparseMatrix.newInstance(m - 1, _zero, _transposed); 409 for (int k1 = 0; k1 < m; k1++) { 410 if (k1 == i) 411 continue; 412 SparseVector<F> row = _rows.get(k1); 413 SparseVector<F> V = SparseVector.newInstance(_n - 1, _zero); 414 M._rows.add(V); 415 for (FastMap.Entry<Index, F> e = row._elements.head(), end = row._elements 416 .tail(); (e = e.getNext()) != end;) { 417 int index = e.getKey().intValue(); 418 if (index < j) { 419 V._elements.put(e.getKey(), e.getValue()); 420 } else if (index > j) { // Position shifted (index minus one). 421 V._elements.put(Index.valueOf(index - 1), e.getValue()); 422 } // Else don't copy element at j. 423 } 424 } 425 return M.determinant(); 426 } 427 428 @Override 429 public SparseMatrix<F> adjoint() { 430 SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero, _transposed); 431 int m = _rows.size(); 432 for (int i = 0; i < m; i++) { 433 SparseVector<F> row = SparseVector.newInstance(_n, _zero); 434 M._rows.add(row); 435 for (int j = 0; j < _n; j++) { 436 F cofactor = _transposed ? cofactor(j, i) : cofactor(i, j); 437 if (!_zero.equals(cofactor)) { 438 row._elements 439 .put(Index.valueOf(j), 440 ((i + j) % 2 == 0) ? cofactor : cofactor 441 .opposite()); 442 } 443 } 444 } 445 return M.transpose(); 446 } 447 448 @Override 449 public SparseMatrix<F> tensor(Matrix<F> that) { 450 final int thism = this.getNumberOfRows(); 451 final int thisn = this.getNumberOfColumns(); 452 final int thatm = that.getNumberOfRows(); 453 final int thatn = that.getNumberOfColumns(); 454 int n = thisn * thatn; // Number of columns, 455 int m = thism * thatm; // Number of rows. 456 SparseMatrix<F> M = SparseMatrix.newInstance(n, _zero, false); 457 for (int i=0; i < m; i++) { // Row index. 458 final int i_rem_thatm = i % thatm; 459 final int i_div_thatm = i / thatm; 460 SparseVector<F> row = SparseVector.newInstance(n, _zero); 461 M._rows.add(row); 462 SparseVector<F> thisRow = this.getRow(i_div_thatm); 463 for (FastMap.Entry<Index, F> e = thisRow._elements.head(), 464 end = thisRow._elements.tail(); (e = e.getNext()) != end;) { 465 F a = e.getValue(); 466 int j = e.getKey().intValue(); 467 for (int k=0; k < thatn; k++) { 468 F b = that.get(i_rem_thatm, k); 469 if (!b.equals(_zero)) { 470 row._elements.put(Index.valueOf(j * thatn + k), a.times(b)); 471 } 472 } 473 } 474 } 475 return M; 476 } 477 478 @Override 479 public SparseVector<F> vectorization() { 480 SparseVector<F> V = SparseVector.newInstance(_n 481 * this.getNumberOfRows(), _zero); 482 int offset = 0; 483 for (int j = 0, n = this.getNumberOfColumns(); j < n; j++) { 484 SparseVector<F> column = this.getColumn(j); 485 for (FastMap.Entry<Index, F> e = column._elements.head(), end = column._elements 486 .tail(); (e = e.getNext()) != end;) { 487 V._elements.put(Index.valueOf(e.getKey().intValue() + offset), 488 e.getValue()); 489 } 490 offset += this.getNumberOfRows(); 491 } 492 return V; 493 } 494 495 @SuppressWarnings("unchecked") 496 @Override 497 public SparseMatrix<F> copy() { 498 SparseMatrix<F> M = newInstance(_n, (F)_zero.copy(), _transposed); 499 for (SparseVector<F> row : _rows) { 500 M._rows.add(row.copy()); 501 } 502 return M; 503 } 504 505 /////////////////////// 506 // Factory creation. // 507 /////////////////////// 508 509 @SuppressWarnings("unchecked") 510 static <F extends Field<F>> SparseMatrix<F> newInstance(int n, F zero, 511 boolean transposed) { 512 SparseMatrix<F> M = FACTORY.object(); 513 M._n = n; 514 M._zero = zero; 515 M._transposed = transposed; 516 return M; 517 } 518 519 private static final ObjectFactory<SparseMatrix> FACTORY = new ObjectFactory<SparseMatrix>() { 520 @Override 521 protected SparseMatrix create() { 522 return new SparseMatrix(); 523 } 524 525 @Override 526 protected void cleanup(SparseMatrix matrix) { 527 matrix._rows.reset(); 528 } 529 }; 530 531 private SparseMatrix() { 532 } 533 534 private static final long serialVersionUID = 1L; 535 536}