/* * 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; } } }