ARTS 2.5.0 (git: 9ee3ac6c)
matpackI.cc
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
24#include "matpackI.h"
25#include <cmath>
26#include <cstring>
27#include "blas.h"
28#include "exceptions.h"
29
30using std::setw;
31
32// Define the global joker object:
33extern const Joker joker = Joker();
34
35static const Numeric RAD2DEG = 57.295779513082323;
36
37std::ostream& operator<<(std::ostream& os, const Range& r) {
38 os << "Range(" << r.get_start() << ", " << r.get_extent() << ", "
39 << r.get_stride() << ")";
40 return os;
41}
42
43// Functions for ConstVectorView:
44// ------------------------------
45
46bool ConstVectorView::empty() const ARTS_NOEXCEPT { return (nelem() == 0); }
47
49
51
53 Numeric s = 0;
55 const ConstIterator1D e = end();
56
57 for (; i != e; ++i) s += *i;
58
59 return s;
60}
61
63 return ConstVectorView(mdata, mrange, r);
64}
65
68}
69
71 return ConstIterator1D(
74}
75
76ConstVectorView::operator ConstMatrixView() const {
77 // The old version (before 2013-01-18) of this was:
78 // return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
79 // Bus this was a bug! The problem is that the matrix index operator adds
80 // the mstart from both row and columm range object to mdata
81
82 return ConstMatrixView(mdata, mrange, Range(0, 1));
83}
84
85/* This one is a bit tricky: We have to cast away the arguments const
86 qualifier, because mdata is not const. This should be safe, since
87 there are no non-const methods for ConstVectorView. */
89 : mrange(0, 1), mdata(&const_cast<Numeric&>(a)) {
90 // Nothing to do here.
91}
92
94 : mrange(range), mdata(data) {
95 // Nothing to do here.
96}
97
99 : mrange(p, n), mdata(data) {
100 // Nothing to do here.
101}
102
103/* Output operator. This demonstrates how iterators can be used to
104 traverse the vector. The iterators know which part of the vector
105 is `active', and also the stride. */
106std::ostream& operator<<(std::ostream& os, const ConstVectorView& v) {
107 ConstIterator1D i = v.begin();
108 const ConstIterator1D end = v.end();
109
110 if (i != end) {
111 os << setw(3) << *i;
112 ++i;
113 }
114 for (; i != end; ++i) {
115 os << " " << setw(3) << *i;
116 }
117
118 return os;
119}
120
121// Functions for VectorView:
122// ------------------------
123
125 ARTS_ASSERT (false,
126 "Creating a VectorView from a const Vector is not allowed.\n"
127 "This is not really a runtime error, but I don't want to start\n"
128 "producing direct output from inside matpack. And just exiting is\n"
129 "not so nice.\n"
130 "If you see this error, there is a bug in the code, not in the\n"
131 "ARTS input.")
132}
133
135 mdata = v.mdata;
136 mrange = v.mrange;
137}
138
140 return VectorView(mdata, mrange, r);
141}
142
145}
146
150}
151
153 // cout << "Assigning VectorView from ConstVectorView.\n";
154
155 // Check that sizes are compatible:
156 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
157
158 copy(v.begin(), v.end(), begin());
159
160 return *this;
161}
162
164 // cout << "Assigning VectorView from VectorView.\n";
165
166 // Check that sizes are compatible:
167 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
168
169 copy(v.begin(), v.end(), begin());
170
171 return *this;
172}
173
175 // cout << "Assigning VectorView from Vector.\n";
176
177 // Check that sizes are compatible:
178 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
179
180 copy(v.begin(), v.end(), begin());
181
182 return *this;
183}
184
186 copy(x, begin(), end());
187 return *this;
188}
189
191 const Iterator1D e = end();
192 for (Iterator1D i = begin(); i != e; ++i) *i *= x;
193 return *this;
194}
195
197 const Iterator1D e = end();
198 for (Iterator1D i = begin(); i != e; ++i) *i /= x;
199 return *this;
200}
201
203 const Iterator1D e = end();
204 for (Iterator1D i = begin(); i != e; ++i) *i += x;
205 return *this;
206}
207
209 const Iterator1D e = end();
210 for (Iterator1D i = begin(); i != e; ++i) *i -= x;
211 return *this;
212}
213
215 ARTS_ASSERT(nelem() == x.nelem());
216
217 ConstIterator1D s = x.begin();
218
219 Iterator1D i = begin();
220 const Iterator1D e = end();
221
222 for (; i != e; ++i, ++s) *i *= *s;
223 return *this;
224}
225
227 ARTS_ASSERT(nelem() == x.nelem());
228
229 ConstIterator1D s = x.begin();
230
231 Iterator1D i = begin();
232 const Iterator1D e = end();
233
234 for (; i != e; ++i, ++s) *i /= *s;
235 return *this;
236}
237
239 ARTS_ASSERT(nelem() == x.nelem());
240
241 ConstIterator1D s = x.begin();
242
243 Iterator1D i = begin();
244 const Iterator1D e = end();
245
246 for (; i != e; ++i, ++s) *i += *s;
247 return *this;
248}
249
251 ARTS_ASSERT(nelem() == x.nelem());
252
253 ConstIterator1D s = x.begin();
254
255 Iterator1D i = begin();
256 const Iterator1D e = end();
257
258 for (; i != e; ++i, ++s) *i -= *s;
259 return *this;
260}
261
262VectorView::operator MatrixView() ARTS_NOEXCEPT {
263 // The old version (before 2013-01-18) of this was:
264 // return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
265 // Bus this was a bug! The problem is that the matrix index operator adds
266 // the mstart from both row and columm range object to mdata
267
268 return MatrixView(mdata, mrange, Range(0, 1));
269}
270
272 ARTS_ASSERT(not (mrange.mstart != 0 || mrange.mstride != 1),
273 "A VectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");
274
275 return mdata;
276}
277
279 ARTS_ASSERT(not (mrange.mstart != 0 || mrange.mstride != 1),
280 "A VectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");
281
282 return mdata;
283}
284
286 // Nothing to do here.
287}
288
290 : ConstVectorView(data, range) {
291 // Nothing to do here.
292}
293
295 : ConstVectorView(data, p, n) {
296 // Nothing to do here.
297}
298
300 const ConstIterator1D& end,
301 Iterator1D target) {
302 if (origin.mstride == 1 && target.mstride == 1)
303 memcpy((void*)target.mx,
304 (void*)origin.mx,
305 sizeof(Numeric) * (end.mx - origin.mx));
306 else
307 for (; origin != end; ++origin, ++target) *target = *origin;
308}
309
311 for (; target != end; ++target) *target = x;
312}
313
314// Functions for Vector:
315// ---------------------
316
317Vector::Vector(std::initializer_list<Numeric> init)
318 : VectorView(new Numeric[init.size()], Range(0, init.size())) {
319 std::copy(init.begin(), init.end(), begin());
320}
321
322Vector::Vector(const Eigen::VectorXd& init)
323 : VectorView(new Numeric[init.size()], Range(0, init.size()))
324{
325 for (Index i=0; i<size(); i++) operator[](i) = init[i];
326}
327
328
330 // Nothing to do here.
331}
332
334 : VectorView(new Numeric[n], Range(0, n)) {
335 // Here we can access the raw memory directly, for slightly
336 // increased efficiency:
337 std::fill_n(mdata, n, fill);
338}
339
341 : VectorView(new Numeric[extent], Range(0, extent)) {
342 // Fill with values:
343 Numeric x = start;
344 Iterator1D i = begin();
345 const Iterator1D e = end();
346 for (; i != e; ++i) {
347 *i = x;
348 x += stride;
349 }
350}
351
353 : VectorView(new Numeric[v.nelem()], Range(0, v.nelem())) {
354 copy(v.begin(), v.end(), begin());
355}
356
358 : VectorView(new Numeric[v.nelem()], Range(0, v.nelem())) {
359 std::memcpy(mdata, v.mdata, nelem() * sizeof(Numeric));
360}
361
362Vector::Vector(const std::vector<Numeric>& v)
363 : VectorView(new Numeric[v.size()], Range(0, v.size())) {
364 std::vector<Numeric>::const_iterator vec_it_end = v.end();
365 Iterator1D this_it = this->begin();
366 for (std::vector<Numeric>::const_iterator vec_it = v.begin();
367 vec_it != vec_it_end;
368 ++vec_it, ++this_it)
369 *this_it = *vec_it;
370}
371
372Vector& Vector::operator=(std::initializer_list<Numeric> v) {
373 resize(v.size());
374 std::copy(v.begin(), v.end(), begin());
375 return *this;
376}
377
379 if (this != &v) {
380 resize(v.nelem());
381 std::memcpy(mdata, v.mdata, nelem() * sizeof(Numeric));
382 }
383 return *this;
384}
385
387 if (this != &v) {
388 delete[] mdata;
389 mdata = v.mdata;
390 mrange = v.mrange;
391 v.mrange = Range(0, 0);
392 v.mdata = nullptr;
393 }
394 return *this;
395}
396
398 resize(x.nelem());
400 return *this;
401}
402
404 std::fill_n(mdata, nelem(), x);
405 return *this;
406}
407
409 ARTS_ASSERT(0 <= n);
410 if (mrange.mextent != n) {
411 delete[] mdata;
412 mdata = new Numeric[n];
413 mrange.mstart = 0;
414 mrange.mextent = n;
415 mrange.mstride = 1;
416 }
417}
418
419void swap(Vector& v1, Vector& v2) {
420 std::swap(v1.mrange, v2.mrange);
421 std::swap(v1.mdata, v2.mdata);
422}
423
424Vector::~Vector() { delete[] mdata; }
425
426// Functions for ConstMatrixView:
427// ------------------------------
428
430bool ConstMatrixView::empty() const ARTS_NOEXCEPT { return (nrows() == 0 || ncols() == 0); }
431
434
437
442 const Range& c) const ARTS_NOEXCEPT {
443 return ConstMatrixView(mdata, mrr, mcr, r, c);
444}
445
452 // Check that c is valid:
453 ARTS_ASSERT(0 <= c);
454 ARTS_ASSERT(c < mcr.mextent);
455
456 return ConstVectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
457}
458
465 // Check that r is valid:
466 ARTS_ASSERT(0 <= r);
467 ARTS_ASSERT(r < mrr.mextent);
468
469 return ConstVectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
470}
471
475}
476
479 return ConstIterator2D(
481 mrr.mstride);
482}
483
485
496 Range(0, n, mrr.mstride + mcr.mstride));
497}
498
503 const Range& rr,
504 const Range& cr) ARTS_NOEXCEPT
505 : mrr(rr), mcr(cr), mdata(data) {
506 // Nothing to do here.
507}
508
524 const Range& pr,
525 const Range& pc,
526 const Range& nr,
527 const Range& nc) ARTS_NOEXCEPT
528 : mrr(pr, nr), mcr(pc, nc), mdata(data) {
529 // Nothing to do here.
530}
531
539std::ostream& operator<<(std::ostream& os, const ConstMatrixView& v) {
540 // Row iterators:
541 ConstIterator2D ir = v.begin();
542 const ConstIterator2D end_row = v.end();
543
544 if (ir != end_row) {
545 ConstIterator1D ic = ir->begin();
546 const ConstIterator1D end_col = ir->end();
548 if (ic != end_col) {
549 os << setw(3) << *ic;
550 ++ic;
551 }
552 for (; ic != end_col; ++ic) {
553 os << " " << setw(3) << *ic;
554 }
555 ++ir;
556 }
557 for (; ir != end_row; ++ir) {
558 ConstIterator1D ic = ir->begin();
559 const ConstIterator1D end_col = ir->end();
560
561 os << "\n";
562 if (ic != end_col) {
563 os << setw(3) << *ic;
564 ++ic;
565 }
566 for (; ic != end_col; ++ic) {
567 os << " " << setw(3) << *ic;
568 }
569 }
570
571 return os;
572}
573
574// Functions for MatrixView:
575// -------------------------
576
581 return MatrixView(mdata, mrr, mcr, r, c);
582}
583
589 // Check that c is valid:
590 ARTS_ASSERT(0 <= c);
591 ARTS_ASSERT(c < mcr.mextent);
592
593 return VectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
594}
595
601 // Check that r is valid:
602 ARTS_ASSERT(0 <= r);
603 ARTS_ASSERT(r < mrr.mextent);
604
605 return VectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
606}
607
609// hiden by the non-const operator of the derived class.*/
610//ConstIterator2D MatrixView::begin() const
611//{
612// return ConstMatrixView::begin();
613//}
614//
616//ConstIterator2D MatrixView::end() const
617//{
618// return ConstMatrixView::end();
619//}
620
624}
625
628 return Iterator2D(
630 mrr.mstride);
631}
632
638 // Check that sizes are compatible:
641
642 copy(m.begin(), m.end(), begin());
643 return *this;
644}
645
652 // Check that sizes are compatible:
655
656 copy(m.begin(), m.end(), begin());
657 return *this;
658}
659
664 // Check that sizes are compatible:
667
668 copy(m.begin(), m.end(), begin());
669 return *this;
670}
671
677 // Check that sizes are compatible:
678 ARTS_ASSERT(mrr.mextent == v.nelem());
679 ARTS_ASSERT(mcr.mextent == 1);
680 // dummy = ConstMatrixView(v.mdata,v.mrange,Range(v.mrange.mstart,1));;
681 ConstMatrixView dummy(v);
682 copy(dummy.begin(), dummy.end(), begin());
683 return *this;
684}
685
689 copy(x, begin(), end());
690 return *this;
691}
692
695 const Iterator2D er = end();
696 for (Iterator2D r = begin(); r != er; ++r) {
697 const Iterator1D ec = r->end();
698 for (Iterator1D c = r->begin(); c != ec; ++c) *c *= x;
699 }
700 return *this;
701}
702
705 const Iterator2D er = end();
706 for (Iterator2D r = begin(); r != er; ++r) {
707 const Iterator1D ec = r->end();
708 for (Iterator1D c = r->begin(); c != ec; ++c) *c /= x;
709 }
710 return *this;
711}
712
715 const Iterator2D er = end();
716 for (Iterator2D r = begin(); r != er; ++r) {
717 const Iterator1D ec = r->end();
718 for (Iterator1D c = r->begin(); c != ec; ++c) *c += x;
719 }
720 return *this;
721}
722
725 const Iterator2D er = end();
726 for (Iterator2D r = begin(); r != er; ++r) {
727 const Iterator1D ec = r->end();
728 for (Iterator1D c = r->begin(); c != ec; ++c) *c -= x;
729 }
730 return *this;
731}
732
740 ARTS_ASSERT(not (mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 || mcr.mstride != 1),
741 "A MatrixView can only be converted to a plain C-array if it's pointing to a continuous block of data");
742
743 return mdata;
744}
745
753 ARTS_ASSERT(not (mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 || mcr.mstride != 1),
754 "A MatrixView can only be converted to a plain C-array if it's pointing to a continuous block of data");
755
756 return mdata;
757}
758
761 ARTS_ASSERT(nrows() == x.nrows());
762 ARTS_ASSERT(ncols() == x.ncols());
763 ConstIterator2D sr = x.begin();
764 Iterator2D r = begin();
765 const Iterator2D er = end();
766 for (; r != er; ++r, ++sr) {
767 ConstIterator1D sc = sr->begin();
768 Iterator1D c = r->begin();
769 const Iterator1D ec = r->end();
770 for (; c != ec; ++c, ++sc) *c *= *sc;
771 }
772 return *this;
773}
774
777 ARTS_ASSERT(nrows() == x.nrows());
778 ARTS_ASSERT(ncols() == x.ncols());
779 ConstIterator2D sr = x.begin();
780 Iterator2D r = begin();
781 const Iterator2D er = end();
782 for (; r != er; ++r, ++sr) {
783 ConstIterator1D sc = sr->begin();
784 Iterator1D c = r->begin();
785 const Iterator1D ec = r->end();
786 for (; c != ec; ++c, ++sc) *c /= *sc;
787 }
788 return *this;
789}
790
793 ARTS_ASSERT(nrows() == x.nrows());
794 ARTS_ASSERT(ncols() == x.ncols());
795 ConstIterator2D sr = x.begin();
796 Iterator2D r = begin();
797 const Iterator2D er = end();
798 for (; r != er; ++r, ++sr) {
799 ConstIterator1D sc = sr->begin();
800 Iterator1D c = r->begin();
801 const Iterator1D ec = r->end();
802 for (; c != ec; ++c, ++sc) *c += *sc;
803 }
804 return *this;
805}
806
809 ARTS_ASSERT(nrows() == x.nrows());
810 ARTS_ASSERT(ncols() == x.ncols());
811 ConstIterator2D sr = x.begin();
812 Iterator2D r = begin();
813 const Iterator2D er = end();
814 for (; r != er; ++r, ++sr) {
815 ConstIterator1D sc = sr->begin();
816 Iterator1D c = r->begin();
817 const Iterator1D ec = r->end();
818 for (; c != ec; ++c, ++sc) *c -= *sc;
819 }
820 return *this;
821}
822
825 ARTS_ASSERT(nrows() == x.nelem());
826 ARTS_ASSERT(ncols() == 1);
827 ConstIterator1D sc = x.begin();
828 Iterator2D r = begin();
829 const Iterator2D er = end();
830 for (; r != er; ++r, ++sc) {
831 Iterator1D c = r->begin();
832 *c *= *sc;
833 }
834 return *this;
835}
836
839 ARTS_ASSERT(nrows() == x.nelem());
840 ARTS_ASSERT(ncols() == 1);
841 ConstIterator1D sc = x.begin();
842 Iterator2D r = begin();
843 const Iterator2D er = end();
844 for (; r != er; ++r, ++sc) {
845 Iterator1D c = r->begin();
846 *c /= *sc;
847 }
848 return *this;
849}
850
853 ARTS_ASSERT(nrows() == x.nelem());
854 ARTS_ASSERT(ncols() == 1);
855 ConstIterator1D sc = x.begin();
856 Iterator2D r = begin();
857 const Iterator2D er = end();
858 for (; r != er; ++r, ++sc) {
859 Iterator1D c = r->begin();
860 *c += *sc;
861 }
862 return *this;
863}
864
867 ARTS_ASSERT(nrows() == x.nelem());
868 ARTS_ASSERT(ncols() == 1);
869 ConstIterator1D sc = x.begin();
870 Iterator2D r = begin();
871 const Iterator2D er = end();
872 for (; r != er; ++r, ++sc) {
873 Iterator1D c = r->begin();
874 *c -= *sc;
875 }
876 return *this;
877}
878
883 : ConstMatrixView(data, rr, cr) {
884 // Nothing to do here.
885}
886
902 const Range& pr,
903 const Range& pc,
904 const Range& nr,
905 const Range& nc) ARTS_NOEXCEPT
906 : ConstMatrixView(data, pr, pc, nr, nc) {
907 // Nothing to do here.
908}
909
920 const ConstIterator2D& end,
921 Iterator2D target) {
922 for (; origin != end; ++origin, ++target) {
923 ConstIterator1D o = origin->begin();
924 const ConstIterator1D e = origin->end();
925 Iterator1D t = target->begin();
926 for (; o != e; ++o, ++t) *t = *o;
927 }
928}
929
932 for (; target != end; ++target) {
933 Iterator1D t = target->begin();
934 const Iterator1D e = target->end();
935 for (; t != e; ++t) *t = x;
936 }
937}
938
939// Functions for Matrix:
940// ---------------------
941
945 : MatrixView(new Numeric[r * c], Range(0, r, c), Range(0, c)) {
946 // Nothing to do here.
947}
948
951 : MatrixView(new Numeric[r * c], Range(0, r, c), Range(0, c)) {
952 // Here we can access the raw memory directly, for slightly
953 // increased efficiency:
954 std::fill_n(mdata, r * c, fill);
955}
956
960 : MatrixView(new Numeric[m.nrows() * m.ncols()],
961 Range(0, m.nrows(), m.ncols()),
962 Range(0, m.ncols())) {
963 copy(m.begin(), m.end(), begin());
964}
965
969 : MatrixView(new Numeric[m.nrows() * m.ncols()],
970 Range(0, m.nrows(), m.ncols()),
971 Range(0, m.ncols())) {
972 // There is a catch here: If m is an empty matrix, then it will have
973 // 0 colunns. But this is used to initialize the stride of the row
974 // Range! Thus, this method has to be consistent with the behaviour
975 // of Range::Range. For now, Range::Range allows also stride 0.
976 std::memcpy(mdata, m.mdata, nrows() * ncols() * sizeof(Numeric));
977}
978
980
1005 if (this != &m) {
1006 resize(m.nrows(), m.ncols());
1007 std::memcpy(mdata, m.mdata, nrows() * ncols() * sizeof(Numeric));
1008 }
1009 return *this;
1010}
1011
1014 if (this != &m) {
1015 delete[] mdata;
1016 mdata = m.mdata;
1017 mrr = m.mrr;
1018 mcr = m.mcr;
1019 m.mrr = Range(0, 0);
1020 m.mcr = Range(0, 0);
1021 m.mdata = nullptr;
1022 }
1023 return *this;
1024}
1025
1029 std::fill_n(mdata, nrows() * ncols(), x);
1030 return *this;
1031}
1032
1034
1047 resize(v.nelem(), 1);
1048 ConstMatrixView dummy(v);
1049 copy(dummy.begin(), dummy.end(), begin());
1050 return *this;
1057 ARTS_ASSERT(0 <= r);
1058 ARTS_ASSERT(0 <= c);
1059
1060 if (mrr.mextent != r || mcr.mextent != c) {
1061 delete[] mdata;
1062 mdata = new Numeric[r * c];
1063
1064 mrr.mstart = 0;
1065 mrr.mextent = r;
1066 mrr.mstride = c;
1067
1068 mcr.mstart = 0;
1069 mcr.mextent = c;
1071 }
1072}
1075void swap(Matrix& m1, Matrix& m2) {
1076 std::swap(m1.mrr, m2.mrr);
1077 std::swap(m1.mcr, m2.mcr);
1078 std::swap(m1.mdata, m2.mdata);
1080
1084 // cout << "Destroying a Matrix:\n"
1085 // << *this << "\n........................................\n";
1086 delete[] mdata;
1087}
1088
1089// Some general Matrix Vector functions:
1093 // Check dimensions:
1094 ARTS_ASSERT(a.nelem() == b.nelem());
1095
1096 const ConstIterator1D ae = a.end();
1097 ConstIterator1D ai = a.begin();
1098 ConstIterator1D bi = b.begin();
1099
1100 Numeric res = 0;
1101 for (; ai != ae; ++ai, ++bi) res += (*ai) * (*bi);
1102
1103 return res;
1104}
1105
1107
1124 ARTS_ASSERT(y.mrange.get_extent() == M.mrr.get_extent());
1125 ARTS_ASSERT(M.mcr.get_extent() == x.mrange.get_extent());
1126 ARTS_ASSERT((M.mcr.get_extent() != 0) && (M.mrr.get_extent() != 0));
1127
1128 if ((M.mcr.get_stride() == 1) || (M.mrr.get_stride() == 1)) {
1129 char trans;
1130 int m, n;
1131 double zero = 0.0;
1132 double one = 1.0;
1133 int LDA, incx, incy;
1134
1135 if (M.mcr.get_stride() != 1) {
1136 trans = 'n';
1137 m = (int)M.mrr.get_extent();
1138 n = (int)M.mcr.get_extent();
1139 LDA = (int)M.mcr.get_stride();
1140 } else {
1141 trans = 't';
1142 m = (int)M.mcr.get_extent();
1143 n = (int)M.mrr.get_extent();
1144 LDA = (int)M.mrr.get_stride();
1145 if (M.mrr.get_stride() == 1) LDA = m;
1146 }
1147
1148 incx = (int)x.mrange.get_stride();
1149 incy = (int)y.mrange.get_stride();
1150
1151 double* mstart = M.mdata + M.mcr.get_start() + M.mrr.get_start();
1152 double* ystart = y.mdata + y.mrange.get_start();
1153 double* xstart = x.mdata + x.mrange.get_start();
1154
1155 dgemv_(&trans,
1156 &m,
1157 &n,
1158 &one,
1159 mstart,
1160 &LDA,
1161 xstart,
1162 &incx,
1163 &zero,
1164 ystart,
1165 &incy);
1166
1167 } else {
1168 mult_general(y, M, x);
1169 }
1170}
1171
1182 const ConstMatrixView& M,
1183 const ConstVectorView& x) ARTS_NOEXCEPT {
1184 // Check dimensions:
1185 ARTS_ASSERT(y.mrange.mextent == M.mrr.mextent);
1186 ARTS_ASSERT(M.mcr.mextent == x.mrange.mextent);
1187 ARTS_ASSERT(M.mcr.mextent != 0 && M.mrr.mextent != 0);
1188
1189 // Let's first find the pointers to the starting positions
1190 Numeric* mdata = M.mdata + M.mcr.mstart + M.mrr.mstart;
1191 Numeric* xdata = x.mdata + x.mrange.mstart;
1192 Numeric* yelem = y.mdata + y.mrange.mstart;
1193
1194 Index i = M.mrr.mextent;
1195 while (i--) {
1196 Numeric* melem = mdata;
1197 Numeric* xelem = xdata; // Reset xelem to first element of source vector
1198
1199 // Multiply first element of matrix row with first element of source
1200 // vector. This is done outside the loop to avoid initialization of the
1201 // target vector's element with zero (saves one assignment)
1202 *yelem = *melem * *xelem;
1204 Index j = M.mcr.mextent; // --j (instead of j-- like in the outer loop)
1205 while (--j) // is important here because we only want
1206 { // mextent-1 cycles
1207 melem += M.mcr.mstride;
1208 xelem += x.mrange.mstride;
1209 *yelem += *melem * *xelem;
1210 }
1211
1212 mdata += M.mrr.mstride; // Jump to next matrix row
1213 yelem += y.mrange.mstride; // Jump to next element in target vector
1214 }
1215}
1216
1218
1242void mult(MatrixView A, const ConstMatrixView& B, const ConstMatrixView& C) {
1243 // Check dimensions:
1244 ARTS_ASSERT(A.nrows() == B.nrows());
1245 ARTS_ASSERT(A.ncols() == C.ncols());
1246 ARTS_ASSERT(B.ncols() == C.nrows());
1247
1248 // Catch trivial case if one of the matrices is empty.
1249 if ((B.nrows() == 0) || (B.ncols() == 0) || (C.ncols() == 0)) return;
1250
1251 // Matrices B and C must be continuous in at least on dimension, C
1252 // must be continuous along the second dimension.
1253 if (((B.mrr.get_stride() == 1) || (B.mcr.get_stride() == 1)) &&
1254 ((C.mrr.get_stride() == 1) || (C.mcr.get_stride() == 1)) &&
1255 (A.mcr.get_stride() == 1)) {
1256 // BLAS uses column-major order while arts uses row-major order.
1257 // Hence instead of C = A * B we compute C^T = A^T * B^T!
1258
1259 int k, m, n;
1260
1261 k = (int)B.ncols();
1262 m = (int)C.ncols();
1263 n = (int)B.nrows();
1264
1265 // Note also the clash in nomenclature: BLAS uses C = A * B while
1266 // arts uses A = B * C. Taking into accout this and the difference in
1267 // memory layouts, we need to map the MatrixViews A, B and C to the BLAS
1268 // arguments as follows:
1269 // A (arts) -> C (BLAS)
1270 // B (arts) -> B (BLAS)
1271 // C (arts) -> A (BLAS)
1272
1273 // Char indicating whether A (BLAS) or B (BLAS) should be transposed.
1274 char transa, transb;
1275 // Sizes of the matrices along the direction in which they are
1276 // traversed.
1277 int lda, ldb, ldc;
1278
1279 // Check if C (arts) is transposed.
1280 if (C.mrr.get_stride() == 1) {
1281 transa = 'T';
1282 lda = (int)C.mcr.get_stride();
1283 } else {
1284 transa = 'N';
1285 lda = (int)C.mrr.get_stride();
1286 }
1287
1288 // Check if B (arts) is transposed.
1289 if (B.mrr.get_stride() == 1) {
1290 transb = 'T';
1291 ldb = (int)B.mcr.get_stride();
1292 } else {
1293 transb = 'N';
1294 ldb = (int)B.mrr.get_stride();
1295 }
1296
1297 // In case B (arts) has only one column, column and row stride are 1.
1298 // We therefore need to set ldb to k, since dgemm_ requires lda to be at
1299 // least k / m if A is non-transposed / transposed.
1300 if ((B.mcr.get_stride() == 1) && (B.mrr.get_stride() == 1)) {
1301 transb = 'N';
1302 ldb = k;
1303 }
1304
1305 // The same holds for C (arts).
1306 if ((C.mcr.get_stride() == 1) && (C.mrr.get_stride() == 1)) {
1307 transa = 'N';
1308 lda = m;
1309 }
1310
1311 ldc = (int)A.mrr.get_stride();
1312 // The same holds for A (arts).
1313 if ((A.mcr.get_stride() == 1) && (A.mrr.get_stride() == 1)) {
1314 ldc = m;
1315 }
1316 double alpha = 1.0, beta = 0.0;
1317
1318 dgemm_(&transa,
1319 &transb,
1320 &m,
1321 &n,
1322 &k,
1323 &alpha,
1324 C.mdata + C.mrr.get_start() + C.mcr.get_start(),
1325 &lda,
1326 B.mdata + B.mrr.get_start() + B.mcr.get_start(),
1327 &ldb,
1328 &beta,
1329 A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1330 &ldc);
1331
1332 } else {
1333 mult_general(A, B, C);
1334 }
1335}
1336
1338
1347 const ConstMatrixView& B,
1348 const ConstMatrixView& C) ARTS_NOEXCEPT {
1349 // Check dimensions:
1350 ARTS_ASSERT(A.nrows() == B.nrows());
1351 ARTS_ASSERT(A.ncols() == C.ncols());
1352 ARTS_ASSERT(B.ncols() == C.nrows());
1353
1354 // Let's get the transpose of C, so that we can use 2D iterators to
1355 // access the columns (= rows of the transpose).
1356 ConstMatrixView CT = transpose(C);
1357
1358 const Iterator2D ae = A.end();
1359 Iterator2D ai = A.begin();
1360 ConstIterator2D bi = B.begin();
1361
1362 // This walks through the rows of A and B:
1363 for (; ai != ae; ++ai, ++bi) {
1364 const Iterator1D ace = ai->end();
1365 Iterator1D aci = ai->begin();
1366 ConstIterator2D cti = CT.begin();
1367
1368 // This walks through the columns of A with a 1D iterator, and
1369 // at the same time through the rows of CT, which are the columns of
1370 // C, with a 2D iterator:
1371 for (; aci != ace; ++aci, ++cti) {
1372 // The operator * is used to compute the scalar product
1373 // between rows of B and rows of C.transpose().
1374 *aci = (*bi) * (*cti);
1375 }
1376 }
1377}
1378
1380
1394 ARTS_ASSERT(a.nelem() == 3);
1395 ARTS_ASSERT(b.nelem() == 3);
1396 ARTS_ASSERT(c.nelem() == 3);
1397
1398 c[0] = a[1] * b[2] - a[2] * b[1];
1399 c[1] = a[2] * b[0] - a[0] * b[2];
1400 c[2] = a[0] * b[1] - a[1] * b[0];
1401}
1402
1413 ARTS_ASSERT(a.nelem() == b.nelem());
1414 Numeric arg = (a * b) / sqrt(a * a) / sqrt(b * b);
1415
1416 // Due to numerical inaccuracy, arg might be slightly larger than 1.
1417 // We catch those cases to avoid spurious returns of NaNs
1418 return fabs(arg) > 1. ? 0. : acos(arg) * RAD2DEG;
1419}
1420
1435 ARTS_ASSERT(a.nelem() == b.nelem());
1436 ARTS_ASSERT(a.nelem() == c.nelem());
1437
1438 const Numeric C = (a * b) / (a * a);
1439 c = a;
1440 c *= C;
1441};
1442
1445 return ConstMatrixView(m.mdata, m.mcr, m.mrr);
1446}
1447
1450MatrixView transpose(MatrixView m) ARTS_NOEXCEPT { return MatrixView(m.mdata, m.mcr, m.mrr); }
1451
1472void transform(VectorView y, double (&my_func)(double), ConstVectorView x) {
1473 // Check dimensions:
1474 ARTS_ASSERT(y.nelem() == x.nelem());
1475
1476 const ConstIterator1D xe = x.end();
1477 ConstIterator1D xi = x.begin();
1478 Iterator1D yi = y.begin();
1479 for (; xi != xe; ++xi, ++yi) *yi = my_func(*xi);
1480}
1481
1500void transform(MatrixView y, double (&my_func)(double), ConstMatrixView x) {
1501 // Check dimensions:
1502 ARTS_ASSERT(y.nrows() == x.nrows());
1503 ARTS_ASSERT(y.ncols() == x.ncols());
1504
1505 const ConstIterator2D rxe = x.end();
1506 ConstIterator2D rx = x.begin();
1507 Iterator2D ry = y.begin();
1508 for (; rx != rxe; ++rx, ++ry) {
1509 const ConstIterator1D cxe = rx->end();
1510 ConstIterator1D cx = rx->begin();
1511 Iterator1D cy = ry->begin();
1512 for (; cx != cxe; ++cx, ++cy) *cy = my_func(*cx);
1513 }
1514}
1515
1518 // Initial value for max:
1519 Numeric max = x[0];
1520
1521 const ConstIterator1D xe = x.end();
1522 ConstIterator1D xi = x.begin();
1523
1524 for (; xi != xe; ++xi) {
1525 if (*xi > max) max = *xi;
1526 }
1527
1528 return max;
1529}
1530
1533 // Initial value for max:
1534 Numeric max = x(0, 0);
1535
1536 const ConstIterator2D rxe = x.end();
1537 ConstIterator2D rx = x.begin();
1538
1539 for (; rx != rxe; ++rx) {
1540 const ConstIterator1D cxe = rx->end();
1541 ConstIterator1D cx = rx->begin();
1542
1543 for (; cx != cxe; ++cx)
1544 if (*cx > max) max = *cx;
1545 }
1546
1547 return max;
1548}
1549
1552 // Initial value for min:
1553 Numeric min = x[0];
1554
1555 const ConstIterator1D xe = x.end();
1556 ConstIterator1D xi = x.begin();
1557
1558 for (; xi != xe; ++xi) {
1559 if (*xi < min) min = *xi;
1560 }
1561
1562 return min;
1563}
1564
1567 // Initial value for min:
1568 Numeric min = x(0, 0);
1569
1570 const ConstIterator2D rxe = x.end();
1571 ConstIterator2D rx = x.begin();
1572
1573 for (; rx != rxe; ++rx) {
1574 const ConstIterator1D cxe = rx->end();
1575 ConstIterator1D cx = rx->begin();
1576
1577 for (; cx != cxe; ++cx)
1578 if (*cx < min) min = *cx;
1579 }
1580
1581 return min;
1582}
1583
1586 // Initial value for mean:
1587 Numeric mean = 0;
1588
1589 const ConstIterator1D xe = x.end();
1590 ConstIterator1D xi = x.begin();
1591
1592 for (; xi != xe; ++xi) mean += *xi;
1593
1594 mean /= (Numeric)x.nelem();
1595
1596 return mean;
1597}
1598
1601 // Initial value for mean:
1602 Numeric nanmean = 0;
1603
1604 // Counter for the number of normal values
1605 Index numnormal = 0;
1606
1607 const ConstIterator1D xe = x.end();
1608 ConstIterator1D xi = x.begin();
1609
1610 for (; xi != xe; ++xi) if (std::isnormal(*xi) or (*xi == 0)) {nanmean += *xi; numnormal++;}
1611
1612 if (numnormal)
1613 nanmean /= Numeric(numnormal);
1614 else
1615 nanmean = std::numeric_limits<Numeric>::quiet_NaN();
1616
1617 return nanmean;
1618}
1619
1622 // Initial value for mean:
1623 Numeric mean = 0;
1624
1625 const ConstIterator2D rxe = x.end();
1626 ConstIterator2D rx = x.begin();
1627
1628 for (; rx != rxe; ++rx) {
1629 const ConstIterator1D cxe = rx->end();
1630 ConstIterator1D cx = rx->begin();
1631
1632 for (; cx != cxe; ++cx) mean += *cx;
1633 }
1634
1635 mean /= (Numeric)(x.nrows() * x.ncols());
1636
1637 return mean;
1638}
1639
1650 // cout << "Assigning VectorView from Array<Numeric>.\n";
1651
1652 // Check that sizes are compatible:
1653 ARTS_ASSERT(mrange.mextent == v.nelem());
1654
1655 // Iterators for Array:
1657 const Array<Numeric>::const_iterator e = v.end();
1658 // Iterator for Vector:
1659 Iterator1D target = begin();
1660
1661 for (; i != e; ++i, ++target) *target = *i;
1662
1663 return *this;
1664}
1665
1666// Const
1667
1668// Converts constant matrix to constant eigen map
1670 return ConstMatrixViewMap(A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1671 A.nrows(),
1672 A.ncols(),
1674}
1675
1676// Converts constant vector to constant eigen row-view
1679 A.nelem(),
1680 1,
1681 StrideType(A.mrange.get_stride(), 1));
1682}
1683
1684// Converts constant vector to constant eigen row-view
1686 return MapToEigen(A);
1687}
1688
1689// Converts constant vector to constant eigen column-view
1692 1,
1693 A.nelem(),
1694 StrideType(1, A.mrange.get_stride()));
1695}
1696
1697// Non- const
1698
1699// Converts matrix to eigen map
1701 return MatrixViewMap(A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1702 A.nrows(),
1703 A.ncols(),
1705}
1706
1707// Converts vector to eigen map row-view
1709 return MatrixViewMap(A.mdata + A.mrange.get_start(),
1710 A.nelem(),
1711 1,
1712 StrideType(A.mrange.get_stride(), 1));
1713}
1714
1715// Converts vector to eigen map row-view
1717
1718// Converts vector to eigen map column-view
1720 return MatrixViewMap(A.mdata + A.mrange.get_start(),
1721 1,
1722 A.nelem(),
1723 StrideType(1, A.mrange.get_stride()));
1724}
1725
1726// Special 4x4
1727
1728// Converts matrix to eigen map
1730 return Matrix4x4ViewMap(A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1731 4,
1732 4,
1734}
1735
1736// Converts constant matrix to constant eigen map
1738 return ConstMatrix4x4ViewMap(
1739 A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1740 4,
1741 4,
1743}
1744
1746// Helper function for debugging
1747#ifndef NDEBUG
1748
1763 return mv(r, c);
1764}
1765
1766#endif
Index nrows
void * data
Index ncols
Interface for BLAS library.
void dgemv_(char *trans, int *m, int *n, double *alpha, double *A, int *LDA, double *x, int *incx, double *beta, double *y, int *incy)
Matrix-Vector Multiplication.
void dgemm_(char *transa, char *transb, int *m, int *n, int *k, double *alpha, double *A, int *lda, double *B, int *ldb, double *beta, double *C, int *ldc)
BLAS matrix multiplication.
This can be used to make arrays out of anything.
Definition: array.h:107
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
The constant iterator class for sub vectors.
Definition: matpackI.h:425
const Numeric * mx
Current position.
Definition: matpackI.h:473
Index mstride
Stride.
Definition: matpackI.h:475
The const row iterator class for sub matrices.
Definition: matpackI.h:825
A constant view of a Matrix.
Definition: matpackI.h:1014
Range mrr
The row range of mdata that is actually used.
Definition: matpackI.h:1109
bool empty() const ARTS_NOEXCEPT
Returns true if variable size is zero.
Definition: matpackI.cc:430
friend class ConstVectorView
Definition: matpackI.h:1063
ConstIterator2D begin() const ARTS_NOEXCEPT
Return const iterator to first row.
Definition: matpackI.cc:473
ConstIterator2D end() const ARTS_NOEXCEPT
Return const iterator behind last row.
Definition: matpackI.cc:478
Numeric operator()(Index r, Index c) const ARTS_NOEXCEPT
Plain const index operator.
Definition: matpackI.h:1031
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1113
Range mcr
The column range of mdata that is actually used.
Definition: matpackI.h:1111
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
ConstVectorView diagonal() const ARTS_NOEXCEPT
Matrix diagonal as vector.
Definition: matpackI.cc:493
ConstMatrixView()=default
A constant view of a Vector.
Definition: matpackI.h:489
ConstIterator1D end() const ARTS_NOEXCEPT
Return const iterator behind last element.
Definition: matpackI.cc:70
bool empty() const ARTS_NOEXCEPT
Returns true if variable size is zero.
Definition: matpackI.cc:46
ConstIterator1D begin() const ARTS_NOEXCEPT
Return const iterator to first element.
Definition: matpackI.cc:66
ConstVectorView()=default
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:612
Numeric sum() const ARTS_NOEXCEPT
The sum of all elements of a Vector.
Definition: matpackI.cc:52
Numeric operator[](Index n) const ARTS_NOEXCEPT
Plain const index operator.
Definition: matpackI.h:523
Index size() const ARTS_NOEXCEPT
Definition: matpackI.cc:50
Range mrange
The range of mdata that is actually used.
Definition: matpackI.h:610
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The iterator class for sub vectors.
Definition: matpackI.h:365
Index mstride
Stride.
Definition: matpackI.h:420
Numeric * mx
Current position.
Definition: matpackI.h:418
The row iterator class for sub matrices.
Definition: matpackI.h:781
The Joker class.
Definition: matpackI.h:131
The MatrixView class.
Definition: matpackI.h:1125
MatrixView & operator=(const ConstMatrixView &m)
Assignment operator.
Definition: matpackI.cc:637
friend class VectorView
Definition: matpackI.h:1194
MatrixView & operator/=(Numeric x) ARTS_NOEXCEPT
Division by scalar.
Definition: matpackI.cc:704
Numeric & operator()(Index r, Index c) ARTS_NOEXCEPT
Plain index operator.
Definition: matpackI.h:1140
Iterator2D end() ARTS_NOEXCEPT
Return iterator behind last row.
Definition: matpackI.cc:627
MatrixView & operator+=(Numeric x) ARTS_NOEXCEPT
Addition of scalar.
Definition: matpackI.cc:714
Iterator2D begin() ARTS_NOEXCEPT
‍** Return const iterator to first row. Has to be redefined here, since it is
Definition: matpackI.cc:622
MatrixView & operator*=(Numeric x) ARTS_NOEXCEPT
Multiplication by scalar.
Definition: matpackI.cc:694
MatrixView & operator-=(Numeric x) ARTS_NOEXCEPT
Subtraction of scalar.
Definition: matpackI.cc:724
MatrixView()=default
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array.
Definition: matpackI.cc:739
The Matrix class.
Definition: matpackI.h:1225
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
Matrix()=default
virtual ~Matrix()
Destructor for Matrix.
Definition: matpackI.cc:1083
Matrix & operator=(const Matrix &m)
Assignment operator from another matrix.
Definition: matpackI.cc:1004
The range class.
Definition: matpackI.h:165
constexpr Index get_stride() const ARTS_NOEXCEPT
Returns the stride of the range.
Definition: matpackI.h:336
Index mstart
The start index.
Definition: matpackI.h:351
constexpr Index get_start() const ARTS_NOEXCEPT
Returns the start index of the range.
Definition: matpackI.h:332
Index mstride
The stride.
Definition: matpackI.h:360
constexpr Index get_extent() const ARTS_NOEXCEPT
Returns the extent of the range.
Definition: matpackI.h:334
Index mextent
The number of elements.
Definition: matpackI.h:358
The VectorView class.
Definition: matpackI.h:626
VectorView operator+=(Numeric x) ARTS_NOEXCEPT
Addition of scalar.
Definition: matpackI.cc:202
VectorView()=default
VectorView operator-=(Numeric x) ARTS_NOEXCEPT
Subtraction of scalar.
Definition: matpackI.cc:208
VectorView & operator=(const ConstVectorView &v)
Assignment operator.
Definition: matpackI.cc:152
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array, const-version.
Definition: matpackI.cc:271
Numeric & operator[](Index n) ARTS_NOEXCEPT
Plain Index operator.
Definition: matpackI.h:647
Iterator1D begin() ARTS_NOEXCEPT
Return iterator to first element.
Definition: matpackI.cc:143
VectorView operator*=(Numeric x) ARTS_NOEXCEPT
Multiplication by scalar.
Definition: matpackI.cc:190
Iterator1D end() ARTS_NOEXCEPT
Return iterator behind last element.
Definition: matpackI.cc:147
VectorView operator/=(Numeric x) ARTS_NOEXCEPT
Division by scalar.
Definition: matpackI.cc:196
The Vector class.
Definition: matpackI.h:876
virtual ~Vector()
Destructor for Vector.
Definition: matpackI.cc:424
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
Vector()=default
Vector & operator=(const Vector &v)
Assignment from another Vector.
Definition: matpackI.cc:378
const Numeric RAD2DEG
Global constant, conversion from radians to degrees.
#define ARTS_NOEXCEPT
Definition: debug.h:80
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
The declarations of all the exception classes.
#define beta
Numeric mean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version.
Definition: matpackI.cc:1585
ConstMatrixViewMap MapToEigenRow(const ConstVectorView &A)
Definition: matpackI.cc:1685
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix-Vector Multiplication.
Definition: matpackI.cc:1123
Numeric vector_angle(ConstVectorView a, ConstVectorView b)
Definition: matpackI.cc:1412
Numeric max(const ConstVectorView &x) ARTS_NOEXCEPT
Max function, vector version.
Definition: matpackI.cc:1517
ConstMatrixViewMap MapToEigen(const ConstMatrixView &A)
Definition: matpackI.cc:1669
Numeric nanmean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version ignoring nans and infs.
Definition: matpackI.cc:1600
ConstMatrixViewMap MapToEigenCol(const ConstVectorView &A)
Definition: matpackI.cc:1690
void proj(Vector &c, ConstVectorView a, ConstVectorView b) ARTS_NOEXCEPT
Definition: matpackI.cc:1434
void mult_general(VectorView y, const ConstMatrixView &M, const ConstVectorView &x) ARTS_NOEXCEPT
Matrix Vector multiplication.
Definition: matpackI.cc:1181
const Joker joker
std::ostream & operator<<(std::ostream &os, const Range &r)
Definition: matpackI.cc:37
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:1472
Matrix4x4ViewMap MapToEigen4x4(MatrixView &A)
Definition: matpackI.cc:1729
void swap(Vector &v1, Vector &v2)
Definition: matpackI.cc:419
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
cross3
Definition: matpackI.cc:1393
Numeric min(const ConstVectorView &x) ARTS_NOEXCEPT
Min function, vector version.
Definition: matpackI.cc:1551
void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Definition: matpackI.cc:299
Numeric debug_matrixview_get_elem(MatrixView &mv, Index r, Index c)
Helper function to access matrix elements.
Definition: matpackI.cc:1762
Numeric operator*(const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
Scalar product.
Definition: matpackI.cc:1092
ConstMatrixView transpose(ConstMatrixView m) ARTS_NOEXCEPT
Const version of transpose.
Definition: matpackI.cc:1444
Implementation of Matrix, Vector, and such stuff.
Eigen::Map< MatrixType, 0, StrideType > MatrixViewMap
Definition: matpackI.h:115
Eigen::Map< const Matrix4x4Type, 0, StrideType > ConstMatrix4x4ViewMap
Definition: matpackI.h:119
Eigen::Map< const MatrixType, 0, StrideType > ConstMatrixViewMap
Definition: matpackI.h:116
Eigen::Map< Matrix4x4Type, 0, StrideType > Matrix4x4ViewMap
Definition: matpackI.h:118
Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > StrideType
Definition: matpackI.h:109
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Workspace & init(Workspace &ws)
Index nelem(const Lines &l)
Number of lines.
constexpr std::string_view Joker
Definition: isotopologues.h:11
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
Definition: zeemandata.h:78
constexpr Rational end(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the largest M for a polarization type of this transition.
Definition: zeemandata.h:109
#define v
#define a
#define c
#define b
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739
#define M
Definition: rng.cc:165