MatrixUtils.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      https://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.linear;

  25. import java.util.ArrayList;
  26. import java.util.Arrays;
  27. import java.util.List;

  29. import org.hipparchus.CalculusFieldElement;
  30. import org.hipparchus.Field;
  31. import org.hipparchus.FieldElement;
  32. import org.hipparchus.exception.LocalizedCoreFormats;
  33. import org.hipparchus.exception.MathIllegalArgumentException;
  34. import org.hipparchus.exception.MathRuntimeException;
  35. import org.hipparchus.exception.NullArgumentException;
  36. import org.hipparchus.fraction.BigFraction;
  37. import org.hipparchus.fraction.Fraction;
  38. import org.hipparchus.util.FastMath;
  39. import org.hipparchus.util.MathArrays;
  40. import org.hipparchus.util.MathUtils;
  41. import org.hipparchus.util.Precision;

  43. /**
  44.  * A collection of static methods that operate on or return matrices.
  45.  *
  46.  */
  47. public class MatrixUtils {

  49.     /**
  50.      * The default format for {@link RealMatrix} objects.
  51.      */
  52.     public static final RealMatrixFormat DEFAULT_FORMAT = RealMatrixFormat.getRealMatrixFormat();

  54.     /**
  55.      * A format for {@link RealMatrix} objects compatible with octave.
  56.      */
  57.     public static final RealMatrixFormat OCTAVE_FORMAT = new RealMatrixFormat("[", "]", "", "", "; ", ", ");

  59.     /** Pade coefficients required for the matrix exponential calculation. */
  60.     private static final double[] PADE_COEFFICIENTS_3 = {
  61.             120.0,
  62.             60.0,
  63.             12.0,
  64.             1.0
  65.     };

  67.     /** Pade coefficients required for the matrix exponential calculation. */
  68.     private static final double[] PADE_COEFFICIENTS_5 = {
  69.             30240.0,
  70.             15120.0,
  71.             3360.0,
  72.             420.0,
  73.             30.0,
  74.             1
  75.     };

  77.     /** Pade coefficients required for the matrix exponential calculation. */
  78.     private static final double[] PADE_COEFFICIENTS_7 = {
  79.             17297280.0,
  80.             8648640.0,
  81.             1995840.0,
  82.             277200.0,
  83.             25200.0,
  84.             1512.0,
  85.             56.0,
  86.             1.0
  87.     };

  89.     /** Pade coefficients required for the matrix exponential calculation. */
  90.     private static final double[] PADE_COEFFICIENTS_9 = {
  91.             17643225600.0,
  92.             8821612800.0,
  93.             2075673600.0,
  94.             302702400.0,
  95.             30270240.0,
  96.             2162160.0,
  97.             110880.0,
  98.             3960.0,
  99.             90.0,
  100.             1.0
  101.     };

  103.     /** Pade coefficients required for the matrix exponential calculation. */
  104.     private static final double[] PADE_COEFFICIENTS_13 = {
  105.             6.476475253248e+16,
  106.             3.238237626624e+16,
  107.             7.7717703038976e+15,
  108.             1.1873537964288e+15,
  109.             129060195264000.0,
  110.             10559470521600.0,
  111.             670442572800.0,
  112.             33522128640.0,
  113.             1323241920.0,
  114.             40840800.0,
  115.             960960.0,
  116.             16380.0,
  117.             182.0,
  118.             1.0
  119.     };

  121.     /**
  122.      * Private constructor.
  123.      */
  124.     private MatrixUtils() {
  125.         super();
  126.     }

  128.     /**
  129.      * Returns a {@link RealMatrix} with specified dimensions.
  130.      * <p>The type of matrix returned depends on the dimension. Below
  131.      * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
  132.      * square matrix) which can be stored in a 32kB array, a {@link
  133.      * Array2DRowRealMatrix} instance is built. Above this threshold a {@link
  134.      * BlockRealMatrix} instance is built.</p>
  135.      * <p>The matrix elements are all set to 0.0.</p>
  136.      * @param rows number of rows of the matrix
  137.      * @param columns number of columns of the matrix
  138.      * @return  RealMatrix with specified dimensions
  139.      * @see #createRealMatrix(double[][])
  140.      */
  141.     public static RealMatrix createRealMatrix(final int rows, final int columns) {
  142.         return (rows * columns <= 4096) ?
  143.                 new Array2DRowRealMatrix(rows, columns) : new BlockRealMatrix(rows, columns);
  144.     }

  146.     /**
  147.      * Returns a {@link FieldMatrix} with specified dimensions.
  148.      * <p>The type of matrix returned depends on the dimension. Below
  149.      * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
  150.      * square matrix), a {@link FieldMatrix} instance is built. Above
  151.      * this threshold a {@link BlockFieldMatrix} instance is built.</p>
  152.      * <p>The matrix elements are all set to field.getZero().</p>
  153.      * @param <T> the type of the field elements
  154.      * @param field field to which the matrix elements belong
  155.      * @param rows number of rows of the matrix
  156.      * @param columns number of columns of the matrix
  157.      * @return  FieldMatrix with specified dimensions
  158.      * @see #createFieldMatrix(FieldElement[][])
  159.      */
  160.     public static <T extends FieldElement<T>> FieldMatrix<T> createFieldMatrix(final Field<T> field,
  161.                                                                                final int rows,
  162.                                                                                final int columns) {
  163.         return (rows * columns <= 4096) ?
  164.                 new Array2DRowFieldMatrix<>(field, rows, columns) : new BlockFieldMatrix<>(field, rows, columns);
  165.     }

  167.     /**
  168.      * Returns a {@link RealMatrix} whose entries are the the values in the
  169.      * the input array.
  170.      * <p>The type of matrix returned depends on the dimension. Below
  171.      * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
  172.      * square matrix) which can be stored in a 32kB array, a {@link
  173.      * Array2DRowRealMatrix} instance is built. Above this threshold a {@link
  174.      * BlockRealMatrix} instance is built.</p>
  175.      * <p>The input array is copied, not referenced.</p>
  176.      *
  177.      * @param data input array
  178.      * @return  RealMatrix containing the values of the array
  179.      * @throws org.hipparchus.exception.MathIllegalArgumentException
  180.      * if {@code data} is not rectangular (not all rows have the same length).
  181.      * @throws MathIllegalArgumentException if a row or column is empty.
  182.      * @throws NullArgumentException if either {@code data} or {@code data[0]}
  183.      * is {@code null}.
  184.      * @throws MathIllegalArgumentException if {@code data} is not rectangular.
  185.      * @see #createRealMatrix(int, int)
  186.      */
  187.     public static RealMatrix createRealMatrix(double[][] data)
  188.         throws MathIllegalArgumentException, NullArgumentException {
  189.         if (data == null ||
  190.             data[0] == null) {
  191.             throw new NullArgumentException();
  192.         }
  193.         return (data.length * data[0].length <= 4096) ?
  194.                 new Array2DRowRealMatrix(data) : new BlockRealMatrix(data);
  195.     }

  197.     /**
  198.      * Returns a {@link FieldMatrix} whose entries are the the values in the
  199.      * the input array.
  200.      * <p>The type of matrix returned depends on the dimension. Below
  201.      * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
  202.      * square matrix), a {@link FieldMatrix} instance is built. Above
  203.      * this threshold a {@link BlockFieldMatrix} instance is built.</p>
  204.      * <p>The input array is copied, not referenced.</p>
  205.      * @param <T> the type of the field elements
  206.      * @param data input array
  207.      * @return a matrix containing the values of the array.
  208.      * @throws org.hipparchus.exception.MathIllegalArgumentException
  209.      * if {@code data} is not rectangular (not all rows have the same length).
  210.      * @throws MathIllegalArgumentException if a row or column is empty.
  211.      * @throws NullArgumentException if either {@code data} or {@code data[0]}
  212.      * is {@code null}.
  213.      * @see #createFieldMatrix(Field, int, int)
  214.      */
  215.     public static <T extends FieldElement<T>> FieldMatrix<T> createFieldMatrix(T[][] data)
  216.         throws MathIllegalArgumentException, NullArgumentException {
  217.         if (data == null ||
  218.             data[0] == null) {
  219.             throw new NullArgumentException();
  220.         }
  221.         return (data.length * data[0].length <= 4096) ?
  222.                 new Array2DRowFieldMatrix<>(data) : new BlockFieldMatrix<>(data);
  223.     }

  225.     /**
  226.      * Returns <code>dimension x dimension</code> identity matrix.
  227.      *
  228.      * @param dimension dimension of identity matrix to generate
  229.      * @return identity matrix
  230.      * @throws IllegalArgumentException if dimension is not positive
  231.      */
  232.     public static RealMatrix createRealIdentityMatrix(int dimension) {
  233.         final RealMatrix m = createRealMatrix(dimension, dimension);
  234.         for (int i = 0; i < dimension; ++i) {
  235.             m.setEntry(i, i, 1.0);
  236.         }
  237.         return m;
  238.     }

  240.     /**
  241.      * Returns <code>dimension x dimension</code> identity matrix.
  242.      *
  243.      * @param <T> the type of the field elements
  244.      * @param field field to which the elements belong
  245.      * @param dimension dimension of identity matrix to generate
  246.      * @return identity matrix
  247.      * @throws IllegalArgumentException if dimension is not positive
  248.      */
  249.     public static <T extends FieldElement<T>> FieldMatrix<T>
  250.         createFieldIdentityMatrix(final Field<T> field, final int dimension) {
  251.         final T zero = field.getZero();
  252.         final T one  = field.getOne();
  253.         final T[][] d = MathArrays.buildArray(field, dimension, dimension);
  254.         for (int row = 0; row < dimension; row++) {
  255.             final T[] dRow = d[row];
  256.             Arrays.fill(dRow, zero);
  257.             dRow[row] = one;
  258.         }
  259.         return new Array2DRowFieldMatrix<>(field, d, false);
  260.     }

  262.     /**
  263.      * Returns a diagonal matrix with specified elements.
  264.      *
  265.      * @param diagonal diagonal elements of the matrix (the array elements
  266.      * will be copied)
  267.      * @return diagonal matrix
  268.      */
  269.     public static RealMatrix createRealDiagonalMatrix(final double[] diagonal) {
  270.         final RealMatrix m = createRealMatrix(diagonal.length, diagonal.length);
  271.         for (int i = 0; i < diagonal.length; ++i) {
  272.             m.setEntry(i, i, diagonal[i]);
  273.         }
  274.         return m;
  275.     }

  277.     /**
  278.      * Returns a diagonal matrix with specified elements.
  279.      *
  280.      * @param <T> the type of the field elements
  281.      * @param diagonal diagonal elements of the matrix (the array elements
  282.      * will be copied)
  283.      * @return diagonal matrix
  284.      */
  285.     public static <T extends FieldElement<T>> FieldMatrix<T>
  286.         createFieldDiagonalMatrix(final T[] diagonal) {
  287.         final FieldMatrix<T> m =
  288.             createFieldMatrix(diagonal[0].getField(), diagonal.length, diagonal.length);
  289.         for (int i = 0; i < diagonal.length; ++i) {
  290.             m.setEntry(i, i, diagonal[i]);
  291.         }
  292.         return m;
  293.     }

  295.     /**
  296.      * Creates a {@link RealVector} using the data from the input array.
  297.      *
  298.      * @param data the input data
  299.      * @return a data.length RealVector
  300.      * @throws MathIllegalArgumentException if {@code data} is empty.
  301.      * @throws NullArgumentException if {@code data} is {@code null}.
  302.      */
  303.     public static RealVector createRealVector(double[] data)
  304.         throws MathIllegalArgumentException, NullArgumentException {
  305.         if (data == null) {
  306.             throw new NullArgumentException();
  307.         }
  308.         return new ArrayRealVector(data, true);
  309.     }

  311.     /**
  312.      * Creates a {@link RealVector} with specified dimensions.
  313.      *
  314.      * @param dimension dimension of the vector
  315.      * @return a new vector
  316.      * @since 1.3
  317.      */
  318.     public static RealVector createRealVector(final int dimension) {
  319.         return new ArrayRealVector(new double[dimension]);
  320.     }

  322.     /**
  323.      * Creates a {@link FieldVector} using the data from the input array.
  324.      *
  325.      * @param <T> the type of the field elements
  326.      * @param data the input data
  327.      * @return a data.length FieldVector
  328.      * @throws MathIllegalArgumentException if {@code data} is empty.
  329.      * @throws NullArgumentException if {@code data} is {@code null}.
  330.      * @throws MathIllegalArgumentException if {@code data} has 0 elements
  331.      */
  332.     public static <T extends FieldElement<T>> FieldVector<T> createFieldVector(final T[] data)
  333.         throws MathIllegalArgumentException, NullArgumentException {
  334.         if (data == null) {
  335.             throw new NullArgumentException();
  336.         }
  337.         if (data.length == 0) {
  338.             throw new MathIllegalArgumentException(LocalizedCoreFormats.VECTOR_MUST_HAVE_AT_LEAST_ONE_ELEMENT);
  339.         }
  340.         return new ArrayFieldVector<>(data[0].getField(), data, true);
  341.     }

  343.     /**
  344.      * Creates a {@link FieldVector} with specified dimensions.
  345.      *
  346.      * @param <T> the type of the field elements
  347.      * @param field field to which array elements belong
  348.      * @param dimension dimension of the vector
  349.      * @return a new vector
  350.      * @since 1.3
  351.      */
  352.     public static <T extends FieldElement<T>> FieldVector<T> createFieldVector(final Field<T> field, final int dimension) {
  353.         return new ArrayFieldVector<>(MathArrays.buildArray(field, dimension));
  354.     }

  356.     /**
  357.      * Create a row {@link RealMatrix} using the data from the input
  358.      * array.
  359.      *
  360.      * @param rowData the input row data
  361.      * @return a 1 x rowData.length RealMatrix
  362.      * @throws MathIllegalArgumentException if {@code rowData} is empty.
  363.      * @throws NullArgumentException if {@code rowData} is {@code null}.
  364.      */
  365.     public static RealMatrix createRowRealMatrix(double[] rowData)
  366.         throws MathIllegalArgumentException, NullArgumentException {
  367.         if (rowData == null) {
  368.             throw new NullArgumentException();
  369.         }
  370.         final int nCols = rowData.length;
  371.         final RealMatrix m = createRealMatrix(1, nCols);
  372.         for (int i = 0; i < nCols; ++i) {
  373.             m.setEntry(0, i, rowData[i]);
  374.         }
  375.         return m;
  376.     }

  378.     /**
  379.      * Create a row {@link FieldMatrix} using the data from the input
  380.      * array.
  381.      *
  382.      * @param <T> the type of the field elements
  383.      * @param rowData the input row data
  384.      * @return a 1 x rowData.length FieldMatrix
  385.      * @throws MathIllegalArgumentException if {@code rowData} is empty.
  386.      * @throws NullArgumentException if {@code rowData} is {@code null}.
  387.      */
  388.     public static <T extends FieldElement<T>> FieldMatrix<T>
  389.         createRowFieldMatrix(final T[] rowData)
  390.         throws MathIllegalArgumentException, NullArgumentException {
  391.         if (rowData == null) {
  392.             throw new NullArgumentException();
  393.         }
  394.         final int nCols = rowData.length;
  395.         if (nCols == 0) {
  396.             throw new MathIllegalArgumentException(LocalizedCoreFormats.AT_LEAST_ONE_COLUMN);
  397.         }
  398.         final FieldMatrix<T> m = createFieldMatrix(rowData[0].getField(), 1, nCols);
  399.         for (int i = 0; i < nCols; ++i) {
  400.             m.setEntry(0, i, rowData[i]);
  401.         }
  402.         return m;
  403.     }

  405.     /**
  406.      * Creates a column {@link RealMatrix} using the data from the input
  407.      * array.
  408.      *
  409.      * @param columnData  the input column data
  410.      * @return a columnData x 1 RealMatrix
  411.      * @throws MathIllegalArgumentException if {@code columnData} is empty.
  412.      * @throws NullArgumentException if {@code columnData} is {@code null}.
  413.      */
  414.     public static RealMatrix createColumnRealMatrix(double[] columnData)
  415.         throws MathIllegalArgumentException, NullArgumentException {
  416.         if (columnData == null) {
  417.             throw new NullArgumentException();
  418.         }
  419.         final int nRows = columnData.length;
  420.         final RealMatrix m = createRealMatrix(nRows, 1);
  421.         for (int i = 0; i < nRows; ++i) {
  422.             m.setEntry(i, 0, columnData[i]);
  423.         }
  424.         return m;
  425.     }

  427.     /**
  428.      * Creates a column {@link FieldMatrix} using the data from the input
  429.      * array.
  430.      *
  431.      * @param <T> the type of the field elements
  432.      * @param columnData  the input column data
  433.      * @return a columnData x 1 FieldMatrix
  434.      * @throws MathIllegalArgumentException if {@code data} is empty.
  435.      * @throws NullArgumentException if {@code columnData} is {@code null}.
  436.      */
  437.     public static <T extends FieldElement<T>> FieldMatrix<T>
  438.         createColumnFieldMatrix(final T[] columnData)
  439.         throws MathIllegalArgumentException, NullArgumentException {
  440.         if (columnData == null) {
  441.             throw new NullArgumentException();
  442.         }
  443.         final int nRows = columnData.length;
  444.         if (nRows == 0) {
  445.             throw new MathIllegalArgumentException(LocalizedCoreFormats.AT_LEAST_ONE_ROW);
  446.         }
  447.         final FieldMatrix<T> m = createFieldMatrix(columnData[0].getField(), nRows, 1);
  448.         for (int i = 0; i < nRows; ++i) {
  449.             m.setEntry(i, 0, columnData[i]);
  450.         }
  451.         return m;
  452.     }

  454.     /**
  455.      * Checks whether a matrix is symmetric, within a given relative tolerance.
  456.      *
  457.      * @param matrix Matrix to check.
  458.      * @param relativeTolerance Tolerance of the symmetry check.
  459.      * @param raiseException If {@code true}, an exception will be raised if
  460.      * the matrix is not symmetric.
  461.      * @return {@code true} if {@code matrix} is symmetric.
  462.      * @throws MathIllegalArgumentException if the matrix is not square.
  463.      * @throws MathIllegalArgumentException if the matrix is not symmetric.
  464.      */
  465.     private static boolean isSymmetricInternal(RealMatrix matrix,
  466.                                                double relativeTolerance,
  467.                                                boolean raiseException) {
  468.         final int rows = matrix.getRowDimension();
  469.         if (rows != matrix.getColumnDimension()) {
  470.             if (raiseException) {
  471.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
  472.                                                        rows, matrix.getColumnDimension());
  473.             } else {
  474.                 return false;
  475.             }
  476.         }
  477.         for (int i = 0; i < rows; i++) {
  478.             for (int j = i + 1; j < rows; j++) {
  479.                 final double mij = matrix.getEntry(i, j);
  480.                 final double mji = matrix.getEntry(j, i);
  481.                 if (FastMath.abs(mij - mji) >
  482.                     FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * relativeTolerance) {
  483.                     if (raiseException) {
  484.                         throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SYMMETRIC_MATRIX,
  485.                                                                i, j, relativeTolerance);
  486.                     } else {
  487.                         return false;
  488.                     }
  489.                 }
  490.             }
  491.         }
  492.         return true;
  493.     }

  495.     /**
  496.      * Checks whether a matrix is symmetric.
  497.      *
  498.      * @param matrix Matrix to check.
  499.      * @param eps Relative tolerance.
  500.      * @throws MathIllegalArgumentException if the matrix is not square.
  501.      * @throws MathIllegalArgumentException if the matrix is not symmetric.
  502.      */
  503.     public static void checkSymmetric(RealMatrix matrix,
  504.                                       double eps) {
  505.         isSymmetricInternal(matrix, eps, true);
  506.     }

  508.     /**
  509.      * Checks whether a matrix is symmetric.
  510.      *
  511.      * @param matrix Matrix to check.
  512.      * @param eps Relative tolerance.
  513.      * @return {@code true} if {@code matrix} is symmetric.
  514.      */
  515.     public static boolean isSymmetric(RealMatrix matrix,
  516.                                       double eps) {
  517.         return isSymmetricInternal(matrix, eps, false);
  518.     }

  520.     /**
  521.      * Check if matrix indices are valid.
  522.      *
  523.      * @param m Matrix.
  524.      * @param row Row index to check.
  525.      * @param column Column index to check.
  526.      * @throws MathIllegalArgumentException if {@code row} or {@code column} is not
  527.      * a valid index.
  528.      */
  529.     public static void checkMatrixIndex(final AnyMatrix m,
  530.                                         final int row, final int column)
  531.         throws MathIllegalArgumentException {
  532.         checkRowIndex(m, row);
  533.         checkColumnIndex(m, column);
  534.     }

  536.     /**
  537.      * Check if a row index is valid.
  538.      *
  539.      * @param m Matrix.
  540.      * @param row Row index to check.
  541.      * @throws MathIllegalArgumentException if {@code row} is not a valid index.
  542.      */
  543.     public static void checkRowIndex(final AnyMatrix m, final int row)
  544.         throws MathIllegalArgumentException {
  545.         if (row < 0 ||
  546.             row >= m.getRowDimension()) {
  547.             throw new MathIllegalArgumentException(LocalizedCoreFormats.ROW_INDEX,
  548.                                           row, 0, m.getRowDimension() - 1);
  549.         }
  550.     }

  552.     /**
  553.      * Check if a column index is valid.
  554.      *
  555.      * @param m Matrix.
  556.      * @param column Column index to check.
  557.      * @throws MathIllegalArgumentException if {@code column} is not a valid index.
  558.      */
  559.     public static void checkColumnIndex(final AnyMatrix m, final int column)
  560.         throws MathIllegalArgumentException {
  561.         if (column < 0 || column >= m.getColumnDimension()) {
  562.             throw new MathIllegalArgumentException(LocalizedCoreFormats.COLUMN_INDEX,
  563.                                            column, 0, m.getColumnDimension() - 1);
  564.         }
  565.     }

  567.     /**
  568.      * Check if submatrix ranges indices are valid.
  569.      * Rows and columns are indicated counting from 0 to {@code n - 1}.
  570.      *
  571.      * @param m Matrix.
  572.      * @param startRow Initial row index.
  573.      * @param endRow Final row index.
  574.      * @param startColumn Initial column index.
  575.      * @param endColumn Final column index.
  576.      * @throws MathIllegalArgumentException if the indices are invalid.
  577.      * @throws MathIllegalArgumentException if {@code endRow < startRow} or
  578.      * {@code endColumn < startColumn}.
  579.      */
  580.     public static void checkSubMatrixIndex(final AnyMatrix m,
  581.                                            final int startRow, final int endRow,
  582.                                            final int startColumn, final int endColumn)
  583.         throws MathIllegalArgumentException {
  584.         checkRowIndex(m, startRow);
  585.         checkRowIndex(m, endRow);
  586.         if (endRow < startRow) {
  587.             throw new MathIllegalArgumentException(LocalizedCoreFormats.INITIAL_ROW_AFTER_FINAL_ROW,
  588.                                                 endRow, startRow, false);
  589.         }

  591.         checkColumnIndex(m, startColumn);
  592.         checkColumnIndex(m, endColumn);
  593.         if (endColumn < startColumn) {
  594.             throw new MathIllegalArgumentException(LocalizedCoreFormats.INITIAL_COLUMN_AFTER_FINAL_COLUMN,
  595.                                                 endColumn, startColumn, false);
  596.         }


  599.     }

  601.     /**
  602.      * Check if submatrix ranges indices are valid.
  603.      * Rows and columns are indicated counting from 0 to n-1.
  604.      *
  605.      * @param m Matrix.
  606.      * @param selectedRows Array of row indices.
  607.      * @param selectedColumns Array of column indices.
  608.      * @throws NullArgumentException if {@code selectedRows} or
  609.      * {@code selectedColumns} are {@code null}.
  610.      * @throws MathIllegalArgumentException if the row or column selections are empty (zero
  611.      * length).
  612.      * @throws MathIllegalArgumentException if row or column selections are not valid.
  613.      */
  614.     public static void checkSubMatrixIndex(final AnyMatrix m,
  615.                                            final int[] selectedRows,
  616.                                            final int[] selectedColumns)
  617.         throws MathIllegalArgumentException, NullArgumentException {
  618.         if (selectedRows == null) {
  619.             throw new NullArgumentException();
  620.         }
  621.         if (selectedColumns == null) {
  622.             throw new NullArgumentException();
  623.         }
  624.         if (selectedRows.length == 0) {
  625.             throw new MathIllegalArgumentException(LocalizedCoreFormats.EMPTY_SELECTED_ROW_INDEX_ARRAY);
  626.         }
  627.         if (selectedColumns.length == 0) {
  628.             throw new MathIllegalArgumentException(LocalizedCoreFormats.EMPTY_SELECTED_COLUMN_INDEX_ARRAY);
  629.         }

  631.         for (final int row : selectedRows) {
  632.             checkRowIndex(m, row);
  633.         }
  634.         for (final int column : selectedColumns) {
  635.             checkColumnIndex(m, column);
  636.         }
  637.     }

  639.     /**
  640.      * Check if matrices are addition compatible.
  641.      *
  642.      * @param left Left hand side matrix.
  643.      * @param right Right hand side matrix.
  644.      * @throws MathIllegalArgumentException if the matrices are not addition
  645.      * compatible.
  646.      */
  647.     public static void checkAdditionCompatible(final AnyMatrix left, final AnyMatrix right)
  648.         throws MathIllegalArgumentException {
  649.         if ((left.getRowDimension()    != right.getRowDimension()) ||
  650.             (left.getColumnDimension() != right.getColumnDimension())) {
  651.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH_2x2,
  652.                                                    left.getRowDimension(), left.getColumnDimension(),
  653.                                                    right.getRowDimension(), right.getColumnDimension());
  654.         }
  655.     }

  657.     /**
  658.      * Check if matrices are subtraction compatible
  659.      *
  660.      * @param left Left hand side matrix.
  661.      * @param right Right hand side matrix.
  662.      * @throws MathIllegalArgumentException if the matrices are not addition
  663.      * compatible.
  664.      */
  665.     public static void checkSubtractionCompatible(final AnyMatrix left, final AnyMatrix right)
  666.         throws MathIllegalArgumentException {
  667.         if ((left.getRowDimension()    != right.getRowDimension()) ||
  668.             (left.getColumnDimension() != right.getColumnDimension())) {
  669.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH_2x2,
  670.                                                    left.getRowDimension(), left.getColumnDimension(),
  671.                                                    right.getRowDimension(), right.getColumnDimension());
  672.         }
  673.     }

  675.     /**
  676.      * Check if matrices are multiplication compatible
  677.      *
  678.      * @param left Left hand side matrix.
  679.      * @param right Right hand side matrix.
  680.      * @throws MathIllegalArgumentException if matrices are not multiplication
  681.      * compatible.
  682.      */
  683.     public static void checkMultiplicationCompatible(final AnyMatrix left, final AnyMatrix right)
  684.         throws MathIllegalArgumentException {
  685.         if (left.getColumnDimension() != right.getRowDimension()) {
  686.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  687.                                                    left.getColumnDimension(), right.getRowDimension());
  688.         }
  689.     }

  691.     /**
  692.      * Check if matrices have the same number of columns.
  693.      *
  694.      * @param left Left hand side matrix.
  695.      * @param right Right hand side matrix.
  696.      * @throws MathIllegalArgumentException if matrices don't have the same number of columns.
  697.      * @since 1.3
  698.      */
  699.     public static void checkSameColumnDimension(final AnyMatrix left, final AnyMatrix right)
  700.         throws MathIllegalArgumentException {
  701.         if (left.getColumnDimension() != right.getColumnDimension()) {
  702.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  703.                                                    left.getColumnDimension(), right.getColumnDimension());
  704.         }
  705.     }

  707.     /**
  708.      * Check if matrices have the same number of rows.
  709.      *
  710.      * @param left Left hand side matrix.
  711.      * @param right Right hand side matrix.
  712.      * @throws MathIllegalArgumentException if matrices don't have the same number of rows.
  713.      * @since 1.3
  714.      */
  715.     public static void checkSameRowDimension(final AnyMatrix left, final AnyMatrix right)
  716.         throws MathIllegalArgumentException {
  717.         if (left.getRowDimension() != right.getRowDimension()) {
  718.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  719.                                                    left.getRowDimension(), right.getRowDimension());
  720.         }
  721.     }

  723.     /**
  724.      * Convert a {@link FieldMatrix}/{@link Fraction} matrix to a {@link RealMatrix}.
  725.      * @param m Matrix to convert.
  726.      * @return the converted matrix.
  727.      */
  728.     public static Array2DRowRealMatrix fractionMatrixToRealMatrix(final FieldMatrix<Fraction> m) {
  729.         final FractionMatrixConverter converter = new FractionMatrixConverter();
  730.         m.walkInOptimizedOrder(converter);
  731.         return converter.getConvertedMatrix();
  732.     }

  734.     /** Converter for {@link FieldMatrix}/{@link Fraction}. */
  735.     private static class FractionMatrixConverter extends DefaultFieldMatrixPreservingVisitor<Fraction> {
  736.         /** Converted array. */
  737.         private double[][] data;
  738.         /** Simple constructor. */
  739.         FractionMatrixConverter() {
  740.             super(Fraction.ZERO);
  741.         }

  743.         /** {@inheritDoc} */
  744.         @Override
  745.         public void start(int rows, int columns,
  746.                           int startRow, int endRow, int startColumn, int endColumn) {
  747.             data = new double[rows][columns];
  748.         }

  750.         /** {@inheritDoc} */
  751.         @Override
  752.         public void visit(int row, int column, Fraction value) {
  753.             data[row][column] = value.doubleValue();
  754.         }

  756.         /**
  757.          * Get the converted matrix.
  758.          *
  759.          * @return the converted matrix.
  760.          */
  761.         Array2DRowRealMatrix getConvertedMatrix() {
  762.             return new Array2DRowRealMatrix(data, false);
  763.         }

  765.     }

  767.     /**
  768.      * Convert a {@link FieldMatrix}/{@link BigFraction} matrix to a {@link RealMatrix}.
  769.      *
  770.      * @param m Matrix to convert.
  771.      * @return the converted matrix.
  772.      */
  773.     public static Array2DRowRealMatrix bigFractionMatrixToRealMatrix(final FieldMatrix<BigFraction> m) {
  774.         final BigFractionMatrixConverter converter = new BigFractionMatrixConverter();
  775.         m.walkInOptimizedOrder(converter);
  776.         return converter.getConvertedMatrix();
  777.     }

  779.     /** Converter for {@link FieldMatrix}/{@link BigFraction}. */
  780.     private static class BigFractionMatrixConverter extends DefaultFieldMatrixPreservingVisitor<BigFraction> {
  781.         /** Converted array. */
  782.         private double[][] data;
  783.         /** Simple constructor. */
  784.         BigFractionMatrixConverter() {
  785.             super(BigFraction.ZERO);
  786.         }

  788.         /** {@inheritDoc} */
  789.         @Override
  790.         public void start(int rows, int columns,
  791.                           int startRow, int endRow, int startColumn, int endColumn) {
  792.             data = new double[rows][columns];
  793.         }

  795.         /** {@inheritDoc} */
  796.         @Override
  797.         public void visit(int row, int column, BigFraction value) {
  798.             data[row][column] = value.doubleValue();
  799.         }

  801.         /**
  802.          * Get the converted matrix.
  803.          *
  804.          * @return the converted matrix.
  805.          */
  806.         Array2DRowRealMatrix getConvertedMatrix() {
  807.             return new Array2DRowRealMatrix(data, false);
  808.         }
  809.     }

  811.     /**Solve  a  system of composed of a Lower Triangular Matrix
  812.      * {@link RealMatrix}.
  813.      * <p>
  814.      * This method is called to solve systems of equations which are
  815.      * of the lower triangular form. The matrix {@link RealMatrix}
  816.      * is assumed, though not checked, to be in lower triangular form.
  817.      * The vector {@link RealVector} is overwritten with the solution.
  818.      * The matrix is checked that it is square and its dimensions match
  819.      * the length of the vector.
  820.      * </p>
  821.      * @param rm RealMatrix which is lower triangular
  822.      * @param b  RealVector this is overwritten
  823.      * @throws MathIllegalArgumentException if the matrix and vector are not
  824.      * conformable
  825.      * @throws MathIllegalArgumentException if the matrix {@code rm} is not square
  826.      * @throws MathRuntimeException if the absolute value of one of the diagonal
  827.      * coefficient of {@code rm} is lower than {@link Precision#SAFE_MIN}
  828.      */
  829.     public static void solveLowerTriangularSystem(RealMatrix rm, RealVector b)
  830.         throws MathIllegalArgumentException, MathRuntimeException {
  831.         if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) {
  832.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  833.                                                    (rm == null) ? 0 : rm.getRowDimension(),
  834.                                                    (b  == null) ? 0 : b.getDimension());
  835.         }
  836.         if( rm.getColumnDimension() != rm.getRowDimension() ){
  837.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
  838.                                                    rm.getRowDimension(), rm.getColumnDimension());
  839.         }
  840.         int rows = rm.getRowDimension();
  841.         for( int i = 0 ; i < rows ; i++ ){
  842.             double diag = rm.getEntry(i, i);
  843.             if( FastMath.abs(diag) < Precision.SAFE_MIN ){
  844.                 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
  845.             }
  846.             double bi = b.getEntry(i)/diag;
  847.             b.setEntry(i,  bi );
  848.             for( int j = i+1; j< rows; j++ ){
  849.                 b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i)  );
  850.             }
  851.         }
  852.     }

  854.     /** Solver a  system composed  of an Upper Triangular Matrix
  855.      * {@link RealMatrix}.
  856.      * <p>
  857.      * This method is called to solve systems of equations which are
  858.      * of the lower triangular form. The matrix {@link RealMatrix}
  859.      * is assumed, though not checked, to be in upper triangular form.
  860.      * The vector {@link RealVector} is overwritten with the solution.
  861.      * The matrix is checked that it is square and its dimensions match
  862.      * the length of the vector.
  863.      * </p>
  864.      * @param rm RealMatrix which is upper triangular
  865.      * @param b  RealVector this is overwritten
  866.      * @throws MathIllegalArgumentException if the matrix and vector are not
  867.      * conformable
  868.      * @throws MathIllegalArgumentException if the matrix {@code rm} is not
  869.      * square
  870.      * @throws MathRuntimeException if the absolute value of one of the diagonal
  871.      * coefficient of {@code rm} is lower than {@link Precision#SAFE_MIN}
  872.      */
  873.     public static void solveUpperTriangularSystem(RealMatrix rm, RealVector b)
  874.         throws MathIllegalArgumentException, MathRuntimeException {
  875.         if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) {
  876.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  877.                                                    (rm == null) ? 0 : rm.getRowDimension(),
  878.                                                    (b  == null) ? 0 : b.getDimension());
  879.         }
  880.         if( rm.getColumnDimension() != rm.getRowDimension() ){
  881.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
  882.                                                    rm.getRowDimension(), rm.getColumnDimension());
  883.         }
  884.         int rows = rm.getRowDimension();
  885.         for( int i = rows-1 ; i >-1 ; i-- ){
  886.             double diag = rm.getEntry(i, i);
  887.             if( FastMath.abs(diag) < Precision.SAFE_MIN ){
  888.                 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
  889.             }
  890.             double bi = b.getEntry(i)/diag;
  891.             b.setEntry(i,  bi );
  892.             for( int j = i-1; j>-1; j-- ){
  893.                 b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i)  );
  894.             }
  895.         }
  896.     }

  898.     /**
  899.      * Computes the inverse of the given matrix by splitting it into
  900.      * 4 sub-matrices.
  901.      *
  902.      * @param m Matrix whose inverse must be computed.
  903.      * @param splitIndex Index that determines the "split" line and
  904.      * column.
  905.      * The element corresponding to this index will part of the
  906.      * upper-left sub-matrix.
  907.      * @return the inverse of {@code m}.
  908.      * @throws MathIllegalArgumentException if {@code m} is not square.
  909.      */
  910.     public static RealMatrix blockInverse(RealMatrix m,
  911.                                           int splitIndex) {
  912.         final int n = m.getRowDimension();
  913.         if (m.getColumnDimension() != n) {
  914.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
  915.                                                    m.getRowDimension(), m.getColumnDimension());
  916.         }

  918.         final int splitIndex1 = splitIndex + 1;

  920.         final RealMatrix a = m.getSubMatrix(0, splitIndex, 0, splitIndex);
  921.         final RealMatrix b = m.getSubMatrix(0, splitIndex, splitIndex1, n - 1);
  922.         final RealMatrix c = m.getSubMatrix(splitIndex1, n - 1, 0, splitIndex);
  923.         final RealMatrix d = m.getSubMatrix(splitIndex1, n - 1, splitIndex1, n - 1);

  925.         final SingularValueDecomposition aDec = new SingularValueDecomposition(a);
  926.         final DecompositionSolver aSolver = aDec.getSolver();
  927.         if (!aSolver.isNonSingular()) {
  928.             throw new MathIllegalArgumentException(LocalizedCoreFormats.SINGULAR_MATRIX);
  929.         }
  930.         final RealMatrix aInv = aSolver.getInverse();

  932.         final SingularValueDecomposition dDec = new SingularValueDecomposition(d);
  933.         final DecompositionSolver dSolver = dDec.getSolver();
  934.         if (!dSolver.isNonSingular()) {
  935.             throw new MathIllegalArgumentException(LocalizedCoreFormats.SINGULAR_MATRIX);
  936.         }
  937.         final RealMatrix dInv = dSolver.getInverse();

  939.         final RealMatrix tmp1 = a.subtract(b.multiply(dInv).multiply(c));
  940.         final SingularValueDecomposition tmp1Dec = new SingularValueDecomposition(tmp1);
  941.         final DecompositionSolver tmp1Solver = tmp1Dec.getSolver();
  942.         if (!tmp1Solver.isNonSingular()) {
  943.             throw new MathIllegalArgumentException(LocalizedCoreFormats.SINGULAR_MATRIX);
  944.         }
  945.         final RealMatrix result00 = tmp1Solver.getInverse();

  947.         final RealMatrix tmp2 = d.subtract(c.multiply(aInv).multiply(b));
  948.         final SingularValueDecomposition tmp2Dec = new SingularValueDecomposition(tmp2);
  949.         final DecompositionSolver tmp2Solver = tmp2Dec.getSolver();
  950.         if (!tmp2Solver.isNonSingular()) {
  951.             throw new MathIllegalArgumentException(LocalizedCoreFormats.SINGULAR_MATRIX);
  952.         }
  953.         final RealMatrix result11 = tmp2Solver.getInverse();

  955.         final RealMatrix result01 = aInv.multiply(b).multiply(result11).scalarMultiply(-1);
  956.         final RealMatrix result10 = dInv.multiply(c).multiply(result00).scalarMultiply(-1);

  958.         final RealMatrix result = new Array2DRowRealMatrix(n, n);
  959.         result.setSubMatrix(result00.getData(), 0, 0);
  960.         result.setSubMatrix(result01.getData(), 0, splitIndex1);
  961.         result.setSubMatrix(result10.getData(), splitIndex1, 0);
  962.         result.setSubMatrix(result11.getData(), splitIndex1, splitIndex1);

  964.         return result;
  965.     }

  967.     /**
  968.      * Computes the inverse of the given matrix.
  969.      * <p>
  970.      * By default, the inverse of the matrix is computed using the QR-decomposition,
  971.      * unless a more efficient method can be determined for the input matrix.
  972.      * <p>
  973.      * Note: this method will use a singularity threshold of 0,
  974.      * use {@link #inverse(RealMatrix, double)} if a different threshold is needed.
  975.      *
  976.      * @param matrix Matrix whose inverse shall be computed
  977.      * @return the inverse of {@code matrix}
  978.      * @throws NullArgumentException if {@code matrix} is {@code null}
  979.      * @throws MathIllegalArgumentException if m is singular
  980.      * @throws MathIllegalArgumentException if matrix is not square
  981.      */
  982.     public static RealMatrix inverse(RealMatrix matrix)
  983.             throws MathIllegalArgumentException, NullArgumentException {
  984.         return inverse(matrix, 0);
  985.     }

  987.     /**
  988.      * Computes the inverse of the given matrix.
  989.      * <p>
  990.      * By default, the inverse of the matrix is computed using the QR-decomposition,
  991.      * unless a more efficient method can be determined for the input matrix.
  992.      *
  993.      * @param matrix Matrix whose inverse shall be computed
  994.      * @param threshold Singularity threshold
  995.      * @return the inverse of {@code m}
  996.      * @throws NullArgumentException if {@code matrix} is {@code null}
  997.      * @throws MathIllegalArgumentException if matrix is singular
  998.      * @throws MathIllegalArgumentException if matrix is not square
  999.      */
  1000.     public static RealMatrix inverse(RealMatrix matrix, double threshold)
  1001.             throws MathIllegalArgumentException, NullArgumentException {

  1003.         MathUtils.checkNotNull(matrix);

  1005.         if (!matrix.isSquare()) {
  1006.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
  1007.                                                    matrix.getRowDimension(), matrix.getColumnDimension());
  1008.         }

  1010.         if (matrix instanceof DiagonalMatrix) {
  1011.             return ((DiagonalMatrix) matrix).inverse(threshold);
  1012.         } else {
  1013.             QRDecomposition decomposition = new QRDecomposition(matrix, threshold);
  1014.             return decomposition.getSolver().getInverse();
  1015.         }
  1016.     }

  1018.     /**
  1019.      * Computes the <a href="https://mathworld.wolfram.com/MatrixExponential.html">
  1020.      * matrix exponential</a> of the given matrix.
  1021.      *
  1022.      * The algorithm implementation follows the Pade approximant method of
  1023.      * <p>Higham, Nicholas J. “The Scaling and Squaring Method for the Matrix Exponential
  1024.      * Revisited.” SIAM Journal on Matrix Analysis and Applications 26, no. 4 (January 2005): 1179–93.</p>
  1025.      *
  1026.      * @param rm RealMatrix whose inverse shall be computed
  1027.      * @return The inverse of {@code rm}
  1028.      * @throws MathIllegalArgumentException if matrix is not square
  1029.      */
  1030.     public static RealMatrix matrixExponential(final RealMatrix rm) {

  1032.         // Check that the input matrix is square
  1033.         if (!rm.isSquare()) {
  1034.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
  1035.                                                    rm.getRowDimension(), rm.getColumnDimension());
  1036.         }

  1038.         // Copy input matrix
  1039.         int dim = rm.getRowDimension();
  1040.         final RealMatrix identity = createRealIdentityMatrix(dim);
  1041.         RealMatrix scaledMatrix = rm.copy();

  1043.         // Select pade degree required
  1044.         final double l1Norm = rm.getNorm1();
  1045.         double[] padeCoefficients;
  1046.         int squaringCount = 0;

  1048.         if (l1Norm < 1.495585217958292e-2) {
  1049.             padeCoefficients = PADE_COEFFICIENTS_3;
  1050.         } else if (l1Norm < 2.539398330063230e-1) {
  1051.             padeCoefficients = PADE_COEFFICIENTS_5;
  1052.         } else if (l1Norm < 9.504178996162932e-1) {
  1053.             padeCoefficients = PADE_COEFFICIENTS_7;
  1054.         } else if (l1Norm < 2.097847961257068) {
  1055.             padeCoefficients = PADE_COEFFICIENTS_9;
  1056.         } else {
  1057.             padeCoefficients = PADE_COEFFICIENTS_13;

  1059.             // Calculate scaling factor
  1060.             final double normScale = 5.371920351148152;
  1061.             squaringCount = FastMath.max(0, FastMath.getExponent(l1Norm / normScale));

  1063.             // Scale matrix by power of 2
  1064.             final int finalSquaringCount = squaringCount;
  1065.             scaledMatrix.mapToSelf(x -> FastMath.scalb(x, -finalSquaringCount));
  1066.         }

  1068.         // Calculate U and V using Horner
  1069.         // See Golub, Gene H., and Charles F. van Loan. Matrix Computations. 4th ed.
  1070.         // John Hopkins University Press, 2013.  pages 530/531
  1071.         final RealMatrix scaledMatrix2 = scaledMatrix.multiply(scaledMatrix);
  1072.         final int coeffLength = padeCoefficients.length;

  1074.         // Calculate V
  1075.         RealMatrix padeV = createRealMatrix(dim, dim);
  1076.         for (int i = coeffLength - 1; i > 1; i -= 2) {
  1077.             padeV = scaledMatrix2.multiply(padeV.add(identity.scalarMultiply(padeCoefficients[i])));
  1078.         }
  1079.         padeV = scaledMatrix.multiply(padeV.add(identity.scalarMultiply(padeCoefficients[1])));

  1081.         // Calculate U
  1082.         RealMatrix padeU = createRealMatrix(dim, dim);
  1083.         for (int i = coeffLength - 2; i > 1; i -= 2) {
  1084.             padeU = scaledMatrix2.multiply(padeU.add(identity.scalarMultiply(padeCoefficients[i])));
  1085.         }
  1086.         padeU = padeU.add(identity.scalarMultiply(padeCoefficients[0]));

  1088.         // Calculate pade approximate by solving (U-V) F = (U+V) for F
  1089.         RealMatrix padeNumer = padeU.add(padeV);
  1090.         RealMatrix padeDenom = padeU.subtract(padeV);

  1092.         // Calculate the matrix ratio
  1093.         QRDecomposition decomposition = new QRDecomposition(padeDenom);
  1094.         RealMatrix result = decomposition.getSolver().solve(padeNumer);

  1096.         // Repeated squaring if matrix was scaled
  1097.         for (int i = 0; i < squaringCount; i++) {
  1098.             result = result.multiply(result);
  1099.         }

  1101.         return result;
  1102.     }

  1104.     /** Orthonormalize a list of vectors.
  1105.      * <p>
  1106.      * Orthonormalization is performed by using the Modified Gram-Schmidt process.
  1107.      * </p>
  1108.      * @param independent list of independent vectors
  1109.      * @param threshold projected vectors with a norm less than or equal to this threshold
  1110.      * are considered to have zero norm, hence the vectors they come from are not independent from
  1111.      * previous vectors
  1112.      * @param handler handler for dependent vectors
  1113.      * @return orthonormal basis having the same span as {@code independent}
  1114.      * @since 2.1
  1115.      */
  1116.     public static List<RealVector> orthonormalize(final List<RealVector> independent,
  1117.                                                   final double threshold, final DependentVectorsHandler handler) {

  1119.         // create separate list
  1120.         final List<RealVector> basis = new ArrayList<>(independent);

  1122.         // loop over basis vectors
  1123.         int index = 0;
  1124.         while (index < basis.size()) {

  1126.             // check dependency
  1127.             final RealVector vi = basis.get(index);
  1128.             final double norm = vi.getNorm();
  1129.             if (norm <= threshold) {
  1130.                 // the current vector is dependent from the previous ones
  1131.                 index = handler.manageDependent(index, basis);
  1132.             } else {

  1134.                 // normalize basis vector in place
  1135.                 vi.mapDivideToSelf(vi.getNorm());

  1137.                 // project remaining vectors in place
  1138.                 for (int j = index + 1; j < basis.size(); ++j) {
  1139.                     final RealVector vj  = basis.get(j);
  1140.                     final double     dot = vi.dotProduct(vj);
  1141.                     for (int k = 0; k < vj.getDimension(); ++k) {
  1142.                         vj.setEntry(k, vj.getEntry(k) - dot * vi.getEntry(k));
  1143.                     }
  1144.                 }

  1146.                 ++index;

  1148.             }

  1150.         }

  1152.         return basis;

  1154.     }

  1156.     /** Orthonormalize a list of vectors.
  1157.      * <p>
  1158.      * Orthonormalization is performed by using the Modified Gram-Schmidt process.
  1159.      * </p>
  1160.      * @param <T> type of the field elements
  1161.      * @param independent list of independent vectors
  1162.      * @param threshold projected vectors with a norm less than or equal to this threshold
  1163.      * are considered to have zero norm, hence the vectors they come from are not independent from
  1164.      * previous vectors
  1165.      * @param field type of the files elements
  1166.      * @param handler handler for dependent vectors
  1167.      * @return orthonormal basis having the same span as {@code independent}
  1168.      * @since 2.1
  1169.      */
  1170.     public static <T extends CalculusFieldElement<T>> List<FieldVector<T>> orthonormalize(final Field<T> field,
  1171.                                                                                           final List<FieldVector<T>> independent,
  1172.                                                                                           final T threshold,
  1173.                                                                                           final DependentVectorsHandler handler) {

  1175.         // create separate list
  1176.         final List<FieldVector<T>> basis = new ArrayList<>(independent);

  1178.         // loop over basis vectors
  1179.         int index = 0;
  1180.         while (index < basis.size()) {

  1182.             // check dependency
  1183.             final FieldVector<T> vi = basis.get(index);
  1184.             final T norm = vi.dotProduct(vi).sqrt();
  1185.             if (norm.subtract(threshold).getReal() <= 0) {
  1186.                 // the current vector is dependent from the previous ones
  1187.                 index = handler.manageDependent(field, index, basis);
  1188.             } else {

  1190.                 // normalize basis vector in place
  1191.                 vi.mapDivideToSelf(norm);

  1193.                 // project remaining vectors in place
  1194.                 for (int j = index + 1; j < basis.size(); ++j) {
  1195.                     final FieldVector<T> vj  = basis.get(j);
  1196.                     final T              dot = vi.dotProduct(vj);
  1197.                     for (int k = 0; k < vj.getDimension(); ++k) {
  1198.                         vj.setEntry(k, vj.getEntry(k).subtract(dot.multiply(vi.getEntry(k))));
  1199.                     }
  1200.                 }

  1202.                 ++index;

  1204.             }
  1205.         }

  1207.         return basis;

  1209.     }

  1211. }
Created with JaCoCo 0.8.12.202403310830