abstractmatrix.java

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

/*
 * Copyright (C) 2003-2006 Bjørn-Ove Heimsund
 * 
 * This file is part of MTJ.
 * 
 * This library is free software; you can redistribute it and/or modify it
 * under the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation; either version 2.1 of the License, or (at your
 * option) any later version.
 * 
 * This library is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
 * for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; if not, write to the Free Software Foundation,
 * Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 */

package no.uib.cipr.matrix;

import java.util.Formatter;
import java.util.Iterator;

/**
 * Partial implementation of <code>Matrix</code>. The following methods throw
 * <code>UnsupportedOperationException</code>, and should be overridden by a
 * subclass:
 * <ul>
 * <li><code>get(int,int)</code></li>
 * <li><code>set(int,int,double)</code></li>
 * <li><code>copy</code></li>
 * <li>All the direct solution methods</li>
 * </ul>
 * <p>
 * For the rest of the methods, simple default implementations using a matrix
 * iterator has been provided. There are some kernel operations which the
 * simpler operations forward to, for instance, <code>mult(Matrix,Matrix)</code>
 * forwards to <code>multAdd(double,Matrix,Matrix)</code>. Subclasses can
 * thus focus on overriding the kernel operations, which are:
 * <ul>
 * <li> <code>multAdd(double,Vector,Vector)</code> and
 * <code>transMultAdd(double,Vector,Vector)</code>. </li>
 * <li> <code>rank1(double,Vector,Vector)</code> and
 * <code>rank1(double,Vector,Vector)</code>.</li>
 * <li> <code>multAdd(double,Matrix,Matrix)</code>,
 * <code>transAmultAdd(double,Matrix,Matrix)</code>,
 * <code>transBmultAdd(double,Matrix,Matrix)</code>, and
 * <code>transABmultAdd(double,Matrix,Matrix)</code>. </li>
 * <li> <code>scale(double)</code>. </li>
 * <li> <code>set(double,Matrix)</code> and <code>add(double,Matrix)</code>.
 * </li>
 * <li> <code>transpose</code> and <code>transpose(Matrix)</code>. </li>
 * <li> All the norms.</li>
 * </ul>
 * <p>
 * Finally, a default iterator is provided by this class, which works by calling
 * the <code>get</code> function. A tailored replacement should be used by
 * subclasses.
 * </ul>
 */
public abstract class AbstractMatrix implements Matrix {

    /**
     * Number of rows
     */
    protected int numRows;

    /**
     * Number of columns
     */
    protected int numColumns;

    /**
     * Constructor for AbstractMatrix
     */
    protected AbstractMatrix(int numRows, int numColumns) {
        if (numRows < 0 || numColumns < 0)
            throw new IndexOutOfBoundsException(
                    "Matrix size cannot be negative");
        this.numRows = numRows;
        this.numColumns = numColumns;
    }

    /**
     * Constructor for AbstractMatrix, same size as A. The invoking constructor
     * should set this matrix equal the argument matrix
     */
    protected AbstractMatrix(Matrix A) {
        this(A.numRows(), A.numColumns());
    }

    public int numRows() {
        return numRows;
    }

    public int numColumns() {
        return numColumns;
    }

    public boolean isSquare() {
        return numRows == numColumns;
    }

    public void set(int row, int column, double value) {
        throw new UnsupportedOperationException();
    }

    public void add(int row, int column, double value) {
        set(row, column, value + get(row, column));
    }

    public double get(int row, int column) {
        throw new UnsupportedOperationException();
    }

    /**
     * Checks the passed row and column indices
     */
    protected void check(int row, int column) {
        if (row < 0)
            throw new IndexOutOfBoundsException("row index is negative (" + row
                    + ")");
        if (column < 0)
            throw new IndexOutOfBoundsException("column index is negative ("
                    + column + ")");
        if (row >= numRows)
            throw new IndexOutOfBoundsException("row index >= numRows (" + row
                    + " >= " + numRows + ")");
        if (column >= numColumns)
            throw new IndexOutOfBoundsException("column index >= numColumns ("
                    + column + " >= " + numColumns + ")");
    }

    public Matrix copy() {
        throw new UnsupportedOperationException();
    }

    public Matrix zero() {
        for (MatrixEntry e : this)
            e.set(0);
        return this;
    }

    public Vector mult(Vector x, Vector y) {
        return mult(1, x, y);
    }

    public Vector mult(double alpha, Vector x, Vector y) {
        return multAdd(alpha, x, y.zero());
    }

    public Vector multAdd(Vector x, Vector y) {
        return multAdd(1, x, y);
    }

    public Vector multAdd(double alpha, Vector x, Vector y) {
        checkMultAdd(x, y);

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

        return y;
    }

    /**
     * Checks the arguments to <code>mult</code> and <code>multAdd</code>
     */
    protected void checkMultAdd(Vector x, Vector y) {
        if (numColumns != x.size())
            throw new IndexOutOfBoundsException("A.numColumns != x.size ("
                    + numColumns + " != " + x.size() + ")");
        if (numRows != y.size())
            throw new IndexOutOfBoundsException("A.numRows != y.size ("
                    + numRows + " != " + y.size() + ")");
    }

