ARTS 2.5.10 (git: 2f1c442c)
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
64 return {mdata + mrange.mstart, mrange.mstride};
65}
66
68 return {
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
145 return {mdata + mrange.mstart, mrange.mstride};
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
299void copy(Numeric x, Iterator1D target, const Iterator1D& end) ARTS_NOEXCEPT {
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
322Vector::Vector(Numeric start, Index extent, Numeric stride)
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 auto vec_it_end = v.end();
347 Iterator1D this_it = this->begin();
348 for (auto 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) noexcept {
402 using std::swap;
403 swap(v1.mrange, v2.mrange);
404 swap(v1.mdata, v2.mdata);
405}
406
407Vector::~Vector() noexcept { delete[] mdata; }
408
409// Functions for ConstMatrixView:
410// ------------------------------
411
416 const Range& r, const Range& c) const ARTS_NOEXCEPT {
417 return ConstMatrixView(mdata, mrr, mcr, r, c);
418}
419
426 Index c) const ARTS_NOEXCEPT {
427 // Check that c is valid:
428 ARTS_ASSERT(0 <= c);
429 ARTS_ASSERT(c < mcr.mextent);
430
431 return ConstVectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
432}
433
441 // Check that r is valid:
442 ARTS_ASSERT(0 <= r);
443 ARTS_ASSERT(r < mrr.mextent);
444
445 return ConstVectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
446}
447
451}
452
455 return {
457 mrr.mstride};
458}
459
461
470 Index n = std::min(mrr.mextent, mcr.mextent);
472 Range(0, n, mrr.mstride + mcr.mstride));
473}
474
479 const Range& rr,
480 const Range& cr) ARTS_NOEXCEPT : mrr(rr),
481 mcr(cr),
482 mdata(data) {
483 // Nothing to do here.
484}
485
501 const Range& pr,
502 const Range& pc,
503 const Range& nr,
504 const Range& nc) ARTS_NOEXCEPT : mrr(pr, nr),
505 mcr(pc, nc),
506 mdata(data) {
507 // Nothing to do here.
508}
509
517std::ostream& operator<<(std::ostream& os, const ConstMatrixView& v) {
518 // Row iterators:
519 ConstIterator2D ir = v.begin();
520 const ConstIterator2D end_row = v.end();
521
522 if (ir != end_row) {
523 ConstIterator1D ic = ir->begin();
524 const ConstIterator1D end_col = ir->end();
525
526 if (ic != end_col) {
527 os << setw(3) << *ic;
528 ++ic;
529 }
530 for (; ic != end_col; ++ic) {
531 os << " " << setw(3) << *ic;
532 }
533 ++ir;
534 }
535 for (; ir != end_row; ++ir) {
536 ConstIterator1D ic = ir->begin();
537 const ConstIterator1D end_col = ir->end();
538
539 os << "\n";
540 if (ic != end_col) {
541 os << setw(3) << *ic;
542 ++ic;
543 }
544 for (; ic != end_col; ++ic) {
545 os << " " << setw(3) << *ic;
546 }
547 }
548
549 return os;
550}
551
552// Functions for MatrixView:
553// -------------------------
554
559 const Range& c) ARTS_NOEXCEPT {
560 return MatrixView(mdata, mrr, mcr, r, c);
561}
562
568 // Check that c is valid:
569 ARTS_ASSERT(0 <= c);
570 ARTS_ASSERT(c < mcr.mextent);
571
572 return VectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
573}
574
580 // Check that r is valid:
581 ARTS_ASSERT(0 <= r);
582 ARTS_ASSERT(r < mrr.mextent);
583
584 return VectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
586
588// hiden by the non-const operator of the derived class.*/
589//ConstIterator2D MatrixView::begin() const
590//{
591// return ConstMatrixView::begin();
592//}
593//
595//ConstIterator2D MatrixView::end() const
596//{
597// return ConstMatrixView::end();
598//}
599
603}
607 return Iterator2D(
609 mrr.mstride);
610}
611
617 // Check that sizes are compatible:
620
621 copy(m.begin(), m.end(), begin());
622 return *this;
623}
624
631 // Check that sizes are compatible:
634
635 copy(m.begin(), m.end(), begin());
636 return *this;
637}
638
643 // Check that sizes are compatible:
646
647 copy(m.begin(), m.end(), begin());
648 return *this;
649}
650
656 // Check that sizes are compatible:
657 ARTS_ASSERT(mrr.mextent == v.nelem());
658 ARTS_ASSERT(mcr.mextent == 1);
659 // dummy = ConstMatrixView(v.mdata,v.mrange,Range(v.mrange.mstart,1));;
660 ConstMatrixView dummy(v);
661 copy(dummy.begin(), dummy.end(), begin());
662 return *this;
663}
664
668 copy(x, begin(), end());
669 return *this;
670}
671
674 const Iterator2D er = end();
675 for (Iterator2D r = begin(); r != er; ++r) {
676 const Iterator1D ec = r->end();
677 for (Iterator1D c = r->begin(); c != ec; ++c) *c *= x;
678 }
679 return *this;
680}
681
684 const Iterator2D er = end();
685 for (Iterator2D r = begin(); r != er; ++r) {
686 const Iterator1D ec = r->end();
687 for (Iterator1D c = r->begin(); c != ec; ++c) *c /= x;
688 }
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 ARTS_ASSERT(nrows() == x.nrows());
715 ARTS_ASSERT(ncols() == x.ncols());
716 ConstIterator2D sr = x.begin();
717 Iterator2D r = begin();
718 const Iterator2D er = end();
719 for (; r != er; ++r, ++sr) {
720 ConstIterator1D sc = sr->begin();
721 Iterator1D c = r->begin();
722 const Iterator1D ec = r->end();
723 for (; c != ec; ++c, ++sc) *c *= *sc;
724 }
725 return *this;
726}
727
730 ARTS_ASSERT(nrows() == x.nrows());
731 ARTS_ASSERT(ncols() == x.ncols());
732 ConstIterator2D sr = x.begin();
733 Iterator2D r = begin();
734 const Iterator2D er = end();
735 for (; r != er; ++r, ++sr) {
736 ConstIterator1D sc = sr->begin();
737 Iterator1D c = r->begin();
738 const Iterator1D ec = r->end();
739 for (; c != ec; ++c, ++sc) *c /= *sc;
740 }
741 return *this;
742}
743
746 ARTS_ASSERT(nrows() == x.nrows());
747 ARTS_ASSERT(ncols() == x.ncols());
748 ConstIterator2D sr = x.begin();
749 Iterator2D r = begin();
750 const Iterator2D er = end();
751 for (; r != er; ++r, ++sr) {
752 ConstIterator1D sc = sr->begin();
753 Iterator1D c = r->begin();
754 const Iterator1D ec = r->end();
755 for (; c != ec; ++c, ++sc) *c += *sc;
756 }
757 return *this;
758}
759
762 ARTS_ASSERT(nrows() == x.nrows());
763 ARTS_ASSERT(ncols() == x.ncols());
764 ConstIterator2D sr = x.begin();
765 Iterator2D r = begin();
766 const Iterator2D er = end();
767 for (; r != er; ++r, ++sr) {
768 ConstIterator1D sc = sr->begin();
769 Iterator1D c = r->begin();
770 const Iterator1D ec = r->end();
771 for (; c != ec; ++c, ++sc) *c -= *sc;
772 }
773 return *this;
774}
775
778 ARTS_ASSERT(nrows() == x.nelem());
779 ARTS_ASSERT(ncols() == 1);
780 ConstIterator1D sc = x.begin();
781 Iterator2D r = begin();
782 const Iterator2D er = end();
783 for (; r != er; ++r, ++sc) {
784 Iterator1D c = r->begin();
785 *c *= *sc;
786 }
787 return *this;
788}
789
792 ARTS_ASSERT(nrows() == x.nelem());
793 ARTS_ASSERT(ncols() == 1);
794 ConstIterator1D sc = x.begin();
795 Iterator2D r = begin();
796 const Iterator2D er = end();
797 for (; r != er; ++r, ++sc) {
798 Iterator1D c = r->begin();
799 *c /= *sc;
800 }
801 return *this;
802}
803
806 ARTS_ASSERT(nrows() == x.nelem());
807 ARTS_ASSERT(ncols() == 1);
808 ConstIterator1D sc = x.begin();
809 Iterator2D r = begin();
810 const Iterator2D er = end();
811 for (; r != er; ++r, ++sc) {
812 Iterator1D c = r->begin();
813 *c += *sc;
814 }
815 return *this;
816}
817
820 ARTS_ASSERT(nrows() == x.nelem());
821 ARTS_ASSERT(ncols() == 1);
822 ConstIterator1D sc = x.begin();
823 Iterator2D r = begin();
824 const Iterator2D er = end();
825 for (; r != er; ++r, ++sc) {
826 Iterator1D c = r->begin();
827 *c -= *sc;
828 }
829 return *this;
830}
831
836 const Range& rr,
837 const Range& cr) ARTS_NOEXCEPT
838 : ConstMatrixView(data, rr, cr) {
839 // Nothing to do here.
840}
841
857 const Range& pr,
858 const Range& pc,
859 const Range& nr,
860 const Range& nc) ARTS_NOEXCEPT
861 : ConstMatrixView(data, pr, pc, nr, nc) {
862 // Nothing to do here.
863}
864
875 const ConstIterator2D& end,
876 Iterator2D target) {
877 for (; origin != end; ++origin, ++target) {
878 ConstIterator1D o = origin->begin();
879 const ConstIterator1D e = origin->end();
880 Iterator1D t = target->begin();
881 for (; o != e; ++o, ++t) *t = *o;
882 }
883}
884
886void copy(Numeric x, Iterator2D target, const Iterator2D& end) ARTS_NOEXCEPT {
887 for (; target != end; ++target) {
888 Iterator1D t = target->begin();
889 const Iterator1D e = target->end();
890 for (; t != e; ++t) *t = x;
891 }
892}
893
894// Functions for Matrix:
895// ---------------------
896
900 : MatrixView(new Numeric[r * c], Range(0, r, c), Range(0, c)) {
901 // Nothing to do here.
902}
903
906 : MatrixView(new Numeric[r * c], Range(0, r, c), Range(0, c)) {
907 // Here we can access the raw memory directly, for slightly
908 // increased efficiency:
909 std::fill_n(mdata, r * c, fill);
910}
911
915 : MatrixView(new Numeric[m.nrows() * m.ncols()],
916 Range(0, m.nrows(), m.ncols()),
917 Range(0, m.ncols())) {
918 copy(m.begin(), m.end(), begin());
919}
920
924 : MatrixView(new Numeric[m.nrows() * m.ncols()],
925 Range(0, m.nrows(), m.ncols()),
926 Range(0, m.ncols())) {
927 // There is a catch here: If m is an empty matrix, then it will have
928 // 0 colunns. But this is used to initialize the stride of the row
929 // Range! Thus, this method has to be consistent with the behaviour
930 // of Range::Range. For now, Range::Range allows also stride 0.
931 std::memcpy(mdata, m.mdata, nrows() * ncols() * sizeof(Numeric));
932}
933
935
960 if (this != &m) {
961 resize(m.nrows(), m.ncols());
962 std::memcpy(mdata, m.mdata, nrows() * ncols() * sizeof(Numeric));
963 }
964 return *this;
965}
966
969 if (this != &m) {
970 delete[] mdata;
971 mdata = m.mdata;
972 mrr = m.mrr;
973 mcr = m.mcr;
974 m.mrr = Range(0, 0);
975 m.mcr = Range(0, 0);
976 m.mdata = nullptr;
977 }
978 return *this;
979}
980
984 std::fill_n(mdata, nrows() * ncols(), x);
985 return *this;
986}
987
989
1002 resize(v.nelem(), 1);
1003 ConstMatrixView dummy(v);
1004 copy(dummy.begin(), dummy.end(), begin());
1005 return *this;
1006}
1007
1012 ARTS_ASSERT(0 <= r);
1013 ARTS_ASSERT(0 <= c);
1014
1015 if (mrr.mextent != r || mcr.mextent != c) {
1016 delete[] mdata;
1017 mdata = new Numeric[r * c];
1018
1019 mrr.mstart = 0;
1020 mrr.mextent = r;
1021 mrr.mstride = c;
1022
1023 mcr.mstart = 0;
1024 mcr.mextent = c;
1025 mcr.mstride = 1;
1026 }
1027}
1028
1030void swap(Matrix& m1, Matrix& m2) noexcept {
1031 using std::swap;
1032 swap(m1.mrr, m2.mrr);
1033 swap(m1.mcr, m2.mcr);
1034 swap(m1.mdata, m2.mdata);
1035}
1036
1040 // cout << "Destroying a Matrix:\n"
1041 // << *this << "\n........................................\n";
1042 delete[] mdata;
1043}
1044
1045// Some general Matrix Vector functions:
1046
1050 // Check dimensions:
1051 ARTS_ASSERT(a.nelem() == b.nelem());
1052
1053 const ConstIterator1D ae = a.end();
1054 ConstIterator1D ai = a.begin();
1055 ConstIterator1D bi = b.begin();
1056
1057 Numeric res = 0;
1058 for (; ai != ae; ++ai, ++bi) res += (*ai) * (*bi);
1059
1060 return res;
1061}
1062
1064
1081 ARTS_ASSERT(y.mrange.get_extent() == M.mrr.get_extent());
1082 ARTS_ASSERT(M.mcr.get_extent() == x.mrange.get_extent());
1083 ARTS_ASSERT((M.mcr.get_extent() != 0) && (M.mrr.get_extent() != 0));
1084
1085 if ((M.mcr.get_stride() == 1) || (M.mrr.get_stride() == 1)) {
1086 char trans;
1087 int m, n;
1088 double zero = 0.0;
1089 double one = 1.0;
1090 int LDA, incx, incy;
1091
1092 if (M.mcr.get_stride() != 1) {
1093 trans = 'n';
1094 m = (int)M.mrr.get_extent();
1095 n = (int)M.mcr.get_extent();
1096 LDA = (int)M.mcr.get_stride();
1097 } else {
1098 trans = 't';
1099 m = (int)M.mcr.get_extent();
1100 n = (int)M.mrr.get_extent();
1101 LDA = (int)M.mrr.get_stride();
1102 if (M.mrr.get_stride() == 1) LDA = m;
1103 }
1104
1105 incx = (int)x.mrange.get_stride();
1106 incy = (int)y.mrange.get_stride();
1107
1108 double* mstart = M.mdata + M.mcr.get_start() + M.mrr.get_start();
1109 double* ystart = y.mdata + y.mrange.get_start();
1110 double* xstart = x.mdata + x.mrange.get_start();
1111
1112 dgemv_(&trans,
1113 &m,
1114 &n,
1115 &one,
1116 mstart,
1117 &LDA,
1118 xstart,
1119 &incx,
1120 &zero,
1121 ystart,
1122 &incy);
1123
1124 } else {
1125 mult_general(y, M, x);
1126 }
1127}
1128
1139 const ConstMatrixView& M,
1140 const ConstVectorView& x) ARTS_NOEXCEPT {
1141 // Check dimensions:
1142 ARTS_ASSERT(y.mrange.mextent == M.mrr.mextent);
1143 ARTS_ASSERT(M.mcr.mextent == x.mrange.mextent);
1144 ARTS_ASSERT(M.mcr.mextent != 0 && M.mrr.mextent != 0);
1145
1146 // Let's first find the pointers to the starting positions
1147 Numeric* mdata = M.mdata + M.mcr.mstart + M.mrr.mstart;
1148 Numeric* xdata = x.mdata + x.mrange.mstart;
1149 Numeric* yelem = y.mdata + y.mrange.mstart;
1150
1151 Index i = M.mrr.mextent;
1152 while (i--) {
1153 Numeric* melem = mdata;
1154 Numeric* xelem = xdata; // Reset xelem to first element of source vector
1155
1156 // Multiply first element of matrix row with first element of source
1157 // vector. This is done outside the loop to avoid initialization of the
1158 // target vector's element with zero (saves one assignment)
1159 *yelem = *melem * *xelem;
1160
1161 Index j = M.mcr.mextent; // --j (instead of j-- like in the outer loop)
1162 while (--j) // is important here because we only want
1163 { // mextent-1 cycles
1164 melem += M.mcr.mstride;
1165 xelem += x.mrange.mstride;
1166 *yelem += *melem * *xelem;
1167 }
1168
1169 mdata += M.mrr.mstride; // Jump to next matrix row
1170 yelem += y.mrange.mstride; // Jump to next element in target vector
1171 }
1172}
1173
1175
1199void mult(MatrixView A, const ConstMatrixView& B, const ConstMatrixView& C) {
1200 // Check dimensions:
1201 ARTS_ASSERT(A.nrows() == B.nrows());
1202 ARTS_ASSERT(A.ncols() == C.ncols());
1203 ARTS_ASSERT(B.ncols() == C.nrows());
1204
1205 // Catch trivial case if one of the matrices is empty.
1206 if ((B.nrows() == 0) || (B.ncols() == 0) || (C.ncols() == 0)) return;
1207
1208 // Matrices B and C must be continuous in at least on dimension, C
1209 // must be continuous along the second dimension.
1210 if (((B.mrr.get_stride() == 1) || (B.mcr.get_stride() == 1)) &&
1211 ((C.mrr.get_stride() == 1) || (C.mcr.get_stride() == 1)) &&
1212 (A.mcr.get_stride() == 1)) {
1213 // BLAS uses column-major order while arts uses row-major order.
1214 // Hence instead of C = A * B we compute C^T = A^T * B^T!
1215
1216 int k, m, n;
1217
1218 k = (int)B.ncols();
1219 m = (int)C.ncols();
1220 n = (int)B.nrows();
1221
1222 // Note also the clash in nomenclature: BLAS uses C = A * B while
1223 // arts uses A = B * C. Taking into accout this and the difference in
1224 // memory layouts, we need to map the MatrixViews A, B and C to the BLAS
1225 // arguments as follows:
1226 // A (arts) -> C (BLAS)
1227 // B (arts) -> B (BLAS)
1228 // C (arts) -> A (BLAS)
1229
1230 // Char indicating whether A (BLAS) or B (BLAS) should be transposed.
1231 char transa, transb;
1232 // Sizes of the matrices along the direction in which they are
1233 // traversed.
1234 int lda, ldb, ldc;
1235
1236 // Check if C (arts) is transposed.
1237 if (C.mrr.get_stride() == 1) {
1238 transa = 'T';
1239 lda = (int)C.mcr.get_stride();
1240 } else {
1241 transa = 'N';
1242 lda = (int)C.mrr.get_stride();
1243 }
1244
1245 // Check if B (arts) is transposed.
1246 if (B.mrr.get_stride() == 1) {
1247 transb = 'T';
1248 ldb = (int)B.mcr.get_stride();
1249 } else {
1250 transb = 'N';
1251 ldb = (int)B.mrr.get_stride();
1252 }
1253
1254 // In case B (arts) has only one column, column and row stride are 1.
1255 // We therefore need to set ldb to k, since dgemm_ requires lda to be at
1256 // least k / m if A is non-transposed / transposed.
1257 if ((B.mcr.get_stride() == 1) && (B.mrr.get_stride() == 1)) {
1258 transb = 'N';
1259 ldb = k;
1260 }
1261
1262 // The same holds for C (arts).
1263 if ((C.mcr.get_stride() == 1) && (C.mrr.get_stride() == 1)) {
1264 transa = 'N';
1265 lda = m;
1266 }
1267
1268 ldc = (int)A.mrr.get_stride();
1269 // The same holds for A (arts).
1270 if ((A.mcr.get_stride() == 1) && (A.mrr.get_stride() == 1)) {
1271 ldc = m;
1272 }
1273 double alpha = 1.0, beta = 0.0;
1274
1275 dgemm_(&transa,
1276 &transb,
1277 &m,
1278 &n,
1279 &k,
1280 &alpha,
1281 C.mdata + C.mrr.get_start() + C.mcr.get_start(),
1282 &lda,
1283 B.mdata + B.mrr.get_start() + B.mcr.get_start(),
1284 &ldb,
1285 &beta,
1286 A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1287 &ldc);
1288
1289 } else {
1290 mult_general(A, B, C);
1291 }
1292}
1293
1295
1304 const ConstMatrixView& B,
1305 const ConstMatrixView& C) ARTS_NOEXCEPT {
1306 // Check dimensions:
1307 ARTS_ASSERT(A.nrows() == B.nrows());
1308 ARTS_ASSERT(A.ncols() == C.ncols());
1309 ARTS_ASSERT(B.ncols() == C.nrows());
1310
1311 // Let's get the transpose of C, so that we can use 2D iterators to
1312 // access the columns (= rows of the transpose).
1313 ConstMatrixView CT = transpose(C);
1314
1315 const Iterator2D ae = A.end();
1316 Iterator2D ai = A.begin();
1317 ConstIterator2D bi = B.begin();
1318
1319 // This walks through the rows of A and B:
1320 for (; ai != ae; ++ai, ++bi) {
1321 const Iterator1D ace = ai->end();
1322 Iterator1D aci = ai->begin();
1323 ConstIterator2D cti = CT.begin();
1324
1325 // This walks through the columns of A with a 1D iterator, and
1326 // at the same time through the rows of CT, which are the columns of
1327 // C, with a 2D iterator:
1328 for (; aci != ace; ++aci, ++cti) {
1329 // The operator * is used to compute the scalar product
1330 // between rows of B and rows of C.transpose().
1331 *aci = (*bi) * (*cti);
1332 }
1333 }
1334}
1335
1337
1351 const ConstVectorView& a,
1353 ARTS_ASSERT(a.nelem() == 3);
1354 ARTS_ASSERT(b.nelem() == 3);
1355 ARTS_ASSERT(c.nelem() == 3);
1356
1357 c[0] = a[1] * b[2] - a[2] * b[1];
1358 c[1] = a[2] * b[0] - a[0] * b[2];
1359 c[2] = a[0] * b[1] - a[1] * b[0];
1360}
1361
1372 ARTS_ASSERT(a.nelem() == b.nelem());
1373 Numeric arg = (a * b) / sqrt(a * a) / sqrt(b * b);
1374
1375 // Due to numerical inaccuracy, arg might be slightly larger than 1.
1376 // We catch those cases to avoid spurious returns of NaNs
1377 return fabs(arg) > 1. ? 0. : acos(arg) * RAD2DEG;
1378}
1379
1394 ARTS_ASSERT(a.nelem() == b.nelem());
1395 ARTS_ASSERT(a.nelem() == c.nelem());
1396
1397 const Numeric C = (a * b) / (a * a);
1398 c = a;
1399 c *= C;
1400};
1401
1404 return ConstMatrixView(m.mdata, m.mcr, m.mrr);
1405}
1406
1410 return MatrixView(m.mdata, m.mcr, m.mrr);
1411}
1412
1433void transform(VectorView y, double (&my_func)(double), ConstVectorView x) {
1434 // Check dimensions:
1435 ARTS_ASSERT(y.nelem() == x.nelem());
1436
1437 const ConstIterator1D xe = x.end();
1438 ConstIterator1D xi = x.begin();
1439 Iterator1D yi = y.begin();
1440 for (; xi != xe; ++xi, ++yi) *yi = my_func(*xi);
1441}
1442
1461void transform(MatrixView y, double (&my_func)(double), ConstMatrixView x) {
1462 // Check dimensions:
1463 ARTS_ASSERT(y.nrows() == x.nrows());
1464 ARTS_ASSERT(y.ncols() == x.ncols());
1465
1466 const ConstIterator2D rxe = x.end();
1467 ConstIterator2D rx = x.begin();
1468 Iterator2D ry = y.begin();
1469 for (; rx != rxe; ++rx, ++ry) {
1470 const ConstIterator1D cxe = rx->end();
1471 ConstIterator1D cx = rx->begin();
1472 Iterator1D cy = ry->begin();
1473 for (; cx != cxe; ++cx, ++cy) *cy = my_func(*cx);
1474 }
1475}
1476
1479 // Initial value for max:
1480 Numeric max = x[0];
1481
1482 const ConstIterator1D xe = x.end();
1483 ConstIterator1D xi = x.begin();
1484
1485 for (; xi != xe; ++xi) {
1486 if (*xi > max) max = *xi;
1487 }
1488
1489 return max;
1490}
1491
1494 // Initial value for max:
1495 Numeric max = x(0, 0);
1496
1497 const ConstIterator2D rxe = x.end();
1498 ConstIterator2D rx = x.begin();
1499
1500 for (; rx != rxe; ++rx) {
1501 const ConstIterator1D cxe = rx->end();
1502 ConstIterator1D cx = rx->begin();
1503
1504 for (; cx != cxe; ++cx)
1505 if (*cx > max) max = *cx;
1506 }
1507
1508 return max;
1509}
1510
1513 // Initial value for min:
1514 Numeric min = x[0];
1515
1516 const ConstIterator1D xe = x.end();
1517 ConstIterator1D xi = x.begin();
1518
1519 for (; xi != xe; ++xi) {
1520 if (*xi < min) min = *xi;
1521 }
1522
1523 return min;
1524}
1525
1528 // Initial value for min:
1529 Numeric min = x(0, 0);
1530
1531 const ConstIterator2D rxe = x.end();
1532 ConstIterator2D rx = x.begin();
1533
1534 for (; rx != rxe; ++rx) {
1535 const ConstIterator1D cxe = rx->end();
1536 ConstIterator1D cx = rx->begin();
1537
1538 for (; cx != cxe; ++cx)
1539 if (*cx < min) min = *cx;
1540 }
1541
1542 return min;
1543}
1544
1547 // Initial value for mean:
1548 Numeric mean = 0;
1549
1550 const ConstIterator1D xe = x.end();
1551 ConstIterator1D xi = x.begin();
1552
1553 for (; xi != xe; ++xi) mean += *xi;
1554
1555 mean /= (Numeric)x.nelem();
1556
1557 return mean;
1558}
1559
1562 // Initial value for mean:
1563 Numeric nanmean = 0;
1564
1565 // Counter for the number of normal values
1566 Index numnormal = 0;
1567
1568 const ConstIterator1D xe = x.end();
1569 ConstIterator1D xi = x.begin();
1570
1571 for (; xi != xe; ++xi)
1572 if (std::isnormal(*xi) or (*xi == 0)) {
1573 nanmean += *xi;
1574 numnormal++;
1575 }
1576
1577 if (numnormal)
1578 nanmean /= Numeric(numnormal);
1579 else
1580 nanmean = std::numeric_limits<Numeric>::quiet_NaN();
1581
1582 return nanmean;
1583}
1584
1587 // Initial value for mean:
1588 Numeric mean = 0;
1589
1590 const ConstIterator2D rxe = x.end();
1591 ConstIterator2D rx = x.begin();
1592
1593 for (; rx != rxe; ++rx) {
1594 const ConstIterator1D cxe = rx->end();
1595 ConstIterator1D cx = rx->begin();
1596
1597 for (; cx != cxe; ++cx) mean += *cx;
1598 }
1599
1600 mean /= (Numeric)(x.nrows() * x.ncols());
1601
1602 return mean;
1603}
1604
1615 // cout << "Assigning VectorView from Array<Numeric>.\n";
1616
1617 // Check that sizes are compatible:
1618 ARTS_ASSERT(mrange.mextent == v.nelem());
1619
1620 // Iterators for Array:
1621 auto i = v.begin();
1622 const auto e = v.end();
1623 // Iterator for Vector:
1624 Iterator1D target = begin();
1625
1626 for (; i != e; ++i, ++target) *target = *i;
1627
1628 return *this;
1629}
1630
1632// Helper function for debugging
1633#ifndef NDEBUG
1634
1649 return mv(r, c);
1650}
1651
1652#endif
base max(const Array< base > &x)
Max function.
Definition: array.h:145
base min(const Array< base > &x)
Min function.
Definition: array.h:161
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:48
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
The constant iterator class for sub vectors.
Definition: matpackI.h:458
const Numeric * mx
Current position.
Definition: matpackI.h:505
Index mstride
Stride.
Definition: matpackI.h:507
The const row iterator class for sub matrices.
Definition: matpackI.h:859
A constant view of a Matrix.
Definition: matpackI.h:1065
Range mrr
The row range of mdata that is actually used.
Definition: matpackI.h:1172
friend class ConstVectorView
Definition: matpackI.h:1127
ConstIterator2D begin() const ARTS_NOEXCEPT
Return const iterator to first row.
Definition: matpackI.cc:449
Index nrows() const noexcept
Definition: matpackI.h:1079
ConstIterator2D end() const ARTS_NOEXCEPT
Return const iterator behind last row.
Definition: matpackI.cc:454
Index ncols() const noexcept
Definition: matpackI.h:1080
Numeric operator()(Index r, Index c) const ARTS_NOEXCEPT
Plain const index operator.
Definition: matpackI.h:1093
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1176
Range mcr
The column range of mdata that is actually used.
Definition: matpackI.h:1174
ConstVectorView diagonal() const ARTS_NOEXCEPT
Matrix diagonal as vector.
Definition: matpackI.cc:469
ConstMatrixView()=default
A constant view of a Vector.
Definition: matpackI.h:521
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:660
Numeric sum() const ARTS_NOEXCEPT
The sum of all elements of a Vector.
Definition: matpackI.cc:48
Numeric operator[](Index n) const ARTS_NOEXCEPT
Plain const index operator.
Definition: matpackI.h:560
Range mrange
The range of mdata that is actually used.
Definition: matpackI.h:658
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
The iterator class for sub vectors.
Definition: matpackI.h:399
Index mstride
Stride.
Definition: matpackI.h:453
Numeric * mx
Current position.
Definition: matpackI.h:451
The row iterator class for sub matrices.
Definition: matpackI.h:815
The Joker class.
Definition: matpackI.h:126
The MatrixView class.
Definition: matpackI.h:1188
MatrixView & operator=(const ConstMatrixView &m)
Assignment operator.
Definition: matpackI.cc:616
friend class VectorView
Definition: matpackI.h:1254
MatrixView & operator/=(Numeric x) ARTS_NOEXCEPT
Division by scalar.
Definition: matpackI.cc:683
Numeric & operator()(Index r, Index c) ARTS_NOEXCEPT
Plain index operator.
Definition: matpackI.h:1205
Iterator2D end() ARTS_NOEXCEPT
Return iterator behind last row.
Definition: matpackI.cc:606
MatrixView & operator+=(Numeric x) ARTS_NOEXCEPT
Addition of scalar.
Definition: matpackI.cc:693
Iterator2D begin() ARTS_NOEXCEPT
‍** Return const iterator to first row. Has to be redefined here, since it is
Definition: matpackI.cc:601
MatrixView & operator*=(Numeric x) ARTS_NOEXCEPT
Multiplication by scalar.
Definition: matpackI.cc:673
MatrixView & operator-=(Numeric x) ARTS_NOEXCEPT
Subtraction of scalar.
Definition: matpackI.cc:703
MatrixView()=default
The Matrix class.
Definition: matpackI.h:1285
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1011
Matrix()=default
~Matrix() noexcept override
Destructor for Matrix.
Definition: matpackI.cc:1039
Matrix & operator=(const Matrix &m)
Assignment operator from another matrix.
Definition: matpackI.cc:959
The range class.
Definition: matpackI.h:160
Index mstart
The start index.
Definition: matpackI.h:377
constexpr Index get_extent() const noexcept
Returns the extent of the range.
Definition: matpackI.h:343
Index mstride
The stride.
Definition: matpackI.h:382
constexpr Index get_start() const noexcept
Returns the start index of the range.
Definition: matpackI.h:341
Index mextent
The number of elements.
Definition: matpackI.h:380
constexpr Index get_stride() const noexcept
Returns the stride of the range.
Definition: matpackI.h:345
The VectorView class.
Definition: matpackI.h:674
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:696
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:910
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
~Vector() noexcept override
Destructor for Vector.
Definition: matpackI.cc:407
Vector()=default
Vector & operator=(const Vector &v)
Assignment from another Vector.
Definition: matpackI.cc:360
#define ARTS_NOEXCEPT
Definition: debug.h:99
#define ARTS_ASSERT(condition,...)
Definition: debug.h:102
constexpr Numeric RAD2DEG
Definition: doit.cc:60
The declarations of all the exception classes.
Numeric mean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version.
Definition: matpackI.cc:1546
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix-Vector Multiplication.
Definition: matpackI.cc:1080
void swap(Vector &v1, Vector &v2) noexcept
Definition: matpackI.cc:401
Numeric vector_angle(ConstVectorView a, ConstVectorView b)
Definition: matpackI.cc:1371
Numeric nanmean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version ignoring nans and infs.
Definition: matpackI.cc:1561
void proj(Vector &c, ConstVectorView a, ConstVectorView b) ARTS_NOEXCEPT
Definition: matpackI.cc:1393
void mult_general(VectorView y, const ConstMatrixView &M, const ConstVectorView &x) ARTS_NOEXCEPT
Matrix Vector multiplication.
Definition: matpackI.cc:1138
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:1433
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
cross3
Definition: matpackI.cc:1350
void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Definition: matpackI.cc:288
Numeric debug_matrixview_get_elem(MatrixView &mv, Index r, Index c)
Helper function to access matrix elements.
Definition: matpackI.cc:1648
Numeric operator*(const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
Scalar product.
Definition: matpackI.cc:1048
ConstMatrixView transpose(ConstMatrixView m) ARTS_NOEXCEPT
Const version of transpose.
Definition: matpackI.cc:1403
Implementation of Matrix, Vector, and such stuff.
void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Definition: matpackI.cc:288
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
constexpr Numeric alpha
Fine structure constant convenience name [-].
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