ARTS 2.5.4 (git: 4c0d3b4d)
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
26#include <cmath>
27#include <cstring>
28
29#include "blas.h"
30#include "exceptions.h"
31
32using std::setw;
33
34// Define the global joker object:
35extern const Joker joker = Joker();
36
37static const Numeric RAD2DEG = 57.295779513082323;
38
39std::ostream& operator<<(std::ostream& os, const Range& r) {
40 os << "Range(" << r.get_start() << ", " << r.get_extent() << ", "
41 << r.get_stride() << ")";
42 return os;
43}
44
45// Functions for ConstVectorView:
46// ------------------------------
47
49 Numeric s = 0;
51 const ConstIterator1D e = end();
52
53 for (; i != e; ++i) s += *i;
54
55 return s;
56}
57
60 return ConstVectorView(mdata, mrange, r);
61}
62
65}
66
68 return ConstIterator1D(
71}
72
73ConstVectorView::operator ConstMatrixView() const {
74 // The old version (before 2013-01-18) of this was:
75 // return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
76 // Bus this was a bug! The problem is that the matrix index operator adds
77 // the mstart from both row and columm range object to mdata
78
79 return ConstMatrixView(mdata, mrange, Range(0, 1));
80}
81
82/* This one is a bit tricky: We have to cast away the arguments const
83 qualifier, because mdata is not const. This should be safe, since
84 there are no non-const methods for ConstVectorView. */
86 : mrange(0, 1),
87 mdata(&const_cast<Numeric&>(a)) {
88 // Nothing to do here.
89}
90
92 const Range& range) ARTS_NOEXCEPT
93 : mrange(range),
94 mdata(data) {
95 // Nothing to do here.
96}
97
99 const Range& p,
100 const Range& n) ARTS_NOEXCEPT : mrange(p, n),
101 mdata(data) {
102 // Nothing to do here.
103}
104
105/* Output operator. This demonstrates how iterators can be used to
106 traverse the vector. The iterators know which part of the vector
107 is `active', and also the stride. */
108std::ostream& operator<<(std::ostream& os, const ConstVectorView& v) {
109 ConstIterator1D i = v.begin();
110 const ConstIterator1D end = v.end();
111
112 if (i != end) {
113 os << setw(3) << *i;
114 ++i;
115 }
116 for (; i != end; ++i) {
117 os << " " << setw(3) << *i;
118 }
119
120 return os;
121}
122
123// Functions for VectorView:
124// ------------------------
125
127 false,
128 "Creating a VectorView from a const Vector is not allowed.\n"
129 "This is not really a runtime error, but I don't want to start\n"
130 "producing direct output from inside matpack. And just exiting is\n"
131 "not so nice.\n"
132 "If you see this error, there is a bug in the code, not in the\n"
133 "ARTS input.")}
134
136 mdata = v.mdata;
137 mrange = v.mrange;
138}
139
141 return VectorView(mdata, mrange, r);
142}
143
146}
147
151}
152
154 // cout << "Assigning VectorView from ConstVectorView.\n";
155
156 // Check that sizes are compatible:
157 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
158
159 copy(v.begin(), v.end(), begin());
160
161 return *this;
162}
163
165 // cout << "Assigning VectorView from VectorView.\n";
166
167 // Check that sizes are compatible:
168 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
169
170 copy(v.begin(), v.end(), begin());
171
172 return *this;
173}
174
176 // cout << "Assigning VectorView from Vector.\n";
177
178 // Check that sizes are compatible:
179 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
180
181 copy(v.begin(), v.end(), begin());
182
183 return *this;
184}
185
187 copy(x, begin(), end());
188 return *this;
189}
190
192 const Iterator1D e = end();
193 for (Iterator1D i = begin(); i != e; ++i) *i *= x;
194 return *this;
195}
196
198 const Iterator1D e = end();
199 for (Iterator1D i = begin(); i != e; ++i) *i /= x;
200 return *this;
201}
202
204 const Iterator1D e = end();
205 for (Iterator1D i = begin(); i != e; ++i) *i += x;
206 return *this;
207}
208
210 const Iterator1D e = end();
211 for (Iterator1D i = begin(); i != e; ++i) *i -= x;
212 return *this;
213}
214
216 ARTS_ASSERT(nelem() == x.nelem());
217
218 ConstIterator1D s = x.begin();
219
220 Iterator1D i = begin();
221 const Iterator1D e = end();
222
223 for (; i != e; ++i, ++s) *i *= *s;
224 return *this;
225}
226
228 ARTS_ASSERT(nelem() == x.nelem());
229
230 ConstIterator1D s = x.begin();
231
232 Iterator1D i = begin();
233 const Iterator1D e = end();
234
235 for (; i != e; ++i, ++s) *i /= *s;
236 return *this;
237}
238
240 ARTS_ASSERT(nelem() == x.nelem());
241
242 ConstIterator1D s = x.begin();
243
244 Iterator1D i = begin();
245 const Iterator1D e = end();
246
247 for (; i != e; ++i, ++s) *i += *s;
248 return *this;
249}
250
252 ARTS_ASSERT(nelem() == x.nelem());
253
254 ConstIterator1D s = x.begin();
255
256 Iterator1D i = begin();
257 const Iterator1D e = end();
258
259 for (; i != e; ++i, ++s) *i -= *s;
260 return *this;
261}
262
263VectorView::operator MatrixView() ARTS_NOEXCEPT {
264 // The old version (before 2013-01-18) of this was:
265 // return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
266 // Bus this was a bug! The problem is that the matrix index operator adds
267 // the mstart from both row and columm range object to mdata
268
269 return MatrixView(mdata, mrange, Range(0, 1));
270}
271
274 mrange.mstart == 0 and (mrange.mstride == 1 or mrange.mextent == 0),
275 mrange)
276
277 return mdata;
278}
279
282 mrange.mstart == 0 and (mrange.mstride == 1 or mrange.mextent == 0),
283 mrange)
284
285 return mdata;
286}
287
289 // Nothing to do here.
290}
291
293 : ConstVectorView(data, range) {
294 // Nothing to do here.
295}
296
298 const Range& p,
299 const Range& n) ARTS_NOEXCEPT
300 : ConstVectorView(data, p, n) {
301 // Nothing to do here.
302}
303
305 const ConstIterator1D& end,
306 Iterator1D target) {
307 if (origin.mstride == 1 && target.mstride == 1)
308 memcpy((void*)target.mx,
309 (void*)origin.mx,
310 sizeof(Numeric) * (end.mx - origin.mx));
311 else
312 for (; origin != end; ++origin, ++target) *target = *origin;
313}
314
316 for (; target != end; ++target) *target = x;
317}
318
319// Functions for Vector:
320// ---------------------
321
322Vector::Vector(std::initializer_list<Numeric> init)
323 : VectorView(new Numeric[init.size()], Range(0, init.size())) {
324 std::copy(init.begin(), init.end(), begin());
325}
326
327Vector::Vector(const Eigen::VectorXd& init)
328 : VectorView(new Numeric[init.size()], Range(0, init.size())) {
329 for (Index i = 0; i < size(); i++) operator[](i) = init[i];
330}
331
333 // Nothing to do here.
334}
335
337 : VectorView(new Numeric[n], Range(0, n)) {
338 // Here we can access the raw memory directly, for slightly
339 // increased efficiency:
340 std::fill_n(mdata, n, fill);
341}
342
344 : VectorView(new Numeric[extent], Range(0, extent)) {
345 // Fill with values:
346 Numeric x = start;
347 Iterator1D i = begin();
348 const Iterator1D e = end();
349 for (; i != e; ++i) {
350 *i = x;
351 x += stride;
352 }
353}
354
356 : VectorView(new Numeric[v.nelem()], Range(0, v.nelem())) {
357 copy(v.begin(), v.end(), begin());
358}
359
361 : VectorView(new Numeric[v.nelem()], Range(0, v.nelem())) {
362 std::memcpy(mdata, v.mdata, nelem() * sizeof(Numeric));
363}
364
365Vector::Vector(const std::vector<Numeric>& v)
366 : VectorView(new Numeric[v.size()], Range(0, v.size())) {
367 std::vector<Numeric>::const_iterator vec_it_end = v.end();
368 Iterator1D this_it = this->begin();
369 for (std::vector<Numeric>::const_iterator vec_it = v.begin();
370 vec_it != vec_it_end;
371 ++vec_it, ++this_it)
372 *this_it = *vec_it;
373}
374
375Vector& Vector::operator=(std::initializer_list<Numeric> v) {
376 resize(v.size());
377 std::copy(v.begin(), v.end(), begin());
378 return *this;
379}
380
382 if (this != &v) {
383 resize(v.nelem());
384 std::memcpy(mdata, v.mdata, nelem() * sizeof(Numeric));
385 }
386 return *this;
387}
388
390 if (this != &v) {
391 delete[] mdata;
392 mdata = v.mdata;
393 mrange = v.mrange;
394 v.mrange = Range(0, 0);
395 v.mdata = nullptr;
396 }
397 return *this;
398}
399
401 resize(x.nelem());
403 return *this;
404}
405
407 std::fill_n(mdata, nelem(), x);
408 return *this;
409}
410
412 ARTS_ASSERT(0 <= n);
413 if (mrange.mextent != n) {
414 delete[] mdata;
415 mdata = new Numeric[n];
416 mrange.mstart = 0;
417 mrange.mextent = n;
418 mrange.mstride = 1;
419 }
420}
421
422void swap(Vector& v1, Vector& v2) {
423 std::swap(v1.mrange, v2.mrange);
424 std::swap(v1.mdata, v2.mdata);
425}
426
427Vector::~Vector() { delete[] mdata; }
428
429// Functions for ConstMatrixView:
430// ------------------------------
431
436 const Range& r, const Range& c) const ARTS_NOEXCEPT {
437 return ConstMatrixView(mdata, mrr, mcr, r, c);
438}
439
446 Index c) const ARTS_NOEXCEPT {
447 // Check that c is valid:
448 ARTS_ASSERT(0 <= c);
449 ARTS_ASSERT(c < mcr.mextent);
450
451 return ConstVectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
452}
453
461 // Check that r is valid:
462 ARTS_ASSERT(0 <= r);
463 ARTS_ASSERT(r < mrr.mextent);
464
465 return ConstVectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
466}
467
471}
472
475 return ConstIterator2D(
477 mrr.mstride);
478}
479
481
492 Range(0, n, mrr.mstride + mcr.mstride));
493}
494
499 const Range& rr,
500 const Range& cr) ARTS_NOEXCEPT : mrr(rr),
501 mcr(cr),
502 mdata(data) {
503 // Nothing to do here.
504}
505
521 const Range& pr,
522 const Range& pc,
523 const Range& nr,
524 const Range& nc) ARTS_NOEXCEPT : mrr(pr, nr),
525 mcr(pc, nc),
526 mdata(data) {
527 // Nothing to do here.
528}
529
537std::ostream& operator<<(std::ostream& os, const ConstMatrixView& v) {
538 // Row iterators:
539 ConstIterator2D ir = v.begin();
540 const ConstIterator2D end_row = v.end();
541
542 if (ir != end_row) {
543 ConstIterator1D ic = ir->begin();
544 const ConstIterator1D end_col = ir->end();
545
546 if (ic != end_col) {
547 os << setw(3) << *ic;
548 ++ic;
549 }
550 for (; ic != end_col; ++ic) {
551 os << " " << setw(3) << *ic;
552 }
553 ++ir;
554 }
555 for (; ir != end_row; ++ir) {
556 ConstIterator1D ic = ir->begin();
557 const ConstIterator1D end_col = ir->end();
558
559 os << "\n";
560 if (ic != end_col) {
561 os << setw(3) << *ic;
562 ++ic;
563 }
564 for (; ic != end_col; ++ic) {
565 os << " " << setw(3) << *ic;
566 }
567 }
568
569 return os;
571
572// Functions for MatrixView:
573// -------------------------
574
580 return MatrixView(mdata, mrr, mcr, r, c);
581}
582
588 // Check that c is valid:
589 ARTS_ASSERT(0 <= c);
590 ARTS_ASSERT(c < mcr.mextent);
591
592 return VectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
593}
594
600 // Check that r is valid:
601 ARTS_ASSERT(0 <= r);
602 ARTS_ASSERT(r < mrr.mextent);
603
604 return VectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
605}
606
608// hiden by the non-const operator of the derived class.*/
609//ConstIterator2D MatrixView::begin() const
610//{
611// return ConstMatrixView::begin();
612//}
615//ConstIterator2D MatrixView::end() const
616//{
617// return ConstMatrixView::end();
618//}
619
623}
624
627 return Iterator2D(
629 mrr.mstride);
630}
631
637 // Check that sizes are compatible:
640
641 copy(m.begin(), m.end(), begin());
642 return *this;
643}
644
651 // Check that sizes are compatible:
654
655 copy(m.begin(), m.end(), begin());
656 return *this;
657}
658
663 // Check that sizes are compatible:
666
667 copy(m.begin(), m.end(), begin());
668 return *this;
669}
670
676 // Check that sizes are compatible:
677 ARTS_ASSERT(mrr.mextent == v.nelem());
678 ARTS_ASSERT(mcr.mextent == 1);
679 // dummy = ConstMatrixView(v.mdata,v.mrange,Range(v.mrange.mstart,1));;
680 ConstMatrixView dummy(v);
681 copy(dummy.begin(), dummy.end(), begin());
682 return *this;
683}
684
688 copy(x, begin(), end());
689 return *this;
690}
691
694 const Iterator2D er = end();
695 for (Iterator2D r = begin(); r != er; ++r) {
696 const Iterator1D ec = r->end();
697 for (Iterator1D c = r->begin(); c != ec; ++c) *c *= x;
698 }
699 return *this;
700}
701
704 const Iterator2D er = end();
705 for (Iterator2D r = begin(); r != er; ++r) {
706 const Iterator1D ec = r->end();
707 for (Iterator1D c = r->begin(); c != ec; ++c) *c /= x;
708 }
709 return *this;
710}
711
714 const Iterator2D er = end();
715 for (Iterator2D r = begin(); r != er; ++r) {
716 const Iterator1D ec = r->end();
717 for (Iterator1D c = r->begin(); c != ec; ++c) *c += x;
718 }
719 return *this;
720}
721
724 const Iterator2D er = end();
725 for (Iterator2D r = begin(); r != er; ++r) {
726 const Iterator1D ec = r->end();
727 for (Iterator1D c = r->begin(); c != ec; ++c) *c -= x;
728 }
729 return *this;
730}
731
739 ARTS_ASSERT(mrr.mstart == 0 and (mrr.mstride == mcr.mextent or size() == 0),
740 "Row ",
741 mrr)
743 mcr.mstart == 0 and (mcr.mstride == 1 or size() == 0), "Column ", mcr)
744
745 return mdata;
746}
747
755 ARTS_ASSERT(mrr.mstart == 0 and (mrr.mstride == mcr.mextent or size() == 0),
756 "Row ",
757 mrr)
759 mcr.mstart == 0 and (mcr.mstride == 1 or size() == 0), "Column ", mcr)
760
761 return mdata;
762}
763
766 ARTS_ASSERT(nrows() == x.nrows());
767 ARTS_ASSERT(ncols() == x.ncols());
768 ConstIterator2D sr = x.begin();
769 Iterator2D r = begin();
770 const Iterator2D er = end();
771 for (; r != er; ++r, ++sr) {
772 ConstIterator1D sc = sr->begin();
773 Iterator1D c = r->begin();
774 const Iterator1D ec = r->end();
775 for (; c != ec; ++c, ++sc) *c *= *sc;
776 }
777 return *this;
778}
779
782 ARTS_ASSERT(nrows() == x.nrows());
783 ARTS_ASSERT(ncols() == x.ncols());
784 ConstIterator2D sr = x.begin();
785 Iterator2D r = begin();
786 const Iterator2D er = end();
787 for (; r != er; ++r, ++sr) {
788 ConstIterator1D sc = sr->begin();
789 Iterator1D c = r->begin();
790 const Iterator1D ec = r->end();
791 for (; c != ec; ++c, ++sc) *c /= *sc;
792 }
793 return *this;
794}
795
798 ARTS_ASSERT(nrows() == x.nrows());
799 ARTS_ASSERT(ncols() == x.ncols());
800 ConstIterator2D sr = x.begin();
801 Iterator2D r = begin();
802 const Iterator2D er = end();
803 for (; r != er; ++r, ++sr) {
804 ConstIterator1D sc = sr->begin();
805 Iterator1D c = r->begin();
806 const Iterator1D ec = r->end();
807 for (; c != ec; ++c, ++sc) *c += *sc;
808 }
809 return *this;
810}
811
814 ARTS_ASSERT(nrows() == x.nrows());
815 ARTS_ASSERT(ncols() == x.ncols());
816 ConstIterator2D sr = x.begin();
817 Iterator2D r = begin();
818 const Iterator2D er = end();
819 for (; r != er; ++r, ++sr) {
820 ConstIterator1D sc = sr->begin();
821 Iterator1D c = r->begin();
822 const Iterator1D ec = r->end();
823 for (; c != ec; ++c, ++sc) *c -= *sc;
824 }
825 return *this;
826}
827
830 ARTS_ASSERT(nrows() == x.nelem());
831 ARTS_ASSERT(ncols() == 1);
832 ConstIterator1D sc = x.begin();
833 Iterator2D r = begin();
834 const Iterator2D er = end();
835 for (; r != er; ++r, ++sc) {
836 Iterator1D c = r->begin();
837 *c *= *sc;
838 }
839 return *this;
840}
841
844 ARTS_ASSERT(nrows() == x.nelem());
845 ARTS_ASSERT(ncols() == 1);
846 ConstIterator1D sc = x.begin();
847 Iterator2D r = begin();
848 const Iterator2D er = end();
849 for (; r != er; ++r, ++sc) {
850 Iterator1D c = r->begin();
851 *c /= *sc;
852 }
853 return *this;
854}
855
858 ARTS_ASSERT(nrows() == x.nelem());
859 ARTS_ASSERT(ncols() == 1);
860 ConstIterator1D sc = x.begin();
861 Iterator2D r = begin();
862 const Iterator2D er = end();
863 for (; r != er; ++r, ++sc) {
864 Iterator1D c = r->begin();
865 *c += *sc;
866 }
867 return *this;
868}
869
872 ARTS_ASSERT(nrows() == x.nelem());
873 ARTS_ASSERT(ncols() == 1);
874 ConstIterator1D sc = x.begin();
875 Iterator2D r = begin();
876 const Iterator2D er = end();
877 for (; r != er; ++r, ++sc) {
878 Iterator1D c = r->begin();
879 *c -= *sc;
880 }
881 return *this;
882}
883
888 const Range& rr,
889 const Range& cr) ARTS_NOEXCEPT
890 : ConstMatrixView(data, rr, cr) {
891 // Nothing to do here.
892}
893
909 const Range& pr,
910 const Range& pc,
911 const Range& nr,
912 const Range& nc) ARTS_NOEXCEPT
913 : ConstMatrixView(data, pr, pc, nr, nc) {
914 // Nothing to do here.
915}
916
927 const ConstIterator2D& end,
928 Iterator2D target) {
929 for (; origin != end; ++origin, ++target) {
930 ConstIterator1D o = origin->begin();
931 const ConstIterator1D e = origin->end();
932 Iterator1D t = target->begin();
933 for (; o != e; ++o, ++t) *t = *o;
934 }
935}
936
939 for (; target != end; ++target) {
940 Iterator1D t = target->begin();
941 const Iterator1D e = target->end();
942 for (; t != e; ++t) *t = x;
943 }
944}
945
946// Functions for Matrix:
947// ---------------------
948
952 : MatrixView(new Numeric[r * c], Range(0, r, c), Range(0, c)) {
953 // Nothing to do here.
954}
955
958 : MatrixView(new Numeric[r * c], Range(0, r, c), Range(0, c)) {
959 // Here we can access the raw memory directly, for slightly
960 // increased efficiency:
961 std::fill_n(mdata, r * c, fill);
962}
963
967 : MatrixView(new Numeric[m.nrows() * m.ncols()],
968 Range(0, m.nrows(), m.ncols()),
969 Range(0, m.ncols())) {
970 copy(m.begin(), m.end(), begin());
971}
972
976 : MatrixView(new Numeric[m.nrows() * m.ncols()],
977 Range(0, m.nrows(), m.ncols()),
978 Range(0, m.ncols())) {
979 // There is a catch here: If m is an empty matrix, then it will have
980 // 0 colunns. But this is used to initialize the stride of the row
981 // Range! Thus, this method has to be consistent with the behaviour
982 // of Range::Range. For now, Range::Range allows also stride 0.
983 std::memcpy(mdata, m.mdata, nrows() * ncols() * sizeof(Numeric));
984}
985
987
1012 if (this != &m) {
1013 resize(m.nrows(), m.ncols());
1014 std::memcpy(mdata, m.mdata, nrows() * ncols() * sizeof(Numeric));
1015 }
1016 return *this;
1017}
1018
1021 if (this != &m) {
1022 delete[] mdata;
1023 mdata = m.mdata;
1024 mrr = m.mrr;
1025 mcr = m.mcr;
1026 m.mrr = Range(0, 0);
1027 m.mcr = Range(0, 0);
1028 m.mdata = nullptr;
1029 }
1030 return *this;
1031}
1032
1036 std::fill_n(mdata, nrows() * ncols(), x);
1037 return *this;
1038}
1039
1041
1054 resize(v.nelem(), 1);
1055 ConstMatrixView dummy(v);
1056 copy(dummy.begin(), dummy.end(), begin());
1057 return *this;
1058}
1059
1064 ARTS_ASSERT(0 <= r);
1065 ARTS_ASSERT(0 <= c);
1066
1067 if (mrr.mextent != r || mcr.mextent != c) {
1068 delete[] mdata;
1069 mdata = new Numeric[r * c];
1070
1071 mrr.mstart = 0;
1072 mrr.mextent = r;
1073 mrr.mstride = c;
1074
1075 mcr.mstart = 0;
1076 mcr.mextent = c;
1077 mcr.mstride = 1;
1078 }
1079}
1080
1082void swap(Matrix& m1, Matrix& m2) {
1083 std::swap(m1.mrr, m2.mrr);
1084 std::swap(m1.mcr, m2.mcr);
1085 std::swap(m1.mdata, m2.mdata);
1086}
1087
1091 // cout << "Destroying a Matrix:\n"
1092 // << *this << "\n........................................\n";
1093 delete[] mdata;
1094}
1096// Some general Matrix Vector functions:
1097
1101 // Check dimensions:
1102 ARTS_ASSERT(a.nelem() == b.nelem());
1103
1104 const ConstIterator1D ae = a.end();
1105 ConstIterator1D ai = a.begin();
1106 ConstIterator1D bi = b.begin();
1107
1108 Numeric res = 0;
1109 for (; ai != ae; ++ai, ++bi) res += (*ai) * (*bi);
1110
1111 return res;
1112}
1113
1115
1132 ARTS_ASSERT(y.mrange.get_extent() == M.mrr.get_extent());
1133 ARTS_ASSERT(M.mcr.get_extent() == x.mrange.get_extent());
1134 ARTS_ASSERT((M.mcr.get_extent() != 0) && (M.mrr.get_extent() != 0));
1136 if ((M.mcr.get_stride() == 1) || (M.mrr.get_stride() == 1)) {
1137 char trans;
1138 int m, n;
1139 double zero = 0.0;
1140 double one = 1.0;
1141 int LDA, incx, incy;
1142
1143 if (M.mcr.get_stride() != 1) {
1144 trans = 'n';
1145 m = (int)M.mrr.get_extent();
1146 n = (int)M.mcr.get_extent();
1147 LDA = (int)M.mcr.get_stride();
1148 } else {
1149 trans = 't';
1150 m = (int)M.mcr.get_extent();
1151 n = (int)M.mrr.get_extent();
1152 LDA = (int)M.mrr.get_stride();
1153 if (M.mrr.get_stride() == 1) LDA = m;
1154 }
1155
1156 incx = (int)x.mrange.get_stride();
1157 incy = (int)y.mrange.get_stride();
1158
1159 double* mstart = M.mdata + M.mcr.get_start() + M.mrr.get_start();
1160 double* ystart = y.mdata + y.mrange.get_start();
1161 double* xstart = x.mdata + x.mrange.get_start();
1162
1163 dgemv_(&trans,
1164 &m,
1165 &n,
1166 &one,
1167 mstart,
1168 &LDA,
1169 xstart,
1170 &incx,
1171 &zero,
1172 ystart,
1173 &incy);
1174
1175 } else {
1176 mult_general(y, M, x);
1177 }
1178}
1179
1190 const ConstMatrixView& M,
1191 const ConstVectorView& x) ARTS_NOEXCEPT {
1192 // Check dimensions:
1193 ARTS_ASSERT(y.mrange.mextent == M.mrr.mextent);
1194 ARTS_ASSERT(M.mcr.mextent == x.mrange.mextent);
1195 ARTS_ASSERT(M.mcr.mextent != 0 && M.mrr.mextent != 0);
1196
1197 // Let's first find the pointers to the starting positions
1198 Numeric* mdata = M.mdata + M.mcr.mstart + M.mrr.mstart;
1199 Numeric* xdata = x.mdata + x.mrange.mstart;
1200 Numeric* yelem = y.mdata + y.mrange.mstart;
1201
1202 Index i = M.mrr.mextent;
1203 while (i--) {
1204 Numeric* melem = mdata;
1205 Numeric* xelem = xdata; // Reset xelem to first element of source vector
1206
1207 // Multiply first element of matrix row with first element of source
1208 // vector. This is done outside the loop to avoid initialization of the
1209 // target vector's element with zero (saves one assignment)
1210 *yelem = *melem * *xelem;
1211
1212 Index j = M.mcr.mextent; // --j (instead of j-- like in the outer loop)
1213 while (--j) // is important here because we only want
1214 { // mextent-1 cycles
1215 melem += M.mcr.mstride;
1216 xelem += x.mrange.mstride;
1217 *yelem += *melem * *xelem;
1218 }
1219
1220 mdata += M.mrr.mstride; // Jump to next matrix row
1221 yelem += y.mrange.mstride; // Jump to next element in target vector
1222 }
1223}
1224
1226
1250void mult(MatrixView A, const ConstMatrixView& B, const ConstMatrixView& C) {
1251 // Check dimensions:
1252 ARTS_ASSERT(A.nrows() == B.nrows());
1253 ARTS_ASSERT(A.ncols() == C.ncols());
1254 ARTS_ASSERT(B.ncols() == C.nrows());
1255
1256 // Catch trivial case if one of the matrices is empty.
1257 if ((B.nrows() == 0) || (B.ncols() == 0) || (C.ncols() == 0)) return;
1258
1259 // Matrices B and C must be continuous in at least on dimension, C
1260 // must be continuous along the second dimension.
1261 if (((B.mrr.get_stride() == 1) || (B.mcr.get_stride() == 1)) &&
1262 ((C.mrr.get_stride() == 1) || (C.mcr.get_stride() == 1)) &&
1263 (A.mcr.get_stride() == 1)) {
1264 // BLAS uses column-major order while arts uses row-major order.
1265 // Hence instead of C = A * B we compute C^T = A^T * B^T!
1266
1267 int k, m, n;
1268
1269 k = (int)B.ncols();
1270 m = (int)C.ncols();
1271 n = (int)B.nrows();
1272
1273 // Note also the clash in nomenclature: BLAS uses C = A * B while
1274 // arts uses A = B * C. Taking into accout this and the difference in
1275 // memory layouts, we need to map the MatrixViews A, B and C to the BLAS
1276 // arguments as follows:
1277 // A (arts) -> C (BLAS)
1278 // B (arts) -> B (BLAS)
1279 // C (arts) -> A (BLAS)
1280
1281 // Char indicating whether A (BLAS) or B (BLAS) should be transposed.
1282 char transa, transb;
1283 // Sizes of the matrices along the direction in which they are
1284 // traversed.
1285 int lda, ldb, ldc;
1286
1287 // Check if C (arts) is transposed.
1288 if (C.mrr.get_stride() == 1) {
1289 transa = 'T';
1290 lda = (int)C.mcr.get_stride();
1291 } else {
1292 transa = 'N';
1293 lda = (int)C.mrr.get_stride();
1294 }
1295
1296 // Check if B (arts) is transposed.
1297 if (B.mrr.get_stride() == 1) {
1298 transb = 'T';
1299 ldb = (int)B.mcr.get_stride();
1300 } else {
1301 transb = 'N';
1302 ldb = (int)B.mrr.get_stride();
1303 }
1304
1305 // In case B (arts) has only one column, column and row stride are 1.
1306 // We therefore need to set ldb to k, since dgemm_ requires lda to be at
1307 // least k / m if A is non-transposed / transposed.
1308 if ((B.mcr.get_stride() == 1) && (B.mrr.get_stride() == 1)) {
1309 transb = 'N';
1310 ldb = k;
1311 }
1312
1313 // The same holds for C (arts).
1314 if ((C.mcr.get_stride() == 1) && (C.mrr.get_stride() == 1)) {
1315 transa = 'N';
1316 lda = m;
1317 }
1318
1319 ldc = (int)A.mrr.get_stride();
1320 // The same holds for A (arts).
1321 if ((A.mcr.get_stride() == 1) && (A.mrr.get_stride() == 1)) {
1322 ldc = m;
1323 }
1324 double alpha = 1.0, beta = 0.0;
1325
1326 dgemm_(&transa,
1327 &transb,
1328 &m,
1329 &n,
1330 &k,
1331 &alpha,
1332 C.mdata + C.mrr.get_start() + C.mcr.get_start(),
1333 &lda,
1334 B.mdata + B.mrr.get_start() + B.mcr.get_start(),
1335 &ldb,
1336 &beta,
1337 A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1338 &ldc);
1339
1340 } else {
1341 mult_general(A, B, C);
1342 }
1343}
1344
1346
1355 const ConstMatrixView& B,
1356 const ConstMatrixView& C) ARTS_NOEXCEPT {
1357 // Check dimensions:
1358 ARTS_ASSERT(A.nrows() == B.nrows());
1359 ARTS_ASSERT(A.ncols() == C.ncols());
1360 ARTS_ASSERT(B.ncols() == C.nrows());
1361
1362 // Let's get the transpose of C, so that we can use 2D iterators to
1363 // access the columns (= rows of the transpose).
1364 ConstMatrixView CT = transpose(C);
1365
1366 const Iterator2D ae = A.end();
1367 Iterator2D ai = A.begin();
1368 ConstIterator2D bi = B.begin();
1369
1370 // This walks through the rows of A and B:
1371 for (; ai != ae; ++ai, ++bi) {
1372 const Iterator1D ace = ai->end();
1373 Iterator1D aci = ai->begin();
1374 ConstIterator2D cti = CT.begin();
1375
1376 // This walks through the columns of A with a 1D iterator, and
1377 // at the same time through the rows of CT, which are the columns of
1378 // C, with a 2D iterator:
1379 for (; aci != ace; ++aci, ++cti) {
1380 // The operator * is used to compute the scalar product
1381 // between rows of B and rows of C.transpose().
1382 *aci = (*bi) * (*cti);
1383 }
1384 }
1385}
1386
1388
1402 const ConstVectorView& a,
1404 ARTS_ASSERT(a.nelem() == 3);
1405 ARTS_ASSERT(b.nelem() == 3);
1406 ARTS_ASSERT(c.nelem() == 3);
1407
1408 c[0] = a[1] * b[2] - a[2] * b[1];
1409 c[1] = a[2] * b[0] - a[0] * b[2];
1410 c[2] = a[0] * b[1] - a[1] * b[0];
1411}
1412
1423 ARTS_ASSERT(a.nelem() == b.nelem());
1424 Numeric arg = (a * b) / sqrt(a * a) / sqrt(b * b);
1425
1426 // Due to numerical inaccuracy, arg might be slightly larger than 1.
1427 // We catch those cases to avoid spurious returns of NaNs
1428 return fabs(arg) > 1. ? 0. : acos(arg) * RAD2DEG;
1429}
1430
1445 ARTS_ASSERT(a.nelem() == b.nelem());
1446 ARTS_ASSERT(a.nelem() == c.nelem());
1447
1448 const Numeric C = (a * b) / (a * a);
1449 c = a;
1450 c *= C;
1451};
1452
1455 return ConstMatrixView(m.mdata, m.mcr, m.mrr);
1456}
1457
1461 return MatrixView(m.mdata, m.mcr, m.mrr);
1462}
1463
1484void transform(VectorView y, double (&my_func)(double), ConstVectorView x) {
1485 // Check dimensions:
1486 ARTS_ASSERT(y.nelem() == x.nelem());
1487
1488 const ConstIterator1D xe = x.end();
1489 ConstIterator1D xi = x.begin();
1490 Iterator1D yi = y.begin();
1491 for (; xi != xe; ++xi, ++yi) *yi = my_func(*xi);
1492}
1493
1512void transform(MatrixView y, double (&my_func)(double), ConstMatrixView x) {
1513 // Check dimensions:
1514 ARTS_ASSERT(y.nrows() == x.nrows());
1515 ARTS_ASSERT(y.ncols() == x.ncols());
1516
1517 const ConstIterator2D rxe = x.end();
1518 ConstIterator2D rx = x.begin();
1519 Iterator2D ry = y.begin();
1520 for (; rx != rxe; ++rx, ++ry) {
1521 const ConstIterator1D cxe = rx->end();
1522 ConstIterator1D cx = rx->begin();
1523 Iterator1D cy = ry->begin();
1524 for (; cx != cxe; ++cx, ++cy) *cy = my_func(*cx);
1525 }
1526}
1527
1530 // Initial value for max:
1531 Numeric max = x[0];
1532
1533 const ConstIterator1D xe = x.end();
1534 ConstIterator1D xi = x.begin();
1535
1536 for (; xi != xe; ++xi) {
1537 if (*xi > max) max = *xi;
1538 }
1539
1540 return max;
1541}
1542
1545 // Initial value for max:
1546 Numeric max = x(0, 0);
1547
1548 const ConstIterator2D rxe = x.end();
1549 ConstIterator2D rx = x.begin();
1550
1551 for (; rx != rxe; ++rx) {
1552 const ConstIterator1D cxe = rx->end();
1553 ConstIterator1D cx = rx->begin();
1554
1555 for (; cx != cxe; ++cx)
1556 if (*cx > max) max = *cx;
1557 }
1558
1559 return max;
1560}
1561
1564 // Initial value for min:
1565 Numeric min = x[0];
1566
1567 const ConstIterator1D xe = x.end();
1568 ConstIterator1D xi = x.begin();
1569
1570 for (; xi != xe; ++xi) {
1571 if (*xi < min) min = *xi;
1572 }
1573
1574 return min;
1575}
1576
1579 // Initial value for min:
1580 Numeric min = x(0, 0);
1581
1582 const ConstIterator2D rxe = x.end();
1583 ConstIterator2D rx = x.begin();
1584
1585 for (; rx != rxe; ++rx) {
1586 const ConstIterator1D cxe = rx->end();
1587 ConstIterator1D cx = rx->begin();
1588
1589 for (; cx != cxe; ++cx)
1590 if (*cx < min) min = *cx;
1591 }
1592
1593 return min;
1594}
1595
1598 // Initial value for mean:
1599 Numeric mean = 0;
1600
1601 const ConstIterator1D xe = x.end();
1602 ConstIterator1D xi = x.begin();
1603
1604 for (; xi != xe; ++xi) mean += *xi;
1605
1606 mean /= (Numeric)x.nelem();
1607
1608 return mean;
1609}
1610
1613 // Initial value for mean:
1614 Numeric nanmean = 0;
1615
1616 // Counter for the number of normal values
1617 Index numnormal = 0;
1618
1619 const ConstIterator1D xe = x.end();
1620 ConstIterator1D xi = x.begin();
1621
1622 for (; xi != xe; ++xi)
1623 if (std::isnormal(*xi) or (*xi == 0)) {
1624 nanmean += *xi;
1625 numnormal++;
1626 }
1627
1628 if (numnormal)
1629 nanmean /= Numeric(numnormal);
1630 else
1631 nanmean = std::numeric_limits<Numeric>::quiet_NaN();
1632
1633 return nanmean;
1634}
1635
1638 // Initial value for mean:
1639 Numeric mean = 0;
1640
1641 const ConstIterator2D rxe = x.end();
1642 ConstIterator2D rx = x.begin();
1643
1644 for (; rx != rxe; ++rx) {
1645 const ConstIterator1D cxe = rx->end();
1646 ConstIterator1D cx = rx->begin();
1647
1648 for (; cx != cxe; ++cx) mean += *cx;
1649 }
1650
1651 mean /= (Numeric)(x.nrows() * x.ncols());
1652
1653 return mean;
1654}
1655
1666 // cout << "Assigning VectorView from Array<Numeric>.\n";
1667
1668 // Check that sizes are compatible:
1669 ARTS_ASSERT(mrange.mextent == v.nelem());
1670
1671 // Iterators for Array:
1673 const Array<Numeric>::const_iterator e = v.end();
1674 // Iterator for Vector:
1675 Iterator1D target = begin();
1676
1677 for (; i != e; ++i, ++target) *target = *i;
1678
1679 return *this;
1680}
1681
1682// Const
1683
1684// Converts constant matrix to constant eigen map
1686 return ConstMatrixViewMap(A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1687 A.nrows(),
1688 A.ncols(),
1690}
1691
1692// Converts constant vector to constant eigen row-view
1695 A.nelem(),
1696 1,
1697 StrideType(A.mrange.get_stride(), 1));
1698}
1699
1700// Converts constant vector to constant eigen row-view
1702 return MapToEigen(A);
1703}
1704
1705// Converts constant vector to constant eigen column-view
1708 1,
1709 A.nelem(),
1710 StrideType(1, A.mrange.get_stride()));
1711}
1712
1713// Non- const
1714
1715// Converts matrix to eigen map
1717 return MatrixViewMap(A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1718 A.nrows(),
1719 A.ncols(),
1721}
1722
1723// Converts vector to eigen map row-view
1725 return MatrixViewMap(A.mdata + A.mrange.get_start(),
1726 A.nelem(),
1727 1,
1728 StrideType(A.mrange.get_stride(), 1));
1729}
1730
1731// Converts vector to eigen map row-view
1733
1734// Converts vector to eigen map column-view
1736 return MatrixViewMap(A.mdata + A.mrange.get_start(),
1737 1,
1738 A.nelem(),
1739 StrideType(1, A.mrange.get_stride()));
1740}
1741
1742// Special 4x4
1743
1744// Converts matrix to eigen map
1746 return Matrix4x4ViewMap(A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1747 4,
1748 4,
1750}
1751
1752// Converts constant matrix to constant eigen map
1754 return ConstMatrix4x4ViewMap(
1755 A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1756 4,
1757 4,
1759}
1760
1762// Helper function for debugging
1763#ifndef NDEBUG
1764
1779 return mv(r, c);
1780}
1781
1782#endif
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.
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:197
The constant iterator class for sub vectors.
Definition: matpackI.h:454
const Numeric * mx
Current position.
Definition: matpackI.h:501
Index mstride
Stride.
Definition: matpackI.h:503
The const row iterator class for sub matrices.
Definition: matpackI.h:857
A constant view of a Matrix.
Definition: matpackI.h:1050
Range mrr
The row range of mdata that is actually used.
Definition: matpackI.h:1153
friend class ConstVectorView
Definition: matpackI.h:1107
ConstIterator2D begin() const ARTS_NOEXCEPT
Return const iterator to first row.
Definition: matpackI.cc:469
Index size() const noexcept
Definition: matpackI.h:1065
Index nrows() const noexcept
Definition: matpackI.h:1061
ConstIterator2D end() const ARTS_NOEXCEPT
Return const iterator behind last row.
Definition: matpackI.cc:474
Index ncols() const noexcept
Definition: matpackI.h:1062
Numeric operator()(Index r, Index c) const ARTS_NOEXCEPT
Plain const index operator.
Definition: matpackI.h:1073
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1157
Range mcr
The column range of mdata that is actually used.
Definition: matpackI.h:1155
ConstVectorView diagonal() const ARTS_NOEXCEPT
Matrix diagonal as vector.
Definition: matpackI.cc:489
ConstMatrixView()=default
A constant view of a Vector.
Definition: matpackI.h:517
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
ConstVectorView()=default
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:644
Numeric sum() const ARTS_NOEXCEPT
The sum of all elements of a Vector.
Definition: matpackI.cc:48
Numeric operator[](Index n) const ARTS_NOEXCEPT
Plain const index operator.
Definition: matpackI.h:554
Range mrange
The range of mdata that is actually used.
Definition: matpackI.h:642
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:541
Index size() const noexcept
Definition: matpackI.h:544
The iterator class for sub vectors.
Definition: matpackI.h:395
Index mstride
Stride.
Definition: matpackI.h:449
Numeric * mx
Current position.
Definition: matpackI.h:447
The row iterator class for sub matrices.
Definition: matpackI.h:813
The Joker class.
Definition: matpackI.h:132
The MatrixView class.
Definition: matpackI.h:1169
MatrixView & operator=(const ConstMatrixView &m)
Assignment operator.
Definition: matpackI.cc:636
friend class VectorView
Definition: matpackI.h:1239
MatrixView & operator/=(Numeric x) ARTS_NOEXCEPT
Division by scalar.
Definition: matpackI.cc:703
Numeric & operator()(Index r, Index c) ARTS_NOEXCEPT
Plain index operator.
Definition: matpackI.h:1184
Iterator2D end() ARTS_NOEXCEPT
Return iterator behind last row.
Definition: matpackI.cc:626
MatrixView & operator+=(Numeric x) ARTS_NOEXCEPT
Addition of scalar.
Definition: matpackI.cc:713
Iterator2D begin() ARTS_NOEXCEPT
‍** Return const iterator to first row. Has to be redefined here, since it is
Definition: matpackI.cc:621
MatrixView & operator*=(Numeric x) ARTS_NOEXCEPT
Multiplication by scalar.
Definition: matpackI.cc:693
MatrixView & operator-=(Numeric x) ARTS_NOEXCEPT
Subtraction of scalar.
Definition: matpackI.cc:723
MatrixView()=default
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array.
Definition: matpackI.cc:738
The Matrix class.
Definition: matpackI.h:1270
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1063
Matrix()=default
virtual ~Matrix()
Destructor for Matrix.
Definition: matpackI.cc:1090
Matrix & operator=(const Matrix &m)
Assignment operator from another matrix.
Definition: matpackI.cc:1011
The range class.
Definition: matpackI.h:166
Index mstart
The start index.
Definition: matpackI.h:368
constexpr Index get_extent() const noexcept
Returns the extent of the range.
Definition: matpackI.h:349
Index mstride
The stride.
Definition: matpackI.h:377
constexpr Index get_start() const noexcept
Returns the start index of the range.
Definition: matpackI.h:347
Index mextent
The number of elements.
Definition: matpackI.h:375
constexpr Index get_stride() const noexcept
Returns the stride of the range.
Definition: matpackI.h:351
The VectorView class.
Definition: matpackI.h:658
VectorView operator+=(Numeric x) ARTS_NOEXCEPT
Addition of scalar.
Definition: matpackI.cc:203
VectorView()=default
VectorView operator-=(Numeric x) ARTS_NOEXCEPT
Subtraction of scalar.
Definition: matpackI.cc:209
VectorView & operator=(const ConstVectorView &v)
Assignment operator.
Definition: matpackI.cc:153
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array, const-version.
Definition: matpackI.cc:272
Numeric & operator[](Index n) ARTS_NOEXCEPT
Plain Index operator.
Definition: matpackI.h:679
Iterator1D begin() ARTS_NOEXCEPT
Return iterator to first element.
Definition: matpackI.cc:144
VectorView operator*=(Numeric x) ARTS_NOEXCEPT
Multiplication by scalar.
Definition: matpackI.cc:191
Iterator1D end() ARTS_NOEXCEPT
Return iterator behind last element.
Definition: matpackI.cc:148
VectorView operator/=(Numeric x) ARTS_NOEXCEPT
Division by scalar.
Definition: matpackI.cc:197
The Vector class.
Definition: matpackI.h:908
virtual ~Vector()
Destructor for Vector.
Definition: matpackI.cc:427
void resize(Index n)
Resize function.
Definition: matpackI.cc:411
Vector()=default
Vector & operator=(const Vector &v)
Assignment from another Vector.
Definition: matpackI.cc:381
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:1597
ConstMatrixViewMap MapToEigenRow(const ConstVectorView &A)
Definition: matpackI.cc:1701
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix-Vector Multiplication.
Definition: matpackI.cc:1131
Numeric vector_angle(ConstVectorView a, ConstVectorView b)
Definition: matpackI.cc:1422
Numeric max(const ConstVectorView &x) ARTS_NOEXCEPT
Max function, vector version.
Definition: matpackI.cc:1529
ConstMatrixViewMap MapToEigen(const ConstMatrixView &A)
Definition: matpackI.cc:1685
Numeric nanmean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version ignoring nans and infs.
Definition: matpackI.cc:1612
ConstMatrixViewMap MapToEigenCol(const ConstVectorView &A)
Definition: matpackI.cc:1706
void proj(Vector &c, ConstVectorView a, ConstVectorView b) ARTS_NOEXCEPT
Definition: matpackI.cc:1444
void mult_general(VectorView y, const ConstMatrixView &M, const ConstVectorView &x) ARTS_NOEXCEPT
Matrix Vector multiplication.
Definition: matpackI.cc:1189
const Joker joker
std::ostream & operator<<(std::ostream &os, const Range &r)
Definition: matpackI.cc:39
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1484
Matrix4x4ViewMap MapToEigen4x4(MatrixView &A)
Definition: matpackI.cc:1745
void swap(Vector &v1, Vector &v2)
Definition: matpackI.cc:422
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
cross3
Definition: matpackI.cc:1401
Numeric min(const ConstVectorView &x) ARTS_NOEXCEPT
Min function, vector version.
Definition: matpackI.cc:1563
void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Definition: matpackI.cc:304
Numeric debug_matrixview_get_elem(MatrixView &mv, Index r, Index c)
Helper function to access matrix elements.
Definition: matpackI.cc:1778
Numeric operator*(const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
Scalar product.
Definition: matpackI.cc:1099
ConstMatrixView transpose(ConstMatrixView m) ARTS_NOEXCEPT
Const version of transpose.
Definition: matpackI.cc:1454
Implementation of Matrix, Vector, and such stuff.
Eigen::Map< const Matrix4x4Type, 0, StrideType > ConstMatrix4x4ViewMap
Definition: matpackI.h:120
Eigen::Map< const MatrixType, 0, StrideType > ConstMatrixViewMap
Definition: matpackI.h:117
Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > StrideType
Definition: matpackI.h:113
Eigen::Map< MatrixType, 0, StrideType > MatrixViewMap
Definition: matpackI.h:116
Eigen::Map< Matrix4x4Type, 0, StrideType > Matrix4x4ViewMap
Definition: matpackI.h:119
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
Workspace & init(Workspace &ws)
Index nelem(const Lines &l)
Number of lines.
constexpr std::string_view Joker
Definition: isotopologues.h:13
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
#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