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 "matpackII.h"
 
#include <algorithm>
#include <cmath>
#include <iostream>  // For debugging.
#include <iterator>
#include <set>
#include <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:
  ARTS_ASSERT(0 <= r);
  ARTS_ASSERT(0 <= c);
  ARTS_ASSERT(r < nrows());
  ARTS_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:
  ARTS_ASSERT(0 <= r);
  ARTS_ASSERT(0 <= c);
  ARTS_ASSERT(r < nrows());
  ARTS_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++;
    }
  }
}
 
Numeric* Sparse::get_element_pointer() { return matrix.valuePtr(); }
int* Sparse::get_column_index_pointer() { return matrix.innerIndexPtr(); }
int* Sparse::get_row_start_pointer() { return matrix.outerIndexPtr(); }
 
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
  ARTS_ASSERT(0 <= r);
  ARTS_ASSERT(r < nrows());
  ARTS_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) {
  ARTS_ASSERT(0 <= r);
  ARTS_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
  ARTS_ASSERT(A.nrows() == B.nrows());
  ARTS_ASSERT(A.ncols() == B.ncols());
 
  A.matrix = B.matrix.cwiseAbs();
}
 
 
void mult(VectorView y, const Sparse& M, ConstVectorView x) {
  // Check dimensions:
  ARTS_ASSERT(y.nelem() == M.nrows());
  ARTS_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:
  ARTS_ASSERT(y.nelem() == M.ncols());
  ARTS_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:
  ARTS_ASSERT(A.nrows() == B.nrows());
  ARTS_ASSERT(A.ncols() == C.ncols());
  ARTS_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:
  ARTS_ASSERT(A.nrows() == B.nrows());
  ARTS_ASSERT(A.ncols() == C.ncols());
  ARTS_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
  ARTS_ASSERT(A.nrows() == B.ncols());
  ARTS_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
  ARTS_ASSERT(A.nrows() == B.nrows());
  ARTS_ASSERT(A.ncols() == C.ncols());
  ARTS_ASSERT(B.ncols() == C.nrows());
 
  A.matrix = B.matrix * C.matrix;
}
 
 
void add(Sparse& A, const Sparse& B, const Sparse& C) {
  // Check dimensions
  ARTS_ASSERT(B.ncols() == C.ncols());
  ARTS_ASSERT(B.nrows() == C.nrows());
 
  A.resize(B.nrows(), B.ncols());
  A.matrix = B.matrix + C.matrix;
}
 
 
void id_mat(Sparse& A) {
  ARTS_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
  ARTS_ASSERT(B.ncols() == C.ncols());
  ARTS_ASSERT(B.nrows() == C.nrows());
 
  A.matrix = B.matrix - C.matrix;
}