ARTS 2.5.4 (git: 31ce4f0e)
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
273 // Nothing to do here.
274}
275
277 : ConstVectorView(data, range) {
278 // Nothing to do here.
279}
280
282 const Range& p,
283 const Range& n) ARTS_NOEXCEPT
284 : ConstVectorView(data, p, n) {
285 // Nothing to do here.
286}
287
289 const ConstIterator1D& end,
290 Iterator1D target) {
291 if (origin.mstride == 1 && target.mstride == 1)
292 memcpy((void*)target.mx,
293 (void*)origin.mx,
294 sizeof(Numeric) * (end.mx - origin.mx));
295 else
296 for (; origin != end; ++origin, ++target) *target = *origin;
297}
298
300 for (; target != end; ++target) *target = x;
301}
302
303// Functions for Vector:
304// ---------------------
305
306Vector::Vector(std::initializer_list<Numeric> init)
307 : VectorView(new Numeric[init.size()], Range(0, init.size())) {
308 std::copy(init.begin(), init.end(), begin());
309}
310
312 // Nothing to do here.
313}
314
316 : VectorView(new Numeric[n], Range(0, n)) {
317 // Here we can access the raw memory directly, for slightly
318 // increased efficiency:
319 std::fill_n(mdata, n, fill);
320}
321
323 : VectorView(new Numeric[extent], Range(0, extent)) {
324 // Fill with values:
325 Numeric x = start;
326 Iterator1D i = begin();
327 const Iterator1D e = end();
328 for (; i != e; ++i) {
329 *i = x;
330 x += stride;
331 }
332}
333
335 : VectorView(new Numeric[v.nelem()], Range(0, v.nelem())) {
336 copy(v.begin(), v.end(), begin());
337}
338
340 : VectorView(new Numeric[v.nelem()], Range(0, v.nelem())) {
341 std::memcpy(mdata, v.mdata, nelem() * sizeof(Numeric));
342}
343
344Vector::Vector(const std::vector<Numeric>& v)
345 : VectorView(new Numeric[v.size()], Range(0, v.size())) {
346 std::vector<Numeric>::const_iterator vec_it_end = v.end();
347 Iterator1D this_it = this->begin();
348 for (std::vector<Numeric>::const_iterator vec_it = v.begin();
349 vec_it != vec_it_end;
350 ++vec_it, ++this_it)
351 *this_it = *vec_it;
352}
353
354Vector& Vector::operator=(std::initializer_list<Numeric> v) {
355 resize(v.size());
356 std::copy(v.begin(), v.end(), begin());
357 return *this;
358}
359
361 if (this != &v) {
362 resize(v.nelem());
363 std::memcpy(mdata, v.mdata, nelem() * sizeof(Numeric));
364 }
365 return *this;
366}
367
369 if (this != &v) {
370 delete[] mdata;
371 mdata = v.mdata;
372 mrange = v.mrange;
373 v.mrange = Range(0, 0);
374 v.mdata = nullptr;
375 }
376 return *this;
377}
378
380 resize(x.nelem());
382 return *this;
383}
384
386 std::fill_n(mdata, nelem(), x);
387 return *this;
388}
389
391 ARTS_ASSERT(0 <= n);
392 if (mrange.mextent != n) {
393 delete[] mdata;
394 mdata = new Numeric[n];
395 mrange.mstart = 0;
396 mrange.mextent = n;
397 mrange.mstride = 1;
398 }
399}
400
401void swap(Vector& v1, Vector& v2) {
402 std::swap(v1.mrange, v2.mrange);
403 std::swap(v1.mdata, v2.mdata);
404}
405
406Vector::~Vector() { delete[] mdata; }
407
408// Functions for ConstMatrixView:
409// ------------------------------
410
415 const Range& r, const Range& c) const ARTS_NOEXCEPT {
416 return ConstMatrixView(mdata, mrr, mcr, r, c);
417}
418
425 Index c) const ARTS_NOEXCEPT {
426 // Check that c is valid:
427 ARTS_ASSERT(0 <= c);
428 ARTS_ASSERT(c < mcr.mextent);
429
430 return ConstVectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
431}
432
440 // Check that r is valid:
441 ARTS_ASSERT(0 <= r);
442 ARTS_ASSERT(r < mrr.mextent);
443
444 return ConstVectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
445}
446
450}
451
454 return ConstIterator2D(
456 mrr.mstride);
457}
458
460
471 Range(0, n, mrr.mstride + mcr.mstride));
472}
473
478 const Range& rr,
479 const Range& cr) ARTS_NOEXCEPT : mrr(rr),
480 mcr(cr),
481 mdata(data) {
482 // Nothing to do here.
483}
484
500 const Range& pr,
501 const Range& pc,
502 const Range& nr,
503 const Range& nc) ARTS_NOEXCEPT : mrr(pr, nr),
504 mcr(pc, nc),
505 mdata(data) {
506 // Nothing to do here.
507}
508
516std::ostream& operator<<(std::ostream& os, const ConstMatrixView& v) {
517 // Row iterators:
518 ConstIterator2D ir = v.begin();
519 const ConstIterator2D end_row = v.end();
520
521 if (ir != end_row) {
522 ConstIterator1D ic = ir->begin();
523 const ConstIterator1D end_col = ir->end();
524
525 if (ic != end_col) {
526 os << setw(3) << *ic;
527 ++ic;
528 }
529 for (; ic != end_col; ++ic) {
530 os << " " << setw(3) << *ic;
531 }
532 ++ir;
533 }
534 for (; ir != end_row; ++ir) {
535 ConstIterator1D ic = ir->begin();
536 const ConstIterator1D end_col = ir->end();
537
538 os << "\n";
539 if (ic != end_col) {
540 os << setw(3) << *ic;
541 ++ic;
542 }
543 for (; ic != end_col; ++ic) {
544 os << " " << setw(3) << *ic;
545 }
546 }
547
548 return os;
549}
550
551// Functions for MatrixView:
552// -------------------------
553
558 const Range& c) ARTS_NOEXCEPT {
559 return MatrixView(mdata, mrr, mcr, r, c);
560}
561
567 // Check that c is valid:
568 ARTS_ASSERT(0 <= c);
569 ARTS_ASSERT(c < mcr.mextent);
570
571 return VectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
572}
573
579 // Check that r is valid:
580 ARTS_ASSERT(0 <= r);
581 ARTS_ASSERT(r < mrr.mextent);
582
583 return VectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
584}
585
587// hiden by the non-const operator of the derived class.*/
588//ConstIterator2D MatrixView::begin() const
589//{
590// return ConstMatrixView::begin();
591//}
592//
594//ConstIterator2D MatrixView::end() const
595//{
596// return ConstMatrixView::end();
597//}
598
602}
603
606 return Iterator2D(
608 mrr.mstride);
609}
610
616 // Check that sizes are compatible:
619
620 copy(m.begin(), m.end(), begin());
621 return *this;
622}
623
630 // Check that sizes are compatible:
633
634 copy(m.begin(), m.end(), begin());
635 return *this;
636}
637
642 // Check that sizes are compatible:
645
646 copy(m.begin(), m.end(), begin());
647 return *this;
648}
649
655 // Check that sizes are compatible:
656 ARTS_ASSERT(mrr.mextent == v.nelem());
657 ARTS_ASSERT(mcr.mextent == 1);
658 // dummy = ConstMatrixView(v.mdata,v.mrange,Range(v.mrange.mstart,1));;
659 ConstMatrixView dummy(v);
660 copy(dummy.begin(), dummy.end(), begin());
661 return *this;
662}
663
667 copy(x, begin(), end());
668 return *this;
669}
670
673 const Iterator2D er = end();
674 for (Iterator2D r = begin(); r != er; ++r) {
675 const Iterator1D ec = r->end();
676 for (Iterator1D c = r->begin(); c != ec; ++c) *c *= x;
677 }
678 return *this;
679}
680
683 const Iterator2D er = end();
684 for (Iterator2D r = begin(); r != er; ++r) {
685 const Iterator1D ec = r->end();
686 for (Iterator1D c = r->begin(); c != ec; ++c) *c /= x;
687 }
688 return *this;
689}
690
693 const Iterator2D er = end();
694 for (Iterator2D r = begin(); r != er; ++r) {
695 const Iterator1D ec = r->end();
696 for (Iterator1D c = r->begin(); c != ec; ++c) *c += x;
697 }
698 return *this;
699}
700
703 const Iterator2D er = end();
704 for (Iterator2D r = begin(); r != er; ++r) {
705 const Iterator1D ec = r->end();
706 for (Iterator1D c = r->begin(); c != ec; ++c) *c -= x;
707 }
708 return *this;
709}
710
713 ARTS_ASSERT(nrows() == x.nrows());
714 ARTS_ASSERT(ncols() == x.ncols());
715 ConstIterator2D sr = x.begin();
716 Iterator2D r = begin();
717 const Iterator2D er = end();
718 for (; r != er; ++r, ++sr) {
719 ConstIterator1D sc = sr->begin();
720 Iterator1D c = r->begin();
721 const Iterator1D ec = r->end();
722 for (; c != ec; ++c, ++sc) *c *= *sc;
723 }
724 return *this;
725}
726
729 ARTS_ASSERT(nrows() == x.nrows());
730 ARTS_ASSERT(ncols() == x.ncols());
731 ConstIterator2D sr = x.begin();
732 Iterator2D r = begin();
733 const Iterator2D er = end();
734 for (; r != er; ++r, ++sr) {
735 ConstIterator1D sc = sr->begin();
736 Iterator1D c = r->begin();
737 const Iterator1D ec = r->end();
738 for (; c != ec; ++c, ++sc) *c /= *sc;
739 }
740 return *this;
741}
742
745 ARTS_ASSERT(nrows() == x.nrows());
746 ARTS_ASSERT(ncols() == x.ncols());
747 ConstIterator2D sr = x.begin();
748 Iterator2D r = begin();
749 const Iterator2D er = end();
750 for (; r != er; ++r, ++sr) {
751 ConstIterator1D sc = sr->begin();
752 Iterator1D c = r->begin();
753 const Iterator1D ec = r->end();
754 for (; c != ec; ++c, ++sc) *c += *sc;
755 }
756 return *this;
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.nelem());
778 ARTS_ASSERT(ncols() == 1);
779 ConstIterator1D sc = x.begin();
780 Iterator2D r = begin();
781 const Iterator2D er = end();
782 for (; r != er; ++r, ++sc) {
783 Iterator1D c = r->begin();
784 *c *= *sc;
785 }
786 return *this;
787}
788
791 ARTS_ASSERT(nrows() == x.nelem());
792 ARTS_ASSERT(ncols() == 1);
793 ConstIterator1D sc = x.begin();
794 Iterator2D r = begin();
795 const Iterator2D er = end();
796 for (; r != er; ++r, ++sc) {
797 Iterator1D c = r->begin();
798 *c /= *sc;
799 }
800 return *this;
801}
802
805 ARTS_ASSERT(nrows() == x.nelem());
806 ARTS_ASSERT(ncols() == 1);
807 ConstIterator1D sc = x.begin();
808 Iterator2D r = begin();
809 const Iterator2D er = end();
810 for (; r != er; ++r, ++sc) {
811 Iterator1D c = r->begin();
812 *c += *sc;
813 }
814 return *this;
815}
816
819 ARTS_ASSERT(nrows() == x.nelem());
820 ARTS_ASSERT(ncols() == 1);
821 ConstIterator1D sc = x.begin();
822 Iterator2D r = begin();
823 const Iterator2D er = end();
824 for (; r != er; ++r, ++sc) {
825 Iterator1D c = r->begin();
826 *c -= *sc;
827 }
828 return *this;
829}
830
835 const Range& rr,
836 const Range& cr) ARTS_NOEXCEPT
837 : ConstMatrixView(data, rr, cr) {
838 // Nothing to do here.
839}
840
856 const Range& pr,
857 const Range& pc,
858 const Range& nr,
859 const Range& nc) ARTS_NOEXCEPT
860 : ConstMatrixView(data, pr, pc, nr, nc) {
861 // Nothing to do here.
862}
863
874 const ConstIterator2D& end,
875 Iterator2D target) {
876 for (; origin != end; ++origin, ++target) {
877 ConstIterator1D o = origin->begin();
878 const ConstIterator1D e = origin->end();
879 Iterator1D t = target->begin();
880 for (; o != e; ++o, ++t) *t = *o;
881 }
882}
883
886 for (; target != end; ++target) {
887 Iterator1D t = target->begin();
888 const Iterator1D e = target->end();
889 for (; t != e; ++t) *t = x;
890 }
891}
892
893// Functions for Matrix:
894// ---------------------
895
899 : MatrixView(new Numeric[r * c], Range(0, r, c), Range(0, c)) {
900 // Nothing to do here.
901}
902
905 : MatrixView(new Numeric[r * c], Range(0, r, c), Range(0, c)) {
906 // Here we can access the raw memory directly, for slightly
907 // increased efficiency:
908 std::fill_n(mdata, r * c, fill);
909}
910
914 : MatrixView(new Numeric[m.nrows() * m.ncols()],
915 Range(0, m.nrows(), m.ncols()),
916 Range(0, m.ncols())) {
917 copy(m.begin(), m.end(), begin());
918}
919
923 : MatrixView(new Numeric[m.nrows() * m.ncols()],
924 Range(0, m.nrows(), m.ncols()),
925 Range(0, m.ncols())) {
926 // There is a catch here: If m is an empty matrix, then it will have
927 // 0 colunns. But this is used to initialize the stride of the row
928 // Range! Thus, this method has to be consistent with the behaviour
929 // of Range::Range. For now, Range::Range allows also stride 0.
930 std::memcpy(mdata, m.mdata, nrows() * ncols() * sizeof(Numeric));
931}
932
934
959 if (this != &m) {
960 resize(m.nrows(), m.ncols());
961 std::memcpy(mdata, m.mdata, nrows() * ncols() * sizeof(Numeric));
962 }
963 return *this;
964}
965
968 if (this != &m) {
969 delete[] mdata;
970 mdata = m.mdata;
971 mrr = m.mrr;
972 mcr = m.mcr;
973 m.mrr = Range(0, 0);
974 m.mcr = Range(0, 0);
975 m.mdata = nullptr;
976 }
977 return *this;
978}
979
983 std::fill_n(mdata, nrows() * ncols(), x);
984 return *this;
985}
986
988
1001 resize(v.nelem(), 1);
1002 ConstMatrixView dummy(v);
1003 copy(dummy.begin(), dummy.end(), begin());
1004 return *this;
1005}
1006
1011 ARTS_ASSERT(0 <= r);
1012 ARTS_ASSERT(0 <= c);
1013
1014 if (mrr.mextent != r || mcr.mextent != c) {
1015 delete[] mdata;
1016 mdata = new Numeric[r * c];
1017
1018 mrr.mstart = 0;
1019 mrr.mextent = r;
1020 mrr.mstride = c;
1021
1022 mcr.mstart = 0;
1023 mcr.mextent = c;
1024 mcr.mstride = 1;
1025 }
1026}
1027
1029void swap(Matrix& m1, Matrix& m2) {
1030 std::swap(m1.mrr, m2.mrr);
1031 std::swap(m1.mcr, m2.mcr);
1032 std::swap(m1.mdata, m2.mdata);
1033}
1034
1038 // cout << "Destroying a Matrix:\n"
1039 // << *this << "\n........................................\n";
1040 delete[] mdata;
1041}
1042
1043// Some general Matrix Vector functions:
1044
1048 // Check dimensions:
1049 ARTS_ASSERT(a.nelem() == b.nelem());
1050
1051 const ConstIterator1D ae = a.end();
1052 ConstIterator1D ai = a.begin();
1053 ConstIterator1D bi = b.begin();
1054
1055 Numeric res = 0;
1056 for (; ai != ae; ++ai, ++bi) res += (*ai) * (*bi);
1057
1058 return res;
1059}
1060
1062
1079 ARTS_ASSERT(y.mrange.get_extent() == M.mrr.get_extent());
1080 ARTS_ASSERT(M.mcr.get_extent() == x.mrange.get_extent());
1081 ARTS_ASSERT((M.mcr.get_extent() != 0) && (M.mrr.get_extent() != 0));
1082
1083 if ((M.mcr.get_stride() == 1) || (M.mrr.get_stride() == 1)) {
1084 char trans;
1085 int m, n;
1086 double zero = 0.0;
1087 double one = 1.0;
1088 int LDA, incx, incy;
1089
1090 if (M.mcr.get_stride() != 1) {
1091 trans = 'n';
1092 m = (int)M.mrr.get_extent();
1093 n = (int)M.mcr.get_extent();
1094 LDA = (int)M.mcr.get_stride();
1095 } else {
1096 trans = 't';
1097 m = (int)M.mcr.get_extent();
1098 n = (int)M.mrr.get_extent();
1099 LDA = (int)M.mrr.get_stride();
1100 if (M.mrr.get_stride() == 1) LDA = m;
1101 }
1102
1103 incx = (int)x.mrange.get_stride();
1104 incy = (int)y.mrange.get_stride();
1105
1106 double* mstart = M.mdata + M.mcr.get_start() + M.mrr.get_start();
1107 double* ystart = y.mdata + y.mrange.get_start();
1108 double* xstart = x.mdata + x.mrange.get_start();
1109
1110 dgemv_(&trans,
1111 &m,
1112 &n,
1113 &one,
1114 mstart,
1115 &LDA,
1116 xstart,
1117 &incx,
1118 &zero,
1119 ystart,
1120 &incy);
1121
1122 } else {
1123 mult_general(y, M, x);
1124 }
1125}
1126
1137 const ConstMatrixView& M,
1138 const ConstVectorView& x) ARTS_NOEXCEPT {
1139 // Check dimensions:
1140 ARTS_ASSERT(y.mrange.mextent == M.mrr.mextent);
1141 ARTS_ASSERT(M.mcr.mextent == x.mrange.mextent);
1142 ARTS_ASSERT(M.mcr.mextent != 0 && M.mrr.mextent != 0);
1143
1144 // Let's first find the pointers to the starting positions
1145 Numeric* mdata = M.mdata + M.mcr.mstart + M.mrr.mstart;
1146 Numeric* xdata = x.mdata + x.mrange.mstart;
1147 Numeric* yelem = y.mdata + y.mrange.mstart;
1148
1149 Index i = M.mrr.mextent;
1150 while (i--) {
1151 Numeric* melem = mdata;
1152 Numeric* xelem = xdata; // Reset xelem to first element of source vector
1153
1154 // Multiply first element of matrix row with first element of source
1155 // vector. This is done outside the loop to avoid initialization of the
1156 // target vector's element with zero (saves one assignment)
1157 *yelem = *melem * *xelem;
1158
1159 Index j = M.mcr.mextent; // --j (instead of j-- like in the outer loop)
1160 while (--j) // is important here because we only want
1161 { // mextent-1 cycles
1162 melem += M.mcr.mstride;
1163 xelem += x.mrange.mstride;
1164 *yelem += *melem * *xelem;
1165 }
1166
1167 mdata += M.mrr.mstride; // Jump to next matrix row
1168 yelem += y.mrange.mstride; // Jump to next element in target vector
1169 }
1170}
1171
1173
1197void mult(MatrixView A, const ConstMatrixView& B, const ConstMatrixView& C) {
1198 // Check dimensions:
1199 ARTS_ASSERT(A.nrows() == B.nrows());
1200 ARTS_ASSERT(A.ncols() == C.ncols());
1201 ARTS_ASSERT(B.ncols() == C.nrows());
1202
1203 // Catch trivial case if one of the matrices is empty.
1204 if ((B.nrows() == 0) || (B.ncols() == 0) || (C.ncols() == 0)) return;
1205
1206 // Matrices B and C must be continuous in at least on dimension, C
1207 // must be continuous along the second dimension.
1208 if (((B.mrr.get_stride() == 1) || (B.mcr.get_stride() == 1)) &&
1209 ((C.mrr.get_stride() == 1) || (C.mcr.get_stride() == 1)) &&
1210 (A.mcr.get_stride() == 1)) {
1211 // BLAS uses column-major order while arts uses row-major order.
1212 // Hence instead of C = A * B we compute C^T = A^T * B^T!
1213
1214 int k, m, n;
1215
1216 k = (int)B.ncols();
1217 m = (int)C.ncols();
1218 n = (int)B.nrows();
1219
1220 // Note also the clash in nomenclature: BLAS uses C = A * B while
1221 // arts uses A = B * C. Taking into accout this and the difference in
1222 // memory layouts, we need to map the MatrixViews A, B and C to the BLAS
1223 // arguments as follows:
1224 // A (arts) -> C (BLAS)
1225 // B (arts) -> B (BLAS)
1226 // C (arts) -> A (BLAS)
1227
1228 // Char indicating whether A (BLAS) or B (BLAS) should be transposed.
1229 char transa, transb;
1230 // Sizes of the matrices along the direction in which they are
1231 // traversed.
1232 int lda, ldb, ldc;
1233
1234 // Check if C (arts) is transposed.
1235 if (C.mrr.get_stride() == 1) {
1236 transa = 'T';
1237 lda = (int)C.mcr.get_stride();
1238 } else {
1239 transa = 'N';
1240 lda = (int)C.mrr.get_stride();
1241 }
1242
1243 // Check if B (arts) is transposed.
1244 if (B.mrr.get_stride() == 1) {
1245 transb = 'T';
1246 ldb = (int)B.mcr.get_stride();
1247 } else {
1248 transb = 'N';
1249 ldb = (int)B.mrr.get_stride();
1250 }
1251
1252 // In case B (arts) has only one column, column and row stride are 1.
1253 // We therefore need to set ldb to k, since dgemm_ requires lda to be at
1254 // least k / m if A is non-transposed / transposed.
1255 if ((B.mcr.get_stride() == 1) && (B.mrr.get_stride() == 1)) {
1256 transb = 'N';
1257 ldb = k;
1258 }
1259
1260 // The same holds for C (arts).
1261 if ((C.mcr.get_stride() == 1) && (C.mrr.get_stride() == 1)) {
1262 transa = 'N';
1263 lda = m;
1264 }
1265
1266 ldc = (int)A.mrr.get_stride();
1267 // The same holds for A (arts).
1268 if ((A.mcr.get_stride() == 1) && (A.mrr.get_stride() == 1)) {
1269 ldc = m;
1270 }
1271 double alpha = 1.0, beta = 0.0;
1272
1273 dgemm_(&transa,
1274 &transb,
1275 &m,
1276 &n,
1277 &k,
1278 &alpha,
1279 C.mdata + C.mrr.get_start() + C.mcr.get_start(),
1280 &lda,
1281 B.mdata + B.mrr.get_start() + B.mcr.get_start(),
1282 &ldb,
1283 &beta,
1284 A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1285 &ldc);
1286
1287 } else {
1288 mult_general(A, B, C);
1289 }
1290}
1291
1293
1302 const ConstMatrixView& B,
1303 const ConstMatrixView& C) ARTS_NOEXCEPT {
1304 // Check dimensions:
1305 ARTS_ASSERT(A.nrows() == B.nrows());
1306 ARTS_ASSERT(A.ncols() == C.ncols());
1307 ARTS_ASSERT(B.ncols() == C.nrows());
1308
1309 // Let's get the transpose of C, so that we can use 2D iterators to
1310 // access the columns (= rows of the transpose).
1311 ConstMatrixView CT = transpose(C);
1312
1313 const Iterator2D ae = A.end();
1314 Iterator2D ai = A.begin();
1315 ConstIterator2D bi = B.begin();
1316
1317 // This walks through the rows of A and B:
1318 for (; ai != ae; ++ai, ++bi) {
1319 const Iterator1D ace = ai->end();
1320 Iterator1D aci = ai->begin();
1321 ConstIterator2D cti = CT.begin();
1322
1323 // This walks through the columns of A with a 1D iterator, and
1324 // at the same time through the rows of CT, which are the columns of
1325 // C, with a 2D iterator:
1326 for (; aci != ace; ++aci, ++cti) {
1327 // The operator * is used to compute the scalar product
1328 // between rows of B and rows of C.transpose().
1329 *aci = (*bi) * (*cti);
1330 }
1331 }
1332}
1333
1335
1349 const ConstVectorView& a,
1351 ARTS_ASSERT(a.nelem() == 3);
1352 ARTS_ASSERT(b.nelem() == 3);
1353 ARTS_ASSERT(c.nelem() == 3);
1354
1355 c[0] = a[1] * b[2] - a[2] * b[1];
1356 c[1] = a[2] * b[0] - a[0] * b[2];
1357 c[2] = a[0] * b[1] - a[1] * b[0];
1358}
1359
1370 ARTS_ASSERT(a.nelem() == b.nelem());
1371 Numeric arg = (a * b) / sqrt(a * a) / sqrt(b * b);
1372
1373 // Due to numerical inaccuracy, arg might be slightly larger than 1.
1374 // We catch those cases to avoid spurious returns of NaNs
1375 return fabs(arg) > 1. ? 0. : acos(arg) * RAD2DEG;
1376}
1377
1392 ARTS_ASSERT(a.nelem() == b.nelem());
1393 ARTS_ASSERT(a.nelem() == c.nelem());
1394
1395 const Numeric C = (a * b) / (a * a);
1396 c = a;
1397 c *= C;
1398};
1399
1402 return ConstMatrixView(m.mdata, m.mcr, m.mrr);
1403}
1404
1408 return MatrixView(m.mdata, m.mcr, m.mrr);
1409}
1410
1431void transform(VectorView y, double (&my_func)(double), ConstVectorView x) {
1432 // Check dimensions:
1433 ARTS_ASSERT(y.nelem() == x.nelem());
1434
1435 const ConstIterator1D xe = x.end();
1436 ConstIterator1D xi = x.begin();
1437 Iterator1D yi = y.begin();
1438 for (; xi != xe; ++xi, ++yi) *yi = my_func(*xi);
1439}
1440
1459void transform(MatrixView y, double (&my_func)(double), ConstMatrixView x) {
1460 // Check dimensions:
1461 ARTS_ASSERT(y.nrows() == x.nrows());
1462 ARTS_ASSERT(y.ncols() == x.ncols());
1463
1464 const ConstIterator2D rxe = x.end();
1465 ConstIterator2D rx = x.begin();
1466 Iterator2D ry = y.begin();
1467 for (; rx != rxe; ++rx, ++ry) {
1468 const ConstIterator1D cxe = rx->end();
1469 ConstIterator1D cx = rx->begin();
1470 Iterator1D cy = ry->begin();
1471 for (; cx != cxe; ++cx, ++cy) *cy = my_func(*cx);
1472 }
1473}
1474
1477 // Initial value for max:
1478 Numeric max = x[0];
1479
1480 const ConstIterator1D xe = x.end();
1481 ConstIterator1D xi = x.begin();
1482
1483 for (; xi != xe; ++xi) {
1484 if (*xi > max) max = *xi;
1485 }
1486
1487 return max;
1488}
1489
1492 // Initial value for max:
1493 Numeric max = x(0, 0);
1494
1495 const ConstIterator2D rxe = x.end();
1496 ConstIterator2D rx = x.begin();
1497
1498 for (; rx != rxe; ++rx) {
1499 const ConstIterator1D cxe = rx->end();
1500 ConstIterator1D cx = rx->begin();
1501
1502 for (; cx != cxe; ++cx)
1503 if (*cx > max) max = *cx;
1504 }
1505
1506 return max;
1507}
1508
1511 // Initial value for min:
1512 Numeric min = x[0];
1513
1514 const ConstIterator1D xe = x.end();
1515 ConstIterator1D xi = x.begin();
1516
1517 for (; xi != xe; ++xi) {
1518 if (*xi < min) min = *xi;
1519 }
1520
1521 return min;
1522}
1523
1526 // Initial value for min:
1527 Numeric min = x(0, 0);
1528
1529 const ConstIterator2D rxe = x.end();
1530 ConstIterator2D rx = x.begin();
1531
1532 for (; rx != rxe; ++rx) {
1533 const ConstIterator1D cxe = rx->end();
1534 ConstIterator1D cx = rx->begin();
1535
1536 for (; cx != cxe; ++cx)
1537 if (*cx < min) min = *cx;
1538 }
1539
1540 return min;
1541}
1542
1545 // Initial value for mean:
1546 Numeric mean = 0;
1547
1548 const ConstIterator1D xe = x.end();
1549 ConstIterator1D xi = x.begin();
1550
1551 for (; xi != xe; ++xi) mean += *xi;
1552
1553 mean /= (Numeric)x.nelem();
1554
1555 return mean;
1556}
1557
1560 // Initial value for mean:
1561 Numeric nanmean = 0;
1562
1563 // Counter for the number of normal values
1564 Index numnormal = 0;
1565
1566 const ConstIterator1D xe = x.end();
1567 ConstIterator1D xi = x.begin();
1568
1569 for (; xi != xe; ++xi)
1570 if (std::isnormal(*xi) or (*xi == 0)) {
1571 nanmean += *xi;
1572 numnormal++;
1573 }
1574
1575 if (numnormal)
1576 nanmean /= Numeric(numnormal);
1577 else
1578 nanmean = std::numeric_limits<Numeric>::quiet_NaN();
1579
1580 return nanmean;
1581}
1582
1585 // Initial value for mean:
1586 Numeric mean = 0;
1587
1588 const ConstIterator2D rxe = x.end();
1589 ConstIterator2D rx = x.begin();
1590
1591 for (; rx != rxe; ++rx) {
1592 const ConstIterator1D cxe = rx->end();
1593 ConstIterator1D cx = rx->begin();
1594
1595 for (; cx != cxe; ++cx) mean += *cx;
1596 }
1597
1598 mean /= (Numeric)(x.nrows() * x.ncols());
1599
1600 return mean;
1601}
1602
1613 // cout << "Assigning VectorView from Array<Numeric>.\n";
1614
1615 // Check that sizes are compatible:
1616 ARTS_ASSERT(mrange.mextent == v.nelem());
1617
1618 // Iterators for Array:
1620 const Array<Numeric>::const_iterator e = v.end();
1621 // Iterator for Vector:
1622 Iterator1D target = begin();
1623
1624 for (; i != e; ++i, ++target) *target = *i;
1625
1626 return *this;
1627}
1628
1630// Helper function for debugging
1631#ifndef NDEBUG
1632
1647 return mv(r, c);
1648}
1649
1650#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
Definition: array.h:92
The constant iterator class for sub vectors.
Definition: matpackI.h:449
const Numeric * mx
Current position.
Definition: matpackI.h:496
Index mstride
Stride.
Definition: matpackI.h:498
The const row iterator class for sub matrices.
Definition: matpackI.h:848
A constant view of a Matrix.
Definition: matpackI.h:1043
Range mrr
The row range of mdata that is actually used.
Definition: matpackI.h:1148
friend class ConstVectorView
Definition: matpackI.h:1103
ConstIterator2D begin() const ARTS_NOEXCEPT
Return const iterator to first row.
Definition: matpackI.cc:448
Index nrows() const noexcept
Definition: matpackI.h:1055
ConstIterator2D end() const ARTS_NOEXCEPT
Return const iterator behind last row.
Definition: matpackI.cc:453
Index ncols() const noexcept
Definition: matpackI.h:1056
Numeric operator()(Index r, Index c) const ARTS_NOEXCEPT
Plain const index operator.
Definition: matpackI.h:1069
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1152
Range mcr
The column range of mdata that is actually used.
Definition: matpackI.h:1150
ConstVectorView diagonal() const ARTS_NOEXCEPT
Matrix diagonal as vector.
Definition: matpackI.cc:468
ConstMatrixView()=default
A constant view of a Vector.
Definition: matpackI.h:512
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:649
Numeric sum() const ARTS_NOEXCEPT
The sum of all elements of a Vector.
Definition: matpackI.cc:48
Numeric operator[](Index n) const ARTS_NOEXCEPT
Plain const index operator.
Definition: matpackI.h:549
Range mrange
The range of mdata that is actually used.
Definition: matpackI.h:647
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The iterator class for sub vectors.
Definition: matpackI.h:390
Index mstride
Stride.
Definition: matpackI.h:444
Numeric * mx
Current position.
Definition: matpackI.h:442
The row iterator class for sub matrices.
Definition: matpackI.h:804
The Joker class.
Definition: matpackI.h:125
The MatrixView class.
Definition: matpackI.h:1164
MatrixView & operator=(const ConstMatrixView &m)
Assignment operator.
Definition: matpackI.cc:615
friend class VectorView
Definition: matpackI.h:1230
MatrixView & operator/=(Numeric x) ARTS_NOEXCEPT
Division by scalar.
Definition: matpackI.cc:682
Numeric & operator()(Index r, Index c) ARTS_NOEXCEPT
Plain index operator.
Definition: matpackI.h:1179
Iterator2D end() ARTS_NOEXCEPT
Return iterator behind last row.
Definition: matpackI.cc:605
MatrixView & operator+=(Numeric x) ARTS_NOEXCEPT
Addition of scalar.
Definition: matpackI.cc:692
Iterator2D begin() ARTS_NOEXCEPT
‍** Return const iterator to first row. Has to be redefined here, since it is
Definition: matpackI.cc:600
MatrixView & operator*=(Numeric x) ARTS_NOEXCEPT
Multiplication by scalar.
Definition: matpackI.cc:672
MatrixView & operator-=(Numeric x) ARTS_NOEXCEPT
Subtraction of scalar.
Definition: matpackI.cc:702
MatrixView()=default
The Matrix class.
Definition: matpackI.h:1261
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1010
Matrix()=default
virtual ~Matrix()
Destructor for Matrix.
Definition: matpackI.cc:1037
Matrix & operator=(const Matrix &m)
Assignment operator from another matrix.
Definition: matpackI.cc:958
The range class.
Definition: matpackI.h:159
Index mstart
The start index.
Definition: matpackI.h:367
constexpr Index get_extent() const noexcept
Returns the extent of the range.
Definition: matpackI.h:342
Index mstride
The stride.
Definition: matpackI.h:372
constexpr Index get_start() const noexcept
Returns the start index of the range.
Definition: matpackI.h:340
Index mextent
The number of elements.
Definition: matpackI.h:370
constexpr Index get_stride() const noexcept
Returns the stride of the range.
Definition: matpackI.h:344
The VectorView class.
Definition: matpackI.h:663
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
Numeric & operator[](Index n) ARTS_NOEXCEPT
Plain Index operator.
Definition: matpackI.h:684
Iterator1D begin() ARTS_NOEXCEPT
Return iterator to first element.
Definition: matpackI.cc:144
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:899
virtual ~Vector()
Destructor for Vector.
Definition: matpackI.cc:406
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
Vector()=default
Vector & operator=(const Vector &v)
Assignment from another Vector.
Definition: matpackI.cc:360
#define ARTS_NOEXCEPT
Definition: debug.h:80
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
constexpr Numeric RAD2DEG
Definition: doit.cc:60
The declarations of all the exception classes.
#define beta
Numeric mean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version.
Definition: matpackI.cc:1544
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix-Vector Multiplication.
Definition: matpackI.cc:1078
Numeric vector_angle(ConstVectorView a, ConstVectorView b)
Definition: matpackI.cc:1369
Numeric max(const ConstVectorView &x) ARTS_NOEXCEPT
Max function, vector version.
Definition: matpackI.cc:1476
Numeric nanmean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version ignoring nans and infs.
Definition: matpackI.cc:1559
void proj(Vector &c, ConstVectorView a, ConstVectorView b) ARTS_NOEXCEPT
Definition: matpackI.cc:1391
void mult_general(VectorView y, const ConstMatrixView &M, const ConstVectorView &x) ARTS_NOEXCEPT
Matrix Vector multiplication.
Definition: matpackI.cc:1136
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:1431
void swap(Vector &v1, Vector &v2)
Definition: matpackI.cc:401
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
cross3
Definition: matpackI.cc:1348
Numeric min(const ConstVectorView &x) ARTS_NOEXCEPT
Min function, vector version.
Definition: matpackI.cc:1510
void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Definition: matpackI.cc:288
Numeric debug_matrixview_get_elem(MatrixView &mv, Index r, Index c)
Helper function to access matrix elements.
Definition: matpackI.cc:1646
Numeric operator*(const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
Scalar product.
Definition: matpackI.cc:1046
ConstMatrixView transpose(ConstMatrixView m) ARTS_NOEXCEPT
Const version of transpose.
Definition: matpackI.cc:1401
Implementation of Matrix, Vector, and such stuff.
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
Index nelem(const Lines &l)
Number of lines.
constexpr Numeric k
Boltzmann constant convenience name [J/K].
constexpr Numeric e
Elementary charge convenience name [C].
constexpr Numeric alpha
Fine structure constant convenience name [-].
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
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:705
#define M
Definition: rng.cc:165
#define v
#define a
#define c
#define b