CompColMatrix.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.sparse;

import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashSet;
import java.util.Iterator;
import java.util.List;
import java.util.Set;

import no.uib.cipr.matrix.AbstractMatrix;
import no.uib.cipr.matrix.DenseVector;
import no.uib.cipr.matrix.Matrix;
import no.uib.cipr.matrix.MatrixEntry;
import no.uib.cipr.matrix.Vector;
import no.uib.cipr.matrix.io.MatrixInfo;
import no.uib.cipr.matrix.io.MatrixSize;
import no.uib.cipr.matrix.io.MatrixVectorReader;

/**
 * Compressed column storage (CCS) matrix.
 * 
 * Only use this class if the matrix structure (the location of nonzeros) is
 * known and static (does not change).
 */
public class CompColMatrix extends AbstractMatrix {

    /**
     * Matrix data
     */
    double[] data;

    /**
     * Column indices. These are kept sorted within each row.
     */
    int[] columnPointer;

    /**
     * Indices to the start of each row
     */
    int[] rowIndex;

    /**
     * Constructor for CompColMatrix
     * 
     * @param A
     *            Copies from this matrix
     * @param deep
     *            True if the copy is to be deep. If it is a shallow copy,
     *            <code>A</code> must be a <code>CompColMatrix</code>
     */
    public CompColMatrix(Matrix A, boolean deep) {
        super(A);
        construct(A, deep);
    }

    /**
     * Constructor for CompColMatrix
     * 
     * @param A
     *            Copies from this matrix. The copy will be deep
     */
    public CompColMatrix(Matrix A) {
        this(A, true);
    }

    /**
     * Constructor for CompColMatrix
     * 
     * @param r
     *            Reader to get sparse matrix from
     */
    public CompColMatrix(MatrixVectorReader r) throws IOException {
        // Start with a zero-sized matrix
        super(0, 0);

        // Get matrix information. Use the header if present, else just assume
        // that the matrix stores real numbers without any symmetry
        MatrixInfo info = null;
        if (r.hasInfo())
            info = r.readMatrixInfo();
        else
            info = new MatrixInfo(true, MatrixInfo.MatrixField.Real,
                    MatrixInfo.MatrixSymmetry.General);

        // Check that the matrix is in an acceptable format
        if (info.isPattern())
            throw new UnsupportedOperationException(
                    "Pattern matrices are not supported");
        if (info.isDense())
            throw new UnsupportedOperationException(
                    "Dense matrices are not supported");
        if (info.isComplex())
            throw new UnsupportedOperationException(
                    "Complex matrices are not supported");

        // Resize the matrix to correct size
        MatrixSize size = r.readMatrixSize(info);
        numRows = size.numRows();
        numColumns = size.numColumns();

        // Start reading entries
        int numEntries = size.numEntries();
        int[] row = new int[numEntries];
        int[] column = new int[numEntries];
        double[] entry = new double[numEntries];
        r.readCoordinate(row, column, entry);

        // Shift the indices from 1 based to 0 based
        r.add(-1, row);
        r.add(-1, column);

        // Find the number of entries on each column
        List<Set<Integer>> cnz = new ArrayList<Set<Integer>>(numColumns);
        for (int i = 0; i < numColumns; ++i)
            cnz.add(new HashSet<Integer>());

        for (int i = 0; i < numEntries; ++i)
            cnz.get(column[i]).add(row[i]);

        // Allocate some more in case of symmetry
        if (info.isSymmetric() || info.isSkewSymmetric())
            for (int i = 0; i < numEntries; ++i)
                if (row[i] != column[i])
                    cnz.get(row[i]).add(column[i]);

        int[][] nz = new int[numColumns][];
        for (int i = 0; i < numColumns; ++i) {
            nz[i] = new int[cnz.get(i).size()];
            int j = 0;
            for (Integer rowind : cnz.get(i))
                nz[i][j++] = rowind;
        }

        // Create the sparse matrix structure
        construct(nz);

        // Insert the entries
        for (int i = 0; i < size.numEntries(); ++i)
            set(row[i], column[i], entry[i]);

        // Put in extra entries from symmetry or skew symmetry
        if (info.isSymmetric())
            for (int i = 0; i < numEntries; ++i) {
                if (row[i] != column[i])
                    set(column[i], row[i], entry[i]);
            }
        else if (info.isSkewSymmetric())
            for (int i = 0; i < numEntries; ++i) {
                if (row[i] != column[i])
                    set(column[i], row[i], -entry[i]);
            }
    }

