001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018package org.apache.commons.math4.linear;
019
020import java.io.IOException;
021import java.io.ObjectInputStream;
022import java.io.ObjectOutputStream;
023import java.util.Arrays;
024
025import org.apache.commons.math4.Field;
026import org.apache.commons.math4.FieldElement;
027import org.apache.commons.math4.exception.DimensionMismatchException;
028import org.apache.commons.math4.exception.MathArithmeticException;
029import org.apache.commons.math4.exception.NoDataException;
030import org.apache.commons.math4.exception.NullArgumentException;
031import org.apache.commons.math4.exception.NumberIsTooSmallException;
032import org.apache.commons.math4.exception.OutOfRangeException;
033import org.apache.commons.math4.exception.ZeroException;
034import org.apache.commons.math4.exception.util.LocalizedFormats;
035import org.apache.commons.math4.util.FastMath;
036import org.apache.commons.math4.util.MathArrays;
037import org.apache.commons.math4.util.MathUtils;
038import org.apache.commons.numbers.core.Precision;
039
040/**
041 * A collection of static methods that operate on or return matrices.
042 *
043 */
044public class MatrixUtils {
045
046    /**
047     * The default format for {@link RealMatrix} objects.
048     * @since 3.1
049     */
050    public static final RealMatrixFormat DEFAULT_FORMAT = RealMatrixFormat.getInstance();
051
052    /**
053     * A format for {@link RealMatrix} objects compatible with octave.
054     * @since 3.1
055     */
056    public static final RealMatrixFormat OCTAVE_FORMAT = new RealMatrixFormat("[", "]", "", "", "; ", ", ");
057
058    /**
059     * Private constructor.
060     */
061    private MatrixUtils() {
062        super();
063    }
064
065    /**
066     * Returns a {@link RealMatrix} with specified dimensions.
067     * <p>The type of matrix returned depends on the dimension. Below
068     * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
069     * square matrix) which can be stored in a 32kB array, a {@link
070     * Array2DRowRealMatrix} instance is built. Above this threshold a {@link
071     * BlockRealMatrix} instance is built.</p>
072     * <p>The matrix elements are all set to 0.0.</p>
073     * @param rows number of rows of the matrix
074     * @param columns number of columns of the matrix
075     * @return  RealMatrix with specified dimensions
076     * @see #createRealMatrix(double[][])
077     */
078    public static RealMatrix createRealMatrix(final int rows, final int columns) {
079        return (rows * columns <= 4096) ?
080                new Array2DRowRealMatrix(rows, columns) : new BlockRealMatrix(rows, columns);
081    }
082
083    /**
084     * Returns a {@link FieldMatrix} with specified dimensions.
085     * <p>The type of matrix returned depends on the dimension. Below
086     * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
087     * square matrix), a {@link FieldMatrix} instance is built. Above
088     * this threshold a {@link BlockFieldMatrix} instance is built.</p>
089     * <p>The matrix elements are all set to field.getZero().</p>
090     * @param <T> the type of the field elements
091     * @param field field to which the matrix elements belong
092     * @param rows number of rows of the matrix
093     * @param columns number of columns of the matrix
094     * @return  FieldMatrix with specified dimensions
095     * @see #createFieldMatrix(FieldElement[][])
096     * @since 2.0
097     */
098    public static <T extends FieldElement<T>> FieldMatrix<T> createFieldMatrix(final Field<T> field,
099                                                                               final int rows,
100                                                                               final int columns) {
101        return (rows * columns <= 4096) ?
102                new Array2DRowFieldMatrix<>(field, rows, columns) : new BlockFieldMatrix<>(field, rows, columns);
103    }
104
105    /**
106     * Returns a {@link RealMatrix} whose entries are the values in the
107     * the input array.
108     * <p>The type of matrix returned depends on the dimension. Below
109     * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
110     * square matrix) which can be stored in a 32kB array, a {@link
111     * Array2DRowRealMatrix} instance is built. Above this threshold a {@link
112     * BlockRealMatrix} instance is built.</p>
113     * <p>The input array is copied, not referenced.</p>
114     *
115     * @param data input array
116     * @return  RealMatrix containing the values of the array
117     * @throws org.apache.commons.math4.exception.DimensionMismatchException
118     * if {@code data} is not rectangular (not all rows have the same length).
119     * @throws NoDataException if a row or column is empty.
120     * @throws NullArgumentException if either {@code data} or {@code data[0]}
121     * is {@code null}.
122     * @throws DimensionMismatchException if {@code data} is not rectangular.
123     * @see #createRealMatrix(int, int)
124     */
125    public static RealMatrix createRealMatrix(double[][] data)
126        throws NullArgumentException, DimensionMismatchException,
127        NoDataException {
128        if (data == null ||
129            data[0] == null) {
130            throw new NullArgumentException();
131        }
132        return (data.length * data[0].length <= 4096) ?
133                new Array2DRowRealMatrix(data) : new BlockRealMatrix(data);
134    }
135
136    /**
137     * Returns a {@link FieldMatrix} whose entries are the values in the
138     * the input array.
139     * <p>The type of matrix returned depends on the dimension. Below
140     * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
141     * square matrix), a {@link FieldMatrix} instance is built. Above
142     * this threshold a {@link BlockFieldMatrix} instance is built.</p>
143     * <p>The input array is copied, not referenced.</p>
144     * @param <T> the type of the field elements
145     * @param data input array
146     * @return a matrix containing the values of the array.
147     * @throws org.apache.commons.math4.exception.DimensionMismatchException
148     * if {@code data} is not rectangular (not all rows have the same length).
149     * @throws NoDataException if a row or column is empty.
150     * @throws NullArgumentException if either {@code data} or {@code data[0]}
151     * is {@code null}.
152     * @see #createFieldMatrix(Field, int, int)
153     * @since 2.0
154     */
155    public static <T extends FieldElement<T>> FieldMatrix<T> createFieldMatrix(T[][] data)
156        throws DimensionMismatchException, NoDataException, NullArgumentException {
157        if (data == null ||
158            data[0] == null) {
159            throw new NullArgumentException();
160        }
161        return (data.length * data[0].length <= 4096) ?
162                new Array2DRowFieldMatrix<>(data) : new BlockFieldMatrix<>(data);
163    }
164
165    /**
166     * Returns <code>dimension x dimension</code> identity matrix.
167     *
168     * @param dimension dimension of identity matrix to generate
169     * @return identity matrix
170     * @throws IllegalArgumentException if dimension is not positive
171     * @since 1.1
172     */
173    public static RealMatrix createRealIdentityMatrix(int dimension) {
174        final RealMatrix m = createRealMatrix(dimension, dimension);
175        for (int i = 0; i < dimension; ++i) {
176            m.setEntry(i, i, 1.0);
177        }
178        return m;
179    }
180
181    /**
182     * Returns <code>dimension x dimension</code> identity matrix.
183     *
184     * @param <T> the type of the field elements
185     * @param field field to which the elements belong
186     * @param dimension dimension of identity matrix to generate
187     * @return identity matrix
188     * @throws IllegalArgumentException if dimension is not positive
189     * @since 2.0
190     */
191    public static <T extends FieldElement<T>> FieldMatrix<T>
192        createFieldIdentityMatrix(final Field<T> field, final int dimension) {
193        final T zero = field.getZero();
194        final T one  = field.getOne();
195        final T[][] d = MathArrays.buildArray(field, dimension, dimension);
196        for (int row = 0; row < dimension; row++) {
197            final T[] dRow = d[row];
198            Arrays.fill(dRow, zero);
199            dRow[row] = one;
200        }
201        return new Array2DRowFieldMatrix<>(field, d, false);
202    }
203
204    /**
205     * Creates a diagonal matrix with the specified diagonal elements.
206     *
207     * @param diagonal Diagonal elements of the matrix.
208     * The array elements will be copied.
209     * @return a diagonal matrix instance.
210     *
211     * @see #createRealMatrixWithDiagonal(double[])
212     * @since 2.0
213     */
214    public static DiagonalMatrix createRealDiagonalMatrix(final double[] diagonal) {
215        return new DiagonalMatrix(diagonal, true);
216    }
217
218    /**
219     * Creates a dense matrix with the specified diagonal elements.
220     *
221     * @param diagonal Diagonal elements of the matrix.
222     * @return a matrix instance.
223     *
224     * @see #createRealDiagonalMatrix(double[])
225     * @since 4.0
226     */
227    public static RealMatrix createRealMatrixWithDiagonal(final double[] diagonal) {
228        final int size = diagonal.length;
229        final RealMatrix m = createRealMatrix(size, size);
230        for (int i = 0; i < size; i++) {
231            m.setEntry(i, i, diagonal[i]);
232        }
233        return m;
234    }
235
236    /**
237     * Returns a diagonal matrix with specified elements.
238     *
239     * @param <T> the type of the field elements
240     * @param diagonal diagonal elements of the matrix (the array elements
241     * will be copied)
242     * @return diagonal matrix
243     * @since 2.0
244     */
245    public static <T extends FieldElement<T>> FieldMatrix<T>
246        createFieldDiagonalMatrix(final T[] diagonal) {
247        final FieldMatrix<T> m =
248            createFieldMatrix(diagonal[0].getField(), diagonal.length, diagonal.length);
249        for (int i = 0; i < diagonal.length; ++i) {
250            m.setEntry(i, i, diagonal[i]);
251        }
252        return m;
253    }
254
255    /**
256     * Creates a {@link RealVector} using the data from the input array.
257     *
258     * @param data the input data
259     * @return a data.length RealVector
260     * @throws NoDataException if {@code data} is empty.
261     * @throws NullArgumentException if {@code data} is {@code null}.
262     */
263    public static RealVector createRealVector(double[] data)
264        throws NoDataException, NullArgumentException {
265        if (data == null) {
266            throw new NullArgumentException();
267        }
268        return new ArrayRealVector(data, true);
269    }
270
271    /**
272     * Creates a {@link FieldVector} using the data from the input array.
273     *
274     * @param <T> the type of the field elements
275     * @param data the input data
276     * @return a data.length FieldVector
277     * @throws NoDataException if {@code data} is empty.
278     * @throws NullArgumentException if {@code data} is {@code null}.
279     * @throws ZeroException if {@code data} has 0 elements
280     */
281    public static <T extends FieldElement<T>> FieldVector<T> createFieldVector(final T[] data)
282        throws NoDataException, NullArgumentException, ZeroException {
283        if (data == null) {
284            throw new NullArgumentException();
285        }
286        if (data.length == 0) {
287            throw new ZeroException(LocalizedFormats.VECTOR_MUST_HAVE_AT_LEAST_ONE_ELEMENT);
288        }
289        return new ArrayFieldVector<>(data[0].getField(), data, true);
290    }
291
292    /**
293     * Create a row {@link RealMatrix} using the data from the input
294     * array.
295     *
296     * @param rowData the input row data
297     * @return a 1 x rowData.length RealMatrix
298     * @throws NoDataException if {@code rowData} is empty.
299     * @throws NullArgumentException if {@code rowData} is {@code null}.
300     */
301    public static RealMatrix createRowRealMatrix(double[] rowData)
302        throws NoDataException, NullArgumentException {
303        if (rowData == null) {
304            throw new NullArgumentException();
305        }
306        final int nCols = rowData.length;
307        final RealMatrix m = createRealMatrix(1, nCols);
308        for (int i = 0; i < nCols; ++i) {
309            m.setEntry(0, i, rowData[i]);
310        }
311        return m;
312    }
313
314    /**
315     * Create a row {@link FieldMatrix} using the data from the input
316     * array.
317     *
318     * @param <T> the type of the field elements
319     * @param rowData the input row data
320     * @return a 1 x rowData.length FieldMatrix
321     * @throws NoDataException if {@code rowData} is empty.
322     * @throws NullArgumentException if {@code rowData} is {@code null}.
323     */
324    public static <T extends FieldElement<T>> FieldMatrix<T>
325        createRowFieldMatrix(final T[] rowData)
326        throws NoDataException, NullArgumentException {
327        if (rowData == null) {
328            throw new NullArgumentException();
329        }
330        final int nCols = rowData.length;
331        if (nCols == 0) {
332            throw new NoDataException(LocalizedFormats.AT_LEAST_ONE_COLUMN);
333        }
334        final FieldMatrix<T> m = createFieldMatrix(rowData[0].getField(), 1, nCols);
335        for (int i = 0; i < nCols; ++i) {
336            m.setEntry(0, i, rowData[i]);
337        }
338        return m;
339    }
340
341    /**
342     * Creates a column {@link RealMatrix} using the data from the input
343     * array.
344     *
345     * @param columnData  the input column data
346     * @return a columnData x 1 RealMatrix
347     * @throws NoDataException if {@code columnData} is empty.
348     * @throws NullArgumentException if {@code columnData} is {@code null}.
349     */
350    public static RealMatrix createColumnRealMatrix(double[] columnData)
351        throws NoDataException, NullArgumentException {
352        if (columnData == null) {
353            throw new NullArgumentException();
354        }
355        final int nRows = columnData.length;
356        final RealMatrix m = createRealMatrix(nRows, 1);
357        for (int i = 0; i < nRows; ++i) {
358            m.setEntry(i, 0, columnData[i]);
359        }
360        return m;
361    }
362
363    /**
364     * Creates a column {@link FieldMatrix} using the data from the input
365     * array.
366     *
367     * @param <T> the type of the field elements
368     * @param columnData  the input column data
369     * @return a columnData x 1 FieldMatrix
370     * @throws NoDataException if {@code data} is empty.
371     * @throws NullArgumentException if {@code columnData} is {@code null}.
372     */
373    public static <T extends FieldElement<T>> FieldMatrix<T>
374        createColumnFieldMatrix(final T[] columnData)
375        throws NoDataException, NullArgumentException {
376        if (columnData == null) {
377            throw new NullArgumentException();
378        }
379        final int nRows = columnData.length;
380        if (nRows == 0) {
381            throw new NoDataException(LocalizedFormats.AT_LEAST_ONE_ROW);
382        }
383        final FieldMatrix<T> m = createFieldMatrix(columnData[0].getField(), nRows, 1);
384        for (int i = 0; i < nRows; ++i) {
385            m.setEntry(i, 0, columnData[i]);
386        }
387        return m;
388    }
389
390    /**
391     * Checks whether a matrix is symmetric, within a given relative tolerance.
392     *
393     * @param matrix Matrix to check.
394     * @param relativeTolerance Tolerance of the symmetry check.
395     * @param raiseException If {@code true}, an exception will be raised if
396     * the matrix is not symmetric.
397     * @return {@code true} if {@code matrix} is symmetric.
398     * @throws NonSquareMatrixException if the matrix is not square.
399     * @throws NonSymmetricMatrixException if the matrix is not symmetric.
400     */
401    private static boolean isSymmetricInternal(RealMatrix matrix,
402                                               double relativeTolerance,
403                                               boolean raiseException) {
404        final int rows = matrix.getRowDimension();
405        if (rows != matrix.getColumnDimension()) {
406            if (raiseException) {
407                throw new NonSquareMatrixException(rows, matrix.getColumnDimension());
408            } else {
409                return false;
410            }
411        }
412        for (int i = 0; i < rows; i++) {
413            for (int j = i + 1; j < rows; j++) {
414                final double mij = matrix.getEntry(i, j);
415                final double mji = matrix.getEntry(j, i);
416                if (FastMath.abs(mij - mji) >
417                    FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * relativeTolerance) {
418                    if (raiseException) {
419                        throw new NonSymmetricMatrixException(i, j, relativeTolerance);
420                    } else {
421                        return false;
422                    }
423                }
424            }
425        }
426        return true;
427    }
428
429    /**
430     * Checks whether a matrix is symmetric.
431     *
432     * @param matrix Matrix to check.
433     * @param eps Relative tolerance.
434     * @throws NonSquareMatrixException if the matrix is not square.
435     * @throws NonSymmetricMatrixException if the matrix is not symmetric.
436     * @since 3.1
437     */
438    public static void checkSymmetric(RealMatrix matrix,
439                                      double eps) {
440        isSymmetricInternal(matrix, eps, true);
441    }
442
443    /**
444     * Checks whether a matrix is symmetric.
445     *
446     * @param matrix Matrix to check.
447     * @param eps Relative tolerance.
448     * @return {@code true} if {@code matrix} is symmetric.
449     * @since 3.1
450     */
451    public static boolean isSymmetric(RealMatrix matrix,
452                                      double eps) {
453        return isSymmetricInternal(matrix, eps, false);
454    }
455
456    /**
457     * Check if matrix indices are valid.
458     *
459     * @param m Matrix.
460     * @param row Row index to check.
461     * @param column Column index to check.
462     * @throws OutOfRangeException if {@code row} or {@code column} is not
463     * a valid index.
464     */
465    public static void checkMatrixIndex(final AnyMatrix m,
466                                        final int row, final int column)
467        throws OutOfRangeException {
468        checkRowIndex(m, row);
469        checkColumnIndex(m, column);
470    }
471
472    /**
473     * Check if a row index is valid.
474     *
475     * @param m Matrix.
476     * @param row Row index to check.
477     * @throws OutOfRangeException if {@code row} is not a valid index.
478     */
479    public static void checkRowIndex(final AnyMatrix m, final int row)
480        throws OutOfRangeException {
481        if (row < 0 ||
482            row >= m.getRowDimension()) {
483            throw new OutOfRangeException(LocalizedFormats.ROW_INDEX,
484                                          row, 0, m.getRowDimension() - 1);
485        }
486    }
487
488    /**
489     * Check if a column index is valid.
490     *
491     * @param m Matrix.
492     * @param column Column index to check.
493     * @throws OutOfRangeException if {@code column} is not a valid index.
494     */
495    public static void checkColumnIndex(final AnyMatrix m, final int column)
496        throws OutOfRangeException {
497        if (column < 0 || column >= m.getColumnDimension()) {
498            throw new OutOfRangeException(LocalizedFormats.COLUMN_INDEX,
499                                           column, 0, m.getColumnDimension() - 1);
500        }
501    }
502
503    /**
504     * Check if submatrix ranges indices are valid.
505     * Rows and columns are indicated counting from 0 to {@code n - 1}.
506     *
507     * @param m Matrix.
508     * @param startRow Initial row index.
509     * @param endRow Final row index.
510     * @param startColumn Initial column index.
511     * @param endColumn Final column index.
512     * @throws OutOfRangeException if the indices are invalid.
513     * @throws NumberIsTooSmallException if {@code endRow < startRow} or
514     * {@code endColumn < startColumn}.
515     */
516    public static void checkSubMatrixIndex(final AnyMatrix m,
517                                           final int startRow, final int endRow,
518                                           final int startColumn, final int endColumn)
519        throws NumberIsTooSmallException, OutOfRangeException {
520        checkRowIndex(m, startRow);
521        checkRowIndex(m, endRow);
522        if (endRow < startRow) {
523            throw new NumberIsTooSmallException(LocalizedFormats.INITIAL_ROW_AFTER_FINAL_ROW,
524                                                endRow, startRow, false);
525        }
526
527        checkColumnIndex(m, startColumn);
528        checkColumnIndex(m, endColumn);
529        if (endColumn < startColumn) {
530            throw new NumberIsTooSmallException(LocalizedFormats.INITIAL_COLUMN_AFTER_FINAL_COLUMN,
531                                                endColumn, startColumn, false);
532        }
533    }
534
535    /**
536     * Check if submatrix ranges indices are valid.
537     * Rows and columns are indicated counting from 0 to n-1.
538     *
539     * @param m Matrix.
540     * @param selectedRows Array of row indices.
541     * @param selectedColumns Array of column indices.
542     * @throws NullArgumentException if {@code selectedRows} or
543     * {@code selectedColumns} are {@code null}.
544     * @throws NoDataException if the row or column selections are empty (zero
545     * length).
546     * @throws OutOfRangeException if row or column selections are not valid.
547     */
548    public static void checkSubMatrixIndex(final AnyMatrix m,
549                                           final int[] selectedRows,
550                                           final int[] selectedColumns)
551        throws NoDataException, NullArgumentException, OutOfRangeException {
552        if (selectedRows == null) {
553            throw new NullArgumentException();
554        }
555        if (selectedColumns == null) {
556            throw new NullArgumentException();
557        }
558        if (selectedRows.length == 0) {
559            throw new NoDataException(LocalizedFormats.EMPTY_SELECTED_ROW_INDEX_ARRAY);
560        }
561        if (selectedColumns.length == 0) {
562            throw new NoDataException(LocalizedFormats.EMPTY_SELECTED_COLUMN_INDEX_ARRAY);
563        }
564
565        for (final int row : selectedRows) {
566            checkRowIndex(m, row);
567        }
568        for (final int column : selectedColumns) {
569            checkColumnIndex(m, column);
570        }
571    }
572
573    /**
574     * Check if matrices are addition compatible.
575     *
576     * @param left Left hand side matrix.
577     * @param right Right hand side matrix.
578     * @throws MatrixDimensionMismatchException if the matrices are not addition
579     * compatible.
580     */
581    public static void checkAdditionCompatible(final AnyMatrix left, final AnyMatrix right) {
582        left.checkAdd(right);
583    }
584
585    /**
586     * Check if matrices are subtraction compatible
587     *
588     * @param left Left hand side matrix.
589     * @param right Right hand side matrix.
590     * @throws MatrixDimensionMismatchException if the matrices are not addition
591     * compatible.
592     */
593    public static void checkSubtractionCompatible(final AnyMatrix left, final AnyMatrix right) {
594        left.checkAdd(right);
595    }
596
597    /**
598     * Check if matrices are multiplication compatible
599     *
600     * @param left Left hand side matrix.
601     * @param right Right hand side matrix.
602     * @throws DimensionMismatchException if matrices are not multiplication
603     * compatible.
604     */
605    public static void checkMultiplicationCompatible(final AnyMatrix left, final AnyMatrix right) {
606        left.checkMultiply(right);
607    }
608
609    /** Serialize a {@link RealVector}.
610     * <p>
611     * This method is intended to be called from within a private
612     * <code>writeObject</code> method (after a call to
613     * <code>oos.defaultWriteObject()</code>) in a class that has a
614     * {@link RealVector} field, which should be declared <code>transient</code>.
615     * This way, the default handling does not serialize the vector (the {@link
616     * RealVector} interface is not serializable by default) but this method does
617     * serialize it specifically.
618     * </p>
619     * <p>
620     * The following example shows how a simple class with a name and a real vector
621     * should be written:
622     * <pre><code>
623     * public class NamedVector implements Serializable {
624     *
625     *     private final String name;
626     *     private final transient RealVector coefficients;
627     *
628     *     // omitted constructors, getters ...
629     *
630     *     private void writeObject(ObjectOutputStream oos) throws IOException {
631     *         oos.defaultWriteObject();  // takes care of name field
632     *         MatrixUtils.serializeRealVector(coefficients, oos);
633     *     }
634     *
635     *     private void readObject(ObjectInputStream ois) throws ClassNotFoundException, IOException {
636     *         ois.defaultReadObject();  // takes care of name field
637     *         MatrixUtils.deserializeRealVector(this, "coefficients", ois);
638     *     }
639     *
640     * }
641     * </code></pre>
642     *
643     * @param vector real vector to serialize
644     * @param oos stream where the real vector should be written
645     * @exception IOException if object cannot be written to stream
646     * @see #deserializeRealVector(Object, String, ObjectInputStream)
647     */
648    public static void serializeRealVector(final RealVector vector,
649                                           final ObjectOutputStream oos)
650        throws IOException {
651        final int n = vector.getDimension();
652        oos.writeInt(n);
653        for (int i = 0; i < n; ++i) {
654            oos.writeDouble(vector.getEntry(i));
655        }
656    }
657
658    /** Deserialize  a {@link RealVector} field in a class.
659     * <p>
660     * This method is intended to be called from within a private
661     * <code>readObject</code> method (after a call to
662     * <code>ois.defaultReadObject()</code>) in a class that has a
663     * {@link RealVector} field, which should be declared <code>transient</code>.
664     * This way, the default handling does not deserialize the vector (the {@link
665     * RealVector} interface is not serializable by default) but this method does
666     * deserialize it specifically.
667     * </p>
668     * @param instance instance in which the field must be set up
669     * @param fieldName name of the field within the class (may be private and final)
670     * @param ois stream from which the real vector should be read
671     * @exception ClassNotFoundException if a class in the stream cannot be found
672     * @exception IOException if object cannot be read from the stream
673     * @see #serializeRealVector(RealVector, ObjectOutputStream)
674     */
675    public static void deserializeRealVector(final Object instance,
676                                             final String fieldName,
677                                             final ObjectInputStream ois)
678      throws ClassNotFoundException, IOException {
679        try {
680
681            // read the vector data
682            final int n = ois.readInt();
683            final double[] data = new double[n];
684            for (int i = 0; i < n; ++i) {
685                data[i] = ois.readDouble();
686            }
687
688            // create the instance
689            final RealVector vector = new ArrayRealVector(data, false);
690
691            // set up the field
692            final java.lang.reflect.Field f =
693                instance.getClass().getDeclaredField(fieldName);
694            f.setAccessible(true);
695            f.set(instance, vector);
696
697        } catch (NoSuchFieldException nsfe) {
698            IOException ioe = new IOException();
699            ioe.initCause(nsfe);
700            throw ioe;
701        } catch (IllegalAccessException iae) {
702            IOException ioe = new IOException();
703            ioe.initCause(iae);
704            throw ioe;
705        }
706
707    }
708
709    /** Serialize a {@link RealMatrix}.
710     * <p>
711     * This method is intended to be called from within a private
712     * <code>writeObject</code> method (after a call to
713     * <code>oos.defaultWriteObject()</code>) in a class that has a
714     * {@link RealMatrix} field, which should be declared <code>transient</code>.
715     * This way, the default handling does not serialize the matrix (the {@link
716     * RealMatrix} interface is not serializable by default) but this method does
717     * serialize it specifically.
718     * </p>
719     * <p>
720     * The following example shows how a simple class with a name and a real matrix
721     * should be written:
722     * <pre><code>
723     * public class NamedMatrix implements Serializable {
724     *
725     *     private final String name;
726     *     private final transient RealMatrix coefficients;
727     *
728     *     // omitted constructors, getters ...
729     *
730     *     private void writeObject(ObjectOutputStream oos) throws IOException {
731     *         oos.defaultWriteObject();  // takes care of name field
732     *         MatrixUtils.serializeRealMatrix(coefficients, oos);
733     *     }
734     *
735     *     private void readObject(ObjectInputStream ois) throws ClassNotFoundException, IOException {
736     *         ois.defaultReadObject();  // takes care of name field
737     *         MatrixUtils.deserializeRealMatrix(this, "coefficients", ois);
738     *     }
739     *
740     * }
741     * </code></pre>
742     *
743     * @param matrix real matrix to serialize
744     * @param oos stream where the real matrix should be written
745     * @exception IOException if object cannot be written to stream
746     * @see #deserializeRealMatrix(Object, String, ObjectInputStream)
747     */
748    public static void serializeRealMatrix(final RealMatrix matrix,
749                                           final ObjectOutputStream oos)
750        throws IOException {
751        final int n = matrix.getRowDimension();
752        final int m = matrix.getColumnDimension();
753        oos.writeInt(n);
754        oos.writeInt(m);
755        for (int i = 0; i < n; ++i) {
756            for (int j = 0; j < m; ++j) {
757                oos.writeDouble(matrix.getEntry(i, j));
758            }
759        }
760    }
761
762    /** Deserialize  a {@link RealMatrix} field in a class.
763     * <p>
764     * This method is intended to be called from within a private
765     * <code>readObject</code> method (after a call to
766     * <code>ois.defaultReadObject()</code>) in a class that has a
767     * {@link RealMatrix} field, which should be declared <code>transient</code>.
768     * This way, the default handling does not deserialize the matrix (the {@link
769     * RealMatrix} interface is not serializable by default) but this method does
770     * deserialize it specifically.
771     * </p>
772     * @param instance instance in which the field must be set up
773     * @param fieldName name of the field within the class (may be private and final)
774     * @param ois stream from which the real matrix should be read
775     * @exception ClassNotFoundException if a class in the stream cannot be found
776     * @exception IOException if object cannot be read from the stream
777     * @see #serializeRealMatrix(RealMatrix, ObjectOutputStream)
778     */
779    public static void deserializeRealMatrix(final Object instance,
780                                             final String fieldName,
781                                             final ObjectInputStream ois)
782      throws ClassNotFoundException, IOException {
783        try {
784
785            // read the matrix data
786            final int n = ois.readInt();
787            final int m = ois.readInt();
788            final double[][] data = new double[n][m];
789            for (int i = 0; i < n; ++i) {
790                final double[] dataI = data[i];
791                for (int j = 0; j < m; ++j) {
792                    dataI[j] = ois.readDouble();
793                }
794            }
795
796            // create the instance
797            final RealMatrix matrix = new Array2DRowRealMatrix(data, false);
798
799            // set up the field
800            final java.lang.reflect.Field f =
801                instance.getClass().getDeclaredField(fieldName);
802            f.setAccessible(true);
803            f.set(instance, matrix);
804
805        } catch (NoSuchFieldException nsfe) {
806            IOException ioe = new IOException();
807            ioe.initCause(nsfe);
808            throw ioe;
809        } catch (IllegalAccessException iae) {
810            IOException ioe = new IOException();
811            ioe.initCause(iae);
812            throw ioe;
813        }
814    }
815
816    /**Solve  a  system of composed of a Lower Triangular Matrix
817     * {@link RealMatrix}.
818     * <p>
819     * This method is called to solve systems of equations which are
820     * of the lower triangular form. The matrix {@link RealMatrix}
821     * is assumed, though not checked, to be in lower triangular form.
822     * The vector {@link RealVector} is overwritten with the solution.
823     * The matrix is checked that it is square and its dimensions match
824     * the length of the vector.
825     * </p>
826     * @param rm RealMatrix which is lower triangular
827     * @param b  RealVector this is overwritten
828     * @throws DimensionMismatchException if the matrix and vector are not
829     * conformable
830     * @throws NonSquareMatrixException if the matrix {@code rm} is not square
831     * @throws MathArithmeticException if the absolute value of one of the diagonal
832     * coefficient of {@code rm} is lower than {@link Precision#SAFE_MIN}
833     */
834    public static void solveLowerTriangularSystem(RealMatrix rm, RealVector b)
835        throws DimensionMismatchException, MathArithmeticException,
836        NonSquareMatrixException {
837        if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) {
838            throw new DimensionMismatchException(
839                    (rm == null) ? 0 : rm.getRowDimension(),
840                    (b == null) ? 0 : b.getDimension());
841        }
842        if( rm.getColumnDimension() != rm.getRowDimension() ){
843            throw new NonSquareMatrixException(rm.getRowDimension(),
844                                               rm.getColumnDimension());
845        }
846        int rows = rm.getRowDimension();
847        for( int i = 0 ; i < rows ; i++ ){
848            double diag = rm.getEntry(i, i);
849            if( FastMath.abs(diag) < Precision.SAFE_MIN ){
850                throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
851            }
852            double bi = b.getEntry(i)/diag;
853            b.setEntry(i,  bi );
854            for( int j = i+1; j< rows; j++ ){
855                b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i)  );
856            }
857        }
858    }
859
860    /** Solver a  system composed  of an Upper Triangular Matrix
861     * {@link RealMatrix}.
862     * <p>
863     * This method is called to solve systems of equations which are
864     * of the lower triangular form. The matrix {@link RealMatrix}
865     * is assumed, though not checked, to be in upper triangular form.
866     * The vector {@link RealVector} is overwritten with the solution.
867     * The matrix is checked that it is square and its dimensions match
868     * the length of the vector.
869     * </p>
870     * @param rm RealMatrix which is upper triangular
871     * @param b  RealVector this is overwritten
872     * @throws DimensionMismatchException if the matrix and vector are not
873     * conformable
874     * @throws NonSquareMatrixException if the matrix {@code rm} is not
875     * square
876     * @throws MathArithmeticException if the absolute value of one of the diagonal
877     * coefficient of {@code rm} is lower than {@link Precision#SAFE_MIN}
878     */
879    public static void solveUpperTriangularSystem(RealMatrix rm, RealVector b)
880        throws DimensionMismatchException, MathArithmeticException,
881        NonSquareMatrixException {
882        if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) {
883            throw new DimensionMismatchException(
884                    (rm == null) ? 0 : rm.getRowDimension(),
885                    (b == null) ? 0 : b.getDimension());
886        }
887        if( rm.getColumnDimension() != rm.getRowDimension() ){
888            throw new NonSquareMatrixException(rm.getRowDimension(),
889                                               rm.getColumnDimension());
890        }
891        int rows = rm.getRowDimension();
892        for( int i = rows-1 ; i >-1 ; i-- ){
893            double diag = rm.getEntry(i, i);
894            if( FastMath.abs(diag) < Precision.SAFE_MIN ){
895                throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
896            }
897            double bi = b.getEntry(i)/diag;
898            b.setEntry(i,  bi );
899            for( int j = i-1; j>-1; j-- ){
900                b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i)  );
901            }
902        }
903    }
904
905    /**
906     * Computes the inverse of the given matrix by splitting it into
907     * 4 sub-matrices.
908     *
909     * @param m Matrix whose inverse must be computed.
910     * @param splitIndex Index that determines the "split" line and
911     * column.
912     * The element corresponding to this index will part of the
913     * upper-left sub-matrix.
914     * @return the inverse of {@code m}.
915     * @throws NonSquareMatrixException if {@code m} is not square.
916     */
917    public static RealMatrix blockInverse(RealMatrix m,
918                                          int splitIndex) {
919        final int n = m.getRowDimension();
920        if (m.getColumnDimension() != n) {
921            throw new NonSquareMatrixException(m.getRowDimension(),
922                                               m.getColumnDimension());
923        }
924
925        final int splitIndex1 = splitIndex + 1;
926
927        final RealMatrix a = m.getSubMatrix(0, splitIndex, 0, splitIndex);
928        final RealMatrix b = m.getSubMatrix(0, splitIndex, splitIndex1, n - 1);
929        final RealMatrix c = m.getSubMatrix(splitIndex1, n - 1, 0, splitIndex);
930        final RealMatrix d = m.getSubMatrix(splitIndex1, n - 1, splitIndex1, n - 1);
931
932        final SingularValueDecomposition aDec = new SingularValueDecomposition(a);
933        final DecompositionSolver aSolver = aDec.getSolver();
934        if (!aSolver.isNonSingular()) {
935            throw new SingularMatrixException();
936        }
937        final RealMatrix aInv = aSolver.getInverse();
938
939        final SingularValueDecomposition dDec = new SingularValueDecomposition(d);
940        final DecompositionSolver dSolver = dDec.getSolver();
941        if (!dSolver.isNonSingular()) {
942            throw new SingularMatrixException();
943        }
944        final RealMatrix dInv = dSolver.getInverse();
945
946        final RealMatrix tmp1 = a.subtract(b.multiply(dInv).multiply(c));
947        final SingularValueDecomposition tmp1Dec = new SingularValueDecomposition(tmp1);
948        final DecompositionSolver tmp1Solver = tmp1Dec.getSolver();
949        if (!tmp1Solver.isNonSingular()) {
950            throw new SingularMatrixException();
951        }
952        final RealMatrix result00 = tmp1Solver.getInverse();
953
954        final RealMatrix tmp2 = d.subtract(c.multiply(aInv).multiply(b));
955        final SingularValueDecomposition tmp2Dec = new SingularValueDecomposition(tmp2);
956        final DecompositionSolver tmp2Solver = tmp2Dec.getSolver();
957        if (!tmp2Solver.isNonSingular()) {
958            throw new SingularMatrixException();
959        }
960        final RealMatrix result11 = tmp2Solver.getInverse();
961
962        final RealMatrix result01 = aInv.multiply(b).multiply(result11).scalarMultiply(-1);
963        final RealMatrix result10 = dInv.multiply(c).multiply(result00).scalarMultiply(-1);
964
965        final RealMatrix result = new Array2DRowRealMatrix(n, n);
966        result.setSubMatrix(result00.getData(), 0, 0);
967        result.setSubMatrix(result01.getData(), 0, splitIndex1);
968        result.setSubMatrix(result10.getData(), splitIndex1, 0);
969        result.setSubMatrix(result11.getData(), splitIndex1, splitIndex1);
970
971        return result;
972    }
973
974    /**
975     * Computes the inverse of the given matrix.
976     * <p>
977     * By default, the inverse of the matrix is computed using the QR-decomposition,
978     * unless a more efficient method can be determined for the input matrix.
979     * <p>
980     * Note: this method will use a singularity threshold of 0,
981     * use {@link #inverse(RealMatrix, double)} if a different threshold is needed.
982     *
983     * @param matrix Matrix whose inverse shall be computed
984     * @return the inverse of {@code matrix}
985     * @throws NullArgumentException if {@code matrix} is {@code null}
986     * @throws SingularMatrixException if m is singular
987     * @throws NonSquareMatrixException if matrix is not square
988     * @since 3.3
989     */
990    public static RealMatrix inverse(RealMatrix matrix)
991            throws NullArgumentException, SingularMatrixException, NonSquareMatrixException {
992        return inverse(matrix, 0);
993    }
994
995    /**
996     * Computes the inverse of the given matrix.
997     * <p>
998     * By default, the inverse of the matrix is computed using the QR-decomposition,
999     * unless a more efficient method can be determined for the input matrix.
1000     *
1001     * @param matrix Matrix whose inverse shall be computed
1002     * @param threshold Singularity threshold
1003     * @return the inverse of {@code m}
1004     * @throws NullArgumentException if {@code matrix} is {@code null}
1005     * @throws SingularMatrixException if matrix is singular
1006     * @throws NonSquareMatrixException if matrix is not square
1007     * @since 3.3
1008     */
1009    public static RealMatrix inverse(RealMatrix matrix, double threshold)
1010            throws NullArgumentException, SingularMatrixException, NonSquareMatrixException {
1011
1012        MathUtils.checkNotNull(matrix);
1013
1014        if (!matrix.isSquare()) {
1015            throw new NonSquareMatrixException(matrix.getRowDimension(),
1016                                               matrix.getColumnDimension());
1017        }
1018
1019        if (matrix instanceof DiagonalMatrix) {
1020            return ((DiagonalMatrix) matrix).inverse(threshold);
1021        } else {
1022            QRDecomposition decomposition = new QRDecomposition(matrix, threshold);
1023            return decomposition.getSolver().getInverse();
1024        }
1025    }
1026}