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