All Downloads are FREE. Search and download functionalities are using the official Maven repository.

no.uib.cipr.matrix.Matrices Maven / Gradle / Ivy

Go to download

Matrix data structures, linear solvers, least squares methods, eigenvalue, and singular value decompositions. For larger random dense matrices (above ~ 350 x 350) matrix-matrix multiplication C = A.B is about 50% faster than MTJ.

The newest version!
/*
 * 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.Arrays;

/**
 * Static utility methods for matrices and vectors.
 */
public final class Matrices {

    private Matrices() {
        // No need to instantiate
    }

    /**
     * max(1, M) provided as a convenience for 'leading dimension'
     * calculations.
     *
     * @param n
     */
    static int ld(int n) {
        return Math.max(1, n);
    }

    /**
     * max(1, max(M, N)) provided as a convenience for 'leading
     * dimension' calculations.
     *
     * @param m
     * @param n
     */
    static int ld(int m, int n) {
        return Math.max(1, Math.max(m, n));
    }

    /**
     * Returns the number of non-zero entries in the given vector
     */
    public static int cardinality(Vector x) {
        int nz = 0;
        for (VectorEntry e : x)
            if (e.get() != 0)
                nz++;
        return nz;
    }

    /**
     * Returns the number of non-zero entries in the given matrix
     */
    public static int cardinality(Matrix A) {
        int nz = 0;
        for (MatrixEntry e : A)
            if (e.get() != 0)
                nz++;
        return nz;
    }

    /**
     * Returns an array of arrays containing a copy of the given matrix. Each
     * array contains one row.
     */
    public static double[][] getArray(Matrix A) {
        double[][] Ad = new double[A.numRows()][A.numColumns()];
        for (MatrixEntry e : A)
            Ad[e.row()][e.column()] = e.get();
        return Ad;
    }

    /**
     * Returns a dense array containing a copy of the given vector
     */
    public static double[] getArray(Vector x) {
        double[] xd = new double[x.size()];
        for (VectorEntry e : x)
            xd[e.index()] = e.get();
        return xd;
    }

    /**
     * Returns the identity matrix of the given size
     *
     * @param size
     *            Number of rows/columns of the matrix
     * @return Matrix of the given size, with ones on the main diagonal
     */
    public static DenseMatrix identity(int size) {
        DenseMatrix A = new DenseMatrix(size, size);
        for (int i = 0; i < size; ++i)
            A.set(i, i, 1);
        return A;
    }

    /**
     * Creates a random vector. Numbers are drawn from a uniform distribution
     * between 0 and 1
     *
     * @param size
     *            Size of the vector
     */
    public static Vector random(int size) {
        return random(new DenseVector(size));
    }

    /**
     * Populates a vector with random numbers drawn from a uniform distribution
     * between 0 and 1
     *
     * @param x
     *            Vector to populate
     */
    public static Vector random(Vector x) {
        for (int i = 0; i < x.size(); ++i)
            x.set(i, Math.random());
        return x;
    }

    /**
     * Creates a random matrix. Numbers are drawn from a uniform distribution
     * between 0 and 1
     *
     * @param numRows
     *            Number of rows
     * @param numColumns
     *            Number of columns
     */
    public static Matrix random(int numRows, int numColumns) {
        return random(new DenseMatrix(numRows, numColumns));
    }

    /**
     * Populates a matrix with random numbers drawn from a uniform distribution
     * between 0 and 1
     *
     * @param A
     *            Matrix to populate
     */
    public static Matrix random(Matrix A) {
        for (int j = 0; j < A.numColumns(); ++j)
            for (int i = 0; i < A.numRows(); ++i)
                A.set(i, j, Math.random());
        return A;
    }

