realmatriximpl.java

Apache的common math数学软件包

            out[row] = solution[row][0];
        }
        return out;
    }

    /**
     * Returns a matrix of (column) solution vectors for linear systems with
     * coefficient matrix = this and constant vectors = columns of
     * <code>b</code>.
     *
     * @param b  matrix of constant vectors forming RHS of linear systems to
     * to solve
     * @return matrix of solution vectors
     * @throws IllegalArgumentException if this.rowDimension != row dimension
     * @throws InvalidMatrixException if this matrix is not square or is singular
     */
    public RealMatrix solve(RealMatrix b) throws IllegalArgumentException, InvalidMatrixException  {
        if (b.getRowDimension() != this.getRowDimension()) {
            throw new IllegalArgumentException("Incorrect row dimension");
        }
        if (!this.isSquare()) {
            throw new InvalidMatrixException("coefficient matrix is not square");
        }
        if (this.isSingular()) { // side effect: compute LU decomp
            throw new InvalidMatrixException("Matrix is singular.");
        }

        int nCol = this.getColumnDimension();
        int nColB = b.getColumnDimension();
        int nRowB = b.getRowDimension();

        // Apply permutations to b
        double[][] bp = new double[nRowB][nColB];
        for (int row = 0; row < nRowB; row++) {
            for (int col = 0; col < nColB; col++) {
                bp[row][col] = b.getEntry(permutation[row], col);
            }
        }

        // Solve LY = b
        for (int col = 0; col < nCol; col++) {
            for (int i = col + 1; i < nCol; i++) {
                for (int j = 0; j < nColB; j++) {
                    bp[i][j] -= bp[col][j] * lu[i][col];
                }
            }
        }

        // Solve UX = Y
        for (int col = nCol - 1; col >= 0; col--) {
            for (int j = 0; j < nColB; j++) {
                bp[col][j] /= lu[col][col];
            }
            for (int i = 0; i < col; i++) {
                for (int j = 0; j < nColB; j++) {
                    bp[i][j] -= bp[col][j] * lu[i][col];
                }
            }
        }

        RealMatrixImpl outMat = new RealMatrixImpl(bp);
        return outMat;
    }

    /**
     * Computes a new
     * <a href="http://www.math.gatech.edu/~bourbaki/math2601/Web-notes/2num.pdf">
     * LU decomposition</a> for this matrix, storing the result for use by other methods.
     * <p>
     * <strong>Implementation Note</strong>:<br>
     * Uses <a href="http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/linear.htm">
     * Crout's algorithm</a>, with partial pivoting.</p>
     * <p>
     * <strong>Usage Note</strong>:<br>
     * This method should rarely be invoked directly. Its only use is
     * to force recomputation of the LU decomposition when changes have been
     * made to the underlying data using direct array references. Changes
     * made using setXxx methods will trigger recomputation when needed
     * automatically.</p>
     *
     * @throws InvalidMatrixException if the matrix is non-square or singular.
     */
    public void luDecompose() throws InvalidMatrixException {

        int nRows = this.getRowDimension();
        int nCols = this.getColumnDimension();
        if (nRows != nCols) {
            throw new InvalidMatrixException("LU decomposition requires that the matrix be square.");
        }
        lu = this.getData();

        // Initialize permutation array and parity
        permutation = new int[nRows];
        for (int row = 0; row < nRows; row++) {
            permutation[row] = row;
        }
        parity = 1;

        // Loop over columns
        for (int col = 0; col < nCols; col++) {

            double sum = 0;

            // upper
            for (int row = 0; row < col; row++) {
                sum = lu[row][col];
                for (int i = 0; i < row; i++) {
                    sum -= lu[row][i] * lu[i][col];
                }
                lu[row][col] = sum;
            }

            // lower
            int max = col; // permutation row
            double largest = 0d;
            for (int row = col; row < nRows; row++) {
                sum = lu[row][col];
                for (int i = 0; i < col; i++) {
                    sum -= lu[row][i] * lu[i][col];
                }
                lu[row][col] = sum;

                // maintain best permutation choice
                if (Math.abs(sum) > largest) {
                    largest = Math.abs(sum);
                    max = row;
                }
            }

            // Singularity check
            if (Math.abs(lu[max][col]) < TOO_SMALL) {
                lu = null;
                throw new InvalidMatrixException("matrix is singular");
            }

            // Pivot if necessary
            if (max != col) {
                double tmp = 0;
                for (int i = 0; i < nCols; i++) {
                    tmp = lu[max][i];
                    lu[max][i] = lu[col][i];
                    lu[col][i] = tmp;
                }
                int temp = permutation[max];
                permutation[max] = permutation[col];
                permutation[col] = temp;
                parity = -parity;
            }

            //Divide the lower elements by the "winning" diagonal elt.
            for (int row = col + 1; row < nRows; row++) {
                lu[row][col] /= lu[col][col];
            }
        }
    }

    /**
     * Get a string representation for this matrix.
     * @return a string representation for this matrix
     */
    public String toString() {
        StringBuffer res = new StringBuffer();
        res.append("RealMatrixImpl{");
        if (data != null) {
            for (int i = 0; i < data.length; i++) {
                if (i > 0)
                    res.append(",");
                res.append("{");
                for (int j = 0; j < data[0].length; j++) {
                    if (j > 0)
                        res.append(",");
                    res.append(data[i][j]);
                } 
                res.append("}");
            } 
        }
        res.append("}");
        return res.toString();
    } 
    
    /**
     * Returns true iff <code>object</code> is a 
     * <code>RealMatrixImpl</code> instance with the same dimensions as this
     * and all corresponding matrix entries are equal.  Corresponding entries
     * are compared using {@link java.lang.Double#doubleToLongBits(double)}
     * 
     * @param object the object to test equality against.
     * @return true if object equals this
     */
    public boolean equals(Object object) {
        if (object == this ) {
            return true;
        }
        if (object instanceof RealMatrixImpl == false) {
            return false;
        }
        RealMatrix m = (RealMatrix) object;
        int nRows = getRowDimension();
        int nCols = getColumnDimension();
        if (m.getColumnDimension() != nCols || m.getRowDimension() != nRows) {
            return false;
        }
        for (int row = 0; row < nRows; row++) {
            for (int col = 0; col < nCols; col++) {
                if (Double.doubleToLongBits(data[row][col]) != 
                    Double.doubleToLongBits(m.getEntry(row, col))) {
                    return false;
                }
            }
        }
        return true;
    }
    
    /**
     * Computes a hashcode for the matrix.
     * 
     * @return hashcode for matrix
     */
    public int hashCode() {
        int ret = 7;
        int nRows = getRowDimension();
        int nCols = getColumnDimension();
        ret = ret * 31 + nRows;
        ret = ret * 31 + nCols;
        for (int row = 0; row < nRows; row++) {
           for (int col = 0; col < nCols; col++) {
               ret = ret * 31 + (11 * (row+1) + 17 * (col+1)) * 
                   MathUtils.hash(data[row][col]);
           }
        }
        return ret;
    }

    //------------------------ Protected methods

    /**
     * Returns <code>dimension x dimension</code> identity matrix.
     *
     * @param dimension dimension of identity matrix to generate
     * @return identity matrix
     * @throws IllegalArgumentException  if dimension is not positive
     * @deprecated use {@link MatrixUtils#createRealIdentityMatrix}
     */
    protected RealMatrix getIdentity(int dimension) {
        return MatrixUtils.createRealIdentityMatrix(dimension);
    }

    /**
     *  Returns the LU decomposition as a RealMatrix.
     *  Returns a fresh copy of the cached LU matrix if this has been computed;
     *  otherwise the composition is computed and cached for use by other methods.
     *  Since a copy is returned in either case, changes to the returned matrix do not
     *  affect the LU decomposition property.
     * <p>
     * The matrix returned is a compact representation of the LU decomposition.
     * Elements below the main diagonal correspond to entries of the "L" matrix;
     * elements on and above the main diagonal correspond to entries of the "U"
     * matrix.</p>
     * <p>
     * Example: <pre>
     *
     *     Returned matrix                L                  U
     *         2  3  1                   1  0  0            2  3  1
     *         5  4  6                   5  1  0            0  4  6
     *         1  7  8                   1  7  1            0  0  8
     * </pre>
     *
     * The L and U matrices satisfy the matrix equation LU = permuteRows(this), <br>
     *  where permuteRows reorders the rows of the matrix to follow the order determined
     *  by the <a href=#getPermutation()>permutation</a> property.</p>
     *
     * @return LU decomposition matrix
     * @throws InvalidMatrixException if the matrix is non-square or singular.
     */
    protected RealMatrix getLUMatrix() throws InvalidMatrixException {
        if (lu == null) {
            luDecompose();
        }
        return new RealMatrixImpl(lu);
    }

    /**
     * Returns the permutation associated with the lu decomposition.
     * The entries of the array represent a permutation of the numbers 0, ... , nRows - 1.
     * <p>
     * Example:
     * permutation = [1, 2, 0] means current 2nd row is first, current third row is second
     * and current first row is last.</p>
     * <p>
     * Returns a fresh copy of the array.</p>
     *
     * @return the permutation
     */
    protected int[] getPermutation() {
        int[] out = new int[permutation.length];
        System.arraycopy(permutation, 0, out, 0, permutation.length);
        return out;
    }

    //------------------------ Private methods

    /**
     * Returns a fresh copy of the underlying data array.
     *
     * @return a copy of the underlying data array.
     */
    private double[][] copyOut() {
        int nRows = this.getRowDimension();
        double[][] out = new double[nRows][this.getColumnDimension()];
        // can't copy 2-d array in one shot, otherwise get row references
        for (int i = 0; i < nRows; i++) {
            System.arraycopy(data[i], 0, out[i], 0, data[i].length);
        }
        return out;
    }

    /**
     * Replaces data with a fresh copy of the input array.
     * <p>
     * Verifies that the input array is rectangular and non-empty.</p>
     *
     * @param in data to copy in
     * @throws IllegalArgumentException if input array is empty or not
     *    rectangular
     * @throws NullPointerException if input array is null
     */
    private void copyIn(double[][] in) {
        setSubMatrix(in,0,0);
    }

    /**
     * Tests a given coordinate as being valid or invalid
     *
     * @param row the row index.
     * @param col the column index.
     * @return true if the coordinate is with the current dimensions
     */
    private boolean isValidCoordinate(int row, int col) {
        int nRows = this.getRowDimension();
        int nCols = this.getColumnDimension();

        return !(row < 0 || row > nRows - 1 || col < 0 || col > nCols -1);
    }

}