ARTS 2.5.4 (git: 31ce4f0e)
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
100#include "array.h"
101#include "matpack.h"
102#include "matpack_concepts.h"
103
104// Declare existance of some classes
105class bofstream;
106class Vector;
107class Matrix;
108struct Sparse;
109class Tensor3;
110class Tensor4;
111class Tensor5;
112class Tensor6;
113class Tensor7;
114
125class Joker {
126 // Nothing here.
127};
128
129// Declare existence of the global joker object:
130extern const Joker joker;
131
132// Declare the existence of class ConstMatrixView:
133class ConstIterator1D;
134
135// Declare the existence of class VectorView:
136class VectorView;
137
138// Declare the existence of class ConstVectorView:
139class ConstVectorView;
140
141// Declare the existence of class ConstMatrixView:
142class ConstMatrixView;
143
159class Range {
160 public:
161 // Constructors:
162
173 constexpr Range(Index start, Index extent, Index stride = 1) ARTS_NOEXCEPT
174 : mstart(start),
175 mextent(extent),
176 mstride(stride) {
177 // Start must be >= 0:
178 ARTS_ASSERT(0 <= mstart);
179 // Extent also. Although internally negative extent means "to the end",
180 // this can not be created this way, only with the joker. Zero
181 // extent is allowed, though, which corresponds to an empty range.
182 ARTS_ASSERT(0 <= mextent);
183 // Stride can be anything except 0.
184 // SAB 2001-09-21: Allow 0 stride.
185 // ARTS_ASSERT( 0!=mstride);
186 }
187
190 constexpr Range(Index start, Joker, Index stride = 1) ARTS_NOEXCEPT
191 : mstart(start),
192 mextent(-1),
193 mstride(stride) {
194 // Start must be >= 0:
195 ARTS_ASSERT(0 <= mstart);
196 }
197
201 constexpr Range(Joker, Index stride = 1) noexcept
202 : mstart(0),
203 mextent(-1),
204 mstride(stride){
205 // Nothing to do here.
206 };
207
213 constexpr Range(Index max_size, const Range& r) ARTS_NOEXCEPT
214 : mstart(r.mstart),
215 mextent(r.mextent),
216 mstride(r.mstride) {
217 // Start must be >= 0:
218 ARTS_ASSERT(0 <= mstart);
219 // ... and < max_size:
220 ARTS_ASSERT(mstart < max_size);
221
222 // Stride must be != 0:
223 ARTS_ASSERT(0 != mstride);
224
225 // Convert negative extent (joker) to explicit extent
226 if (mextent < 0) {
227 if (0 < mstride)
228 mextent = 1 + (max_size - 1 - mstart) / mstride;
229 else
230 mextent = 1 + (0 - mstart) / mstride;
231 } else {
232#ifndef NDEBUG
233 // Check that extent is ok:
234 Index fin = mstart + (mextent - 1) * mstride;
235 ARTS_ASSERT(0 <= fin);
236 ARTS_ASSERT(fin < max_size);
237#endif
238 }
239 }
240
246 constexpr Range(const Range& p, const Range& n) ARTS_NOEXCEPT
247 : mstart(p.mstart + n.mstart * p.mstride),
248 mextent(n.mextent),
249 mstride(p.mstride* n.mstride) {
250 // We have to juggle here a bit with previous, new, and resulting
251 // quantities. I.e.;
252 // p.mstride: Previous stride
253 // n.mstride: New stride (as specified)
254 // mstride: Resulting stride (old*new)
255
256 // Get the previous final element:
257 Index prev_fin = p.mstart + (p.mextent - 1) * p.mstride;
258
259 // Resulting start must be >= previous start:
260 ARTS_ASSERT(p.mstart <= mstart);
261 // and <= prev_fin, except for Joker:
262 ARTS_ASSERT(mstart <= prev_fin || mextent == -1);
263
264 // Resulting stride must be != 0:
265 ARTS_ASSERT(0 != mstride,
266 "Problem: mstride==",
267 mstride,
268 " for mstart==",
269 mstart,
270 " and mextent==",
271 mextent);
272
273 // Convert negative extent (joker) to explicit extent
274 if (mextent < 0) {
275 if (0 < mstride)
276 mextent = 1 + (prev_fin - mstart) / mstride;
277 else
278 mextent = 1 + (p.mstart - mstart) / mstride;
279 } else {
280#ifndef NDEBUG
281 // Check that extent is ok:
282 Index fin = mstart + (mextent - 1) * mstride;
283 ARTS_ASSERT(p.mstart <= fin);
284 ARTS_ASSERT(fin <= prev_fin);
285#endif
286 }
287 };
288
289 // Friends:
290 friend class ConstVectorView;
291 friend class VectorView;
292 friend class Vector;
293 friend class ConstMatrixView;
294 friend class MatrixView;
295 friend class Matrix;
296 friend class Iterator2D;
297 friend class Iterator3D;
298 friend class Iterator4D;
299 friend class Iterator5D;
300 friend class Iterator6D;
301 friend class Iterator7D;
302 friend class ConstIterator2D;
303 friend class ConstIterator3D;
304 friend class ConstIterator4D;
305 friend class ConstIterator5D;
306 friend class ConstIterator6D;
307 friend class ConstIterator7D;
308 friend class ConstTensor3View;
309 friend class Tensor3View;
310 friend class Tensor3;
311 friend class ConstTensor4View;
312 friend class Tensor4View;
313 friend class Tensor4;
314 friend class ConstTensor5View;
315 friend class Tensor5View;
316 friend class Tensor5;
317 friend class ConstTensor6View;
318 friend class Tensor6View;
319 friend class Tensor6;
320 friend class ConstTensor7View;
321 friend class Tensor7View;
322 friend class Tensor7;
323 friend struct Sparse;
325 friend class ComplexVectorView;
326 friend class ComplexVector;
328 friend class ComplexMatrixView;
329 friend class ComplexMatrix;
330 friend class ComplexIterator2D;
332
333 friend void mult_general(VectorView,
334 const ConstMatrixView&,
336
337 // Member functions:
338
340 [[nodiscard]] constexpr Index get_start() const noexcept { return mstart; }
342 [[nodiscard]] constexpr Index get_extent() const noexcept { return mextent; }
344 [[nodiscard]] constexpr Index get_stride() const noexcept { return mstride; }
345
347 constexpr Range operator()(const Range r) const ARTS_NOEXCEPT {
348 return (r.mextent < 0) ? (mextent < 0) ? Range(mstart + r.mstart * mstride,
349 joker,
350 r.mstride * mstride)
351 : Range(mstart + r.mstart * mstride,
352 mextent,
353 r.mstride * mstride)
354 : Range(mstart + r.mstart * mstride,
355 r.mextent,
356 r.mstride * mstride);
357 }
358
359 constexpr Index operator()(const Index i) const ARTS_NOEXCEPT {
360 return mstart + i * mstride;
361 };
362
363 friend std::ostream& operator<<(std::ostream& os, const Range& r);
364
365 private:
373};
374
376template <size_t N>
377struct Shape {
378 std::array<Index, N> data;
379 bool operator==(const Shape& other) { return data == other.data; }
380 bool operator!=(const Shape& other) { return data not_eq other.data; }
381 friend std::ostream& operator<<(std::ostream& os, const Shape& shape) {
382 os << shape.data[0];
383 for (size_t i = 1; i < N; i++) os << 'x' << shape.data[i];
384 return os;
385 }
386};
387
391 public:
394 using pointer = Numeric*;
396 using iterator_category = std::random_access_iterator_tag;
397
399 Iterator1D() = default;
400
403 : mx(x),
404 mstride(stride) { /* Nothing to do here. */
405 }
406
407 // Operators:
408
411 mx += mstride;
412 return *this;
413 }
414
416 Numeric& operator*() const ARTS_NOEXCEPT { return *mx; }
417
419 bool operator!=(const Iterator1D& other) const ARTS_NOEXCEPT {
420 if (mx != other.mx) return true;
421 return false;
422 }
423
424 bool operator==(const Iterator1D& other) const ARTS_NOEXCEPT {
425 return !operator!=(other);
426 }
427
429 return (Index)(mx - other.mx) / mstride;
430 }
431
436 friend void copy(ConstIterator1D origin,
437 const ConstIterator1D& end,
438 Iterator1D target);
439
440 private:
442 Numeric* mx{nullptr};
445};
446
450 public:
452 using value_type = const Numeric;
453 using pointer = const Numeric*;
454 using reference = const Numeric&;
455 using iterator_category = std::random_access_iterator_tag;
456
458 ConstIterator1D() = default;
459
462 : mx(x),
463 mstride(stride) { /* Nothing to do here. */
464 }
465
466 // Operators:
469 mx += mstride;
470 return *this;
471 }
472
474 const Numeric& operator*() const ARTS_NOEXCEPT { return *mx; }
475
477 bool operator!=(const ConstIterator1D& other) const ARTS_NOEXCEPT {
478 if (mx != other.mx) return true;
479 return false;
480 }
481
482 bool operator==(const ConstIterator1D& other) const ARTS_NOEXCEPT {
483 return !operator!=(other);
484 }
485
487 return (Index)(mx - other.mx) / mstride;
488 }
489
490 friend void copy(ConstIterator1D origin,
491 const ConstIterator1D& end,
492 Iterator1D target);
493
494 private:
496 const Numeric* mx{nullptr};
499};
500
501// Declare the vector class:
502class Vector;
503
504// Declare the MatrixView class
505class MatrixView;
506
513 public:
514 constexpr ConstVectorView(const ConstVectorView&) = default;
515 constexpr ConstVectorView(ConstVectorView&&) = default;
518
519 // Typedef for compatibility with STL
521
522 // Member functions:
523
525 [[nodiscard]] bool empty() const noexcept { return mrange.mextent == 0; };
526
536 [[nodiscard]] Index nelem() const noexcept { return mrange.mextent; }
537
539 [[nodiscard]] Index size() const noexcept { return mrange.mextent; }
540
542 [[nodiscard]] Shape<1> shape() const { return {nelem()}; }
543
545 [[nodiscard]] Numeric sum() const ARTS_NOEXCEPT;
546
547 // Const index operators:
549 Numeric operator[](Index n) const ARTS_NOEXCEPT { // Check if index is valid:
550 ARTS_ASSERT(0 <= n);
552 return get(n);
553 }
554
556 [[nodiscard]] Numeric get(Index n) const ARTS_NOEXCEPT {
557 return *(mdata + mrange.mstart + n * mrange.mstride);
558 }
559
564
565 friend Numeric operator*(const ConstVectorView& a,
567
568 // Functions returning iterators:
569
571 [[nodiscard]] ConstIterator1D begin() const ARTS_NOEXCEPT;
572
574 [[nodiscard]] ConstIterator1D end() const ARTS_NOEXCEPT;
575
577 operator ConstMatrixView() const;
578
580 virtual ~ConstVectorView() = default;
581
582 // Friends:
583 friend class VectorView;
584 friend class ConstIterator2D;
585 friend class ConstMatrixView;
586 friend class ConstTensor3View;
587 friend class ConstTensor4View;
588 friend class ConstTensor5View;
589 friend class ConstTensor6View;
590 friend class ConstTensor7View;
592 friend int poly_root_solve(Matrix& roots, Vector& coeffs);
593 friend void mult(VectorView, const ConstMatrixView&, const ConstVectorView&);
594 friend void mult(VectorView, const Sparse&, ConstVectorView);
595 friend void transpose_mult(VectorView, const Sparse&, ConstVectorView);
596 friend void mult_general(VectorView,
597 const ConstMatrixView&,
599 friend void lubacksub(VectorView,
602 const ArrayOfIndex&);
604
605 friend std::ostream& operator<<(std::ostream& os, const ConstVectorView& v);
606
610
612 [[nodiscard]] Index selem() const noexcept {return mrange.mstart;}
613
615 [[nodiscard]] Index delem() const noexcept {return mrange.mstride;}
616
622 [[nodiscard]] Numeric* get_c_array() const noexcept {return mdata;}
623
624 protected:
625 // Constructors:
626 ConstVectorView() = default;
627
630 ConstVectorView(Numeric* data, const Range& range) ARTS_NOEXCEPT;
631
642 ConstVectorView(Numeric* data, const Range& p, const Range& n) ARTS_NOEXCEPT;
643
644 // Data members:
645 // -------------
649 Numeric* mdata{nullptr};
650};
651
664 public:
665 // Make const methods visible from base class
668 using ConstVectorView::operator[];
670
671 constexpr VectorView(const VectorView&) = default;
672
675 VectorView(const Vector&);
676
679
680 // Typedef for compatibility with STL
682
684 Numeric& operator[](Index n) ARTS_NOEXCEPT { // Check if index is valid:
685 ARTS_ASSERT(0 <= n);
687 return get(n);
688 }
689
692 return *(mdata + mrange.mstart + n * mrange.mstride);
693 }
694
699
700 // Iterators:
701
704
707
708 // Assignment operators:
709
714 VectorView& operator=(const ConstVectorView& v);
715
721 VectorView& operator=(const VectorView& v);
722
725 VectorView& operator=(const Vector& v);
726
727 VectorView& operator=(const Array<Numeric>& v);
728
731 VectorView& operator=(Numeric x);
732
733 // Other operators:
734
736 VectorView operator*=(Numeric x) ARTS_NOEXCEPT;
737
739 VectorView operator/=(Numeric x) ARTS_NOEXCEPT;
740
742 VectorView operator+=(Numeric x) ARTS_NOEXCEPT;
743
745 VectorView operator-=(Numeric x) ARTS_NOEXCEPT;
746
748 VectorView operator*=(const ConstVectorView& x) ARTS_NOEXCEPT;
749
751 VectorView operator/=(const ConstVectorView& x) ARTS_NOEXCEPT;
752
754 VectorView operator+=(const ConstVectorView& x) ARTS_NOEXCEPT;
755
757 VectorView operator-=(const ConstVectorView& x) ARTS_NOEXCEPT;
758
760 operator MatrixView() ARTS_NOEXCEPT;
761
763 virtual ~VectorView() = default;
764
765 // Friends:
766 friend class ConstIterator2D;
767 friend class Iterator2D;
768 friend class MatrixView;
769 friend class Tensor3View;
770 friend class Tensor4View;
771 friend class Tensor5View;
772 friend class Tensor6View;
773 friend class Tensor7View;
774 friend class ComplexVectorView;
775
779
780 protected:
781 // Constructors:
782 VectorView() = default;
783
786 VectorView(Numeric* data, const Range& range) ARTS_NOEXCEPT;
787
798 VectorView(Numeric* data, const Range& p, const Range& n) ARTS_NOEXCEPT;
799};
800
805 public:
806 // Constructors:
808 Iterator2D() = default;
809
812 : msv(x),
813 mstride(stride) { /* Nothing to do here. */
814 }
815
816 // Operators:
819 msv.mdata += mstride;
820 return *this;
821 }
822
824 bool operator!=(const Iterator2D& other) const ARTS_NOEXCEPT {
825 if (msv.mdata + msv.mrange.mstart !=
826 other.msv.mdata + other.msv.mrange.mstart)
827 return true;
828 return false;
829 }
830
834
837
838 private:
842 Index mstride{0};
843};
844
849 public:
850 // Constructors:
852 ConstIterator2D() = default;
853
856 : msv(x),
857 mstride(stride) { /* Nothing to do here. */
858 }
859
860 // Operators:
863 msv.mdata += mstride;
864 return *this;
865 }
866
868 bool operator!=(const ConstIterator2D& other) const ARTS_NOEXCEPT {
869 if (msv.mdata + msv.mrange.mstart !=
870 other.msv.mdata + other.msv.mrange.mstart)
871 return true;
872 return false;
873 }
874
877 const ConstVectorView* operator->() const ARTS_NOEXCEPT { return &msv; }
878
880 const ConstVectorView& operator*() const ARTS_NOEXCEPT { return msv; }
881
882 private:
886 Index mstride{0};
887};
888
899class Vector : public VectorView {
900 public:
901 // Constructors:
902 Vector() = default;
903
905 Vector(std::initializer_list<Numeric> init);
906
908 explicit Vector(const matpack::vector_like auto& init) : Vector(init.size()) {
909 for (Index i=0; i<size(); i++) operator[](i) = init[i];
910 }
911
913 explicit Vector(Index n);
914
916 Vector(Index n, Numeric fill);
917
925 Vector(Numeric start, Index extent, Numeric stride);
926
932 Vector(const ConstVectorView& v);
933
936 Vector(const Vector& v);
937 Vector(Vector&& v) noexcept : VectorView(std::forward<VectorView>(v)) {
938 v.mdata = nullptr;
939 }
940
942 Vector(const std::vector<Numeric>&);
943
953 ARTS_ASSERT(r0.get_extent() >= 0, "Must have size. Has: ", r0.get_extent());
954 }
955
956 // Assignment operators:
957
982 Vector& operator=(const Vector& v);
983
985 Vector& operator=(Vector&& v) noexcept;
986
988 Vector& operator=(std::initializer_list<Numeric> v);
989
1006 Vector& operator=(const Array<Numeric>& x);
1007
1011
1015 void resize(Index n);
1016
1018 friend void swap(Vector& v1, Vector& v2);
1019
1022 virtual ~Vector();
1023
1024 template <class F>
1025 void transform_elementwise(F&& func) {
1026 std::transform(mdata, mdata + size(), mdata, func);
1027 }
1028};
1029
1030// Declare class Matrix:
1031class Matrix;
1032
1044 public:
1045 constexpr ConstMatrixView(const ConstMatrixView&) = default;
1046 constexpr ConstMatrixView(ConstMatrixView&&) = default;
1049
1050 // Typedef for compatibility with STL
1052
1053 // Member functions:
1054 [[nodiscard]] Index selem() const noexcept { return mrr.mstart + mcr.mstart; }
1055 [[nodiscard]] Index nrows() const noexcept { return mrr.mextent; }
1056 [[nodiscard]] Index ncols() const noexcept { return mcr.mextent; }
1057 [[nodiscard]] Index drows() const noexcept { return mrr.mstride; }
1058 [[nodiscard]] Index dcols() const noexcept { return mcr.mstride; }
1059
1060 // Total size
1061 [[nodiscard]] Index size() const noexcept { return nrows() * ncols(); }
1062 [[nodiscard]] bool empty() const noexcept { return size() == 0; }
1063
1065 [[nodiscard]] Shape<2> shape() const { return {nrows(), ncols()}; }
1066
1067 // Const index operators:
1070 ARTS_NOEXCEPT { // Check if indices are valid:
1071 ARTS_ASSERT(0 <= r);
1072 ARTS_ASSERT(0 <= c);
1073 ARTS_ASSERT(r < mrr.mextent);
1074 ARTS_ASSERT(c < mcr.mextent);
1075
1076 return get(r, c);
1077 }
1078
1080 [[nodiscard]] Numeric get(Index r, Index c) const ARTS_NOEXCEPT {
1081 return *(mdata + mrr.mstart + r * mrr.mstride + mcr.mstart +
1082 c * mcr.mstride);
1083 }
1084
1085 ConstMatrixView operator()(const Range& r,
1086 const Range& c) const ARTS_NOEXCEPT;
1087 ConstVectorView operator()(const Range& r, Index c) const ARTS_NOEXCEPT;
1088 ConstVectorView operator()(Index r, const Range& c) const ARTS_NOEXCEPT;
1089
1090 // Functions returning iterators:
1091 [[nodiscard]] ConstIterator2D begin() const ARTS_NOEXCEPT;
1092 [[nodiscard]] ConstIterator2D end() const ARTS_NOEXCEPT;
1093
1094 // View on diagonal vector
1095 [[nodiscard]] ConstVectorView diagonal() const ARTS_NOEXCEPT;
1096
1098 virtual ~ConstMatrixView() = default;
1099
1100 // Friends:
1101 friend class MatrixView;
1102 friend class ConstIterator3D;
1103 friend class ConstVectorView;
1104 friend class ConstTensor3View;
1105 friend class ConstTensor4View;
1106 friend class ConstTensor5View;
1107 friend class ConstTensor6View;
1108 friend class ConstTensor7View;
1111 friend int poly_root_solve(Matrix& roots, Vector& coeffs);
1112 friend void mult(VectorView, const ConstMatrixView&, const ConstVectorView&);
1113 friend void mult(MatrixView, const ConstMatrixView&, const ConstMatrixView&);
1114 friend void mult(MatrixView, const Sparse&, const ConstMatrixView&);
1115 friend void mult(MatrixView, const ConstMatrixView&, const Sparse&);
1116 friend void mult_general(VectorView,
1117 const ConstMatrixView&,
1119 friend void mult_general(MatrixView,
1120 const ConstMatrixView&,
1122 friend void ludcmp(Matrix&, ArrayOfIndex&, ConstMatrixView);
1123 friend void lubacksub(VectorView,
1126 const ArrayOfIndex&);
1127 friend void inv(MatrixView, ConstMatrixView);
1129
1130 friend std::ostream& operator<<(std::ostream& os, const ConstMatrixView& v);
1131
1132 // Conversion to a plain C-array
1133 [[nodiscard]] Numeric* get_c_array() const noexcept {return mdata;}
1134
1135 protected:
1136 // Constructors:
1137 ConstMatrixView() = default;
1138 ConstMatrixView(Numeric* data, const Range& r, const Range& c) ARTS_NOEXCEPT;
1140 const Range& pr,
1141 const Range& pc,
1142 const Range& nr,
1143 const Range& nc) ARTS_NOEXCEPT;
1144
1145 // Data members:
1146 // -------------
1148 Range mrr{0, 0, 1};
1150 Range mcr{0, 0, 1};
1152 Numeric* mdata{nullptr};
1153};
1154
1165 public:
1166 // Make const methods visible from base class
1169 using ConstMatrixView::operator();
1171
1172 constexpr MatrixView(const MatrixView&) = default;
1173
1174 // Typedef for compatibility with STL
1176
1177 // Index Operators:
1180 Index c) ARTS_NOEXCEPT { // Check if indices are valid:
1181 ARTS_ASSERT(0 <= r);
1182 ARTS_ASSERT(0 <= c);
1183 ARTS_ASSERT(r < mrr.mextent);
1184 ARTS_ASSERT(c < mcr.mextent);
1185
1186 return get(r, c);
1187 }
1188
1191 return *(mdata + mrr.mstart + r * mrr.mstride + mcr.mstart +
1192 c * mcr.mstride);
1193 }
1194
1195 MatrixView operator()(const Range& r, const Range& c) ARTS_NOEXCEPT;
1196 VectorView operator()(const Range& r, Index c) ARTS_NOEXCEPT;
1197 VectorView operator()(Index r, const Range& c) ARTS_NOEXCEPT;
1198
1199 // Functions returning iterators:
1202
1203 // Assignment operators:
1204 MatrixView& operator=(const ConstMatrixView& m);
1205 MatrixView& operator=(const MatrixView& m);
1206 MatrixView& operator=(const Matrix& m);
1207 MatrixView& operator=(const ConstVectorView& m);
1208 MatrixView& operator=(Numeric x);
1209
1210 // Other operators:
1211 MatrixView& operator*=(Numeric x) ARTS_NOEXCEPT;
1212 MatrixView& operator/=(Numeric x) ARTS_NOEXCEPT;
1213 MatrixView& operator+=(Numeric x) ARTS_NOEXCEPT;
1214 MatrixView& operator-=(Numeric x) ARTS_NOEXCEPT;
1215
1216 MatrixView& operator*=(const ConstMatrixView& x) ARTS_NOEXCEPT;
1217 MatrixView& operator/=(const ConstMatrixView& x) ARTS_NOEXCEPT;
1218 MatrixView& operator+=(const ConstMatrixView& x) ARTS_NOEXCEPT;
1219 MatrixView& operator-=(const ConstMatrixView& x) ARTS_NOEXCEPT;
1220
1221 MatrixView& operator*=(const ConstVectorView& x) ARTS_NOEXCEPT;
1222 MatrixView& operator/=(const ConstVectorView& x) ARTS_NOEXCEPT;
1223 MatrixView& operator+=(const ConstVectorView& x) ARTS_NOEXCEPT;
1224 MatrixView& operator-=(const ConstVectorView& x) ARTS_NOEXCEPT;
1225
1227 virtual ~MatrixView() = default;
1228
1229 // Friends:
1230 friend class VectorView;
1231 friend class Iterator3D;
1232 friend class Tensor3View;
1233 friend class Tensor4View;
1234 friend class Tensor5View;
1235 friend class Tensor6View;
1236 friend class Tensor7View;
1237 friend class ComplexMatrixView;
1240 friend void mult(MatrixView, const ConstMatrixView&, const ConstMatrixView&);
1241
1242 protected:
1243 // Constructors:
1244 MatrixView() = default;
1245 MatrixView(Numeric* data, const Range& rr, const Range& cr) ARTS_NOEXCEPT;
1246 MatrixView(Numeric* data,
1247 const Range& pr,
1248 const Range& pc,
1249 const Range& nr,
1250 const Range& nc) ARTS_NOEXCEPT;
1251};
1252
1261class Matrix : public MatrixView {
1262 public:
1263 // Constructors:
1264 Matrix() = default;
1265 Matrix(Index r, Index c);
1266 Matrix(Index r, Index c, Numeric fill);
1267 Matrix(const ConstMatrixView& m);
1268 Matrix(const Matrix& m);
1269 Matrix(Matrix&& v) noexcept : MatrixView(std::forward<MatrixView>(v)) {
1270 v.mdata = nullptr;
1271 }
1272
1274 explicit Matrix(const matpack::matrix_like auto& init) : Matrix(matpack::row_size(init), matpack::column_size(init)) {
1275 for (Index i=0; i<nrows(); i++) for (Index j=0; j<ncols(); j++) operator()(i, j) = init(i, j);
1276 }
1277
1288 : MatrixView(d, r0, r1) {
1289 ARTS_ASSERT(r0.get_extent() >= 0, "Must have size. Has: ", r0.get_extent());
1290 ARTS_ASSERT(r1.get_extent() >= 0, "Must have size. Has: ", r1.get_extent());
1291 }
1292
1293 // Assignment operators:
1294 Matrix& operator=(const Matrix& m);
1295 Matrix& operator=(Matrix&& m) noexcept;
1298
1299 // Resize function:
1300 void resize(Index r, Index c);
1301
1302 // Swap function:
1303 friend void swap(Matrix& m1, Matrix& m2);
1304
1305 // Destructor:
1306 virtual ~Matrix();
1307
1309 template <std::size_t dim0>
1311 static_assert(dim0 < 2, "Bad Dimension, Out-of-Bounds");
1312
1313 Range r0(0, dim0 == 0 ? nrows() : ncols());
1314
1315 Vector out(mdata, r0);
1316 ARTS_ASSERT(size() == out.size(),
1317 "Can only reduce size on same size input. "
1318 "The sizes are (input v. output): ",
1319 size(),
1320 " v. ",
1321 out.size());
1322 mdata = nullptr;
1323 return out;
1324 }
1325
1327
1328 template <class F>
1329 void transform_elementwise(F&& func) {
1330 std::transform(mdata, mdata + size(), mdata, func);
1331 }
1332};
1333
1334// Function declarations:
1335// ----------------------
1336
1337void copy(ConstIterator1D origin,
1338 const ConstIterator1D& end,
1339 Iterator1D target);
1340
1342void copy(Numeric x, Iterator1D target, const Iterator1D& end) ARTS_NOEXCEPT;
1343
1344void copy(ConstIterator2D origin,
1345 const ConstIterator2D& end,
1346 Iterator2D target);
1347
1348void copy(Numeric x, Iterator2D target, const Iterator2D& end) ARTS_NOEXCEPT;
1349
1350void mult(VectorView y, const ConstMatrixView& M, const ConstVectorView& x);
1351
1353 const ConstMatrixView& B,
1355
1356void mult(MatrixView A, const ConstMatrixView& B, const ConstMatrixView& C);
1357
1359 const ConstMatrixView& B,
1361
1362void cross3(VectorView c,
1363 const ConstVectorView& a,
1365
1367
1369
1371
1373
1374void transform(VectorView y, double (&my_func)(double), ConstVectorView x);
1375
1376void transform(MatrixView y, double (&my_func)(double), ConstMatrixView x);
1377
1379
1381
1383
1385
1387
1389
1391
1394
1396// Helper function for debugging
1397#ifndef NDEBUG
1398
1400
1401#endif
1403
1406
1408
1411
1413
1414#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:449
bool operator!=(const ConstIterator1D &other) const ARTS_NOEXCEPT
Not equal operator, needed for algorithms like copy.
Definition: matpackI.h:477
const Numeric * pointer
Definition: matpackI.h:453
std::random_access_iterator_tag iterator_category
Definition: matpackI.h:455
const Numeric * mx
Current position.
Definition: matpackI.h:496
Index mstride
Stride.
Definition: matpackI.h:498
ConstIterator1D()=default
Default constructor.
const Numeric value_type
Definition: matpackI.h:452
const Numeric & reference
Definition: matpackI.h:454
bool operator==(const ConstIterator1D &other) const ARTS_NOEXCEPT
Definition: matpackI.h:482
ConstIterator1D & operator++() ARTS_NOEXCEPT
Prefix increment operator.
Definition: matpackI.h:468
Index difference_type
Definition: matpackI.h:451
ConstIterator1D(const Numeric *x, Index stride) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:461
const Numeric & operator*() const ARTS_NOEXCEPT
Dereferencing.
Definition: matpackI.h:474
Index operator-(const ConstIterator1D &other) const ARTS_NOEXCEPT
Definition: matpackI.h:486
friend void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Definition: matpackI.cc:288
The const row iterator class for sub matrices.
Definition: matpackI.h:848
const ConstVectorView & operator*() const ARTS_NOEXCEPT
Dereferencing.
Definition: matpackI.h:880
ConstVectorView msv
Current position.
Definition: matpackI.h:884
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:877
bool operator!=(const ConstIterator2D &other) const ARTS_NOEXCEPT
Not equal operator, needed for algorithms like copy.
Definition: matpackI.h:868
ConstIterator2D(const ConstVectorView &x, Index stride) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:855
ConstIterator2D & operator++() ARTS_NOEXCEPT
Prefix increment operator.
Definition: matpackI.h:862
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:1043
Index drows() const noexcept
Definition: matpackI.h:1057
constexpr ConstMatrixView(ConstMatrixView &&)=default
ConstMatrixView & operator=(const ConstMatrixView &)=default
bool empty() const noexcept
Definition: matpackI.h:1062
Index dcols() const noexcept
Definition: matpackI.h:1058
ConstIterator2D begin() const ARTS_NOEXCEPT
Return const iterator to first row.
Definition: matpackI.cc:448
Index size() const noexcept
Definition: matpackI.h:1061
Index nrows() const noexcept
Definition: matpackI.h:1055
constexpr ConstMatrixView(const ConstMatrixView &)=default
Shape< 2 > shape() const
Definition: matpackI.h:1065
ConstIterator2D end() const ARTS_NOEXCEPT
Return const iterator behind last row.
Definition: matpackI.cc:453
Index ncols() const noexcept
Definition: matpackI.h:1056
Numeric operator()(Index r, Index c) const ARTS_NOEXCEPT
Plain const index operator.
Definition: matpackI.h:1069
ConstMatrixView & operator=(ConstMatrixView &&)=default
ConstMatrixView()=default
Index selem() const noexcept
Definition: matpackI.h:1054
Numeric get(Index r, Index c) const ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:1080
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:512
Numeric get(Index n) const ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:556
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:92
ConstVectorView()=default
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:649
Index delem() const noexcept
Steps in memory between elements.
Definition: matpackI.h:615
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:1136
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:549
constexpr ConstVectorView(ConstVectorView &&)=default
constexpr ConstVectorView(const ConstVectorView &)=default
friend void mult(VectorView, const ConstMatrixView &, const ConstVectorView &)
Matrix-Vector Multiplication.
Definition: matpackI.cc:1078
Index selem() const noexcept
Start element in memory.
Definition: matpackI.h:612
Shape< 1 > shape() const
Definition: matpackI.h:542
Range mrange
The range of mdata that is actually used.
Definition: matpackI.h:647
friend void transpose_mult(VectorView, const Sparse &, ConstVectorView)
Sparse matrix - Vector multiplication.
Definition: matpackII.cc:452
Numeric * get_c_array() const noexcept
Conversion to plain C-array, const-version.
Definition: matpackI.h:622
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
friend class ConstMatrixView
Definition: matpackI.h:585
friend void diagonalize(MatrixView, VectorView, VectorView, ConstMatrixView)
Matrix Diagonalization.
Definition: lin_alg.cc:250
Index size() const noexcept
Definition: matpackI.h:539
ConstVectorView & operator=(const ConstVectorView &)=default
ConstVectorView & operator=(ConstVectorView &&)=default
friend Numeric operator*(const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
Scalar product.
Definition: matpackI.cc:1046
bool empty() const noexcept
Returns true if variable size is zero.
Definition: matpackI.h:525
The iterator class for sub vectors.
Definition: matpackI.h:390
Iterator1D & operator++() ARTS_NOEXCEPT
Prefix increment operator.
Definition: matpackI.h:410
bool operator==(const Iterator1D &other) const ARTS_NOEXCEPT
Definition: matpackI.h:424
Index operator-(const Iterator1D &other) const ARTS_NOEXCEPT
Definition: matpackI.h:428
Iterator1D()=default
Default constructor.
Index difference_type
Definition: matpackI.h:392
Numeric & operator*() const ARTS_NOEXCEPT
Dereferencing.
Definition: matpackI.h:416
Index mstride
Stride.
Definition: matpackI.h:444
std::random_access_iterator_tag iterator_category
Definition: matpackI.h:396
Numeric * mx
Current position.
Definition: matpackI.h:442
friend void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Copy data between begin and end to target.
Definition: matpackI.cc:288
Numeric * pointer
Definition: matpackI.h:394
bool operator!=(const Iterator1D &other) const ARTS_NOEXCEPT
Not equal operator, needed for algorithms like copy.
Definition: matpackI.h:419
Numeric & reference
Definition: matpackI.h:395
Numeric value_type
Definition: matpackI.h:393
Iterator1D(Numeric *x, Index stride) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:402
The row iterator class for sub matrices.
Definition: matpackI.h:804
Iterator2D(const VectorView &x, Index stride) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:811
VectorView & operator*() ARTS_NOEXCEPT
Dereferencing.
Definition: matpackI.h:836
bool operator!=(const Iterator2D &other) const ARTS_NOEXCEPT
Not equal operator, needed for algorithms like copy.
Definition: matpackI.h:824
Iterator2D()=default
Default constructor.
VectorView msv
Current position.
Definition: matpackI.h:840
Iterator2D & operator++() ARTS_NOEXCEPT
Prefix increment operator.
Definition: matpackI.h:818
VectorView * operator->() ARTS_NOEXCEPT
The -> operator is needed, so that we can write i->begin() to get the 1D iterators.
Definition: matpackI.h:833
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:125
The MatrixView class.
Definition: matpackI.h:1164
Numeric & operator()(Index r, Index c) ARTS_NOEXCEPT
Plain index operator.
Definition: matpackI.h:1179
Numeric & get(Index r, Index c) ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:1190
constexpr MatrixView(const MatrixView &)=default
The Matrix class.
Definition: matpackI.h:1261
Matrix(Matrix &&v) noexcept
Definition: matpackI.h:1269
Matrix()=default
Vector reduce_rank() &&ARTS_NOEXCEPT
Definition: matpackI.h:1310
Matrix(const matpack::matrix_like auto &init)
Initialization from a vector type.
Definition: matpackI.h:1274
Numeric * get_raw_data() ARTS_NOEXCEPT
Definition: matpackI.h:1326
Matrix(Numeric *d, const Range &r0, const Range &r1) ARTS_NOEXCEPT
Definition: matpackI.h:1287
void transform_elementwise(F &&func)
Definition: matpackI.h:1329
The range class.
Definition: matpackI.h:159
constexpr Range(Index start, Joker, Index stride=1) ARTS_NOEXCEPT
Constructor with joker extent.
Definition: matpackI.h:190
friend void mult_general(VectorView, const ConstMatrixView &, const ConstVectorView &) ARTS_NOEXCEPT
Matrix Vector multiplication.
Definition: matpackI.cc:1136
constexpr Index operator()(const Index i) const ARTS_NOEXCEPT
Definition: matpackI.h:359
constexpr Range(const Range &p, const Range &n) ARTS_NOEXCEPT
Constructor of a new range relative to an old range.
Definition: matpackI.h:246
constexpr Range(Index start, Index extent, Index stride=1) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:173
Index mstart
The start index.
Definition: matpackI.h:367
constexpr Index get_extent() const noexcept
Returns the extent of the range.
Definition: matpackI.h:342
constexpr Range operator()(const Range r) const ARTS_NOEXCEPT
Range of range.
Definition: matpackI.h:347
friend std::ostream & operator<<(std::ostream &os, const Range &r)
Definition: matpackI.cc:39
Index mstride
The stride.
Definition: matpackI.h:372
constexpr Index get_start() const noexcept
Returns the start index of the range.
Definition: matpackI.h:340
Index mextent
The number of elements.
Definition: matpackI.h:370
constexpr Range(Index max_size, const Range &r) ARTS_NOEXCEPT
Constructor which converts a range with joker to an explicit range.
Definition: matpackI.h:213
constexpr Index get_stride() const noexcept
Returns the stride of the range.
Definition: matpackI.h:344
constexpr Range(Joker, Index stride=1) noexcept
Constructor with just a joker.
Definition: matpackI.h:201
The Tensor3View class.
Definition: matpackIII.h:246
The Tensor3 class.
Definition: matpackIII.h:346
The Tensor4View class.
Definition: matpackIV.h:292
The Tensor4 class.
Definition: matpackIV.h:429
The Tensor5View class.
Definition: matpackV.h:343
The Tensor5 class.
Definition: matpackV.h:516
The Tensor6View class.
Definition: matpackVI.h:632
The Tensor6 class.
Definition: matpackVI.h:1099
The Tensor7View class.
Definition: matpackVII.h:1303
The Tensor7 class.
Definition: matpackVII.h:2399
The VectorView class.
Definition: matpackI.h:663
Numeric & get(Index n) ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:691
VectorView()=default
friend class MatrixView
Definition: matpackI.h:768
friend class ConstIterator2D
Definition: matpackI.h:766
constexpr VectorView(const VectorView &)=default
VectorView & operator=(const ConstVectorView &v)
Assignment operator.
Definition: matpackI.cc:153
friend class Iterator2D
Definition: matpackI.h:767
Numeric & operator[](Index n) ARTS_NOEXCEPT
Plain Index operator.
Definition: matpackI.h:684
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:899
Vector(Numeric *d, const Range &r0) ARTS_NOEXCEPT
Definition: matpackI.h:952
Vector(Vector &&v) noexcept
Definition: matpackI.h:937
Vector()=default
Vector(const matpack::vector_like auto &init)
Initialization from a vector type.
Definition: matpackI.h:908
void transform_elementwise(F &&func)
Definition: matpackI.h:1025
Binary output file stream class.
Definition: bofstream.h:42
A concept for an Arts matrix-like type with access operations.
A concept for an Arts vector-like type with access operations.
#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:57
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:168
Numeric mean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version.
Definition: matpackI.cc:1544
Numeric vector_angle(ConstVectorView a, ConstVectorView b)
Definition: matpackI.cc:1369
Numeric max(const ConstVectorView &x) ARTS_NOEXCEPT
Max function, vector version.
Definition: matpackI.cc:1476
Numeric nanmean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version ignoring nans and infs.
Definition: matpackI.cc:1559
void proj(Vector &c, ConstVectorView a, ConstVectorView b) ARTS_NOEXCEPT
Definition: matpackI.cc:1391
const Joker joker
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:1431
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
cross3
Definition: matpackI.cc:1348
Numeric min(const ConstVectorView &x) ARTS_NOEXCEPT
Min function, vector version.
Definition: matpackI.cc:1510
void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Definition: matpackI.cc:288
Numeric debug_matrixview_get_elem(MatrixView &mv, Index r, Index c)
Helper function to access matrix elements.
Definition: matpackI.cc:1646
ConstMatrixView transpose(ConstMatrixView m) ARTS_NOEXCEPT
Const version of transpose.
Definition: matpackI.cc:1401
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.
constexpr Numeric r0
The reference radius in IGRF13.
Definition: igrf13.cc:203
VectorView std(VectorView std, const Vector &y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the standard deviation of the ranged ys.
Definition: raw.cc:205
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
constexpr auto column_size(column_keeper auto &&x)
Get a column size from x.
constexpr auto row_size(row_keeper auto &&x)
Get a row size from x.
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 N
Definition: rng.cc:164
#define M
Definition: rng.cc:165
Helper shape class.
Definition: matpackI.h:377
std::array< Index, N > data
Definition: matpackI.h:378
bool operator==(const Shape &other)
Definition: matpackI.h:379
bool operator!=(const Shape &other)
Definition: matpackI.h:380
friend std::ostream & operator<<(std::ostream &os, const Shape &shape)
Definition: matpackI.h:381
The Sparse class.
Definition: matpackII.h:75
#define d
#define v
#define a
#define c
#define b