ARTS 2.5.0 (git: 9ee3ac6c)
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
52using std::cout;
53using std::endl;
54using std::setw;
55using std::vector;
56
57// Simple member Functions
58// ----------------
59
61bool Sparse::empty() const {
62 return (matrix.rows() == 0 || matrix.cols() == 0);
63}
64
66Index Sparse::nrows() const { return matrix.rows(); }
67
69Index Sparse::ncols() const { return matrix.cols(); }
70
72Index Sparse::nnz() const { return matrix.nonZeros(); }
73
74// Index Operators
75// ---------------
76
78
92 // Check if indices are valid:
93 ARTS_ASSERT(0 <= r);
94 ARTS_ASSERT(0 <= c);
95 ARTS_ASSERT(r < nrows());
96 ARTS_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 ARTS_ASSERT(0 <= r);
133 ARTS_ASSERT(0 <= c);
134 ARTS_ASSERT(r < nrows());
136
137 return matrix.coeff((int)r, (int)c);
138}
139
140// Constructors
141// ------------
149 // Nothing to do here
150}
151
153
157Sparse::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
189Sparse::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
296void 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 ARTS_ASSERT(0 <= r);
317 ARTS_ASSERT(r < nrows());
318 ARTS_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 ARTS_ASSERT(0 <= r);
363 ARTS_ASSERT(0 <= c);
364
365 matrix.resize((int)r, (int)c);
366}
367
369
375std::ostream& operator<<(std::ostream& os, const Sparse& M) {
376 os << M.matrix;
377 return os;
378}
379
380// General matrix functions
381
383
394void abs(Sparse& A, const Sparse& B) {
395 // Check dimensions
396 ARTS_ASSERT(A.nrows() == B.nrows());
397 ARTS_ASSERT(A.ncols() == B.ncols());
398
399 A.matrix = B.matrix.cwiseAbs();
400}
401
403
418 // Check dimensions:
419 ARTS_ASSERT(y.nelem() == M.nrows());
420 ARTS_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 ARTS_ASSERT(y.nelem() == M.ncols());
455 ARTS_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
490void mult(MatrixView A, const Sparse& B, const ConstMatrixView& C) {
491 // Check dimensions:
492 ARTS_ASSERT(A.nrows() == B.nrows());
493 ARTS_ASSERT(A.ncols() == C.ncols());
494 ARTS_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
539void mult(MatrixView A, const ConstMatrixView& B, const Sparse& C) {
540 // Check dimensions:
541 ARTS_ASSERT(A.nrows() == B.nrows());
542 ARTS_ASSERT(A.ncols() == C.ncols());
543 ARTS_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;
569
571
582void transpose(Sparse& A, const Sparse& B) {
583 // Check dimensions
584 ARTS_ASSERT(A.nrows() == B.ncols());
585 ARTS_ASSERT(A.ncols() == B.nrows());
586
587 A.matrix = B.matrix.transpose();
588}
589
591
606void mult(Sparse& A, const Sparse& B, const Sparse& C) {
607 // Check dimensions and make sure that A is empty
608 ARTS_ASSERT(A.nrows() == B.nrows());
609 ARTS_ASSERT(A.ncols() == C.ncols());
610 ARTS_ASSERT(B.ncols() == C.nrows());
611
612 A.matrix = B.matrix * C.matrix;
613}
614
616
630void add(Sparse& A, const Sparse& B, const Sparse& C) {
631 // Check dimensions
632 ARTS_ASSERT(B.ncols() == C.ncols());
633 ARTS_ASSERT(B.nrows() == C.nrows());
634
635 A.resize(B.nrows(), B.ncols());
636 A.matrix = B.matrix + C.matrix;
637}
638
640
647void id_mat(Sparse& A) {
648 ARTS_ASSERT(A.ncols() == A.nrows());
649 A.matrix.setIdentity();
650}
651
653
667void sub(Sparse& A, const Sparse& B, const Sparse& C) {
668 A.resize(B.nrows(), B.ncols());
669
670 // Check dimensions
671 ARTS_ASSERT(B.ncols() == C.ncols());
672 ARTS_ASSERT(B.nrows() == C.nrows());
673
674 A.matrix = B.matrix - C.matrix;
675}
Index nrows
void * data
Index ncols
A constant view of a Matrix.
Definition: matpackI.h:1014
Range mrr
The row range of mdata that is actually used.
Definition: matpackI.h:1109
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1113
Range mcr
The column range of mdata that is actually used.
Definition: matpackI.h:1111
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
A constant view of a Vector.
Definition: matpackI.h:489
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:612
Range mrange
The range of mdata that is actually used.
Definition: matpackI.h:610
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The MatrixView class.
Definition: matpackI.h:1125
The Matrix class.
Definition: matpackI.h:1225
constexpr Index get_stride() const ARTS_NOEXCEPT
Returns the stride of the range.
Definition: matpackI.h:336
constexpr Index get_start() const ARTS_NOEXCEPT
Returns the start index of the range.
Definition: matpackI.h:332
The Sparse class.
Definition: matpackII.h:67
Numeric operator()(Index r, Index c) const
Plain index operator.
Definition: matpackII.cc:111
Eigen::SparseMatrix< Numeric, Eigen::RowMajor > matrix
The actual matrix.
Definition: matpackII.h:148
Index nnz() const
Returns the number of nonzero elements.
Definition: matpackII.cc:72
Sparse()
Default constructor.
Definition: matpackII.cc:148
void split(Index offset, Index nrows)
Reduce matrix to the row range [offset, offset + nrows].
Definition: matpackII.cc:296
Numeric & rw(Index r, Index c)
Read and write index operator.
Definition: matpackII.cc:91
void resize(Index r, Index c)
Resize function.
Definition: matpackII.cc:361
void list_elements(Vector &values, ArrayOfIndex &row_indices, ArrayOfIndex &column_indices) const
List elements in matrix.
Definition: matpackII.cc:270
void insert_row(Index r, Vector v)
Insert row function.
Definition: matpackII.cc:314
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:66
Vector diagonal() const
Diagonal elements as vector.
Definition: matpackII.cc:172
Numeric ro(Index r, Index c) const
Read only index operator.
Definition: matpackII.cc:130
Sparse & operator/=(Numeric x)
Scale matrix by reciprocal.
Definition: matpackII.cc:255
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:69
bool empty() const
Returns true if variable size is zero.
Definition: matpackII.cc:61
Sparse & operator*=(Numeric x)
Scale matrix.
Definition: matpackII.cc:243
Sparse & operator-=(const Sparse &x)
Subtract sparse matrix.
Definition: matpackII.cc:231
Sparse & operator+=(const Sparse &x)
Add sparse matrix.
Definition: matpackII.cc:215
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
The VectorView class.
Definition: matpackI.h:626
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define min(a, b)
void id_mat(Sparse &A)
Sparse identity matrix.
Definition: matpackII.cc:647
std::ostream & operator<<(std::ostream &os, const Sparse &M)
Output operator for Sparse.
Definition: matpackII.cc:375
void transpose(Sparse &A, const Sparse &B)
Transpose of sparse matrix.
Definition: matpackII.cc:582
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
void transpose_mult(VectorView y, const Sparse &M, ConstVectorView x)
Sparse matrix - Vector multiplication.
Definition: matpackII.cc:452
void sub(Sparse &A, const Sparse &B, const Sparse &C)
Sparse - Sparse subtraction.
Definition: matpackII.cc:667
void mult(VectorView y, const Sparse &M, ConstVectorView x)
Sparse matrix - Vector multiplication.
Definition: matpackII.cc:417
void add(Sparse &A, const Sparse &B, const Sparse &C)
Sparse - Sparse addition.
Definition: matpackII.cc:630
Header file for sparse matrices.
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Matrix matrix(const std::array< Numeric, 196 > data)
Create the square matrix from the static data.
Definition: igrf13.cc:206
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:33
#define v
#define c
#define M
Definition: rng.cc:165