Source code

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.ConcurrentContext;
015import javolution.context.ObjectFactory;
016import javolution.lang.MathLib;
017import javolution.util.FastTable;
018import org.jscience.mathematics.structure.Field;
019
020/**
021 * <p> This class represents a matrix made of {@link DenseVector dense
022 *     vectors} (as rows). To create a dense matrix made of column vectors the 
023 *     {@link #transpose} method can be used. 
024 *     For example:[code]
025 *        DenseVector<Rational> column0 = DenseVector.valueOf(...);
026 *        DenseVector<Rational> column1 = DenseVector.valueOf(...);
027 *        DenseMatrix<Rational> M = DenseMatrix.valueOf(column0, column1).transpose();
028 *     [/code]</p>
029 * <p> As for any concrete {@link org.jscience.mathematics.structure.Structure
030 *     structure}, this class is declared <code>final</code> (otherwise most
031 *     operations would have to be overridden to return the appropriate type).
032 *     Specialized dense matrix should sub-class {@link Matrix} directly.
033 *     For example:[code]
034 *        // Extension through composition.
035 *        final class TriangularMatrix <F extends Field<F>> extends Matrix<F> {
036 *             private DenseMatrix<F> _value; // Possible implementation.
037 *             ...
038 *             public TriangularMatrix opposite() { // Returns the right type.
039 *                 return TriangularMatrix.valueOf(_value.opposite());
040 *             }
041 *             ...
042 *        }[/code]
043 *     </p>   
044 * @author <a href="mailto:jean-marie@dautelle.com">Jean-Marie Dautelle</a>
045 * @version 3.3, January 2, 2007
046 */
047public final class DenseMatrix<F extends Field<F>> extends Matrix<F> {
048
049    /**
050     * Holds the number of columns n.
051     */
052    int _n;
053
054    /**
055     * Indicates if this matrix is transposed (the rows are then the columns).
056     */
057    boolean _transposed;
058
059    /**
060     * Holds this matrix rows (or columns when transposed).
061     */
062    final FastTable<DenseVector<F>> _rows = new FastTable<DenseVector<F>>();
063
064    /**
065     * Returns a dense matrix from the specified 2-dimensional array.
066     * The first dimension being the row and the second being the column.
067     * 
068     * @param elements this matrix elements.
069     * @return the matrix having the specified elements.
070     * @throws DimensionException if rows have different length.
071     * @see    DenseMatrix 
072     */
073    public static <F extends Field<F>> DenseMatrix<F> valueOf(F[][] elements) {
074        int m = elements.length;
075        int n = elements[0].length;
076        DenseMatrix<F> M = DenseMatrix.newInstance(n, false);
077        for (int i = 0; i < m; i++) {
078            DenseVector<F> row = DenseVector.valueOf(elements[i]);
079            if (row.getDimension() != n)
080                throw new DimensionException();
081            M._rows.add(row);
082        }
083        return M;
084    }
085
086    /**
087     * Returns a dense matrix holding the specified row vectors 
088     * (column vectors if {@link #transpose transposed}).
089     *
090     * @param rows the row vectors.
091     * @return the matrix having the specified rows.
092     * @throws DimensionException if the rows do not have the same dimension.
093     */
094    public static <F extends Field<F>> DenseMatrix<F> valueOf(
095            DenseVector<F>... rows) {
096        final int n = rows[0].getDimension();
097        DenseMatrix<F> M = DenseMatrix.newInstance(n, false);
098        for (int i = 0, m = rows.length; i < m; i++) {
099            DenseVector<F> rowi = rows[i];
100            if (rowi.getDimension() != n)
101                throw new DimensionException(
102                        "All vectors must have the same dimension.");
103            M._rows.add(rowi);
104        }
105        return M;
106    }
107    
108    /**
109     * Returns a dense matrix holding the row vectors from the specified 
110     * collection (column vectors if {@link #transpose transposed}).
111     *
112     * @param rows the list of row vectors.
113     * @return the matrix having the specified rows.
114     * @throws DimensionException if the rows do not have the same dimension.
115     */
116    public static <F extends Field<F>> DenseMatrix<F> valueOf(
117            List<DenseVector<F>> rows) {
118        final int n = rows.get(0).getDimension();
119        DenseMatrix<F> M = DenseMatrix.newInstance(n, false);
120        Iterator<DenseVector<F>> iterator = rows.iterator();
121        for (int i = 0, m = rows.size(); i < m; i++) {
122            DenseVector<F> rowi = iterator.next();
123            if (rowi.getDimension() != n)
124                throw new DimensionException(
125                        "All vectors must have the same dimension.");
126            M._rows.add(rowi);
127        }
128        return M;
129    }
130
131    /**
132     * Returns a dense matrix equivalent to the specified matrix.
133     *
134     * @param that the matrix to convert.
135     * @return <code>that</code> or a dense matrix holding the same elements
136     *         as the specified matrix.
137     */
138    public static <F extends Field<F>> DenseMatrix<F> valueOf(Matrix<F> that) {
139        if (that instanceof DenseMatrix)
140            return (DenseMatrix<F>) that;
141        int n = that.getNumberOfColumns();
142        int m = that.getNumberOfRows();
143        DenseMatrix<F> M = DenseMatrix.newInstance(n, false);
144        for (int i = 0; i < m; i++) {
145            DenseVector<F> rowi = DenseVector.valueOf(that.getRow(i));
146            M._rows.add(rowi);
147        }
148        return M;
149    }
150
151    @Override
152    public int getNumberOfRows() {
153        return _transposed ? _n : _rows.size();
154    }
155
156    @Override
157    public int getNumberOfColumns() {
158        return _transposed ? _rows.size() : _n;
159    }
160
161    @Override
162    public F get(int i, int j) {
163        return _transposed ? _rows.get(j).get(i) : _rows.get(i).get(j);
164    }
165
166    @Override
167    public DenseVector<F> getRow(int i) {
168        if (!_transposed)
169            return _rows.get(i);
170        // Else transposed.
171        int n = _rows.size();
172        int m = _n;
173        if ((i < 0) || (i >= m))
174            throw new DimensionException();
175        DenseVector<F> V = DenseVector.newInstance();
176        for (int j = 0; j < n; j++) {
177            V._elements.add(_rows.get(j).get(i));
178        }
179        return V;
180    }
181
182    @Override
183    public DenseVector<F> getColumn(int j) {
184        if (_transposed)
185            return _rows.get(j);
186        int m = _rows.size();
187        if ((j < 0) || (j >= _n))
188            throw new DimensionException();
189        DenseVector<F> V = DenseVector.newInstance();
190        for (int i = 0; i < m; i++) {
191            V._elements.add(_rows.get(i).get(j));
192        }
193        return V;
194    }
195
196    @Override
197    public DenseVector<F> getDiagonal() {
198        int m = this.getNumberOfRows();
199        int n = this.getNumberOfColumns();
200        int dimension = MathLib.min(m, n);
201        DenseVector<F> V = DenseVector.newInstance();
202        for (int i = 0; i < dimension; i++) {
203            V._elements.add(this.get(i, i));
204        }
205        return V;
206    }
207
208    @Override
209    public DenseMatrix<F> opposite() {
210        DenseMatrix<F> M = DenseMatrix.newInstance(_n, _transposed);
211        for (int i = 0, p = _rows.size(); i < p; i++) {
212            M._rows.add(_rows.get(i).opposite());
213        }
214        return M;
215    }
216
217    @Override
218    public DenseMatrix<F> plus(Matrix<F> that) {
219        if (this.getNumberOfRows() != that.getNumberOfRows())
220            throw new DimensionException();
221        DenseMatrix<F> M = DenseMatrix.newInstance(_n, _transposed);
222        for (int i = 0, p = _rows.size(); i < p; i++) {
223            M._rows.add(_rows.get(i).plus(
224                    _transposed ? that.getColumn(i) : that.getRow(i)));
225        }
226        return M;
227    }
228
229    @Override
230    public DenseMatrix<F> minus(Matrix<F> that) { // Returns more specialized type.
231        return this.plus(that.opposite());
232    }
233
234    @Override
235    public DenseMatrix<F> times(F k) {
236        DenseMatrix<F> M = DenseMatrix.newInstance(_n, _transposed);
237        for (int i = 0, p = _rows.size(); i < p; i++) {
238            M._rows.add(_rows.get(i).times(k));
239        }
240        return M;
241    }
242
243    @Override
244    public DenseVector<F> times(Vector<F> v) {
245        if (v.getDimension() != this.getNumberOfColumns())
246            throw new DimensionException();
247        final int m = this.getNumberOfRows();
248        DenseVector<F> V = DenseVector.newInstance();
249        for (int i = 0; i < m; i++) {
250            V._elements.add(this.getRow(i).times(v));
251        }
252        return V;
253    }
254
255    @Override
256    public DenseMatrix<F> times(Matrix<F> that) {
257        final int n = this.getNumberOfColumns();
258        final int m = this.getNumberOfRows();
259        final int p = that.getNumberOfColumns();
260        if (that.getNumberOfRows() != n)
261            throw new DimensionException();
262        // Creates a mxp matrix in transposed form (p columns vectors of size m)
263        DenseMatrix<F> M = DenseMatrix.newInstance(m, true); // Transposed.
264        M._rows.setSize(p);
265        Multiply<F> multiply = Multiply.valueOf(this, that, 0, p, M._rows);
266        multiply.run();
267        Multiply.recycle(multiply);
268        return M;
269    }
270
271    // Logic to multiply two matrices. 
272    private static class Multiply<F extends Field<F>> implements Runnable {
273        private static final ObjectFactory<Multiply> FACTORY = new ObjectFactory<Multiply>() {
274
275            @Override
276            protected Multiply create() {
277                return new Multiply();
278            }
279        };
280
281        private DenseMatrix<F> _left;
282
283        private Matrix<F> _right;
284
285        private int _rightColumnStart;
286
287        private int _rightColumnEnd;
288
289        private FastTable<DenseVector<F>> _columnsResult;
290
291        @SuppressWarnings("unchecked")
292        static <F extends Field<F>> Multiply<F> valueOf(DenseMatrix<F> left, Matrix<F> right,
293                int rightColumnStart, int rightColumnEnd,
294                FastTable<DenseVector<F>> columnsResult) {
295            Multiply<F> multiply = Multiply.FACTORY.object();
296            multiply._left = left;
297            multiply._right = right;
298            multiply._rightColumnStart = rightColumnStart;
299            multiply._rightColumnEnd = rightColumnEnd;
300            multiply._columnsResult = columnsResult;
301            return multiply;
302        }
303
304        static <F extends Field<F>> void recycle(Multiply<F> multiply) {
305            multiply._left = null;
306            multiply._right = null;
307            multiply._columnsResult = null;
308            Multiply.FACTORY.recycle(multiply);
309        }
310
311        public void run() {
312            if (_rightColumnEnd - _rightColumnStart < 32) { // Direct calculation.
313                FastTable<DenseVector<F>> rows = _left.getRows();
314                final int m = rows.size();
315                for (int j = _rightColumnStart; j < _rightColumnEnd; j++) {
316                    Vector<F> thatColj = _right.getColumn(j);
317                    DenseVector<F> column = DenseVector.newInstance();
318                    _columnsResult.set(j, column);
319                    for (int i = 0; i < m; i++) {
320                        column._elements.add(rows.get(i).times(thatColj));
321                    }
322                }
323            } else { // Concurrent/Recursive calculation.
324                int halfIndex = (_rightColumnStart + _rightColumnEnd) >> 1;
325                Multiply<F> firstHalf = Multiply.valueOf(_left, _right,
326                        _rightColumnStart, halfIndex, _columnsResult);
327                Multiply<F> secondHalf = Multiply.valueOf(_left, _right,
328                        halfIndex, _rightColumnEnd, _columnsResult);
329                ConcurrentContext.enter();
330                try {
331                    ConcurrentContext.execute(firstHalf);
332                    ConcurrentContext.execute(secondHalf);
333                } finally {
334                    ConcurrentContext.exit();
335                }
336                Multiply.recycle(firstHalf);
337                Multiply.recycle(secondHalf);
338            }
339        }
340    }
341
342    private FastTable<DenseVector<F>> getRows() {
343        if (!_transposed)
344            return _rows;
345        FastTable<DenseVector<F>> rows = FastTable.newInstance();
346        for (int i = 0; i < _n; i++) {
347            rows.add(this.getRow(i));
348        }
349        return rows;
350    }
351
352    @Override
353    public DenseMatrix<F> inverse() {
354        if (!isSquare())
355            throw new DimensionException("Matrix not square");
356        return LUDecomposition.valueOf(this).inverse();
357    }
358
359    @Override
360    public F determinant() {
361        return LUDecomposition.valueOf(this).determinant();
362    }
363
364    @Override
365    public DenseMatrix<F> transpose() {
366        DenseMatrix<F> M = DenseMatrix.newInstance(_n, !_transposed);
367        M._rows.addAll(this._rows);
368        return M;
369    }
370
371    @Override
372    public F cofactor(int i, int j) {
373        if (_transposed) {
374            int k = i;
375            i = j;
376            j = k; // Swaps i,j
377        }
378        int m = _rows.size();
379        DenseMatrix<F> M = DenseMatrix.newInstance(m - 1, _transposed);
380        for (int k1 = 0; k1 < m; k1++) {
381            if (k1 == i)
382                continue;
383            DenseVector<F> row = _rows.get(k1);
384            DenseVector<F> V = DenseVector.newInstance();
385            M._rows.add(V);
386            for (int k2 = 0; k2 < _n; k2++) {
387                if (k2 == j)
388                    continue;
389                V._elements.add(row.get(k2));
390            }
391        }
392        return M.determinant();
393    }
394
395    @Override
396    public DenseMatrix<F> adjoint() {
397        DenseMatrix<F> M = DenseMatrix.newInstance(_n, _transposed);
398        int m = _rows.size();
399        for (int i = 0; i < m; i++) {
400            DenseVector<F> row = DenseVector.newInstance();
401            M._rows.add(row);
402            for (int j = 0; j < _n; j++) {
403                F cofactor = _transposed ? cofactor(j, i) : cofactor(i, j);
404                row._elements.add(((i + j) % 2 == 0) ? cofactor : cofactor
405                        .opposite());
406            }
407        }
408        return M.transpose();
409    }
410
411    @Override
412    public Matrix<F> tensor(Matrix<F> that) {
413        final int thism = this.getNumberOfRows();
414        final int thisn = this.getNumberOfColumns();
415        final int thatm = that.getNumberOfRows();
416        final int thatn = that.getNumberOfColumns();
417        int n = thisn * thatn; // Number of columns,
418        int m = thism * thatm; // Number of rows.
419        DenseMatrix<F> M = DenseMatrix.newInstance(n, false);
420        for (int i = 0; i < m; i++) { // Row index.
421            final int i_rem_thatm = i % thatm;
422            final int i_div_thatm = i / thatm;
423            DenseVector<F> row = DenseVector.newInstance();
424            M._rows.add(row);
425            for (int j = 0; j < thisn; j++) {
426                F a = this.get(i_div_thatm, j);
427                for (int k = 0; k < thatn; k++) {
428                    row._elements.add(a.times(that.get(i_rem_thatm, k)));
429                }
430            }
431        }
432        return M;
433    }
434
435    @Override
436    public Vector<F> vectorization() {
437        DenseVector<F> V = DenseVector.newInstance();
438        for (int j = 0, n = this.getNumberOfColumns(); j < n; j++) {
439            Vector<F> column = this.getColumn(j);
440            for (int i = 0, m = column.getDimension(); i < m; i++) {
441                V._elements.add(column.get(i));
442            }
443        }
444        return V;
445    }
446
447    @Override
448    public DenseMatrix<F> copy() {
449        DenseMatrix<F> M = DenseMatrix.newInstance(_n, _transposed);
450        for (DenseVector<F> row : _rows) {
451            M._rows.add(row.copy());
452        }
453        return M;
454    }
455
456    ///////////////////////////////
457    // Package Private Utilities //
458    ///////////////////////////////
459
460    void set(int i, int j, F e) {
461        if (_transposed) {
462            _rows.get(j)._elements.set(i, e);
463        } else {
464            _rows.get(i)._elements.set(j, e);
465        }
466    }
467
468    ///////////////////////
469    // Factory creation. //
470    ///////////////////////
471
472    @SuppressWarnings("unchecked")
473    static <F extends Field<F>> DenseMatrix<F> newInstance(int n,
474            boolean transposed) {
475        DenseMatrix<F> M = FACTORY.object();
476        M._n = n;
477        M._transposed = transposed;
478        return M;
479    }
480
481    private static ObjectFactory<DenseMatrix> FACTORY = new ObjectFactory<DenseMatrix>() {
482        @Override
483        protected DenseMatrix create() {
484            return new DenseMatrix();
485        }
486
487        @Override
488        protected void cleanup(DenseMatrix matrix) {
489            matrix._rows.reset();
490        }
491    };
492
493    private DenseMatrix() {
494    }
495
496    private static final long serialVersionUID = 1L;
497
498}