matpackII.cc

Go to the documentation of this file.

/* Copyright (C) 2001-2012
   Stefan Buehler <sbuehler@ltu.se>
   Mattias Ekstroem <ekstrom@rss.chalmers.se>

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.

   This program 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 General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */

// #include <vector>
#include "matpackII.h"
#include <Eigen/Core>
#include <algorithm>
#include <cmath>
#include <iostream>  // For debugging.
#include <iterator>
#include <set>

using std::cout;
using std::endl;
using std::setw;
using std::vector;

// Simple member Functions
// ----------------

bool Sparse::empty() const {
  return (matrix.rows() == 0 || matrix.cols() == 0);
}

Index Sparse::nrows() const { return matrix.rows(); }

Index Sparse::ncols() const { return matrix.cols(); }

Index Sparse::nnz() const { return matrix.nonZeros(); }

// Index Operators
// ---------------


Numeric& Sparse::rw(Index r, Index c) {
  // Check if indices are valid:
  assert(0 <= r);
  assert(0 <= c);
  assert(r < nrows());
  assert(c < ncols());

  return matrix.coeffRef((int)r, (int)c);
}


Numeric Sparse::operator()(Index r, Index c) const {
  return matrix.coeff((int)r, (int)c);
}


Numeric Sparse::ro(Index r, Index c) const {
  // Check if indices are valid:
  assert(0 <= r);
  assert(0 <= c);
  assert(r < nrows());
  assert(c < ncols());

  return matrix.coeff((int)r, (int)c);
}

// Constructors
// ------------


Sparse::Sparse() : matrix(0, 0) {
  // Nothing to do here
}


Sparse::Sparse(Index r, Index c) : matrix((int)r, (int)c) {
  // Nothing to do here.
}

Sparse Sparse::diagonal(ConstVectorView v) {
  Index n = v.nelem();
  ArrayOfIndex indices(n);
  for (Index i = 0; i < n; ++i) {
    indices[i] = i;
  }
  Sparse m(n, n);
  m.insert_elements(n, indices, indices, v);
  return m;
}

Vector Sparse::diagonal() const {
  Index m = std::min(nrows(), ncols());
  Vector diag(m);

  for (int i = 0; i < m; i++) {
    diag[i] = matrix.coeff(i, i);
  }
  return diag;
}


Sparse::operator Matrix() const {
  Index m, n;
  m = nrows();
  n = ncols();
  Matrix A(m, n);
  A = 0.0;

  for (int i = 0; i < m; i++) {
    Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(matrix, i);
    for (; it; ++it) {
      A(it.row(), it.col()) = it.value();
    }
  }
  return A;
}


Sparse& Sparse::operator+=(const Sparse& A) {
  matrix += A.matrix;
  return *this;
}


Sparse& Sparse::operator-=(const Sparse& A) {
  matrix -= A.matrix;
  return *this;
}


Sparse& Sparse::operator*=(Numeric x) {
  matrix = matrix * x;
  return *this;
}


Sparse& Sparse::operator/=(Numeric x) {
  matrix = matrix / x;
  return *this;
}


void Sparse::list_elements(Vector& values,
                           ArrayOfIndex& row_indices,
                           ArrayOfIndex& column_indices) const {
  Index m, n;
  m = nrows();
  n = ncols();
  Matrix A(m, n);
  A = 0.0;

  values.resize(nnz());
  row_indices.resize(nnz());
  column_indices.resize(nnz());

  Index j = 0;
  for (int i = 0; i < m; i++) {
    Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(matrix, i);
    for (; it; ++it) {
      values[j] = it.value();
      row_indices[j] = it.row();
      column_indices[j] = it.col();
      j++;
    }
  }
}

void Sparse::split(Index offset, Index nrows_block) {
  Eigen::SparseMatrix<Numeric, Eigen::RowMajor> block_copy(
      matrix.block((int)offset, 0, (int)nrows_block, (int)ncols()));
  matrix.swap(block_copy);
}


void Sparse::insert_row(Index r, Vector v) {
  // Check if the row index and the Vector length are valid
  assert(0 <= r);
  assert(r < nrows());
  assert(v.nelem() == ncols());

  for (int i = 0; i < ncols(); i++) {
    if (v[i] != 0) matrix.coeffRef((int)r, i) = v[i];
  }
}


void Sparse::insert_elements(Index nelems,
                             const ArrayOfIndex& rowind,
                             const ArrayOfIndex& colind,
                             ConstVectorView data) {
  typedef Eigen::Triplet<Numeric> T;
  std::vector<T> tripletList(nelems);

  for (Index i = 0; i < nelems; i++) {
    tripletList[i] = T((int)rowind[i], (int)colind[i], data[i]);
  }

  matrix.setFromTriplets(tripletList.begin(), tripletList.end());
}


void Sparse::resize(Index r, Index c) {
  assert(0 <= r);
  assert(0 <= c);

  matrix.resize((int)r, (int)c);
}


std::ostream& operator<<(std::ostream& os, const Sparse& M) {
  os << M.matrix;
  return os;
}

// General matrix functions


void abs(Sparse& A, const Sparse& B) {
  // Check dimensions
  assert(A.nrows() == B.nrows());
  assert(A.ncols() == B.ncols());

  A.matrix = B.matrix.cwiseAbs();
}