    /**
     * Constructor for CompColMatrix
     * 
     * @param numRows
     *            Number of rows
     * @param numColumns
     *            Number of columns
     * @param nz
     *            The nonzero row indices on each column
     */
    public CompColMatrix(int numRows, int numColumns, int[][] nz) {
        super(numRows, numColumns);
        construct(nz);
    }

    private void construct(int[][] nz) {
        int nnz = 0;
        for (int i = 0; i < nz.length; ++i)
            nnz += nz[i].length;

        columnPointer = new int[numColumns + 1];
        rowIndex = new int[nnz];
        data = new double[nnz];

        if (nz.length != numColumns)
            throw new IllegalArgumentException("nz.length != numColumns");

        for (int i = 1; i <= numColumns; ++i) {
            columnPointer[i] = columnPointer[i - 1] + nz[i - 1].length;

            for (int j = columnPointer[i - 1], k = 0; j < columnPointer[i]; ++j, ++k) {
                rowIndex[j] = nz[i - 1][k];
                if (nz[i - 1][k] < 0 || nz[i - 1][k] >= numRows)
                    throw new IllegalArgumentException("nz[" + (i - 1) + "]["
                            + k + "]=" + nz[i - 1][k]
                            + ", which is not a valid row index");
            }

            Arrays.sort(rowIndex, columnPointer[i - 1], columnPointer[i]);
        }
    }

    private void construct(Matrix A, boolean deep) {
        if (deep) {
            if (A instanceof CompColMatrix) {
                CompColMatrix Ac = (CompColMatrix) A;
                data = new double[Ac.data.length];
                columnPointer = new int[Ac.columnPointer.length];
                rowIndex = new int[Ac.rowIndex.length];

                System.arraycopy(Ac.data, 0, data, 0, data.length);
                System.arraycopy(Ac.columnPointer, 0, columnPointer, 0,
                        columnPointer.length);
                System.arraycopy(Ac.rowIndex, 0, rowIndex, 0, rowIndex.length);
            } else {

                List<Set<Integer>> cnz = new ArrayList<Set<Integer>>(numColumns);
                for (int i = 0; i < numColumns; ++i)
                    cnz.add(new HashSet<Integer>());

                for (MatrixEntry e : A)
                    cnz.get(e.column()).add(e.row());

                int[][] nz = new int[numColumns][];
                for (int i = 0; i < numColumns; ++i) {
                    nz[i] = new int[cnz.get(i).size()];
                    int j = 0;
                    for (Integer rowind : cnz.get(i))
                        nz[i][j++] = rowind;
                }

                construct(nz);
                set(A);

            }
        } else {
            CompColMatrix Ac = (CompColMatrix) A;
            columnPointer = Ac.getColumnPointers();
            rowIndex = Ac.getRowIndices();
            data = Ac.getData();
        }
    }

    /**
     * Returns the column pointers
     */
    public int[] getColumnPointers() {
        return columnPointer;
    }

    /**
     * Returns the row indices
     */
    public int[] getRowIndices() {
        return rowIndex;
    }

    /**
     * Returns the internal data storage
     */
    public double[] getData() {
        return data;
    }

    @Override
    public Vector multAdd(double alpha, Vector x, Vector y) {
        if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
            return super.multAdd(alpha, x, y);

        checkMultAdd(x, y);

        double[] xd = ((DenseVector) x).getData(), yd = ((DenseVector) y)
                .getData();

        // y = 1/alpha * y
        y.scale(1 / alpha);

        // y = A*x + y
        for (int i = 0; i < numColumns; ++i)
            for (int j = columnPointer[i]; j < columnPointer[i + 1]; ++j)
                yd[rowIndex[j]] += data[j] * xd[i];

        // y = alpha*y = alpha*A*x + y
        return y.scale(alpha);
    }

