abstractmatrix.java

另一个功能更强大的矩阵运算软件开源代码

    public Matrix transABmult(double alpha, Matrix B, Matrix C) {
        return transABmultAdd(alpha, B, C.zero());
    }

    public Matrix transABmultAdd(Matrix B, Matrix C) {
        return transABmultAdd(1, B, C);
    }

    public Matrix transABmultAdd(double alpha, Matrix B, Matrix C) {
        checkTransABmultAdd(B, C);

        if (alpha != 0)
            for (int i = 0; i < numColumns; ++i)
                for (int j = 0; j < C.numColumns(); ++j) {
                    double dot = 0;
                    for (int k = 0; k < numRows; ++k)
                        dot += get(k, i) * B.get(j, k);
                    C.add(i, j, alpha * dot);
                }

        return C;
    }

    /**
     * Checks the arguments to <code>transABmultAdd</code> and
     * <code>transABmultAdd</code>
     */
    protected void checkTransABmultAdd(Matrix B, Matrix C) {
        if (numRows != B.numColumns())
            throw new IndexOutOfBoundsException("A.numRows != B.numColumns ("
                    + numRows + " != " + B.numColumns() + ")");
        if (numColumns != C.numRows())
            throw new IndexOutOfBoundsException("A.numColumns != C.numRows ("
                    + numColumns + " != " + C.numRows() + ")");
        if (B.numRows() != C.numColumns())
            throw new IndexOutOfBoundsException("B.numRows != C.numColumns ("
                    + B.numRows() + " != " + C.numColumns() + ")");
    }

    public Matrix solve(Matrix B, Matrix X) {    	
        throw new UnsupportedOperationException();
    }

    public Matrix transSolve(Matrix B, Matrix X) {
        throw new UnsupportedOperationException();
    }

    /**
     * Checks that a matrix inversion is legal for the given arguments. This is
     * for the square case, not for least-squares problems
     */
    protected void checkSolve(Matrix B, Matrix X) {
        if (!isSquare())
            throw new IndexOutOfBoundsException("!A.isSquare");
        if (B.numRows() != numRows)
            throw new IndexOutOfBoundsException("B.numRows != A.numRows ("
                    + B.numRows() + " != " + numRows + ")");
        if (B.numColumns() != X.numColumns())
            throw new IndexOutOfBoundsException(
                    "B.numColumns != X.numColumns (" + B.numColumns() + " != "
                            + X.numColumns() + ")");
        if (X.numRows() != numColumns)
            throw new IndexOutOfBoundsException("X.numRows != A.numColumns ("
                    + X.numRows() + " != " + numColumns + ")");
    }

    public Matrix rank1(Matrix C) {
        return rank1(1, C);
    }

    public Matrix rank1(double alpha, Matrix C) {
        checkRank1(C);

        if (alpha == 0)
            return this;

        return C.transBmultAdd(alpha, C, this);
    }

    /**
     * Checks that a matrix rank1 update is possible for the given matrix
     */
    protected void checkRank1(Matrix C) {
        if (!isSquare())
            throw new IndexOutOfBoundsException("!A.isSquare");
        if (numRows != C.numRows())
            throw new IndexOutOfBoundsException("A.numRows != C.numRows ("
                    + numRows + " != " + C.numRows() + ")");
    }

    public Matrix transRank1(Matrix C) {
        return transRank1(1, C);
    }

    public Matrix transRank1(double alpha, Matrix C) {
        checkTransRank1(C);

        if (alpha == 0)
            return this;

        return C.transAmultAdd(alpha, C, this);
    }

    /**
     * Checks that a transposed rank1 update is leagal with the given argument
     */
    protected void checkTransRank1(Matrix C) {
        if (!isSquare())
            throw new IndexOutOfBoundsException("!A.isSquare");
        if (numRows != C.numColumns())
            throw new IndexOutOfBoundsException("A.numRows != C.numColumns ("
                    + numRows + " != " + C.numColumns() + ")");
    }

    public Matrix rank2(Matrix B, Matrix C) {
        return rank2(1, B, C);
    }

    public Matrix rank2(double alpha, Matrix B, Matrix C) {
        checkRank2(B, C);

        if (alpha == 0)
            return this;

        return B.transBmultAdd(alpha, C, C.transBmultAdd(alpha, B, this));
    }

    /**
     * Checks that a rank2 update is legal for the given arguments
     */
    protected void checkRank2(Matrix B, Matrix C) {
        if (!isSquare())
            throw new IndexOutOfBoundsException("!A.isSquare");
        if (B.numRows() != C.numRows())
            throw new IndexOutOfBoundsException("B.numRows != C.numRows ("
                    + B.numRows() + " != " + C.numRows() + ")");
        if (B.numColumns() != C.numColumns())
            throw new IndexOutOfBoundsException(
                    "B.numColumns != C.numColumns (" + B.numColumns() + " != "
                            + C.numColumns() + ")");
    }

    public Matrix transRank2(Matrix B, Matrix C) {
        return transRank2(1, B, C);
    }

    public Matrix transRank2(double alpha, Matrix B, Matrix C) {
        checkTransRank2(B, C);

        if (alpha == 0)
            return this;

        return B.transAmultAdd(alpha, C, C.transAmultAdd(alpha, B, this));
    }

    /**
     * Checks that a transposed rank2 update is leagal with the given arguments
     */
    protected void checkTransRank2(Matrix B, Matrix C) {
        if (!isSquare())
            throw new IndexOutOfBoundsException("!A.isSquare");
        if (numRows != B.numColumns())
            throw new IndexOutOfBoundsException("A.numRows != B.numColumns ("
                    + numRows + " != " + B.numColumns() + ")");
        if (B.numRows() != C.numRows())
            throw new IndexOutOfBoundsException("B.numRows != C.numRows ("
                    + B.numRows() + " != " + C.numRows() + ")");
        if (B.numColumns() != C.numColumns())
            throw new IndexOutOfBoundsException(
                    "B.numColumns != C.numColumns (" + B.numColumns() + " != "
                            + C.numColumns() + ")");
    }

    public Matrix scale(double alpha) {
        if (alpha == 1)
            return this;
        else if (alpha == 0)
            return zero();

        for (MatrixEntry e : this)
            e.set(alpha * e.get());

        return this;
    }

    public Matrix set(Matrix B) {
        return set(1, B);
    }

    public Matrix set(double alpha, Matrix B) {
        checkSize(B);

        if (alpha == 0.)
            return zero();
        if (B == this)
            return scale(alpha);

        zero();
        for (MatrixEntry e : B)
            set(e.row(), e.column(), alpha * e.get());

        return this;
    }

    public Matrix add(Matrix B) {
        return add(1, B);
    }

    public Matrix add(double alpha, Matrix B) {
        checkSize(B);

        if (alpha != 0)
            for (MatrixEntry e : B)
                add(e.row(), e.column(), alpha * e.get());

        return this;
    }

    /**
     * Checks that the sizes of this matrix and the given conform
     */
    protected void checkSize(Matrix B) {
        if (numRows != B.numRows())
            throw new IndexOutOfBoundsException("A.numRows != B.numRows ("
                    + numRows + " != " + B.numRows() + ")");
        if (numColumns != B.numColumns())
            throw new IndexOutOfBoundsException(
                    "A.numColumns != B.numColumns (" + numColumns + " != "
                            + B.numColumns() + ")");
    }

    public Matrix transpose() {
        checkTranspose();

        for (int j = 0; j < numColumns; ++j)
            for (int i = j + 1; i < numRows; ++i) {
                double value = get(i, j);
                set(i, j, get(j, i));
                set(j, i, value);
            }

        return this;
    }

    /**
     * Checks that the matrix may be transposed
     */
    protected void checkTranspose() {
        if (!isSquare())
            throw new IndexOutOfBoundsException("!A.isSquare");
    }

    public Matrix transpose(Matrix B) {
        checkTranspose(B);

        if (B == this)
            return transpose();

        B.zero();
        for (MatrixEntry e : this)
            B.set(e.column(), e.row(), e.get());

        return B;
    }

    /**
     * Checks that this matrix can be transposed into the given matrix
     */
    protected void checkTranspose(Matrix B) {
        if (numRows != B.numColumns())
            throw new IndexOutOfBoundsException("A.numRows != B.numColumns ("
                    + numRows + " != " + B.numColumns() + ")");
        if (numColumns != B.numRows())
            throw new IndexOutOfBoundsException("A.numColumns != B.numRows ("
                    + numColumns + " != " + B.numRows() + ")");
    }

    public double norm(Norm type) {
        if (type == Norm.One)
            return norm1();
        else if (type == Norm.Frobenius)
            return normF();
        else if (type == Norm.Infinity)
            return normInf();
        else
            // Maxvalue
            return max();
    }

    /**
     * Computes the 1 norm
     */
    protected double norm1() {
        double[] rowSum = new double[numRows];
        for (MatrixEntry e : this)
            rowSum[e.row()] += Math.abs(e.get());
        return max(rowSum);
    }

    /**
     * Computes the Frobenius norm. This implementation is overflow resistant
     */
    protected double normF() {
        double scale = 0, ssq = 1;
        for (MatrixEntry e : this) {
            double Aval = e.get();
            if (Aval != 0) {
                double absxi = Math.abs(Aval);
                if (scale < absxi) {
                    ssq = 1 + ssq * Math.pow(scale / absxi, 2);
                    scale = absxi;
                } else
                    ssq = ssq + Math.pow(absxi / scale, 2);
            }
        }
        return scale * Math.sqrt(ssq);
    }

    /**
     * Computes the infinity norm
     */
    protected double normInf() {
        double[] columnSum = new double[numColumns];
        for (MatrixEntry e : this)
            columnSum[e.column()] += Math.abs(e.get());
        return max(columnSum);
    }

    /**
     * Returns the largest absolute value
     */
    protected double max() {
        double max = 0;
        for (MatrixEntry e : this)
            max = Math.max(Math.abs(e.get()), max);
        return max;
    }

    /**
     * Returns the largest element of the passed array
     */
    protected double max(double[] x) {
        double max = 0;
        for (int i = 0; i < x.length; ++i)
            max = Math.max(x[i], max);
        return max;
    }

    @Override
    public String toString() {
        // Output into coordinate format. Indices start from 1 instead of 0
        Formatter out = new Formatter();

        out.format("%10d %10d %19d\n", numRows, numColumns, Matrices
                .cardinality(this));

        for (MatrixEntry e : this)
            if (e.get() != 0)
                out.format("%10d %10d % .12e\n", e.row() + 1, e.column() + 1, e
                        .get());

        return out.toString();
    }

    public Iterator<MatrixEntry> iterator() {
        return new RefMatrixIterator();
    }

    /**
     * Iterator over a general matrix. Uses column-major traversal
     */
    class RefMatrixIterator implements Iterator<MatrixEntry> {

        /**
         * Matrix cursor
         */
        int row, column;

        /**
         * Matrix entry
         */
        final RefMatrixEntry entry = new RefMatrixEntry();

        public boolean hasNext() {
            return (row < numRows) && (column < numColumns);
        }

        public MatrixEntry next() {
            entry.update(row, column);

            // Traversal first down the columns, then the rows
            if (row < numRows - 1)
                row++;
            else {
                column++;
                row = 0;
            }

            return entry;
        }

        public void remove() {
            entry.set(0);
        }

    }

    /**
     * Matrix entry backed by the matrix. May be reused for higher performance
     */
    class RefMatrixEntry implements MatrixEntry {

        /**
         * Matrix position
         */
        private int row, column;

        /**
         * Updates the entry
         */
        public void update(int row, int column) {
            this.row = row;
            this.column = column;
        }

        public int row() {
            return row;
        }

        public int column() {
            return column;
        }

        public double get() {
            return AbstractMatrix.this.get(row, column);
        }

        public void set(double value) {
            AbstractMatrix.this.set(row, column, value);
        }
    }

}