ARTS 2.5.9 (git: 825fa5f2)
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 "matpackII.h"
44
45#include <algorithm>
46#include <cmath>
47#include <iostream> // For debugging.
48#include <iterator>
49#include <set>
50#include <vector>
51
52
53// Simple member Functions
54// ----------------
55
57bool Sparse::empty() const {
58 return (matrix.rows() == 0 || matrix.cols() == 0);
59}
60
62Index Sparse::nrows() const { return matrix.rows(); }
63
65Index Sparse::ncols() const { return matrix.cols(); }
66
68Index Sparse::nnz() const { return matrix.nonZeros(); }
69
70// Index Operators
71// ---------------
72
74
88 // Check if indices are valid:
89 ARTS_ASSERT(0 <= r);
90 ARTS_ASSERT(0 <= c);
91 ARTS_ASSERT(r < nrows());
92 ARTS_ASSERT(c < ncols());
93
94 return matrix.coeffRef((int)r, (int)c);
95}
96
98
108 return matrix.coeff((int)r, (int)c);
109}
110
112
127 // Check if indices are valid:
128 ARTS_ASSERT(0 <= r);
129 ARTS_ASSERT(0 <= c);
130 ARTS_ASSERT(r < nrows());
131 ARTS_ASSERT(c < ncols());
132
133 return matrix.coeff((int)r, (int)c);
134}
135
136// Constructors
137// ------------
138
140
144Sparse::Sparse() : matrix(0, 0) {
145 // Nothing to do here
146}
147
153Sparse::Sparse(Index r, Index c) : matrix((int)r, (int)c) {
154 // Nothing to do here.
155}
156
158 Index n = v.nelem();
159 ArrayOfIndex indices(n);
160 for (Index i = 0; i < n; ++i) {
161 indices[i] = i;
162 }
163 Sparse m(n, n);
164 m.insert_elements(n, indices, indices, v);
165 return m;
166}
167
169 Index m = std::min(nrows(), ncols());
170 Vector diag(m);
171
172 for (int i = 0; i < m; i++) {
173 diag[i] = matrix.coeff(i, i);
174 }
175 return diag;
176}
177
179
185Sparse::operator Matrix() const {
186 Index m, n;
187 m = nrows();
188 n = ncols();
189 Matrix A(m, n);
190 A = 0.0;
191
192 for (int i = 0; i < m; i++) {
193 Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(matrix, i);
194 for (; it; ++it) {
195 A(it.row(), it.col()) = it.value();
196 }
197 }
198 return A;
199}
200
202
212 matrix += A.matrix;
213 return *this;
214}
215
217
228 matrix -= A.matrix;
229 return *this;
230}
231
233
240 matrix = matrix * x;
241 return *this;
242}
243
245
252 matrix = matrix / x;
253 return *this;
254}
255
257
267 ArrayOfIndex& row_indices,
268 ArrayOfIndex& column_indices) const {
269 Index m, n;
270 m = nrows();
271 n = ncols();
272 Matrix A(m, n);
273 A = 0.0;
274
275 values.resize(nnz());
276 row_indices.resize(nnz());
277 column_indices.resize(nnz());
278
279 Index j = 0;
280 for (int i = 0; i < m; i++) {
281 Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(matrix, i);
282 for (; it; ++it) {
283 values[j] = it.value();
284 row_indices[j] = it.row();
285 column_indices[j] = it.col();
286 j++;
287 }
288 }
289}
290
292int* Sparse::get_column_index_pointer() { return matrix.innerIndexPtr(); }
293int* Sparse::get_row_start_pointer() { return matrix.outerIndexPtr(); }
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,
340 ConstVectorView data) {
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;
568}
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}
A constant view of a Matrix.
Definition: matpackI.h:1065
Range mrr
The row range of mdata that is actually used.
Definition: matpackI.h:1172
Index nrows() const noexcept
Definition: matpackI.h:1079
Index ncols() const noexcept
Definition: matpackI.h:1080
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1176
Range mcr
The column range of mdata that is actually used.
Definition: matpackI.h:1174
A constant view of a Vector.
Definition: matpackI.h:521
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:660
Range mrange
The range of mdata that is actually used.
Definition: matpackI.h:658
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
The MatrixView class.
Definition: matpackI.h:1188
The Matrix class.
Definition: matpackI.h:1285
constexpr Index get_start() const noexcept
Returns the start index of the range.
Definition: matpackI.h:341
constexpr Index get_stride() const noexcept
Returns the stride of the range.
Definition: matpackI.h:345
The VectorView class.
Definition: matpackI.h:674
The Vector class.
Definition: matpackI.h:910
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
#define ARTS_ASSERT(condition,...)
Definition: debug.h:82
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.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
#define M
Definition: rng.cc:165
The Sparse class.
Definition: matpackII.h:75
Numeric operator()(Index r, Index c) const
Plain index operator.
Definition: matpackII.cc:107
Eigen::SparseMatrix< Numeric, Eigen::RowMajor > matrix
The actual matrix.
Definition: matpackII.h:156
Index nnz() const
Returns the number of nonzero elements.
Definition: matpackII.cc:68
int * get_row_start_pointer()
Definition: matpackII.cc:293
Sparse()
Default constructor.
Definition: matpackII.cc:144
void split(Index offset, Index nrows)
Reduce matrix to the row range [offset, offset + nrows].
Definition: matpackII.cc:296
Numeric * get_element_pointer()
Definition: matpackII.cc:291
Numeric & rw(Index r, Index c)
Read and write index operator.
Definition: matpackII.cc:87
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:266
void insert_row(Index r, Vector v)
Insert row function.
Definition: matpackII.cc:314
int * get_column_index_pointer()
Definition: matpackII.cc:292
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:62
Vector diagonal() const
Diagonal elements as vector.
Definition: matpackII.cc:168
Numeric ro(Index r, Index c) const
Read only index operator.
Definition: matpackII.cc:126
Sparse & operator/=(Numeric x)
Scale matrix by reciprocal.
Definition: matpackII.cc:251
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:65
bool empty() const
Returns true if variable size is zero.
Definition: matpackII.cc:57
Sparse & operator*=(Numeric x)
Scale matrix.
Definition: matpackII.cc:239
Sparse & operator-=(const Sparse &x)
Subtract sparse matrix.
Definition: matpackII.cc:227
Sparse & operator+=(const Sparse &x)
Add sparse matrix.
Definition: matpackII.cc:211
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
#define v
#define c