ARTS 2.5.10 (git: 2f1c442c)
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#include <utility>
100
101#include "array.h"
102#include "matpack.h"
103#include "matpack_concepts.h"
104
105// Declare existance of some classes
106class bofstream;
107class Vector;
108class Matrix;
109struct Sparse;
110class Tensor3;
111class Tensor4;
112class Tensor5;
113class Tensor6;
114class Tensor7;
115
126class Joker {
127 // Nothing here.
128};
129
130// Declare existence of the global joker object:
131extern const Joker joker;
132
133// Declare the existence of class ConstMatrixView:
134class ConstIterator1D;
135
136// Declare the existence of class VectorView:
137class VectorView;
138
139// Declare the existence of class ConstVectorView:
140class ConstVectorView;
141
142// Declare the existence of class ConstMatrixView:
143class ConstMatrixView;
144
160class Range {
161 public:
162 // Constructors:
163
174 constexpr Range(Index start, Index extent, Index stride = 1) ARTS_NOEXCEPT
175 : mstart(start),
176 mextent(extent),
177 mstride(stride) {
178 // Start must be >= 0:
179 ARTS_ASSERT(0 <= mstart);
180 // Extent also. Although internally negative extent means "to the end",
181 // this can not be created this way, only with the joker. Zero
182 // extent is allowed, though, which corresponds to an empty range.
183 ARTS_ASSERT(0 <= mextent);
184 // Stride can be anything except 0.
185 // SAB 2001-09-21: Allow 0 stride.
186 // ARTS_ASSERT( 0!=mstride);
187 }
188
191 constexpr Range(Index start, Joker, Index stride = 1) ARTS_NOEXCEPT
192 : mstart(start),
193 mextent(-1),
194 mstride(stride) {
195 // Start must be >= 0:
196 ARTS_ASSERT(0 <= mstart);
197 }
198
202 constexpr Range(Joker, Index stride = 1) noexcept
203 : mstart(0),
204 mextent(-1),
205 mstride(stride){
206 // Nothing to do here.
207 };
208
214 constexpr Range(Index max_size, const Range& r) ARTS_NOEXCEPT
215 : mstart(r.mstart),
216 mextent(r.mextent),
217 mstride(r.mstride) {
218 // Start must be >= 0:
219 ARTS_ASSERT(0 <= mstart);
220 // ... and < max_size:
221 ARTS_ASSERT(mstart < max_size);
222
223 // Stride must be != 0:
224 ARTS_ASSERT(0 != mstride);
225
226 // Convert negative extent (joker) to explicit extent
227 if (mextent < 0) {
228 if (0 < mstride)
229 mextent = 1 + (max_size - 1 - mstart) / mstride;
230 else
231 mextent = 1 + (0 - mstart) / mstride;
232 } else {
233#ifndef NDEBUG
234 // Check that extent is ok:
235 Index fin = mstart + (mextent - 1) * mstride;
236 ARTS_ASSERT(0 <= fin);
237 ARTS_ASSERT(fin < max_size);
238#endif
239 }
240 }
241
247 constexpr Range(const Range& p, const Range& n) ARTS_NOEXCEPT
248 : mstart(p.mstart + n.mstart * p.mstride),
249 mextent(n.mextent),
250 mstride(p.mstride* n.mstride) {
251 // We have to juggle here a bit with previous, new, and resulting
252 // quantities. I.e.;
253 // p.mstride: Previous stride
254 // n.mstride: New stride (as specified)
255 // mstride: Resulting stride (old*new)
256
257 // Get the previous final element:
258 Index const prev_fin = p.mstart + (p.mextent - 1) * p.mstride;
259
260 // Resulting start must be >= previous start:
261 ARTS_ASSERT(p.mstart <= mstart);
262 // and <= prev_fin, except for Joker:
263 ARTS_ASSERT(mstart <= prev_fin || mextent == -1);
264
265 // Resulting stride must be != 0:
266 ARTS_ASSERT(0 != mstride,
267 "Problem: mstride==",
268 mstride,
269 " for mstart==",
270 mstart,
271 " and mextent==",
272 mextent);
273
274 // Convert negative extent (joker) to explicit extent
275 if (mextent < 0) {
276 if (0 < mstride)
277 mextent = 1 + (prev_fin - mstart) / mstride;
278 else
279 mextent = 1 + (p.mstart - mstart) / mstride;
280 } else {
281#ifndef NDEBUG
282 // Check that extent is ok:
283 Index fin = mstart + (mextent - 1) * mstride;
284 ARTS_ASSERT(p.mstart <= fin);
285 ARTS_ASSERT(fin <= prev_fin);
286#endif
287 }
288 }
289
290 // Friends:
291 friend class ConstVectorView;
292 friend class VectorView;
293 friend class Vector;
294 friend class ConstMatrixView;
295 friend class MatrixView;
296 friend class Matrix;
297 friend class Iterator2D;
298 friend class Iterator3D;
299 friend class Iterator4D;
300 friend class Iterator5D;
301 friend class Iterator6D;
302 friend class Iterator7D;
303 friend class ConstIterator2D;
304 friend class ConstIterator3D;
305 friend class ConstIterator4D;
306 friend class ConstIterator5D;
307 friend class ConstIterator6D;
308 friend class ConstIterator7D;
309 friend class ConstTensor3View;
310 friend class Tensor3View;
311 friend class Tensor3;
312 friend class ConstTensor4View;
313 friend class Tensor4View;
314 friend class Tensor4;
315 friend class ConstTensor5View;
316 friend class Tensor5View;
317 friend class Tensor5;
318 friend class ConstTensor6View;
319 friend class Tensor6View;
320 friend class Tensor6;
321 friend class ConstTensor7View;
322 friend class Tensor7View;
323 friend class Tensor7;
324 friend struct Sparse;
326 friend class ComplexVectorView;
327 friend class ComplexVector;
329 friend class ComplexMatrixView;
330 friend class ComplexMatrix;
331 friend class ComplexIterator2D;
333
334 friend void mult_general(VectorView,
335 const ConstMatrixView&,
337
338 // Member functions:
339
341 [[nodiscard]] constexpr Index get_start() const noexcept { return mstart; }
343 [[nodiscard]] constexpr Index get_extent() const noexcept { return mextent; }
345 [[nodiscard]] constexpr Index get_stride() const noexcept { return mstride; }
346
348 constexpr Range operator()(const Range r) const ARTS_NOEXCEPT {
349 return (r.mextent < 0) ? (mextent < 0) ? Range(mstart + r.mstart * mstride,
350 joker,
351 r.mstride * mstride)
352 : Range(mstart + r.mstart * mstride,
353 mextent,
354 r.mstride * mstride)
355 : Range(mstart + r.mstart * mstride,
356 r.mextent,
357 r.mstride * mstride);
358 }
359
360 constexpr Index operator()(const Index i) const ARTS_NOEXCEPT {
361 return mstart + i * mstride;
362 };
363
364 friend std::ostream& operator<<(std::ostream& os, const Range& r);
365
366 constexpr void swap(Range& other) noexcept {
367 using std::swap;
368 swap(mstart, other.mstart);
369 swap(mextent, other.mextent);
370 swap(mstride, other.mstride);
371 }
372
373 friend constexpr void swap(Range& a, Range& b) noexcept { a.swap(b); }
374
375 private:
383};
384
386template <size_t N>
387struct Shape {
388 std::array<Index, N> data;
389 bool operator<=>(const Shape& other) const = default;
390 friend std::ostream& operator<<(std::ostream& os, const Shape& shape) {
391 os << shape.data[0];
392 for (size_t i = 1; i < N; i++) os << 'x' << shape.data[i];
393 return os;
394 }
395};
396
400 public:
403 using pointer = Numeric*;
405 using iterator_category = std::random_access_iterator_tag;
406
408 Iterator1D() = default;
409
412 : mx(x),
413 mstride(stride) { /* Nothing to do here. */
414 }
415
416 // Operators:
417
420 mx += mstride;
421 return *this;
422 }
423
425 Numeric& operator*() const ARTS_NOEXCEPT { return *mx; }
426
428 bool operator!=(const Iterator1D& other) const ARTS_NOEXCEPT {
429 if (mx != other.mx) return true;
430 return false;
431 }
432
433 bool operator==(const Iterator1D& other) const ARTS_NOEXCEPT {
434 return !operator!=(other);
435 }
436
438 return (Index)(mx - other.mx) / mstride;
439 }
440
445 friend void copy(ConstIterator1D origin,
446 const ConstIterator1D& end,
447 Iterator1D target);
448
449 private:
451 Numeric* mx{nullptr};
454};
455
459 public:
461 using value_type = const Numeric;
462 using pointer = const Numeric*;
463 using reference = const Numeric&;
464 using iterator_category = std::random_access_iterator_tag;
465
467 ConstIterator1D() = default;
468
471 : mx(x),
472 mstride(stride) { /* Nothing to do here. */
473 }
474
475 // Operators:
478 mx += mstride;
479 return *this;
480 }
481
483 const Numeric& operator*() const ARTS_NOEXCEPT { return *mx; }
484
486 bool operator!=(const ConstIterator1D& other) const ARTS_NOEXCEPT {
487 if (mx != other.mx) return true;
488 return false;
489 }
490
491 bool operator==(const ConstIterator1D& other) const ARTS_NOEXCEPT {
492 return !operator!=(other);
493 }
494
496 return (Index)(mx - other.mx) / mstride;
497 }
498
499 friend void copy(ConstIterator1D origin,
500 const ConstIterator1D& end,
501 Iterator1D target);
502
503 private:
505 const Numeric* mx{nullptr};
508};
509
510// Declare the vector class:
511class Vector;
512
513// Declare the MatrixView class
514class MatrixView;
515
522 public:
523 static constexpr bool matpack_type{true};
524
525 constexpr ConstVectorView(const ConstVectorView&) = default;
526 constexpr ConstVectorView(ConstVectorView&&) = default;
529
530 // Typedef for compatibility with STL
532
533 // Member functions:
534
536 [[nodiscard]] bool empty() const noexcept { return mrange.mextent == 0; };
537
547 [[nodiscard]] Index nelem() const noexcept { return mrange.mextent; }
548
550 [[nodiscard]] Index size() const noexcept { return mrange.mextent; }
551
553 [[nodiscard]] Shape<1> shape() const { return {nelem()}; }
554
556 [[nodiscard]] Numeric sum() const ARTS_NOEXCEPT;
557
558 // Const index operators:
560 Numeric operator[](Index n) const ARTS_NOEXCEPT { // Check if index is valid:
561 ARTS_ASSERT(0 <= n);
563 return get(n);
564 }
565
567 [[nodiscard]] Numeric get(Index n) const ARTS_NOEXCEPT {
568 return *(mdata + mrange.mstart + n * mrange.mstride);
569 }
570
575
576 friend Numeric operator*(const ConstVectorView& a,
578
579 // Functions returning iterators:
580
582 [[nodiscard]] ConstIterator1D begin() const ARTS_NOEXCEPT;
583
585 [[nodiscard]] ConstIterator1D end() const ARTS_NOEXCEPT;
586
588 operator ConstMatrixView() const;
589
591 virtual ~ConstVectorView() = default;
592
593 // Friends:
594 friend class VectorView;
595 friend class ConstIterator2D;
596 friend class ConstMatrixView;
597 friend class ConstTensor3View;
598 friend class ConstTensor4View;
599 friend class ConstTensor5View;
600 friend class ConstTensor6View;
601 friend class ConstTensor7View;
603 friend int poly_root_solve(Matrix& roots, Vector& coeffs);
604 friend void mult(VectorView, const ConstMatrixView&, const ConstVectorView&);
605 friend void mult(VectorView, const Sparse&, ConstVectorView);
606 friend void transpose_mult(VectorView, const Sparse&, ConstVectorView);
607 friend void mult_general(VectorView,
608 const ConstMatrixView&,
610 friend void lubacksub(VectorView,
613 const ArrayOfIndex&);
615
616 friend std::ostream& operator<<(std::ostream& os, const ConstVectorView& v);
617
621
623 [[nodiscard]] Index selem() const noexcept {return mrange.mstart;}
624
626 [[nodiscard]] Index delem() const noexcept {return mrange.mstride;}
627
633 [[nodiscard]] Numeric* get_c_array() const noexcept {return mdata;}
634
635 protected:
636 // Constructors:
637 ConstVectorView() = default;
638
641 ConstVectorView(Numeric* data, const Range& range) ARTS_NOEXCEPT;
642
653 ConstVectorView(Numeric* data, const Range& p, const Range& n) ARTS_NOEXCEPT;
654
655 // Data members:
656 // -------------
660 Numeric* mdata{nullptr};
661};
662
675 public:
676 // Make const methods visible from base class
679 using ConstVectorView::operator[];
681
682 constexpr VectorView(const VectorView&) = default;
683
686 VectorView(const Vector&);
687
690
691 // Typedef for compatibility with STL
693
694#define GETFUN(n) *(mdata + mrange.mstart + n * mrange.mstride)
696 Numeric& operator[](Index n) ARTS_NOEXCEPT { // Check if index is valid:
697 ARTS_ASSERT(0 <= n);
699 return GETFUN(n);
700 }
701
704#undef GETFUN
705
710
711 // Iterators:
712
715
718
719 // Assignment operators:
720
725 VectorView& operator=(const ConstVectorView& v);
726
732 VectorView& operator=(const VectorView& v);
733
736 VectorView& operator=(const Vector& v);
737
738 VectorView& operator=(const Array<Numeric>& v);
739
742 VectorView& operator=(Numeric x);
743
744 // Other operators:
745
747 VectorView operator*=(Numeric x) ARTS_NOEXCEPT;
748
750 VectorView operator/=(Numeric x) ARTS_NOEXCEPT;
751
753 VectorView operator+=(Numeric x) ARTS_NOEXCEPT;
754
756 VectorView operator-=(Numeric x) ARTS_NOEXCEPT;
757
759 VectorView operator*=(const ConstVectorView& x) ARTS_NOEXCEPT;
760
762 VectorView operator/=(const ConstVectorView& x) ARTS_NOEXCEPT;
763
765 VectorView operator+=(const ConstVectorView& x) ARTS_NOEXCEPT;
766
768 VectorView operator-=(const ConstVectorView& x) ARTS_NOEXCEPT;
769
771 operator MatrixView() ARTS_NOEXCEPT;
772
774 ~VectorView() override = default;
775
776 // Friends:
777 friend class ConstIterator2D;
778 friend class Iterator2D;
779 friend class MatrixView;
780 friend class Tensor3View;
781 friend class Tensor4View;
782 friend class Tensor5View;
783 friend class Tensor6View;
784 friend class Tensor7View;
785 friend class ComplexVectorView;
786
790
791 protected:
792 // Constructors:
793 VectorView() = default;
794
797 VectorView(Numeric* data, const Range& range) ARTS_NOEXCEPT;
798
809 VectorView(Numeric* data, const Range& p, const Range& n) ARTS_NOEXCEPT;
810};
811
816 public:
817 // Constructors:
819 Iterator2D() = default;
820
823 : msv(x),
824 mstride(stride) { /* Nothing to do here. */
825 }
826
827 // Operators:
830 msv.mdata += mstride;
831 return *this;
832 }
833
835 bool operator!=(const Iterator2D& other) const ARTS_NOEXCEPT {
836 if (msv.mdata + msv.mrange.mstart !=
837 other.msv.mdata + other.msv.mrange.mstart)
838 return true;
839 return false;
840 }
841
845
848
849 private:
853 Index mstride{0};
854};
855
860 public:
861 // Constructors:
863 ConstIterator2D() = default;
864
867 : msv(std::move(x)),
868 mstride(stride) { /* Nothing to do here. */
869 }
870
871 // Operators:
874 msv.mdata += mstride;
875 return *this;
876 }
877
879 bool operator!=(const ConstIterator2D& other) const ARTS_NOEXCEPT {
880 if (msv.mdata + msv.mrange.mstart !=
881 other.msv.mdata + other.msv.mrange.mstart)
882 return true;
883 return false;
884 }
885
888 const ConstVectorView* operator->() const ARTS_NOEXCEPT { return &msv; }
889
891 const ConstVectorView& operator*() const ARTS_NOEXCEPT { return msv; }
892
893 private:
897 Index mstride{0};
898};
899
910class Vector : public VectorView {
911 public:
912 // Constructors:
913 Vector() = default;
914
916 Vector(std::initializer_list<Numeric> init);
917
919 explicit Vector(const matpack::vector_like_not_vector auto &init)
920 : Vector(matpack::column_size(init)) {
921 *this = init;
922 }
923
925 explicit Vector(Index n);
926
928 Vector(Index n, Numeric fill);
929
937 Vector(Numeric start, Index extent, Numeric stride);
938
944 Vector(const ConstVectorView& v);
945
948 Vector(const Vector& v);
949 Vector(Vector&& v) noexcept : VectorView(std::forward<VectorView>(v)) {
950 v.mdata = nullptr;
951 }
952
954 Vector(const std::vector<Numeric>&);
955
965 ARTS_ASSERT(r0.get_extent() >= 0, "Must have size. Has: ", r0.get_extent());
966 }
967
968 // Assignment operators:
969
994 Vector& operator=(const Vector& v);
995
997 Vector& operator=(Vector&& v) ARTS_NOEXCEPT;
998
1000 Vector& operator=(std::initializer_list<Numeric> v);
1001
1004 if (const auto s = matpack::shape<Index, 1>(init); shape().data not_eq s) resize(s[0]);
1005
1006 for (Index i = 0; i < size(); i++)
1007 operator[](i) = init[i];
1008
1009 return *this;
1010 }
1011
1028 Vector& operator=(const Array<Numeric>& x);
1029
1032 Vector& operator=(Numeric x);
1033
1037 void resize(Index n);
1038
1040 friend void swap(Vector& v1, Vector& v2) noexcept;
1041
1044 ~Vector() noexcept override;
1045
1046 template <class F>
1047 void transform_elementwise(F&& func) {
1048 std::transform(mdata, mdata + size(), mdata, func);
1049 }
1050};
1051
1052// Declare class Matrix:
1053class Matrix;
1054
1066 public:
1067 static constexpr bool matpack_type{true};
1068
1069 constexpr ConstMatrixView(const ConstMatrixView&) = default;
1070 constexpr ConstMatrixView(ConstMatrixView&&) = default;
1073
1074 // Typedef for compatibility with STL
1076
1077 // Member functions:
1078 [[nodiscard]] Index selem() const noexcept { return mrr.mstart + mcr.mstart; }
1079 [[nodiscard]] Index nrows() const noexcept { return mrr.mextent; }
1080 [[nodiscard]] Index ncols() const noexcept { return mcr.mextent; }
1081 [[nodiscard]] Index drows() const noexcept { return mrr.mstride; }
1082 [[nodiscard]] Index dcols() const noexcept { return mcr.mstride; }
1083
1084 // Total size
1085 [[nodiscard]] Index size() const noexcept { return nrows() * ncols(); }
1086 [[nodiscard]] bool empty() const noexcept { return size() == 0; }
1087
1089 [[nodiscard]] Shape<2> shape() const { return {nrows(), ncols()}; }
1090
1091 // Const index operators:
1094 ARTS_NOEXCEPT { // Check if indices are valid:
1095 ARTS_ASSERT(0 <= r);
1096 ARTS_ASSERT(0 <= c);
1097 ARTS_ASSERT(r < mrr.mextent);
1098 ARTS_ASSERT(c < mcr.mextent);
1099
1100 return get(r, c);
1101 }
1102
1104 [[nodiscard]] Numeric get(Index r, Index c) const ARTS_NOEXCEPT {
1105 return *(mdata + mrr.mstart + r * mrr.mstride + mcr.mstart +
1106 c * mcr.mstride);
1107 }
1108
1109 ConstMatrixView operator()(const Range& r,
1110 const Range& c) const ARTS_NOEXCEPT;
1111 ConstVectorView operator()(const Range& r, Index c) const ARTS_NOEXCEPT;
1112 ConstVectorView operator()(Index r, const Range& c) const ARTS_NOEXCEPT;
1113
1114 // Functions returning iterators:
1115 [[nodiscard]] ConstIterator2D begin() const ARTS_NOEXCEPT;
1116 [[nodiscard]] ConstIterator2D end() const ARTS_NOEXCEPT;
1117
1118 // View on diagonal vector
1119 [[nodiscard]] ConstVectorView diagonal() const ARTS_NOEXCEPT;
1120
1122 virtual ~ConstMatrixView() = default;
1123
1124 // Friends:
1125 friend class MatrixView;
1126 friend class ConstIterator3D;
1127 friend class ConstVectorView;
1128 friend class ConstTensor3View;
1129 friend class ConstTensor4View;
1130 friend class ConstTensor5View;
1131 friend class ConstTensor6View;
1132 friend class ConstTensor7View;
1135 friend int poly_root_solve(Matrix& roots, Vector& coeffs);
1136 friend void mult(VectorView, const ConstMatrixView&, const ConstVectorView&);
1137 friend void mult(MatrixView, const ConstMatrixView&, const ConstMatrixView&);
1138 friend void mult(MatrixView, const Sparse&, const ConstMatrixView&);
1139 friend void mult(MatrixView, const ConstMatrixView&, const Sparse&);
1140 friend void mult_general(VectorView,
1141 const ConstMatrixView&,
1143 friend void mult_general(MatrixView,
1144 const ConstMatrixView&,
1146 friend void ludcmp(Matrix&, ArrayOfIndex&, ConstMatrixView);
1147 friend void lubacksub(VectorView,
1150 const ArrayOfIndex&);
1151 friend void inv(MatrixView, ConstMatrixView);
1153
1154 friend std::ostream& operator<<(std::ostream& os, const ConstMatrixView& v);
1155
1156 // Conversion to a plain C-array
1157 [[nodiscard]] Numeric* get_c_array() const noexcept {return mdata;}
1158
1159 protected:
1160 // Constructors:
1161 ConstMatrixView() = default;
1162 ConstMatrixView(Numeric* data, const Range& r, const Range& c) ARTS_NOEXCEPT;
1164 const Range& pr,
1165 const Range& pc,
1166 const Range& nr,
1167 const Range& nc) ARTS_NOEXCEPT;
1168
1169 // Data members:
1170 // -------------
1172 Range mrr{0, 0, 1};
1174 Range mcr{0, 0, 1};
1176 Numeric* mdata{nullptr};
1177};
1178
1189 public:
1190 // Make const methods visible from base class
1193 using ConstMatrixView::operator();
1195
1196 constexpr MatrixView(const MatrixView&) = default;
1197
1198 // Typedef for compatibility with STL
1200
1201#define GETFUN(r, c) \
1202 *(mdata + mrr.mstart + r * mrr.mstride + mcr.mstart + c * mcr.mstride)
1203 // Index Operators:
1206 Index c) ARTS_NOEXCEPT { // Check if indices are valid:
1207 ARTS_ASSERT(0 <= r);
1208 ARTS_ASSERT(0 <= c);
1209 ARTS_ASSERT(r < mrr.mextent);
1210 ARTS_ASSERT(c < mcr.mextent);
1211
1212 return GETFUN(r, c);
1213 }
1214
1217#undef GETFUN
1218
1219 MatrixView operator()(const Range& r, const Range& c) ARTS_NOEXCEPT;
1220 VectorView operator()(const Range& r, Index c) ARTS_NOEXCEPT;
1221 VectorView operator()(Index r, const Range& c) ARTS_NOEXCEPT;
1222
1223 // Functions returning iterators:
1226
1227 // Assignment operators:
1228 MatrixView& operator=(const ConstMatrixView& m);
1229 MatrixView& operator=(const MatrixView& m);
1230 MatrixView& operator=(const Matrix& m);
1231 MatrixView& operator=(const ConstVectorView& m);
1232 MatrixView& operator=(Numeric x);
1233
1234 // Other operators:
1235 MatrixView& operator*=(Numeric x) ARTS_NOEXCEPT;
1236 MatrixView& operator/=(Numeric x) ARTS_NOEXCEPT;
1237 MatrixView& operator+=(Numeric x) ARTS_NOEXCEPT;
1238 MatrixView& operator-=(Numeric x) ARTS_NOEXCEPT;
1239
1240 MatrixView& operator*=(const ConstMatrixView& x) ARTS_NOEXCEPT;
1241 MatrixView& operator/=(const ConstMatrixView& x) ARTS_NOEXCEPT;
1242 MatrixView& operator+=(const ConstMatrixView& x) ARTS_NOEXCEPT;
1243 MatrixView& operator-=(const ConstMatrixView& x) ARTS_NOEXCEPT;
1244
1245 MatrixView& operator*=(const ConstVectorView& x) ARTS_NOEXCEPT;
1246 MatrixView& operator/=(const ConstVectorView& x) ARTS_NOEXCEPT;
1247 MatrixView& operator+=(const ConstVectorView& x) ARTS_NOEXCEPT;
1248 MatrixView& operator-=(const ConstVectorView& x) ARTS_NOEXCEPT;
1249
1251 ~MatrixView() override = default;
1252
1253 // Friends:
1254 friend class VectorView;
1255 friend class Iterator3D;
1256 friend class Tensor3View;
1257 friend class Tensor4View;
1258 friend class Tensor5View;
1259 friend class Tensor6View;
1260 friend class Tensor7View;
1261 friend class ComplexMatrixView;
1264 friend void mult(MatrixView, const ConstMatrixView&, const ConstMatrixView&);
1265
1266 protected:
1267 // Constructors:
1268 MatrixView() = default;
1269 MatrixView(Numeric* data, const Range& rr, const Range& cr) ARTS_NOEXCEPT;
1270 MatrixView(Numeric* data,
1271 const Range& pr,
1272 const Range& pc,
1273 const Range& nr,
1274 const Range& nc) ARTS_NOEXCEPT;
1275};
1276
1285class Matrix : public MatrixView {
1286 public:
1287 // Constructors:
1288 Matrix() = default;
1289 Matrix(Index r, Index c);
1290 Matrix(Index r, Index c, Numeric fill);
1291 Matrix(const ConstMatrixView& m);
1292 Matrix(const Matrix& m);
1293 Matrix(Matrix&& v) noexcept : MatrixView(std::forward<MatrixView>(v)) {
1294 v.mdata = nullptr;
1295 }
1296
1298 explicit Matrix(const matpack::matrix_like_not_matrix auto &init)
1299 : Matrix(matpack::row_size(init), matpack::column_size(init)) {
1300 *this = init;
1301 }
1302
1312 Matrix(Numeric* d, const Range& r0, const Range& r1) ARTS_NOEXCEPT
1313 : MatrixView(d, r0, r1) {
1314 ARTS_ASSERT(r0.get_extent() >= 0, "Must have size. Has: ", r0.get_extent());
1315 ARTS_ASSERT(r1.get_extent() >= 0, "Must have size. Has: ", r1.get_extent());
1316 }
1317
1318 // Assignment operators:
1319 Matrix& operator=(const Matrix& m);
1323
1326 if (const auto s = matpack::shape<Index, 2>(init);
1327 shape().data not_eq s)
1328 resize(s[0], s[1]);
1329
1330 auto [I, J] = shape().data;
1331 for (Index i = 0; i < I; i++)
1332 for (Index j = 0; j < J; j++)
1333 operator()(i, j) = init(i, j);
1334
1335 return *this;
1336 }
1337
1338 // Resize function:
1339 void resize(Index r, Index c);
1340
1341 // Swap function:
1342 friend void swap(Matrix& m1, Matrix& m2) noexcept;
1343
1344 // Destructor:
1345 ~Matrix() noexcept override;
1346
1348 template <std::size_t dim0>
1349 Vector reduce_rank() && ARTS_NOEXCEPT {
1350 static_assert(dim0 < 2, "Bad Dimension, Out-of-Bounds");
1351
1352 Range r0(0, dim0 == 0 ? nrows() : ncols());
1353
1354 Vector out(mdata, r0);
1355 ARTS_ASSERT(size() == out.size(),
1356 "Can only reduce size on same size input. "
1357 "The sizes are (input v. output): ",
1358 size(),
1359 " v. ",
1360 out.size());
1361 mdata = nullptr;
1362 return out;
1363 }
1364
1366
1367 template <class F>
1368 void transform_elementwise(F&& func) {
1369 std::transform(mdata, mdata + size(), mdata, func);
1370 }
1371};
1372
1373// Function declarations:
1374// ----------------------
1375
1376void copy(ConstIterator1D origin,
1377 const ConstIterator1D& end,
1378 Iterator1D target);
1379
1381void copy(Numeric x, Iterator1D target, const Iterator1D& end) ARTS_NOEXCEPT;
1382
1383void copy(ConstIterator2D origin,
1384 const ConstIterator2D& end,
1385 Iterator2D target);
1386
1387void copy(Numeric x, Iterator2D target, const Iterator2D& end) ARTS_NOEXCEPT;
1388
1389void mult(VectorView y, const ConstMatrixView& M, const ConstVectorView& x);
1390
1392 const ConstMatrixView& B,
1394
1395void mult(MatrixView A, const ConstMatrixView& B, const ConstMatrixView& C);
1396
1398 const ConstMatrixView& B,
1400
1401void cross3(VectorView c,
1402 const ConstVectorView& a,
1404
1406
1408
1410
1412
1413void transform(VectorView y, double (&my_func)(double), ConstVectorView x);
1414
1415void transform(MatrixView y, double (&my_func)(double), ConstMatrixView x);
1416
1418
1420
1422
1424
1426
1428
1430
1433
1435// Helper function for debugging
1436#ifndef NDEBUG
1437
1439
1440#endif
1442
1445
1447
1450
1452
1453#endif // matpackI_h
This file contains the definition of Array.
base max(const Array< base > &x)
Max function.
Definition: array.h:145
base min(const Array< base > &x)
Min function.
Definition: array.h:161
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:458
bool operator!=(const ConstIterator1D &other) const ARTS_NOEXCEPT
Not equal operator, needed for algorithms like copy.
Definition: matpackI.h:486
const Numeric * pointer
Definition: matpackI.h:462
std::random_access_iterator_tag iterator_category
Definition: matpackI.h:464
const Numeric * mx
Current position.
Definition: matpackI.h:505
Index mstride
Stride.
Definition: matpackI.h:507
ConstIterator1D()=default
Default constructor.
const Numeric value_type
Definition: matpackI.h:461
const Numeric & reference
Definition: matpackI.h:463
bool operator==(const ConstIterator1D &other) const ARTS_NOEXCEPT
Definition: matpackI.h:491
ConstIterator1D & operator++() ARTS_NOEXCEPT
Prefix increment operator.
Definition: matpackI.h:477
Index difference_type
Definition: matpackI.h:460
ConstIterator1D(const Numeric *x, Index stride) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:470
const Numeric & operator*() const ARTS_NOEXCEPT
Dereferencing.
Definition: matpackI.h:483
Index operator-(const ConstIterator1D &other) const ARTS_NOEXCEPT
Definition: matpackI.h:495
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:859
const ConstVectorView & operator*() const ARTS_NOEXCEPT
Dereferencing.
Definition: matpackI.h:891
ConstIterator2D(ConstVectorView x, Index stride) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:866
ConstVectorView msv
Current position.
Definition: matpackI.h:895
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:888
bool operator!=(const ConstIterator2D &other) const ARTS_NOEXCEPT
Not equal operator, needed for algorithms like copy.
Definition: matpackI.h:879
ConstIterator2D & operator++() ARTS_NOEXCEPT
Prefix increment operator.
Definition: matpackI.h:873
Const version of Iterator3D.
Definition: matpackIII.h:78
Const version of Iterator4D.
Definition: matpackIV.h:78
Const version of Iterator5D.
Definition: matpackV.h:87
Const version of Iterator6D.
Definition: matpackVI.h:89
Const version of Iterator7D.
Definition: matpackVII.h:86
A constant view of a Matrix.
Definition: matpackI.h:1065
Index drows() const noexcept
Definition: matpackI.h:1081
constexpr ConstMatrixView(ConstMatrixView &&)=default
ConstMatrixView & operator=(const ConstMatrixView &)=default
bool empty() const noexcept
Definition: matpackI.h:1086
Index dcols() const noexcept
Definition: matpackI.h:1082
ConstIterator2D begin() const ARTS_NOEXCEPT
Return const iterator to first row.
Definition: matpackI.cc:449
Index size() const noexcept
Definition: matpackI.h:1085
Index nrows() const noexcept
Definition: matpackI.h:1079
constexpr ConstMatrixView(const ConstMatrixView &)=default
Shape< 2 > shape() const
Definition: matpackI.h:1089
ConstIterator2D end() const ARTS_NOEXCEPT
Return const iterator behind last row.
Definition: matpackI.cc:454
Index ncols() const noexcept
Definition: matpackI.h:1080
Numeric operator()(Index r, Index c) const ARTS_NOEXCEPT
Plain const index operator.
Definition: matpackI.h:1093
ConstMatrixView & operator=(ConstMatrixView &&)=default
ConstMatrixView()=default
Index selem() const noexcept
Definition: matpackI.h:1078
Numeric get(Index r, Index c) const ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:1104
A constant view of a Tensor3.
Definition: matpackIII.h:133
A constant view of a Tensor4.
Definition: matpackIV.h:133
A constant view of a Tensor5.
Definition: matpackV.h:144
A constant view of a Tensor6.
Definition: matpackVI.h:150
A constant view of a Tensor7.
Definition: matpackVII.h:148
A constant view of a Vector.
Definition: matpackI.h:521
Numeric get(Index n) const ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:567
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
ConstVectorView()=default
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:660
Index delem() const noexcept
Steps in memory between elements.
Definition: matpackI.h:626
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:1138
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:560
constexpr ConstVectorView(ConstVectorView &&)=default
friend class ConstIterator2D
Definition: matpackI.h:595
constexpr ConstVectorView(const ConstVectorView &)=default
friend void mult(VectorView, const ConstMatrixView &, const ConstVectorView &)
Matrix-Vector Multiplication.
Definition: matpackI.cc:1080
static constexpr bool matpack_type
Definition: matpackI.h:523
Index selem() const noexcept
Start element in memory.
Definition: matpackI.h:623
Shape< 1 > shape() const
Definition: matpackI.h:553
Range mrange
The range of mdata that is actually used.
Definition: matpackI.h:658
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:633
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
friend class ConstMatrixView
Definition: matpackI.h:596
friend void diagonalize(MatrixView, VectorView, VectorView, ConstMatrixView)
Matrix Diagonalization.
Definition: lin_alg.cc:240
Index size() const noexcept
Definition: matpackI.h:550
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:1048
bool empty() const noexcept
Returns true if variable size is zero.
Definition: matpackI.h:536
The iterator class for sub vectors.
Definition: matpackI.h:399
Iterator1D & operator++() ARTS_NOEXCEPT
Prefix increment operator.
Definition: matpackI.h:419
bool operator==(const Iterator1D &other) const ARTS_NOEXCEPT
Definition: matpackI.h:433
Index operator-(const Iterator1D &other) const ARTS_NOEXCEPT
Definition: matpackI.h:437
Iterator1D()=default
Default constructor.
Index difference_type
Definition: matpackI.h:401
Numeric & operator*() const ARTS_NOEXCEPT
Dereferencing.
Definition: matpackI.h:425
Index mstride
Stride.
Definition: matpackI.h:453
std::random_access_iterator_tag iterator_category
Definition: matpackI.h:405
Numeric * mx
Current position.
Definition: matpackI.h:451
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:403
bool operator!=(const Iterator1D &other) const ARTS_NOEXCEPT
Not equal operator, needed for algorithms like copy.
Definition: matpackI.h:428
Numeric & reference
Definition: matpackI.h:404
Numeric value_type
Definition: matpackI.h:402
Iterator1D(Numeric *x, Index stride) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:411
The row iterator class for sub matrices.
Definition: matpackI.h:815
Iterator2D(const VectorView &x, Index stride) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:822
VectorView & operator*() ARTS_NOEXCEPT
Dereferencing.
Definition: matpackI.h:847
bool operator!=(const Iterator2D &other) const ARTS_NOEXCEPT
Not equal operator, needed for algorithms like copy.
Definition: matpackI.h:835
Iterator2D()=default
Default constructor.
VectorView msv
Current position.
Definition: matpackI.h:851
Iterator2D & operator++() ARTS_NOEXCEPT
Prefix increment operator.
Definition: matpackI.h:829
VectorView * operator->() ARTS_NOEXCEPT
The -> operator is needed, so that we can write i->begin() to get the 1D iterators.
Definition: matpackI.h:844
Implementation of Tensors of Rank 3.
Definition: matpackIII.h:37
Implementation of Tensors of Rank 4.
Definition: matpackIV.h:40
Implementation of Tensors of Rank 5.
Definition: matpackV.h:41
The outermost iterator class for rank 6 tensors.
Definition: matpackVI.h:43
Implementation of Tensors of Rank 7.
Definition: matpackVII.h:39
The Joker class.
Definition: matpackI.h:126
The MatrixView class.
Definition: matpackI.h:1188
Numeric & operator()(Index r, Index c) ARTS_NOEXCEPT
Plain index operator.
Definition: matpackI.h:1205
Numeric & get(Index r, Index c) ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:1216
constexpr MatrixView(const MatrixView &)=default
The Matrix class.
Definition: matpackI.h:1285
Matrix(Matrix &&v) noexcept
Definition: matpackI.h:1293
Matrix()=default
Matrix & operator=(const matpack::matrix_like_not_matrix auto &init)
Set from a matrix type.
Definition: matpackI.h:1325
Numeric * get_raw_data() ARTS_NOEXCEPT
Definition: matpackI.h:1365
Matrix(const matpack::matrix_like_not_matrix auto &init)
Initialization from a vector type.
Definition: matpackI.h:1298
Matrix(Numeric *d, const Range &r0, const Range &r1) ARTS_NOEXCEPT
Definition: matpackI.h:1312
void transform_elementwise(F &&func)
Definition: matpackI.h:1368
The range class.
Definition: matpackI.h:160
constexpr void swap(Range &other) noexcept
Definition: matpackI.h:366
constexpr Range(Index start, Joker, Index stride=1) ARTS_NOEXCEPT
Constructor with joker extent.
Definition: matpackI.h:191
friend void mult_general(VectorView, const ConstMatrixView &, const ConstVectorView &) ARTS_NOEXCEPT
Matrix Vector multiplication.
Definition: matpackI.cc:1138
constexpr Index operator()(const Index i) const ARTS_NOEXCEPT
Definition: matpackI.h:360
friend constexpr void swap(Range &a, Range &b) noexcept
Definition: matpackI.h:373
constexpr Range(const Range &p, const Range &n) ARTS_NOEXCEPT
Constructor of a new range relative to an old range.
Definition: matpackI.h:247
constexpr Range(Index start, Index extent, Index stride=1) ARTS_NOEXCEPT
Explicit constructor.
Definition: matpackI.h:174
Index mstart
The start index.
Definition: matpackI.h:377
constexpr Index get_extent() const noexcept
Returns the extent of the range.
Definition: matpackI.h:343
constexpr Range operator()(const Range r) const ARTS_NOEXCEPT
Range of range.
Definition: matpackI.h:348
friend std::ostream & operator<<(std::ostream &os, const Range &r)
Definition: matpackI.cc:39
Index mstride
The stride.
Definition: matpackI.h:382
constexpr Index get_start() const noexcept
Returns the start index of the range.
Definition: matpackI.h:341
Index mextent
The number of elements.
Definition: matpackI.h:380
constexpr Range(Index max_size, const Range &r) ARTS_NOEXCEPT
Constructor which converts a range with joker to an explicit range.
Definition: matpackI.h:214
constexpr Index get_stride() const noexcept
Returns the stride of the range.
Definition: matpackI.h:345
constexpr Range(Joker, Index stride=1) noexcept
Constructor with just a joker.
Definition: matpackI.h:202
The Tensor3View class.
Definition: matpackIII.h:251
The Tensor3 class.
Definition: matpackIII.h:352
The Tensor4View class.
Definition: matpackIV.h:296
The Tensor4 class.
Definition: matpackIV.h:435
The Tensor5View class.
Definition: matpackV.h:348
The Tensor5 class.
Definition: matpackV.h:524
The Tensor6View class.
Definition: matpackVI.h:635
The Tensor6 class.
Definition: matpackVI.h:1105
The Tensor7View class.
Definition: matpackVII.h:1308
The Tensor7 class.
Definition: matpackVII.h:2407
The VectorView class.
Definition: matpackI.h:674
Numeric & get(Index n) ARTS_NOEXCEPT
Get element implementation without assertions.
Definition: matpackI.h:703
VectorView()=default
friend class MatrixView
Definition: matpackI.h:779
constexpr VectorView(const VectorView &)=default
VectorView & operator=(const ConstVectorView &v)
Assignment operator.
Definition: matpackI.cc:153
friend class Iterator2D
Definition: matpackI.h:778
Numeric & operator[](Index n) ARTS_NOEXCEPT
Plain Index operator.
Definition: matpackI.h:696
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:910
Vector(Numeric *d, const Range &r0) ARTS_NOEXCEPT
Definition: matpackI.h:964
Vector(Vector &&v) noexcept
Definition: matpackI.h:949
Vector & operator=(const matpack::vector_like_not_vector auto &init)
Set from a vector type.
Definition: matpackI.h:1003
Vector()=default
Vector(const matpack::vector_like_not_vector auto &init)
Initialization from a vector type.
Definition: matpackI.h:919
Binary output file stream class.
Definition: bofstream.h:42
A concept precluding Arts matrix objects but allowing things similar to matrices.
A concept precluding Arts vector objects but allowing things similar to vectors.
#define ARTS_NOEXCEPT
Definition: debug.h:99
#define ARTS_ASSERT(condition,...)
Definition: debug.h:102
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:167
Numeric mean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version.
Definition: matpackI.cc:1546
Numeric vector_angle(ConstVectorView a, ConstVectorView b)
Definition: matpackI.cc:1371
Numeric nanmean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version ignoring nans and infs.
Definition: matpackI.cc:1561
void proj(Vector &c, ConstVectorView a, ConstVectorView b) ARTS_NOEXCEPT
Definition: matpackI.cc:1393
const Joker joker
#define GETFUN(n)
Definition: matpackI.h:1201
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:1433
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
cross3
Definition: matpackI.cc:1350
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:1648
ConstMatrixView transpose(ConstMatrixView m) ARTS_NOEXCEPT
Const version of transpose.
Definition: matpackI.cc:1403
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) noexcept
Swaps two objects.
#define N
Definition: rng.cc:164
#define M
Definition: rng.cc:165
Helper shape class.
Definition: matpackI.h:387
std::array< Index, N > data
Definition: matpackI.h:388
friend std::ostream & operator<<(std::ostream &os, const Shape &shape)
Definition: matpackI.h:390
bool operator<=>(const Shape &other) const =default
The Sparse class.
Definition: matpackII.h:75
#define d
#define v
#define a
#define c
#define b