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