    public Vector transMult(Vector x, Vector y) {
        return transMult(1, x, y);
    }

    public Vector transMult(double alpha, Vector x, Vector y) {
        return transMultAdd(alpha, x, y.zero());
    }

    public Vector transMultAdd(Vector x, Vector y) {
        return transMultAdd(1, x, y);
    }

    public Vector transMultAdd(double alpha, Vector x, Vector y) {
        checkTransMultAdd(x, y);

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

        return y;
    }

    /**
     * Checks the arguments to <code>transMult</code> and
     * <code>transMultAdd</code>
     */
    protected void checkTransMultAdd(Vector x, Vector y) {
        if (numRows != x.size())
            throw new IndexOutOfBoundsException("A.numRows != x.size ("
                    + numRows + " != " + x.size() + ")");
        if (numColumns != y.size())
            throw new IndexOutOfBoundsException("A.numColumns != y.size ("
                    + numColumns + " != " + y.size() + ")");
    }

    public Vector solve(Vector b, Vector x) {
        throw new UnsupportedOperationException();
    }

    public Vector transSolve(Vector b, Vector 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(Vector b, Vector x) {
        if (!isSquare())
            throw new IndexOutOfBoundsException("!A.isSquare");
        if (numRows != b.size())
            throw new IndexOutOfBoundsException("numRows != b.size (" + numRows
                    + " != " + b.size() + ")");
        if (numColumns != x.size())
            throw new IndexOutOfBoundsException("numColumns != x.size ("
                    + numColumns + " != " + x.size() + ")");
    }

    public Matrix rank1(Vector x) {
        return rank1(1, x);
    }

    public Matrix rank1(double alpha, Vector x) {
        return rank1(1, x, x);
    }

    public Matrix rank1(Vector x, Vector y) {
        return rank1(1, x, y);
    }

    public Matrix rank1(double alpha, Vector x, Vector y) {
        checkRank1(x, y);

        if (alpha == 0)
            return this;

        for (VectorEntry ei : x)
            if (ei.get() != 0)
                for (VectorEntry ej : y)
                    if (ej.get() != 0)
                        add(ei.index(), ej.index(), alpha * ei.get() * ej.get());

        return this;
    }

    /**
     * Checks that a vector rank1 update is possible for the given vectors
     */
    protected void checkRank1(Vector x, Vector y) {
        if (!isSquare())
            throw new IndexOutOfBoundsException("!A.isSquare");
        if (x.size() != numRows)
            throw new IndexOutOfBoundsException("x.size != A.numRows ("
                    + x.size() + " != " + numRows + ")");
        if (y.size() != numColumns)
            throw new IndexOutOfBoundsException("y.size != A.numColumns ("
                    + y.size() + " != " + numColumns + ")");
    }

    public Matrix rank2(Vector x, Vector y) {
        return rank2(1, x, y);
    }

    public Matrix rank2(double alpha, Vector x, Vector y) {
        checkRank2(x, y);

        if (alpha == 0)
            return this;

        for (VectorEntry ei : x)
            for (VectorEntry ej : y) {
                add(ei.index(), ej.index(), alpha * ei.get() * ej.get());
                add(ej.index(), ei.index(), alpha * ei.get() * ej.get());
            }

        return this;
    }

    /**
     * Checks that a vector rank2 update is legal with the given vectors
     */
    protected void checkRank2(Vector x, Vector y) {
        if (!isSquare())
            throw new IndexOutOfBoundsException("!A.isSquare");
        if (x.size() != numRows)
            throw new IndexOutOfBoundsException("x.size != A.numRows ("
                    + x.size() + " != " + numRows + ")");
        if (y.size() != numRows)
            throw new IndexOutOfBoundsException("y.size != A.numRows ("
                    + y.size() + " != " + numRows + ")");
    }

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

    public Matrix mult(double alpha, Matrix B, Matrix C) {
        return multAdd(alpha, B, C.zero());
    }

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

    public Matrix multAdd(double alpha, Matrix B, Matrix C) {
        checkMultAdd(B, C);

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

        return C;
    }

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

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

    public Matrix transAmult(double alpha, Matrix B, Matrix C) {
        return transAmultAdd(alpha, B, C.zero());
    }

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

    public Matrix transAmultAdd(double alpha, Matrix B, Matrix C) {
        checkTransAmultAdd(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(k, j);
                    C.add(i, j, alpha * dot);
                }

        return C;
    }

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

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

    public Matrix transBmult(double alpha, Matrix B, Matrix C) {
        return transBmultAdd(alpha, B, C.zero());
    }

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

    public Matrix transBmultAdd(double alpha, Matrix B, Matrix C) {
        checkTransBmultAdd(B, C);

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

        return C;
    }

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

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