ARTS  2.3.1285(git:92a29ea9-dirty)
matpackII.cc
Go to the documentation of this file.
1 /* Copyright (C) 2001-2012
2  Stefan Buehler <sbuehler@ltu.se>
3  Mattias Ekstroem <ekstrom@rss.chalmers.se>
4 
5  This program is free software; you can redistribute it and/or modify it
6  under the terms of the GNU General Public License as published by the
7  Free Software Foundation; either version 2, or (at your option) any
8  later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program; if not, write to the Free Software
17  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18  USA. */
19 
43 // #include <vector>
44 #include "matpackII.h"
45 #include <Eigen/Core>
46 #include <algorithm>
47 #include <cmath>
48 #include <iostream> // For debugging.
49 #include <iterator>
50 #include <set>
51 
52 using std::cout;
53 using std::endl;
54 using std::setw;
55 using std::vector;
56 
57 // Simple member Functions
58 // ----------------
59 
61 bool Sparse::empty() const {
62  return (matrix.rows() == 0 || matrix.cols() == 0);
63 }
64 
66 Index Sparse::nrows() const { return matrix.rows(); }
67 
69 Index Sparse::ncols() const { return matrix.cols(); }
70 
72 Index Sparse::nnz() const { return matrix.nonZeros(); }
73 
74 // Index Operators
75 // ---------------
76 
78 
92  // Check if indices are valid:
93  assert(0 <= r);
94  assert(0 <= c);
95  assert(r < nrows());
96  assert(c < ncols());
97 
98  return matrix.coeffRef((int)r, (int)c);
99 }
100 
102 
112  return matrix.coeff((int)r, (int)c);
113 }
114 
116 
131  // Check if indices are valid:
132  assert(0 <= r);
133  assert(0 <= c);
134  assert(r < nrows());
135  assert(c < ncols());
136 
137  return matrix.coeff((int)r, (int)c);
138 }
139 
140 // Constructors
141 // ------------
142 
144 
149  // Nothing to do here
150 }
151 
153 
157 Sparse::Sparse(Index r, Index c) : matrix((int)r, (int)c) {
158  // Nothing to do here.
159 }
160 
162  Index n = v.nelem();
163  ArrayOfIndex indices(n);
164  for (Index i = 0; i < n; ++i) {
165  indices[i] = i;
166  }
167  Sparse m(n, n);
168  m.insert_elements(n, indices, indices, v);
169  return m;
170 }
171 
173  Index m = std::min(nrows(), ncols());
174  Vector diag(m);
175 
176  for (int i = 0; i < m; i++) {
177  diag[i] = matrix.coeff(i, i);
178  }
179  return diag;
180 }
181 
183 
189 Sparse::operator Matrix() const {
190  Index m, n;
191  m = nrows();
192  n = ncols();
193  Matrix A(m, n);
194  A = 0.0;
195 
196  for (int i = 0; i < m; i++) {
197  Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(matrix, i);
198  for (; it; ++it) {
199  A(it.row(), it.col()) = it.value();
200  }
201  }
202  return A;
203 }
204 
206 
216  matrix += A.matrix;
217  return *this;
218 }
219 
221 
232  matrix -= A.matrix;
233  return *this;
234 }
235 
237 
244  matrix = matrix * x;
245  return *this;
246 }
247 
249 
256  matrix = matrix / x;
257  return *this;
258 }
259 
261 
271  ArrayOfIndex& row_indices,
272  ArrayOfIndex& column_indices) const {
273  Index m, n;
274  m = nrows();
275  n = ncols();
276  Matrix A(m, n);
277  A = 0.0;
278 
279  values.resize(nnz());
280  row_indices.resize(nnz());
281  column_indices.resize(nnz());
282 
283  Index j = 0;
284  for (int i = 0; i < m; i++) {
285  Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(matrix, i);
286  for (; it; ++it) {
287  values[j] = it.value();
288  row_indices[j] = it.row();
289  column_indices[j] = it.col();
290  j++;
291  }
292  }
293 }
294 
296 void Sparse::split(Index offset, Index nrows_block) {
297  Eigen::SparseMatrix<Numeric, Eigen::RowMajor> block_copy(
298  matrix.block((int)offset, 0, (int)nrows_block, (int)ncols()));
299  matrix.swap(block_copy);
300 }
301 
303 
315  // Check if the row index and the Vector length are valid
316  assert(0 <= r);
317  assert(r < nrows());
318  assert(v.nelem() == ncols());
319 
320  for (int i = 0; i < ncols(); i++) {
321  if (v[i] != 0) matrix.coeffRef((int)r, i) = v[i];
322  }
323 }
324 
326 
338  const ArrayOfIndex& rowind,
339  const ArrayOfIndex& colind,
341  typedef Eigen::Triplet<Numeric> T;
342  std::vector<T> tripletList(nelems);
343 
344  for (Index i = 0; i < nelems; i++) {
345  tripletList[i] = T((int)rowind[i], (int)colind[i], data[i]);
346  }
347 
348  matrix.setFromTriplets(tripletList.begin(), tripletList.end());
349 }
350 
352 
362  assert(0 <= r);
363  assert(0 <= c);
364 
365  matrix.resize((int)r, (int)c);
366 }
367 
369 
375 std::ostream& operator<<(std::ostream& os, const Sparse& M) {
376  os << M.matrix;
377  return os;
378 }
379 
380 // General matrix functions
381 
383 
394 void abs(Sparse& A, const Sparse& B) {
395  // Check dimensions
396  assert(A.nrows() == B.nrows());
397  assert(A.ncols() == B.ncols());
398 
399  A.matrix = B.matrix.cwiseAbs();
400 }
401 
403 
417 void mult(VectorView y, const Sparse& M, ConstVectorView x) {
418  // Check dimensions:
419  assert(y.nelem() == M.nrows());
420  assert(M.ncols() == x.nelem());
421 
422  // Typedefs for Eigen interface
423  typedef Eigen::Matrix<Numeric, Eigen::Dynamic, 1, Eigen::ColMajor>
424  EigenColumnVector;
425  typedef Eigen::Stride<1, Eigen::Dynamic> Stride;
426  typedef Eigen::Map<EigenColumnVector, 0, Stride> ColumnMap;
427 
428  Numeric* data;
429  data = x.mdata + x.mrange.get_start();
430  ColumnMap x_map(data, x.nelem(), Stride(1, x.mrange.get_stride()));
431  data = y.mdata + y.mrange.get_start();
432  ColumnMap y_map(data, y.nelem(), Stride(1, y.mrange.get_stride()));
433 
434  y_map = M.matrix * x_map;
435 }
436 
438 
453  // Check dimensions:
454  assert(y.nelem() == M.ncols());
455  assert(M.nrows() == x.nelem());
456 
457  // Typedefs for Eigen interface
458  typedef Eigen::Matrix<Numeric, 1, Eigen::Dynamic, Eigen::RowMajor>
459  EigenColumnVector;
460  typedef Eigen::Stride<1, Eigen::Dynamic> Stride;
461  typedef Eigen::Map<EigenColumnVector, 0, Stride> ColumnMap;
462 
463  Numeric* data;
464  data = x.mdata + x.mrange.get_start();
465  ColumnMap x_map(data, x.nelem(), Stride(1, x.mrange.get_stride()));
466  data = y.mdata + y.mrange.get_start();
467  ColumnMap y_map(data, y.nelem(), Stride(1, y.mrange.get_stride()));
468 
469  y_map = x_map * M.matrix;
470 }
471 
473 
490 void mult(MatrixView A, const Sparse& B, const ConstMatrixView& C) {
491  // Check dimensions:
492  assert(A.nrows() == B.nrows());
493  assert(A.ncols() == C.ncols());
494  assert(B.ncols() == C.nrows());
495 
496  // Typedefs for Eigen interface
497  typedef Eigen::
498  Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
499  EigenMatrix;
500  typedef Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic> Stride;
501  typedef Eigen::Map<EigenMatrix, 0, Stride> MatrixMap;
502 
503  Index row_stride, column_stride;
504  row_stride = C.mrr.get_stride();
505  column_stride = C.mcr.get_stride();
506 
507  Numeric* data;
508  data = C.mdata + C.mrr.get_start() + C.mcr.get_start();
509  Stride c_stride(row_stride, column_stride);
510  MatrixMap C_map(data, C.nrows(), C.ncols(), c_stride);
511 
512  row_stride = A.mrr.get_stride();
513  column_stride = A.mcr.get_stride();
514  data = A.mdata + A.mrr.get_start() + A.mcr.get_start();
515  Stride a_stride(row_stride, column_stride);
516  MatrixMap A_map(data, A.nrows(), A.ncols(), a_stride);
517 
518  A_map = B.matrix * C_map;
519 }
520 
522 
539 void mult(MatrixView A, const ConstMatrixView& B, const Sparse& C) {
540  // Check dimensions:
541  assert(A.nrows() == B.nrows());
542  assert(A.ncols() == C.ncols());
543  assert(B.ncols() == C.nrows());
544 
545  // Typedefs for Eigen interface.
546  typedef Eigen::
547  Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
548  EigenMatrix;
549  typedef Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic> Stride;
550  typedef Eigen::Map<EigenMatrix, 0, Stride> MatrixMap;
551 
552  Index row_stride, column_stride;
553  row_stride = B.mrr.get_stride();
554  column_stride = B.mcr.get_stride();
555 
556  Numeric* data;
557  data = B.mdata + B.mrr.get_start() + B.mcr.get_start();
558  Stride b_stride(row_stride, column_stride);
559  MatrixMap B_map(data, B.nrows(), B.ncols(), b_stride);
560 
561  row_stride = A.mrr.get_stride();
562  column_stride = A.mcr.get_stride();
563  data = A.mdata + A.mrr.get_start() + A.mcr.get_start();
564  Stride a_stride(row_stride, column_stride);
565  MatrixMap A_map(data, A.nrows(), A.ncols(), a_stride);
566 
567  A_map = B_map * C.matrix;
568 }
569 
571 
582 void transpose(Sparse& A, const Sparse& B) {
583  // Check dimensions
584  assert(A.nrows() == B.ncols());
585  assert(A.ncols() == B.nrows());
586 
587  A.matrix = B.matrix.transpose();
588 }
589 
591 
606 void mult(Sparse& A, const Sparse& B, const Sparse& C) {
607  // Check dimensions and make sure that A is empty
608  assert(A.nrows() == B.nrows());
609  assert(A.ncols() == C.ncols());
610  assert(B.ncols() == C.nrows());
611 
612  A.matrix = B.matrix * C.matrix;
613 }
614 
616 
630 void add(Sparse& A, const Sparse& B, const Sparse& C) {
631  // Check dimensions
632  assert(B.ncols() == C.ncols());
633  assert(B.nrows() == C.nrows());
634 
635  A.resize(B.nrows(), B.ncols());
636  A.matrix = B.matrix + C.matrix;
637 }
638 
640 
647 void id_mat(Sparse& A) {
648  assert(A.ncols() == A.nrows());
649  A.matrix.setIdentity();
650 }
651 
653 
667 void sub(Sparse& A, const Sparse& B, const Sparse& C) {
668  A.resize(B.nrows(), B.ncols());
669 
670  // Check dimensions
671  assert(B.ncols() == C.ncols());
672  assert(B.nrows() == C.nrows());
673 
674  A.matrix = B.matrix - C.matrix;
675 }
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
friend void id_mat(Sparse &A)
Sparse identity matrix.
Definition: matpackII.cc:647
The VectorView class.
Definition: matpackI.h:610
Index nnz() const
Returns the number of nonzero elements.
Definition: matpackII.cc:72
Numeric operator()(Index r, Index c) const
Plain index operator.
Definition: matpackII.cc:111
Vector diagonal() const
Diagonal elements as vector.
Definition: matpackII.cc:172
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1081
The Vector class.
Definition: matpackI.h:860
The MatrixView class.
Definition: matpackI.h:1093
The Sparse class.
Definition: matpackII.h:60
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:69
Numeric & rw(Index r, Index c)
Read and write index operator.
Definition: matpackII.cc:91
bool empty() const
Returns true if variable size is zero.
Definition: matpackII.cc:61
Eigen::SparseMatrix< Numeric, Eigen::RowMajor > matrix
The actual matrix.
Definition: matpackII.h:141
#define min(a, b)
constexpr Index get_start() const
Returns the start index of the range.
Definition: matpackI.h:327
G0 G2 FVC Y DV Numeric Numeric Numeric Zeeman LowerQuantumNumbers void * data
friend void transpose_mult(VectorView y, const Sparse &M, ConstVectorView x)
Sparse matrix - Vector multiplication.
Definition: matpackII.cc:452
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
Sparse()
Default constructor.
Definition: matpackII.cc:148
Range mrr
The row range of mdata that is actually used.
Definition: matpackI.h:1077
Header file for sparse matrices.
friend void transpose(Sparse &A, const Sparse &B)
Transpose of sparse matrix.
Definition: matpackII.cc:582
void split(Index offset, Index nrows)
Reduce matrix to the row range [offset, offset + nrows].
Definition: matpackII.cc:296
friend void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:66
Sparse & operator/=(Numeric x)
Scale matrix by reciprocal.
Definition: matpackII.cc:255
friend void mult(VectorView y, const Sparse &M, ConstVectorView x)
Sparse matrix - Vector multiplication.
Definition: matpackII.cc:417
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Sparse & operator-=(const Sparse &x)
Subtract sparse matrix.
Definition: matpackII.cc:231
The Matrix class.
Definition: matpackI.h:1193
friend void sub(Sparse &A, const Sparse &B, const Sparse &C)
Sparse - Sparse subtraction.
Definition: matpackII.cc:667
constexpr Index get_stride() const
Returns the stride of the range.
Definition: matpackI.h:331
Range mcr
The column range of mdata that is actually used.
Definition: matpackI.h:1079
Sparse & operator+=(const Sparse &x)
Add sparse matrix.
Definition: matpackII.cc:215
friend std::ostream & operator<<(std::ostream &os, const Sparse &v)
Output operator for Sparse.
Definition: matpackII.cc:375
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
Numeric ro(Index r, Index c) const
Read only index operator.
Definition: matpackII.cc:130
A constant view of a Vector.
Definition: matpackI.h:476
void list_elements(Vector &values, ArrayOfIndex &row_indices, ArrayOfIndex &column_indices) const
List elements in matrix.
Definition: matpackII.cc:270
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:596
A constant view of a Matrix.
Definition: matpackI.h:982
void resize(Index r, Index c)
Resize function.
Definition: matpackII.cc:361
Range mrange
The range of mdata that is actually used.
Definition: matpackI.h:594
friend void add(Sparse &A, const Sparse &B, const Sparse &C)
Sparse - Sparse addition.
Definition: matpackII.cc:630
void insert_row(Index r, Vector v)
Insert row function.
Definition: matpackII.cc:314
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:34
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
void insert_elements(Index nnz, const ArrayOfIndex &rowind, const ArrayOfIndex &colind, ConstVectorView data)
Insert vector of elements with given row and column indices.
Definition: matpackII.cc:337
Sparse & operator*=(Numeric x)
Scale matrix.
Definition: matpackII.cc:243