void mult(VectorView y, const Sparse& M, ConstVectorView x) {
  // Check dimensions:
  assert(y.nelem() == M.nrows());
  assert(M.ncols() == x.nelem());

  // Typedefs for Eigen interface
  typedef Eigen::Matrix<Numeric, Eigen::Dynamic, 1, Eigen::ColMajor>
      EigenColumnVector;
  typedef Eigen::Stride<1, Eigen::Dynamic> Stride;
  typedef Eigen::Map<EigenColumnVector, 0, Stride> ColumnMap;

  Numeric* data;
  data = x.mdata + x.mrange.get_start();
  ColumnMap x_map(data, x.nelem(), Stride(1, x.mrange.get_stride()));
  data = y.mdata + y.mrange.get_start();
  ColumnMap y_map(data, y.nelem(), Stride(1, y.mrange.get_stride()));

  y_map = M.matrix * x_map;
}


void transpose_mult(VectorView y, const Sparse& M, ConstVectorView x) {
  // Check dimensions:
  assert(y.nelem() == M.ncols());
  assert(M.nrows() == x.nelem());

  // Typedefs for Eigen interface
  typedef Eigen::Matrix<Numeric, 1, Eigen::Dynamic, Eigen::RowMajor>
      EigenColumnVector;
  typedef Eigen::Stride<1, Eigen::Dynamic> Stride;
  typedef Eigen::Map<EigenColumnVector, 0, Stride> ColumnMap;

  Numeric* data;
  data = x.mdata + x.mrange.get_start();
  ColumnMap x_map(data, x.nelem(), Stride(1, x.mrange.get_stride()));
  data = y.mdata + y.mrange.get_start();
  ColumnMap y_map(data, y.nelem(), Stride(1, y.mrange.get_stride()));

  y_map = x_map * M.matrix;
}


void mult(MatrixView A, const Sparse& B, const ConstMatrixView& C) {
  // Check dimensions:
  assert(A.nrows() == B.nrows());
  assert(A.ncols() == C.ncols());
  assert(B.ncols() == C.nrows());

  // Typedefs for Eigen interface
  typedef Eigen::
      Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
          EigenMatrix;
  typedef Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic> Stride;
  typedef Eigen::Map<EigenMatrix, 0, Stride> MatrixMap;

  Index row_stride, column_stride;
  row_stride = C.mrr.get_stride();
  column_stride = C.mcr.get_stride();

  Numeric* data;
  data = C.mdata + C.mrr.get_start() + C.mcr.get_start();
  Stride c_stride(row_stride, column_stride);
  MatrixMap C_map(data, C.nrows(), C.ncols(), c_stride);

  row_stride = A.mrr.get_stride();
  column_stride = A.mcr.get_stride();
  data = A.mdata + A.mrr.get_start() + A.mcr.get_start();
  Stride a_stride(row_stride, column_stride);
  MatrixMap A_map(data, A.nrows(), A.ncols(), a_stride);

  A_map = B.matrix * C_map;
}


void mult(MatrixView A, const ConstMatrixView& B, const Sparse& C) {
  // Check dimensions:
  assert(A.nrows() == B.nrows());
  assert(A.ncols() == C.ncols());
  assert(B.ncols() == C.nrows());

  // Typedefs for Eigen interface.
  typedef Eigen::
      Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
          EigenMatrix;
  typedef Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic> Stride;
  typedef Eigen::Map<EigenMatrix, 0, Stride> MatrixMap;

  Index row_stride, column_stride;
  row_stride = B.mrr.get_stride();
  column_stride = B.mcr.get_stride();

  Numeric* data;
  data = B.mdata + B.mrr.get_start() + B.mcr.get_start();
  Stride b_stride(row_stride, column_stride);
  MatrixMap B_map(data, B.nrows(), B.ncols(), b_stride);

  row_stride = A.mrr.get_stride();
  column_stride = A.mcr.get_stride();
  data = A.mdata + A.mrr.get_start() + A.mcr.get_start();
  Stride a_stride(row_stride, column_stride);
  MatrixMap A_map(data, A.nrows(), A.ncols(), a_stride);

  A_map = B_map * C.matrix;
}


void transpose(Sparse& A, const Sparse& B) {
  // Check dimensions
  assert(A.nrows() == B.ncols());
  assert(A.ncols() == B.nrows());

  A.matrix = B.matrix.transpose();
}


void mult(Sparse& A, const Sparse& B, const Sparse& C) {
  // Check dimensions and make sure that A is empty
  assert(A.nrows() == B.nrows());
  assert(A.ncols() == C.ncols());
  assert(B.ncols() == C.nrows());

  A.matrix = B.matrix * C.matrix;
}


void add(Sparse& A, const Sparse& B, const Sparse& C) {
  // Check dimensions
  assert(B.ncols() == C.ncols());
  assert(B.nrows() == C.nrows());

  A.resize(B.nrows(), B.ncols());
  A.matrix = B.matrix + C.matrix;
}


void id_mat(Sparse& A) {
  assert(A.ncols() == A.nrows());
  A.matrix.setIdentity();
}


void sub(Sparse& A, const Sparse& B, const Sparse& C) {
  A.resize(B.nrows(), B.ncols());

  // Check dimensions
  assert(B.ncols() == C.ncols());
  assert(B.nrows() == C.nrows());

  A.matrix = B.matrix - C.matrix;
}