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×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×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×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×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 }