ARTS 2.5.4 (git: 4c0d3b4d)
matpackI.h
Go to the documentation of this file.
1/* Copyright (C) 2001-2019 Stefan Buehler <sbuehler@ltu.se>
2
3 This program is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2, or (at your option) any
6 later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16 USA. */
17
95#ifndef matpackI_h
96#define matpackI_h
97
98#include <algorithm>
99#pragma GCC diagnostic push
100#pragma GCC diagnostic ignored "-Wshadow"
101#pragma GCC diagnostic ignored "-Wconversion"
102#pragma GCC diagnostic ignored "-Wfloat-conversion"
103#include <Eigen/Dense>
104#pragma GCC diagnostic pop
105
106#include "array.h"
107#include "matpack.h"
108
109// Declare existance of some classes
110class bofstream;
111
112// Declaration of Eigen types
113using StrideType = Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>;
115 Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
116using MatrixViewMap = Eigen::Map<MatrixType, 0, StrideType>;
117using ConstMatrixViewMap = Eigen::Map<const MatrixType, 0, StrideType>;
118using Matrix4x4Type = Eigen::Matrix<Numeric, 4, 4, Eigen::RowMajor>;
119using Matrix4x4ViewMap = Eigen::Map<Matrix4x4Type, 0, StrideType>;
120using ConstMatrix4x4ViewMap = Eigen::Map<const Matrix4x4Type, 0, StrideType>;
121
132class Joker {
133 // Nothing here.
134};
135
136// Declare existence of the global joker object:
137extern const Joker joker;
138
139// Declare the existence of class ConstMatrixView:
140class ConstIterator1D;
141
142// Declare the existence of class VectorView:
143class VectorView;
144
145// Declare the existence of class ConstVectorView:
146class ConstVectorView;
147
148// Declare the existence of class ConstMatrixView:
149class ConstMatrixView;
150
166class Range {
167 public:
168 // Constructors:
169
180 constexpr Range(Index start, Index extent, Index stride = 1) ARTS_NOEXCEPT
181 : mstart(start),
182 mextent(extent),
183 mstride(stride) {
184 // Start must be >= 0:
185 ARTS_ASSERT(0 <= mstart);
186 // Extent also. Although internally negative extent means "to the end",
187 // this can not be created this way, only with the joker. Zero
188 // extent is allowed, though, which corresponds to an empty range.
189 ARTS_ASSERT(0 <= mextent);
190 // Stride can be anything except 0.
191 // SAB 2001-09-21: Allow 0 stride.
192 // ARTS_ASSERT( 0!=mstride);
193 }
194
197 constexpr Range(Index start, Joker, Index stride = 1) ARTS_NOEXCEPT
198 : mstart(start),
199 mextent(-1),
200 mstride(stride) {
201 // Start must be >= 0:
202 ARTS_ASSERT(0 <= mstart);
203 }
204
208 constexpr Range(Joker, Index stride = 1) noexcept
209 : mstart(0),
210 mextent(-1),
211 mstride(stride){
212 // Nothing to do here.
213 };
214
220 constexpr Range(Index max_size, const Range& r) ARTS_NOEXCEPT
221 : mstart(r.mstart),
222 mextent(r.mextent),
223 mstride(r.mstride) {
224 // Start must be >= 0:
225 ARTS_ASSERT(0 <= mstart);
226 // ... and < max_size:
227 ARTS_ASSERT(mstart < max_size);
228
229 // Stride must be != 0:
230 ARTS_ASSERT(0 != mstride);
231
232 // Convert negative extent (joker) to explicit extent
233 if (mextent < 0) {
234 if (0 < mstride)
235 mextent = 1 + (max_size - 1 - mstart) / mstride;
236 else
237 mextent = 1 + (0 - mstart) / mstride;
238 } else {
239#ifndef NDEBUG
240 // Check that extent is ok:
241 Index fin = mstart + (mextent - 1) * mstride;
242 ARTS_ASSERT(0 <= fin);
243 ARTS_ASSERT(fin < max_size);
244#endif
245 }
246 }
247
253 constexpr Range(const Range& p, const Range& n) ARTS_NOEXCEPT
254 : mstart(p.mstart + n.mstart * p.mstride),
255 mextent(n.mextent),
256 mstride(p.mstride* n.mstride) {
257 // We have to juggle here a bit with previous, new, and resulting
258 // quantities. I.e.;
259 // p.mstride: Previous stride
260 // n.mstride: New stride (as specified)
261 // mstride: Resulting stride (old*new)
262
263 // Get the previous final element:
264 Index prev_fin = p.mstart + (p.mextent - 1) * p.mstride;
265
266 // Resulting start must be >= previous start:
267 ARTS_ASSERT(p.mstart <= mstart);
268 // and <= prev_fin, except for Joker:
269 ARTS_ASSERT(mstart <= prev_fin || mextent == -1);
270
271 // Resulting stride must be != 0:
272 ARTS_ASSERT(0 != mstride,
273 "Problem: mstride==",
274 mstride,
275 " for mstart==",
276 mstart,
277 " and mextent==",
278 mextent);
279
280 // Convert negative extent (joker) to explicit extent
281 if (mextent < 0) {
282 if (0 < mstride)
283 mextent = 1 + (prev_fin - mstart) / mstride;
284 else
285 mextent = 1 + (p.mstart - mstart) / mstride;
286 } else {
287#ifndef NDEBUG
288 // Check that extent is ok:
289 Index fin = mstart + (mextent - 1) * mstride;
290 ARTS_ASSERT(p.mstart <= fin);
291 ARTS_ASSERT(fin <= prev_fin);
292#endif
293 }
294 };
295
296 // Friends:
297 friend class ConstVectorView;
298 friend class VectorView;
299 friend class Vector;
300 friend class ConstMatrixView;
301 friend class MatrixView;
302 friend class Matrix;
303 friend class Iterator2D;
304 friend class Iterator3D;
305 friend class Iterator4D;
306 friend class Iterator5D;
307 friend class Iterator6D;
308 friend class Iterator7D;
309 friend class ConstIterator2D;
310 friend class ConstIterator3D;
311 friend class ConstIterator4D;
312 friend class ConstIterator5D;
313 friend class ConstIterator6D;
314 friend class ConstIterator7D;
315 friend class ConstTensor3View;
316 friend class Tensor3View;
317 friend class Tensor3;
318 friend class ConstTensor4View;
319 friend class Tensor4View;
320 friend class Tensor4;
321 friend class ConstTensor5View;
322 friend class Tensor5View;
323 friend class Tensor5;
324 friend class ConstTensor6View;
325 friend class Tensor6View;
326 friend class Tensor6;
327 friend class ConstTensor7View;
328 friend class Tensor7View;
329 friend class Tensor7;
330 friend struct Sparse;
332 friend class ComplexVectorView;
333 friend class ComplexVector;
335 friend class ComplexMatrixView;
336 friend class ComplexMatrix;
337 friend class ComplexIterator2D;
339
340 friend void mult_general(VectorView,
341 const ConstMatrixView&,
343
344 // Member functions:
345
347 [[nodiscard]] constexpr Index get_start() const noexcept { return mstart; }
349 [[nodiscard]] constexpr Index get_extent() const noexcept { return mextent; }
351 [[nodiscard]] constexpr Index get_stride() const noexcept { return mstride; }
352
354 constexpr Range operator()(const Range r) const ARTS_NOEXCEPT {
355 return (r.mextent < 0) ? (mextent < 0) ? Range(mstart + r.mstart * mstride,
356 joker,
357 r.mstride * mstride)
358 : Range(mstart + r.mstart * mstride,
359 mextent,
360 r.mstride * mstride)
361 : Range(mstart + r.mstart * mstride,
362 r.mextent,
363 r.mstride * mstride);
364 }
365
366 constexpr Index operator()(const Index i) const ARTS_NOEXCEPT {
367 return mstart + i * mstride;
368 };
369
370 private:
378};
379
381template <size_t N>
382struct Shape {
383 std::array<Index, N> data;
384 bool operator==(const Shape& other) { return data == other.data; }
385 bool operator!=(const Shape& other) { return data not_eq other.data; }
386 friend std::ostream& operator<<(std::ostream& os, const Shape& shape) {
387 os << shape.data[0];
388 for (size_t i = 1; i < N; i++) os << 'x' << shape.data[i];
389 return os;
390 }
391};
392
396 public:
399 using pointer = Numeric*;
401 using iterator_category = std::random_access_iterator_tag;
402
404 Iterator1D() = default;
405
408 : mx(x),
409 mstride(stride) { /* Nothing to do here. */
410 }
411
412 // Operators:
413
416 mx += mstride;
417 return *this;
418 }
419
421 Numeric& operator*() const ARTS_NOEXCEPT { return *mx; }
422
424 bool operator!=(const Iterator1D& other) const ARTS_NOEXCEPT {
425 if (mx != other.mx) return true;
426 return false;
427 }
428
429 bool operator==(const Iterator1D& other) const ARTS_NOEXCEPT {
430 return !operator!=(other);
431 }
432
434 return (Index)(mx - other.mx) / mstride;
435 }
436
441 friend void copy(ConstIterator1D origin,
442 const ConstIterator1D& end,
443 Iterator1D target);
444
445 private:
447 Numeric* mx{nullptr};
450};
451
455 public:
457 using value_type = const Numeric;
458 using pointer = const Numeric*;
459 using reference = const Numeric&;
460 using iterator_category = std::random_access_iterator_tag;
461
463 ConstIterator1D() = default;
464
467 : mx(x),
468 mstride(stride) { /* Nothing to do here. */
469 }
470
471 // Operators:
474 mx += mstride;
475 return *this;
476 }
477
479 const Numeric& operator*() const ARTS_NOEXCEPT { return *mx; }
480
482 bool operator!=(const ConstIterator1D& other) const ARTS_NOEXCEPT {
483 if (mx != other.mx) return true;
484 return false;
485 }
486
487 bool operator==(const ConstIterator1D& other) const ARTS_NOEXCEPT {
488 return !operator!=(other);
489 }
490
492 return (Index)(mx - other.mx) / mstride;
493 }
494
495 friend void copy(ConstIterator1D origin,
496 const ConstIterator1D& end,
497 Iterator1D target);
498
499 private:
501 const Numeric* mx{nullptr};
504};
505
506// Declare the vector class:
507class Vector;
508
509// Declare the MatrixView class
510class MatrixView;
511
518 public:
519 constexpr ConstVectorView(const ConstVectorView&) = default;
520 constexpr ConstVectorView(ConstVectorView&&) = default;
523
524 // Typedef for compatibility with STL
526
527 // Member functions:
528
530 [[nodiscard]] bool empty() const noexcept { return mrange.mextent == 0; };
531
541 [[nodiscard]] Index nelem() const noexcept { return mrange.mextent; }
542
544 [[nodiscard]] Index size() const noexcept { return mrange.mextent; }
545
547 [[nodiscard]] Shape<1> shape() const { return {nelem()}; }
548
550 [[nodiscard]] Numeric sum() const ARTS_NOEXCEPT;
551
552 // Const index operators:
554 Numeric operator[](Index n) const ARTS_NOEXCEPT { // Check if index is valid:
555 ARTS_ASSERT(0 <= n);
557 return get(n);
558 }
559
561 [[nodiscard]] Numeric get(Index n) const ARTS_NOEXCEPT {
562 return *(mdata + mrange.mstart + n * mrange.mstride);
563 }
564
569
570 friend Numeric operator*(const ConstVectorView& a,
572
573 // Functions returning iterators:
574
576 [[nodiscard]] ConstIterator1D begin() const ARTS_NOEXCEPT;
577
579 [[nodiscard]] ConstIterator1D end() const ARTS_NOEXCEPT;
580
582 operator ConstMatrixView() const;
583
585 virtual ~ConstVectorView() = default;
586
587 // Friends:
588 friend class VectorView;
589 friend class ConstIterator2D;
590 friend class ConstMatrixView;
591 friend class ConstTensor3View;
592 friend class ConstTensor4View;
593 friend class ConstTensor5View;
594 friend class ConstTensor6View;
595 friend class ConstTensor7View;
597 friend int poly_root_solve(Matrix& roots, Vector& coeffs);
598 friend void mult(VectorView, const ConstMatrixView&, const ConstVectorView&);
599 friend void mult(VectorView, const Sparse&, ConstVectorView);
600 friend void transpose_mult(VectorView, const Sparse&, ConstVectorView);
601 friend void mult_general(VectorView,
602 const ConstMatrixView&,
604 friend void lubacksub(VectorView,
607 const ArrayOfIndex&);
609
614
618
619 protected:
620 // Constructors:
621 ConstVectorView() = default;
622
625 ConstVectorView(Numeric* data, const Range& range) ARTS_NOEXCEPT;
626
637 ConstVectorView(Numeric* data, const Range& p, const Range& n) ARTS_NOEXCEPT;
638
639 // Data members:
640 // -------------
644 Numeric* mdata{nullptr};
645};
646
659 public:
660 // Make const methods visible from base class
663 using ConstVectorView::operator[];
665
666 constexpr VectorView(const VectorView&) = default;
667
670 VectorView(const Vector&);
671
674
675 // Typedef for compatibility with STL
677
679 Numeric& operator[](Index n) ARTS_NOEXCEPT { // Check if index is valid:
680 ARTS_ASSERT(0 <= n);
682 return get(n);
683 }
684
687 return *(mdata + mrange.mstart + n * mrange.mstride);
688 }
689
694
695 // Iterators:
696
699
702
703 // Assignment operators:
704
709 VectorView& operator=(const ConstVectorView& v);
710
716 VectorView& operator=(const VectorView& v);
717
720 VectorView& operator=(const Vector& v);
721
722 VectorView& operator=(const Array<Numeric>& v);
723
726 VectorView& operator=(Numeric x);
727
728 // Other operators:
729
731 VectorView operator*=(Numeric x) ARTS_NOEXCEPT;
732
734 VectorView operator/=(Numeric x) ARTS_NOEXCEPT;
735
737 VectorView operator+=(Numeric x) ARTS_NOEXCEPT;
738
740 VectorView operator-=(Numeric x) ARTS_NOEXCEPT;
741
743 VectorView operator*=(const ConstVectorView& x) ARTS_NOEXCEPT;
744
746 VectorView operator/=(const ConstVectorView& x) ARTS_NOEXCEPT;
747
749 VectorView operator+=(const ConstVectorView& x) ARTS_NOEXCEPT;
750
752 VectorView operator-=(const ConstVectorView& x) ARTS_NOEXCEPT;
753
755 operator MatrixView() ARTS_NOEXCEPT;
756
762 [[nodiscard]] const Numeric* get_c_array() const ARTS_NOEXCEPT;
763
770
772 virtual ~VectorView() = default;
773
774 // Friends:
775 friend class ConstIterator2D;
776 friend class Iterator2D;
777 friend class MatrixView;
778 friend class Tensor3View;
779 friend class Tensor4View;
780 friend class Tensor5View;
781 friend class Tensor6View;
782 friend class Tensor7View;
783 friend class ComplexVectorView;
784
788
789 protected:
790 // Constructors:
791 VectorView() = default;
792
795 VectorView(Numeric* data, const Range& range) ARTS_NOEXCEPT;
796
807 VectorView(Numeric* data, const Range& p, const Range& n) ARTS_NOEXCEPT;
808};
809
814 public:
815 // Constructors:
817 Iterator2D() = default;
818
821 : msv(x),
822 mstride(stride) { /* Nothing to do here. */
823 }
824
825 // Operators:
828 msv.mdata += mstride;
829 return *this;
830 }
831
833 bool operator!=(const Iterator2D& other) const ARTS_NOEXCEPT {
834 if (msv.mdata + msv.mrange.mstart !=
835 other.msv.mdata + other.msv.mrange.mstart)
836 return true;
837 return false;
838 }
839
843
846
847 private:
851 Index mstride{0};
852};
853
858 public:
859 // Constructors:
861 ConstIterator2D() = default;
862
865 : msv(x),
866 mstride(stride) { /* Nothing to do here. */
867 }
868
869 // Operators:
872 msv.mdata += mstride;
873 return *this;
874 }
875
877 bool operator!=(const ConstIterator2D& other) const ARTS_NOEXCEPT {
878 if (msv.mdata + msv.mrange.mstart !=
879 other.msv.mdata + other.msv.mrange.mstart)
880 return true;
881 return false;
882 }
883
886 const ConstVectorView* operator->() const ARTS_NOEXCEPT { return &msv; }
887
889 const ConstVectorView& operator*() const ARTS_NOEXCEPT { return msv; }
890
891 private:
895 Index mstride{0};
896};
897
908class Vector : public VectorView {
909 public:
910 // Constructors:
911 Vector() = default;
912
914 Vector(std::initializer_list<Numeric> init);
915
917 Vector(const Eigen::VectorXd& init);
918
920 explicit Vector(Index n);
921
923 Vector(Index n, Numeric fill);
924
932 Vector(Numeric start, Index extent, Numeric stride);
933
939 Vector(const ConstVectorView& v);
940
943 Vector(const Vector& v);
944 Vector(Vector&& v) noexcept : VectorView(std::forward<VectorView>(v)) {
945 v.mdata = nullptr;
946 }
947
949 Vector(const std::vector<Numeric>&);
950
960 ARTS_ASSERT(r0.get_extent() >= 0, "Must have size. Has: ", r0.get_extent());
961 }
962
963 // Assignment operators:
964
989 Vector& operator=(const Vector& v);
990
992 Vector& operator=(Vector&& v) noexcept;
993
995 Vector& operator=(std::initializer_list<Numeric> v);
996
1013 Vector& operator=(const Array<Numeric>& x);
1014
1018
1022 void resize(Index n);
1023
1025 friend void swap(Vector& v1, Vector& v2);
1026
1029 virtual ~Vector();
1030
1031 template <class F>
1032 void transform_elementwise(F&& func) {
1033 std::transform(mdata, mdata + size(), mdata, func);
1034 }
1035};
1036
1037// Declare class Matrix:
1038class Matrix;
1039
1051 public:
1052 constexpr ConstMatrixView(const ConstMatrixView&) = default;
1053 constexpr ConstMatrixView(ConstMatrixView&&) = default;
1056
1057 // Typedef for compatibility with STL
1059
1060 // Member functions:
1061 [[nodiscard]] Index nrows() const noexcept { return mrr.mextent; }
1062 [[nodiscard]] Index ncols() const noexcept { return mcr.mextent; }
1063
1064 // Total size
1065 [[nodiscard]] Index size() const noexcept { return nrows() * ncols(); }
1066 [[nodiscard]] bool empty() const noexcept { return size() == 0; }
1067
1069 [[nodiscard]] Shape<2> shape() const { return {nrows(), ncols()}; }
1070
1071 // Const index operators:
1074 ARTS_NOEXCEPT { // Check if indices are valid:
1075 ARTS_ASSERT(0 <= r);
1076 ARTS_ASSERT(0 <= c);
1077 ARTS_ASSERT(r < mrr.mextent);
1078 ARTS_ASSERT(c < mcr.mextent);
1079
1080 return get(r, c);
1081 }
1082
1084 [[nodiscard]] Numeric get(Index r, Index c) const ARTS_NOEXCEPT {
1085 return *(mdata + mrr.mstart + r * mrr.mstride + mcr.mstart +
1086 c * mcr.mstride);
1087 }
1088
1089 ConstMatrixView operator()(const Range& r,
1090 const Range& c) const ARTS_NOEXCEPT;
1091 ConstVectorView operator()(const Range& r, Index c) const ARTS_NOEXCEPT;
1092 ConstVectorView operator()(Index r, const Range& c) const ARTS_NOEXCEPT;
1093
1094 // Functions returning iterators:
1095 [[nodiscard]] ConstIterator2D begin() const ARTS_NOEXCEPT;
1096 [[nodiscard]] ConstIterator2D end() const ARTS_NOEXCEPT;
1097
1098 // View on diagonal vector
1099 [[nodiscard]] ConstVectorView diagonal() const ARTS_NOEXCEPT;
1100
1102 virtual ~ConstMatrixView() = default;
1103
1104 // Friends:
1105 friend class MatrixView;
1106 friend class ConstIterator3D;
1107 friend class ConstVectorView;
1108 friend class ConstTensor3View;
1109 friend class ConstTensor4View;
1110 friend class ConstTensor5View;
1111 friend class ConstTensor6View;
1112 friend class ConstTensor7View;
1115 friend int poly_root_solve(Matrix& roots, Vector& coeffs);
1116 friend void mult(VectorView, const ConstMatrixView&, const ConstVectorView&);
1117 friend void mult(MatrixView, const ConstMatrixView&, const ConstMatrixView&);
1118 friend void mult(MatrixView, const Sparse&, const ConstMatrixView&);
1119 friend void mult(MatrixView, const ConstMatrixView&, const Sparse&);
1120 friend void mult_general(VectorView,
1121 const ConstMatrixView&,
1123 friend void mult_general(MatrixView,
1124 const ConstMatrixView&,
1126 friend void ludcmp(Matrix&, ArrayOfIndex&, ConstMatrixView);
1127 friend void lubacksub(VectorView,
1130 const ArrayOfIndex&);
1131 friend void inv(MatrixView, ConstMatrixView);
1133
1136
1139
1140 protected:
1141 // Constructors:
1142 ConstMatrixView() = default;
1143 ConstMatrixView(Numeric* data, const Range& r, const Range& c) ARTS_NOEXCEPT;
1145 const Range& pr,
1146 const Range& pc,
1147 const Range& nr,
1148 const Range& nc) ARTS_NOEXCEPT;
1149
1150 // Data members:
1151 // -------------
1153 Range mrr{0, 0, 1};
1155 Range mcr{0, 0, 1};
1157 Numeric* mdata{nullptr};
1158};
1159
1170 public:
1171 // Make const methods visible from base class
1174 using ConstMatrixView::operator();
1176
1177 constexpr MatrixView(const MatrixView&) = default;
1178
1179 // Typedef for compatibility with STL
1181
1182 // Index Operators:
1185 Index c) ARTS_NOEXCEPT { // Check if indices are valid:
1186 ARTS_ASSERT(0 <= r);
1187 ARTS_ASSERT(0 <= c);
1188 ARTS_ASSERT(r < mrr.mextent);
1189 ARTS_ASSERT(c < mcr.mextent);
1190
1191 return get(r, c);
1192 }
1193
1196 return *(mdata + mrr.mstart + r * mrr.mstride + mcr.mstart +
1197 c * mcr.mstride);
1198 }
1199
1200 MatrixView operator()(const Range& r, const Range& c) ARTS_NOEXCEPT;
1201 VectorView operator()(const Range& r, Index c) ARTS_NOEXCEPT;
1202 VectorView operator()(Index r, const Range& c) ARTS_NOEXCEPT;
1203
1204 // Functions returning iterators:
1207
1208 // Assignment operators:
1209 MatrixView& operator=(const ConstMatrixView& m);
1210 MatrixView& operator=(const MatrixView& m);
1211 MatrixView& operator=(const Matrix& m);
1212 MatrixView& operator=(const ConstVectorView& m);
1213 MatrixView& operator=(Numeric x);
1214
1215 // Other operators:
1216 MatrixView& operator*=(Numeric x) ARTS_NOEXCEPT;
1217 MatrixView& operator/=(Numeric x) ARTS_NOEXCEPT;
1218 MatrixView& operator+=(Numeric x) ARTS_NOEXCEPT;
1219 MatrixView& operator-=(Numeric x) ARTS_NOEXCEPT;
1220
1221 MatrixView& operator*=(const ConstMatrixView& x) ARTS_NOEXCEPT;
1222 MatrixView& operator/=(const ConstMatrixView& x) ARTS_NOEXCEPT;
1223 MatrixView& operator+=(const ConstMatrixView& x) ARTS_NOEXCEPT;
1224 MatrixView& operator-=(const ConstMatrixView& x) ARTS_NOEXCEPT;
1225
1226 MatrixView& operator*=(const ConstVectorView& x) ARTS_NOEXCEPT;
1227 MatrixView& operator/=(const ConstVectorView& x) ARTS_NOEXCEPT;
1228 MatrixView& operator+=(const ConstVectorView& x) ARTS_NOEXCEPT;
1229 MatrixView& operator-=(const ConstVectorView& x) ARTS_NOEXCEPT;
1230
1231 // Conversion to a plain C-array
1232 [[nodiscard]] const Numeric* get_c_array() const ARTS_NOEXCEPT;
1234
1236 virtual ~MatrixView() = default;
1237
1238 // Friends:
1239 friend class VectorView;
1240 friend class Iterator3D;
1241 friend class Tensor3View;
1242 friend class Tensor4View;
1243 friend class Tensor5View;
1244 friend class Tensor6View;
1245 friend class Tensor7View;
1246 friend class ComplexMatrixView;
1249 friend void mult(MatrixView, const ConstMatrixView&, const ConstMatrixView&);
1250
1251 protected:
1252 // Constructors:
1253 MatrixView() = default;
1254 MatrixView(Numeric* data, const Range& rr, const Range& cr) ARTS_NOEXCEPT;
1255 MatrixView(Numeric* data,
1256 const Range& pr,
1257 const Range& pc,
1258 const Range& nr,
1259 const Range& nc) ARTS_NOEXCEPT;
1260};
1261
1270class Matrix : public MatrixView {
1271 public:
1272 // Constructors:
1273 Matrix() = default;
1274 Matrix(Index r, Index c);
1275 Matrix(Index r, Index c, Numeric fill);
1276 Matrix(const ConstMatrixView& m);
1277 Matrix(const Matrix& m);
1278 Matrix(Matrix&& v) noexcept : MatrixView(std::forward<MatrixView>(v)) {
1279 v.mdata = nullptr;
1280 }
1281
1292 : MatrixView(d, r0, r1) {
1293 ARTS_ASSERT(r0.get_extent() >= 0, "Must have size. Has: ", r0.get_extent());
1294 ARTS_ASSERT(r1.get_extent() >= 0, "Must have size. Has: ", r1.get_extent());
1295 }
1296
1297 // Assignment operators:
1298 Matrix& operator=(const Matrix& m);
1299 Matrix& operator=(Matrix&& m) noexcept;
1302
1303 // Resize function:
1304 void resize(Index r, Index c);
1305
1306 // Swap function:
1307 friend void swap(Matrix& m1, Matrix& m2);
1308
1309 // Destructor:
1310 virtual ~Matrix();
1311
1313 template <std::size_t dim0>
1315 static_assert(dim0 < 2, "Bad Dimension, Out-of-Bounds");
1316
1317 Range r0(0, dim0 == 0 ? nrows() : ncols());
1318
1319 Vector out(mdata, r0);
1320 ARTS_ASSERT(size() == out.size(),
1321 "Can only reduce size on same size input. "
1322 "The sizes are (input v. output): ",
1323 size(),
1324 " v. ",
1325 out.size());
1326 mdata = nullptr;
1327 return out;
1328 }
1329
1331
1332 template <class F>
1333 void transform_elementwise(F&& func) {
1334 std::transform(mdata, mdata + size(), mdata, func);
1335 }
1336};
1337
1338// Function declarations:
1339// ----------------------
1340
1341void copy(ConstIterator1D origin,
1342 const ConstIterator1D& end,
1343 Iterator1D target);
1344
1346void copy(Numeric x, Iterator1D target, const Iterator1D& end) ARTS_NOEXCEPT;
1347
1348void copy(ConstIterator2D origin,
1349 const ConstIterator2D& end,
1350 Iterator2D target);
1351
1352void copy(Numeric x, Iterator2D target, const Iterator2D& end) ARTS_NOEXCEPT;
1353
1354void mult(VectorView y, const ConstMatrixView& M, const ConstVectorView& x);
1355
1357 const ConstMatrixView& B,
1359
1360void mult(MatrixView A, const ConstMatrixView& B, const ConstMatrixView& C);
1361
1363 const ConstMatrixView& B,
1365
1366void cross3(VectorView c,
1367 const ConstVectorView& a,
1369
1371
1373
1375
1377
1378void transform(VectorView y, double (&my_func)(double), ConstVectorView x);
1379
1380void transform(MatrixView y, double (&my_func)(double), ConstMatrixView x);
1381
1383
1385
1387
1389
1391
1393
1395
1398
1399std::ostream& operator<<(std::ostream& os, const ConstVectorView& v);
1400
1401std::ostream& operator<<(std::ostream& os, const ConstMatrixView& v);
1402
1403std::ostream& operator<<(std::ostream& os, const Range& r);
1404
1405// Converts constant matrix to constant eigen map
1407// Converts constant vector to constant eigen row-view
1409// Converts constant vector to constant eigen row-view
1411// Converts constant vector to constant eigen row-view
1413// Converts constant vector to constant eigen column-view
1415// Converts matrix to eigen map
1417// Converts vector to eigen map row-view
1419// Converts vector to eigen map row-view
1421// Converts vector to eigen map row-view
1423// Converts vector to eigen map column-view
1425
1427// Helper function for debugging
1428#ifndef NDEBUG
1429
1431
1432#endif
1434
1435#endif // matpackI_h
This file contains the definition of Array.
The row iterator class for sub matrices.
The ComplexMatrixView class.
The ComplexMatrix class.
The ComplexVectorView class.
The ComplexVector class.
The const row iterator class for sub matrices.
A constant view of a ComplexMatrix.
A constant view of a ComplexVector.
The constant iterator class for sub vectors.
Definition: matpackI.h:454
bool operator!=(const ConstIterator1D &other) const ARTS_NOEXCEPT
Not equal operator, needed for algorithms like copy.
Definition: matpackI.h:482
const Numeric * pointer
Definition: matpackI.h:458
std::random_access_iterator_tag iterator_category
Definition: matpackI.h:460
const Numeric * mx
Current position.
Definition: matpackI.h:501
Index mstride
Stride.
Definition: matpackI.h:503
ConstIterator1D()=default
Default constructor.
const Numeric value_type
Definition: matpackI.h:457
const Numeric & reference
Definition: matpackI.h:459
bool operator==(const ConstIterator1D &other) const ARTS_NOEXCEPT
Definition: matpackI.h:487
ConstIterator1D & operator++() ARTS_NOEXCEPT
Prefix increment operator.
Definition: matpackI.h:473
Index difference_type
Definition: matpackI.h:456
ConstIterator1D(const Numeric *x, Index stride) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:466
const Numeric & operator*() const ARTS_NOEXCEPT
Dereferencing.
Definition: matpackI.h:479
Index operator-(const ConstIterator1D &other) const ARTS_NOEXCEPT
Definition: matpackI.h:491
friend void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Definition: matpackI.cc:304
The const row iterator class for sub matrices.
Definition: matpackI.h:857
const ConstVectorView & operator*() const ARTS_NOEXCEPT
Dereferencing.
Definition: matpackI.h:889
ConstVectorView msv
Current position.
Definition: matpackI.h:893
ConstIterator2D()=default
Default constructor.
const ConstVectorView * operator->() const ARTS_NOEXCEPT
The -> operator is needed, so that we can write i->begin() to get the 1D iterators.
Definition: matpackI.h:886
bool operator!=(const ConstIterator2D &other) const ARTS_NOEXCEPT
Not equal operator, needed for algorithms like copy.
Definition: matpackI.h:877
ConstIterator2D(const ConstVectorView &x, Index stride) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:864
ConstIterator2D & operator++() ARTS_NOEXCEPT
Prefix increment operator.
Definition: matpackI.h:871
Const version of Iterator3D.
Definition: matpackIII.h:75
Const version of Iterator4D.
Definition: matpackIV.h:76
Const version of Iterator5D.
Definition: matpackV.h:84
Const version of Iterator6D.
Definition: matpackVI.h:86
Const version of Iterator7D.
Definition: matpackVII.h:83
A constant view of a Matrix.
Definition: matpackI.h:1050
constexpr ConstMatrixView(ConstMatrixView &&)=default
ConstMatrixView & operator=(const ConstMatrixView &)=default
bool empty() const noexcept
Definition: matpackI.h:1066
ConstIterator2D begin() const ARTS_NOEXCEPT
Return const iterator to first row.
Definition: matpackI.cc:469
Index size() const noexcept
Definition: matpackI.h:1065
Index nrows() const noexcept
Definition: matpackI.h:1061
constexpr ConstMatrixView(const ConstMatrixView &)=default
Shape< 2 > shape() const
Definition: matpackI.h:1069
ConstIterator2D end() const ARTS_NOEXCEPT
Return const iterator behind last row.
Definition: matpackI.cc:474
Index ncols() const noexcept
Definition: matpackI.h:1062
Numeric operator()(Index r, Index c) const ARTS_NOEXCEPT
Plain const index operator.
Definition: matpackI.h:1073
ConstMatrixView & operator=(ConstMatrixView &&)=default
Numeric get(Index r, Index c) const ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:1084
A constant view of a Tensor3.
Definition: matpackIII.h:130
A constant view of a Tensor4.
Definition: matpackIV.h:131
A constant view of a Tensor5.
Definition: matpackV.h:141
A constant view of a Tensor6.
Definition: matpackVI.h:147
A constant view of a Tensor7.
Definition: matpackVII.h:145
A constant view of a Vector.
Definition: matpackI.h:517
Numeric get(Index n) const ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:561
ConstIterator1D end() const ARTS_NOEXCEPT
Return const iterator behind last element.
Definition: matpackI.cc:67
ConstIterator1D begin() const ARTS_NOEXCEPT
Return const iterator to first element.
Definition: matpackI.cc:63
friend void lubacksub(VectorView, ConstMatrixView, ConstVectorView, const ArrayOfIndex &)
LU backsubstitution.
Definition: lin_alg.cc:91
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:644
friend int poly_root_solve(Matrix &roots, Vector &coeffs)
Definition: poly_roots.cc:90
friend void mult_general(VectorView, const ConstMatrixView &, const ConstVectorView &) ARTS_NOEXCEPT
Matrix Vector multiplication.
Definition: matpackI.cc:1189
Numeric sum() const ARTS_NOEXCEPT
The sum of all elements of a Vector.
Definition: matpackI.cc:48
Numeric operator[](Index n) const ARTS_NOEXCEPT
Plain const index operator.
Definition: matpackI.h:554
constexpr ConstVectorView(ConstVectorView &&)=default
constexpr ConstVectorView(const ConstVectorView &)=default
friend void mult(VectorView, const ConstMatrixView &, const ConstVectorView &)
Matrix-Vector Multiplication.
Definition: matpackI.cc:1131
friend ConstMatrixViewMap MapToEigen(const ConstVectorView &)
Definition: matpackI.cc:1693
Shape< 1 > shape() const
Definition: matpackI.h:547
Range mrange
The range of mdata that is actually used.
Definition: matpackI.h:642
friend void transpose_mult(VectorView, const Sparse &, ConstVectorView)
Sparse matrix - Vector multiplication.
Definition: matpackII.cc:452
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:541
friend void diagonalize(MatrixView, VectorView, VectorView, ConstMatrixView)
Matrix Diagonalization.
Definition: lin_alg.cc:241
Index size() const noexcept
Definition: matpackI.h:544
ConstVectorView & operator=(const ConstVectorView &)=default
friend ConstMatrixViewMap MapToEigenCol(const ConstVectorView &)
Definition: matpackI.cc:1706
ConstVectorView & operator=(ConstVectorView &&)=default
friend Numeric operator*(const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
Scalar product.
Definition: matpackI.cc:1099
bool empty() const noexcept
Returns true if variable size is zero.
Definition: matpackI.h:530
The iterator class for sub vectors.
Definition: matpackI.h:395
Iterator1D & operator++() ARTS_NOEXCEPT
Prefix increment operator.
Definition: matpackI.h:415
bool operator==(const Iterator1D &other) const ARTS_NOEXCEPT
Definition: matpackI.h:429
Index operator-(const Iterator1D &other) const ARTS_NOEXCEPT
Definition: matpackI.h:433
Iterator1D()=default
Default constructor.
Index difference_type
Definition: matpackI.h:397
Numeric & operator*() const ARTS_NOEXCEPT
Dereferencing.
Definition: matpackI.h:421
Index mstride
Stride.
Definition: matpackI.h:449
std::random_access_iterator_tag iterator_category
Definition: matpackI.h:401
Numeric * mx
Current position.
Definition: matpackI.h:447
friend void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Copy data between begin and end to target.
Definition: matpackI.cc:304
Numeric * pointer
Definition: matpackI.h:399
bool operator!=(const Iterator1D &other) const ARTS_NOEXCEPT
Not equal operator, needed for algorithms like copy.
Definition: matpackI.h:424
Numeric & reference
Definition: matpackI.h:400
Numeric value_type
Definition: matpackI.h:398
Iterator1D(Numeric *x, Index stride) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:407
The row iterator class for sub matrices.
Definition: matpackI.h:813
Iterator2D(const VectorView &x, Index stride) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:820
VectorView & operator*() ARTS_NOEXCEPT
Dereferencing.
Definition: matpackI.h:845
bool operator!=(const Iterator2D &other) const ARTS_NOEXCEPT
Not equal operator, needed for algorithms like copy.
Definition: matpackI.h:833
Iterator2D()=default
Default constructor.
VectorView msv
Current position.
Definition: matpackI.h:849
Iterator2D & operator++() ARTS_NOEXCEPT
Prefix increment operator.
Definition: matpackI.h:827
VectorView * operator->() ARTS_NOEXCEPT
The -> operator is needed, so that we can write i->begin() to get the 1D iterators.
Definition: matpackI.h:842
Implementation of Tensors of Rank 3.
Definition: matpackIII.h:34
Implementation of Tensors of Rank 4.
Definition: matpackIV.h:38
Implementation of Tensors of Rank 5.
Definition: matpackV.h:38
The outermost iterator class for rank 6 tensors.
Definition: matpackVI.h:40
Implementation of Tensors of Rank 7.
Definition: matpackVII.h:36
The Joker class.
Definition: matpackI.h:132
The MatrixView class.
Definition: matpackI.h:1169
Numeric & operator()(Index r, Index c) ARTS_NOEXCEPT
Plain index operator.
Definition: matpackI.h:1184
Numeric & get(Index r, Index c) ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:1195
constexpr MatrixView(const MatrixView &)=default
The Matrix class.
Definition: matpackI.h:1270
Matrix(Matrix &&v) noexcept
Definition: matpackI.h:1278
Matrix()=default
Vector reduce_rank() &&ARTS_NOEXCEPT
Definition: matpackI.h:1314
Numeric * get_raw_data() ARTS_NOEXCEPT
Definition: matpackI.h:1330
Matrix(Numeric *d, const Range &r0, const Range &r1) ARTS_NOEXCEPT
Definition: matpackI.h:1291
void transform_elementwise(F &&func)
Definition: matpackI.h:1333
The range class.
Definition: matpackI.h:166
constexpr Range(Index start, Joker, Index stride=1) ARTS_NOEXCEPT
Constructor with joker extent.
Definition: matpackI.h:197
friend void mult_general(VectorView, const ConstMatrixView &, const ConstVectorView &) ARTS_NOEXCEPT
Matrix Vector multiplication.
Definition: matpackI.cc:1189
constexpr Index operator()(const Index i) const ARTS_NOEXCEPT
Definition: matpackI.h:366
constexpr Range(const Range &p, const Range &n) ARTS_NOEXCEPT
Constructor of a new range relative to an old range.
Definition: matpackI.h:253
constexpr Range(Index start, Index extent, Index stride=1) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:180
Index mstart
The start index.
Definition: matpackI.h:368
constexpr Index get_extent() const noexcept
Returns the extent of the range.
Definition: matpackI.h:349
constexpr Range operator()(const Range r) const ARTS_NOEXCEPT
Range of range.
Definition: matpackI.h:354
Index mstride
The stride.
Definition: matpackI.h:377
constexpr Index get_start() const noexcept
Returns the start index of the range.
Definition: matpackI.h:347
Index mextent
The number of elements.
Definition: matpackI.h:375
constexpr Range(Index max_size, const Range &r) ARTS_NOEXCEPT
Constructor which converts a range with joker to an explicit range.
Definition: matpackI.h:220
constexpr Index get_stride() const noexcept
Returns the stride of the range.
Definition: matpackI.h:351
constexpr Range(Joker, Index stride=1) noexcept
Constructor with just a joker.
Definition: matpackI.h:208
The Tensor3View class.
Definition: matpackIII.h:244
The Tensor3 class.
Definition: matpackIII.h:344
The Tensor4View class.
Definition: matpackIV.h:290
The Tensor4 class.
Definition: matpackIV.h:427
The Tensor5View class.
Definition: matpackV.h:341
The Tensor5 class.
Definition: matpackV.h:514
The Tensor6View class.
Definition: matpackVI.h:630
The Tensor6 class.
Definition: matpackVI.h:1097
The Tensor7View class.
Definition: matpackVII.h:1301
The Tensor7 class.
Definition: matpackVII.h:2397
The VectorView class.
Definition: matpackI.h:658
Numeric & get(Index n) ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:686
VectorView()=default
friend class MatrixView
Definition: matpackI.h:777
friend class ConstIterator2D
Definition: matpackI.h:775
constexpr VectorView(const VectorView &)=default
VectorView & operator=(const ConstVectorView &v)
Assignment operator.
Definition: matpackI.cc:153
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array, const-version.
Definition: matpackI.cc:272
friend class Iterator2D
Definition: matpackI.h:776
Numeric & operator[](Index n) ARTS_NOEXCEPT
Plain Index operator.
Definition: matpackI.h:679
Iterator1D begin() ARTS_NOEXCEPT
Return iterator to first element.
Definition: matpackI.cc:144
Iterator1D end() ARTS_NOEXCEPT
Return iterator behind last element.
Definition: matpackI.cc:148
The Vector class.
Definition: matpackI.h:908
Vector(Numeric *d, const Range &r0) ARTS_NOEXCEPT
Definition: matpackI.h:959
Vector(Vector &&v) noexcept
Definition: matpackI.h:944
Vector()=default
void transform_elementwise(F &&func)
Definition: matpackI.h:1032
Binary output file stream class.
Definition: bofstream.h:42
#define ARTS_NOEXCEPT
Definition: debug.h:80
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
void ludcmp(Matrix &LU, ArrayOfIndex &indx, ConstMatrixView A)
LU decomposition.
Definition: lin_alg.cc:56
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:167
Numeric mean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version.
Definition: matpackI.cc:1597
ConstMatrixViewMap MapToEigenRow(const ConstVectorView &A)
Definition: matpackI.cc:1701
Eigen::Map< const Matrix4x4Type, 0, StrideType > ConstMatrix4x4ViewMap
Definition: matpackI.h:120
Numeric vector_angle(ConstVectorView a, ConstVectorView b)
Definition: matpackI.cc:1422
Numeric max(const ConstVectorView &x) ARTS_NOEXCEPT
Max function, vector version.
Definition: matpackI.cc:1529
Eigen::Map< const MatrixType, 0, StrideType > ConstMatrixViewMap
Definition: matpackI.h:117
ConstMatrix4x4ViewMap MapToEigen4x4(const ConstMatrixView &A)
Definition: matpackI.cc:1753
Numeric nanmean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version ignoring nans and infs.
Definition: matpackI.cc:1612
void proj(Vector &c, ConstVectorView a, ConstVectorView b) ARTS_NOEXCEPT
Definition: matpackI.cc:1444
const Joker joker
std::ostream & operator<<(std::ostream &os, const ConstVectorView &v)
Definition: matpackI.cc:108
Eigen::Matrix< Numeric, 4, 4, Eigen::RowMajor > Matrix4x4Type
Definition: matpackI.h:118
Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > StrideType
Definition: matpackI.h:113
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1484
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
cross3
Definition: matpackI.cc:1401
Numeric min(const ConstVectorView &x) ARTS_NOEXCEPT
Min function, vector version.
Definition: matpackI.cc:1563
Eigen::Map< MatrixType, 0, StrideType > MatrixViewMap
Definition: matpackI.h:116
void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Definition: matpackI.cc:304
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixType
Definition: matpackI.h:115
Numeric debug_matrixview_get_elem(MatrixView &mv, Index r, Index c)
Helper function to access matrix elements.
Definition: matpackI.cc:1778
Eigen::Map< Matrix4x4Type, 0, StrideType > Matrix4x4ViewMap
Definition: matpackI.h:119
ConstMatrixView transpose(ConstMatrixView m) ARTS_NOEXCEPT
Const version of transpose.
Definition: matpackI.cc:1454
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
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
Workspace & init(Workspace &ws)
constexpr Numeric r0
The reference radius in IGRF13.
Definition: igrf13.cc:203
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
Definition: zeemandata.h:80
constexpr Rational end(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the largest M for a polarization type of this transition.
Definition: zeemandata.h:113
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:33
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
#define d
#define v
#define a
#define c
#define b
#define N
Definition: rng.cc:164
#define M
Definition: rng.cc:165
Helper shape class.
Definition: matpackI.h:382
std::array< Index, N > data
Definition: matpackI.h:383
bool operator==(const Shape &other)
Definition: matpackI.h:384
bool operator!=(const Shape &other)
Definition: matpackI.h:385
friend std::ostream & operator<<(std::ostream &os, const Shape &shape)
Definition: matpackI.h:386
The Sparse class.
Definition: matpackII.h:67