    /**
     * Returns a synchronized vector which wraps the given vector. Only the
     * set(int, double) and add(int, double) methods
     * and their blocked versions are synchronized.
     * 

* Note: Do not use the wrapped vector for any operations besides * matrix assembly, as these operations may be very slow. * * @param x * Vector to be wrapped * @return A thin wrapper around x */ public static Vector synchronizedVector(Vector x) { return new SynchronizedVector(x); } /** * Returns a synchronized matrix which wraps the given matrix. Only the * set(int, int, double) and add(int, int, double) * methods and their blocked versions are synchronized. *

* Note: Do not use the wrapped matrix for any operations besides * matrix assembly, as these operations may be very slow. * * @param A * Matrix to be wrapped * @return A thin wrapper around A */ public static Matrix synchronizedMatrix(Matrix A) { return new SynchronizedMatrix(A); } /** * Returns a synchronized matrix which wraps the given matrix. Only the * set(int, int, double) and add(int, int, double) * methods and their blocked versions are synchronized. *

* The locking provided is finer than the locking of the whole matrix, as * different threads can access different rows simultaneous, while only one * thread can access a given row at a time. Use this for row-major matrices, * not for column-major matrices. *

* Note: Do not use the wrapped matrix for any operations besides * matrix assembly, as these operations may be very slow. * * @param A * Matrix to be wrapped * @return A thin wrapper around A. Individual rows are locked */ public static Matrix synchronizedMatrixByRows(Matrix A) { return new SynchronizedRowMatrix(A); } /** * Returns a synchronized matrix which wraps the given matrix. Only the * set(int, int, double) and add(int, int, double) * methods and their blocked versions are synchronized. *

* The locking provided is finer than the locking of the whole matrix, as * different threads can access different columns simultaneous, while only * one thread can access a given column at a time. Use this for column-major * matrices, not for row-major matrices. *

* Note: Do not use the wrapped matrix for any operations besides * matrix assembly, as these operations may be very slow. * * @param A * Matrix to be wrapped * @return A thin wrapper around A. Individual columns are * locked */ public static Matrix synchronizedMatrixByColumns(Matrix A) { return new SynchronizedColumnMatrix(A); } /** * Returns a view into the given matrix. This view is only for easing some * matrix-assembly cases, not for general use. To extract a more * higher-performing and general matrix, create a copy of the submatrix. The * result is a {@link no.uib.cipr.matrix.DenseMatrix DenseMatrix}. * * @param A * Matrix to create view on * @param row * Rows to access. Must be within the bounds of A * @param column * Columns to access. Must be within the bounds of A * @return Submatrix of A. Changing it will change the backing * matrix */ public static Matrix getSubMatrix(Matrix A, int[] row, int[] column) { return new RefMatrix(A, row, column); } /** * Returns a view into the given vector. This view is only for easing some * vector-assembly cases, not for general use. To extract a more * higher-performing and general vector, create a copy of the subvector. The * result is a {@link no.uib.cipr.matrix.DenseVector DenseVector}. * * @param x * Vector to create view on * @param index * Indices to access. Must be within the bounds of x * @return Submatrix of x. Changing it will change the backing * matrix */ public static Vector getSubVector(Vector x, int[] index) { return new RefVector(x, index); } /** * Matrix backed by another matrix. Used by getSubMatrix */ private static final class RefMatrix extends AbstractMatrix { private Matrix A; private int[] row, column; public RefMatrix(Matrix A, int[] row, int[] column) { super(row.length, column.length); this.A = A; this.row = row; this.column = column; } @Override public void add(int row, int column, double value) { A.add(this.row[row], this.column[column], value); } @Override public DenseMatrix copy() { return new DenseMatrix(this); } @Override public double get(int row, int column) { return A.get(this.row[row], this.column[column]); } @Override public void set(int row, int column, double value) { A.set(this.row[row], this.column[column], value); } } /** * Vector backed by another vector. Used by getSubVector */ @SuppressWarnings("serial") private static final class RefVector extends AbstractVector { private Vector x; private int[] index; public RefVector(Vector x, int[] index) { super(index.length); this.x = x; this.index = index; } @Override public void add(int index, double value) { x.add(this.index[index], value); } @Override public DenseVector copy() { return new DenseVector(this); } @Override public double get(int index) { return x.get(this.index[index]); } @Override public void set(int index, double value) { x.set(this.index[index], value); } } /** * Ensures correctness in the vector assembly. Since it extends the * AbstractVector class, algebraic operations will be slow. It is not * possible to implement Vector and delegate calls to the imbedded vector, * since casting to the imbedded vector is not possible */ @SuppressWarnings("serial") private static final class SynchronizedVector extends AbstractVector { private Vector x; public SynchronizedVector(Vector x) { super(x); this.x = x; } @Override public synchronized void add(int index, double value) { x.add(index, value); } @Override public synchronized void set(int index, double value) { x.set(index, value); } @Override public synchronized double get(int index) { return x.get(index); } @Override public Vector copy() { return Matrices.synchronizedVector(x.copy()); } } /** * Ensures correctness in the matrix assembly. Since it extends the * AbstractMatrix class, algebraic operations will be slow. It is not * possible to implement Matrix and delegate calls to the imbedded matrix, * since casting to the imbedded matrix is not possible */ private static final class SynchronizedMatrix extends AbstractMatrix { private Matrix A; public SynchronizedMatrix(Matrix A) { super(A); this.A = A; } @Override public synchronized void add(int row, int column, double value) { A.add(row, column, value); } @Override public synchronized void set(int row, int column, double value) { A.set(row, column, value); } @Override public synchronized double get(int row, int column) { return A.get(row, column); } @Override public Matrix copy() { return Matrices.synchronizedMatrix(A.copy()); } } /** * Ensures correctness in the matrix assembly. Since it extends the * AbstractMatrix class, algebraic operations will be slow. It is not * possible to implement Matrix and delegate calls to the imbedded matrix, * since casting to the imbedded matrix is not possible *

* Locks individual rows instead of the whole matrix */ private static final class SynchronizedRowMatrix extends AbstractMatrix { private Matrix A; private Object[] lock; public SynchronizedRowMatrix(Matrix A) { super(A); this.A = A; lock = new Object[A.numRows()]; for (int i = 0; i < lock.length; ++i) lock[i] = new Object(); } @Override public void add(int row, int column, double value) { synchronized (lock[row]) { A.add(row, column, value); } } @Override public void set(int row, int column, double value) { synchronized (lock[row]) { A.set(row, column, value); } } @Override public double get(int row, int column) { return A.get(row, column); } @Override public Matrix copy() { return Matrices.synchronizedMatrixByRows(A.copy()); } } /** * Ensures correctness in the matrix assembly. Implements matrix instead of * subclassing the abstract matrix in order to correctly delegate every * method to possibly overridden method in the encapsulated matrix. *

* Locks individual columns instead of the whole matrix */ private static final class SynchronizedColumnMatrix extends AbstractMatrix { private Matrix A; private Object[] lock; public SynchronizedColumnMatrix(Matrix A) { super(A); this.A = A; lock = new Object[A.numColumns()]; for (int i = 0; i < lock.length; ++i) lock[i] = new Object(); } @Override public void add(int row, int column, double value) { synchronized (lock[column]) { A.add(row, column, value); } } @Override public void set(int row, int column, double value) { synchronized (lock[column]) { A.set(row, column, value); } } @Override public double get(int row, int column) { return A.get(row, column); } @Override public Matrix copy() { return Matrices.synchronizedMatrixByColumns(A.copy()); } } /** * Creates a continuous linear index. * * @param from * Start, inclusive * @param to * Stop, exclusive */ public static int[] index(int from, int to) { int length = to - from; if (length < 0) length = 0; int[] index = new int[length]; for (int i = from, j = 0; j < length; ++i, ++j) index[j] = i; return index; } /** * Creates a strided linear index. * * @param from * Start, inclusive * @param stride * stride=1 for continuous. Negative strides are * allowed * @param to * Stop, exclusive */ public static int[] index(int from, int stride, int to) { if (stride == 1) return index(from, to); else if (stride == 0) return new int[0]; if (to <= from && stride > 0) return new int[0]; if (from <= to && stride < 0) return new int[0]; int length = Math.abs((to - from) / stride); if (Math.abs((to - from) % stride) > 0) length++; if (length < 0) length = 0; int[] index = new int[length]; for (int i = from, j = 0; j < length; i += stride, ++j) index[j] = i; return index; } /** * Finds the number of non-zero entries on each row */ public static int[] rowBandwidth(Matrix A) { int[] nz = new int[A.numRows()]; for (MatrixEntry e : A) nz[e.row()]++; return nz; } /** * Finds the number of non-zero entries on each column */ public static int[] columnBandwidth(Matrix A) { int[] nz = new int[A.numColumns()]; for (MatrixEntry e : A) nz[e.column()]++; return nz; } /** * Finds the number of diagonals below the main diagonal. Useful for * converting a general matrix into a banded matrix */ public static int getNumSubDiagonals(Matrix A) { int kl = 0; for (MatrixEntry e : A) kl = Math.max(kl, e.row() - e.column()); return kl; } /** * Finds the number of diagonals above the main diagonal. Useful for * converting a general matrix into a banded matrix */ public static int getNumSuperDiagonals(Matrix A) { int ku = 0; for (MatrixEntry e : A) ku = Math.max(ku, e.column() - e.row()); return ku; } /** * Sets the selected rows of A equal zero, and puts * diagonal on the diagonal of those rows. Useful for enforcing * boundary conditions */ public static void zeroRows(Matrix A, double diagonal, int... row) { // Sort the rows int[] rowS = row.clone(); Arrays.sort(rowS); for (MatrixEntry e : A) { int j = java.util.Arrays.binarySearch(rowS, e.row()); if (j >= 0) { // Found if (e.row() == e.column()) // Diagonal e.set(diagonal); else // Off diagonal e.set(0); } } // Ensure the diagonal is set. This is necessary in case of missing // rows if (diagonal != 0) for (int rowI : row) A.set(rowI, rowI, diagonal); } /** * Sets the selected columns of A equal zero, and puts * diagonal on the diagonal of those columns. Useful for * enforcing boundary conditions */ public static void zeroColumns(Matrix A, double diagonal, int... column) { // Sort the columns int[] columnS = column.clone(); Arrays.sort(columnS); for (MatrixEntry e : A) { int j = java.util.Arrays.binarySearch(columnS, e.column()); if (j >= 0) { // Found if (e.row() == e.column()) // Diagonal e.set(diagonal); else // Off diagonal e.set(0); } } // Ensure the diagonal is set. This is necessary in case of missing // columns if (diagonal != 0) for (int columnI : column) A.set(columnI, columnI, diagonal); } public static DenseVector getColumn(Matrix m, int j) { DenseVector v = new DenseVector(m.numRows()); for (int i = 0; i < v.size(); i++) { v.set(i, m.get(i, j)); } return v; } }





© 2015 - 2024 Weber Informatics LLC | Privacy Policy