    @Override
    public Vector transMult(Vector x, Vector y) {
        if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
            return super.transMult(x, y);

        checkTransMultAdd(x, y);

        double[] xd = ((DenseVector) x).getData();
        double[] yd = ((DenseVector) y).getData();

        for (int i = 0; i < numColumns; ++i) {
            double dot = 0;
            for (int j = columnPointer[i]; j < columnPointer[i + 1]; ++j)
                dot += data[j] * xd[rowIndex[j]];
            yd[i] = dot;
        }

        return y;
    }

    @Override
    public Vector transMultAdd(double alpha, Vector x, Vector y) {
        if (!(x instanceof DenseVector) || !(y instanceof DenseVector))
            return super.transMultAdd(alpha, x, y);

        checkTransMultAdd(x, y);

        double[] xd = ((DenseVector) x).getData();
        double[] yd = ((DenseVector) y).getData();

        for (int i = 0; i < numColumns; ++i) {
            double dot = 0;
            for (int j = columnPointer[i]; j < columnPointer[i + 1]; ++j)
                dot += data[j] * xd[rowIndex[j]];
            yd[i] += alpha * dot;
        }

        return y;
    }

    @Override
    public void set(int row, int column, double value) {
        check(row, column);

        int index = getIndex(row, column);
        data[index] = value;
    }

    @Override
    public void add(int row, int column, double value) {
        check(row, column);

        int index = getIndex(row, column);
        data[index] += value;
    }

    @Override
    public double get(int row, int column) {
        check(row, column);

        int index = no.uib.cipr.matrix.sparse.Arrays.binarySearch(rowIndex,
                row, columnPointer[column], columnPointer[column + 1]);

        if (index >= 0)
            return data[index];
        else
            return 0;
    }

    /**
     * Finds the insertion index
     */
    private int getIndex(int row, int column) {
        int i = no.uib.cipr.matrix.sparse.Arrays.binarySearch(rowIndex, row,
                columnPointer[column], columnPointer[column + 1]);

        if (i >= 0 && rowIndex[i] == row)
            return i;
        else
            throw new IndexOutOfBoundsException("Entry (" + (row + 1) + ", "
                    + (column + 1) + ") is not in the matrix structure");
    }

    @Override
    public CompColMatrix copy() {
        return new CompColMatrix(this);
    }

    @Override
    public Iterator<MatrixEntry> iterator() {
        return new CompColMatrixIterator();
    }

    @Override
    public CompColMatrix zero() {
        Arrays.fill(data, 0);
        return this;
    }

    /**
     * Iterator over a compressed column matrix
     */
    private class CompColMatrixIterator implements Iterator<MatrixEntry> {

        private int column, cursor;

        private CompColMatrixEntry entry = new CompColMatrixEntry();

        public CompColMatrixIterator() {
            // Find first non-empty column
            nextNonEmptyColumn();
        }

        /**
         * Locates the first non-empty column, starting at the current. After
         * the new column has been found, the cursor is also updated
         */
        private void nextNonEmptyColumn() {
            while (column < numColumns()
                    && columnPointer[column] == columnPointer[column + 1])
                column++;
            cursor = columnPointer[column];
        }

        public boolean hasNext() {
            return cursor < data.length;
        }

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

            // Next position is in the same column
            if (cursor < columnPointer[column + 1] - 1)
                cursor++;

            // Next position is at the following (non-empty) column
            else {
                column++;
                nextNonEmptyColumn();
            }

            return entry;
        }

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

    }

    /**
     * Entry of a compressed column matrix
     */
    private class CompColMatrixEntry implements MatrixEntry {

        private int column, cursor;

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

        public int row() {
            return rowIndex[cursor];
        }

        public int column() {
            return column;
        }

        public double get() {
            return data[cursor];
        }

        public void set(double value) {
            data[cursor] = value;
        }
    }
}