ARTS 2.5.10 (git: 2f1c442c)
matpack_complex.cc
Go to the documentation of this file.
1/* Copyright (C) 2002-2012 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
26#include "blas.h"
27#include "exceptions.h"
28#include "matpack_complex.h"
29#include "matpack_eigen.h"
30#include <algorithm>
31#include <cmath>
32#include <cstring>
33
34// Functions for ConstComplexVectorView:
35// ------------------------------
36
38// bool ConstComplexVectorView::empty() const
39// {
40// return (nelem() == 0);
41// }
42
52// Index ConstComplexVectorView::nelem() const
53// {
54// return mrange.mextent;
55// }
56
59 Complex s = 0;
61 const ConstComplexIterator1D e = end();
62
63 for (; i != e; ++i) s += *i;
64
65 return s;
66}
67
72 const Range& r) const {
74}
75
79}
80
86}
87
89ConstComplexVectorView::operator ConstComplexMatrixView() const {
90 return ConstComplexMatrixView(mdata, mrange, Range(0, 1));
91}
92
101 : mrange(0, 1), mdata(&const_cast<Complex&>(a)) {
102 // Nothing to do here.
103}
104
108 const Range& range)
109 : mrange(range), mdata(data) {
110 // Nothing to do here.
111}
112
124 const Range& p,
125 const Range& n)
126 : mrange(p, n), mdata(data) {
127 // Nothing to do here.
128}
129
133std::ostream& operator<<(std::ostream& os, const ConstComplexVectorView& v) {
134 ConstComplexIterator1D i = v.begin();
135 const ConstComplexIterator1D end = v.end();
136
137 if (i != end) {
138 os << std::setw(3) << *i;
139 ++i;
140 }
141 for (; i != end; ++i) {
142 os << " " << std::setw(3) << *i;
143 }
144
145 return os;
146}
147
148// Functions for ComplexVectorView:
149// ------------------------
150
154 ARTS_ASSERT (false,
155 "Creating a ComplexVectorView from a const ComplexVector is not allowed.\n"
156 "This is not really a runtime error, but I don't want to start\n"
157 "producing direct output from inside matpack. And just exiting is\n"
158 "not so nice.\n"
159 "If you see this error, there is a bug in the code, not in the\n"
160 "ARTS input.")
161}
162
165 mdata = v.mdata;
166 mrange = v.mrange;
167}
168
173 return ComplexVectorView(mdata, mrange, r);
174}
175
179}
180
183 return ComplexIterator1D(
186}
187
193 const ConstComplexVectorView& v) {
194 // cout << "Assigning VectorView from ConstVectorView.\n";
195
196 // Check that sizes are compatible:
197 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
198
199 copy(v.begin(), v.end(), begin());
200
201 return *this;
202}
203
210 // cout << "Assigning ComplexVectorView from ComplexVectorView.\n";
211
212 // Check that sizes are compatible:
213 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
214
215 copy(v.begin(), v.end(), begin());
216
217 return *this;
218}
219
223 // cout << "Assigning ComplexVectorView from Vector.\n";
224
225 // Check that sizes are compatible:
226 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
227
228 copy(v.begin(), v.end(), begin());
229
230 return *this;
231}
232
236 copy(x, begin(), end());
237 return *this;
238}
239
243 copy(x, begin(), end());
244 return *this;
245}
246
249 const ComplexIterator1D e = end();
250 for (ComplexIterator1D i = begin(); i != e; ++i) *i *= x;
251 return *this;
252}
253
256 const ComplexIterator1D e = end();
257 for (ComplexIterator1D i = begin(); i != e; ++i) *i *= x;
258 return *this;
259}
260
263 const ComplexIterator1D e = end();
264 for (ComplexIterator1D i = begin(); i != e; ++i) *i /= x;
265 return *this;
266}
267
270 const ComplexIterator1D e = end();
271 for (ComplexIterator1D i = begin(); i != e; ++i) *i /= x;
272 return *this;
273}
274
277 const ComplexIterator1D e = end();
278 for (ComplexIterator1D i = begin(); i != e; ++i) *i += x;
279 return *this;
280}
281
284 const ComplexIterator1D e = end();
285 for (ComplexIterator1D i = begin(); i != e; ++i) *i += x;
286 return *this;
287}
288
291 const ComplexIterator1D e = end();
292 for (ComplexIterator1D i = begin(); i != e; ++i) *i -= x;
293 return *this;
294}
295
298 const ComplexIterator1D e = end();
299 for (ComplexIterator1D i = begin(); i != e; ++i) *i -= x;
300 return *this;
301}
302
305 const ConstComplexVectorView& x) {
306 ARTS_ASSERT(nelem() == x.nelem());
307
309
311 const ComplexIterator1D e = end();
312
313 for (; i != e; ++i, ++s) *i *= *s;
314 return *this;
315}
316
319 ARTS_ASSERT(nelem() == x.nelem());
320
321 ConstIterator1D s = x.begin();
322
324 const ComplexIterator1D e = end();
325
326 for (; i != e; ++i, ++s) *i *= *s;
327 return *this;
328}
329
332 const ConstComplexVectorView& x) {
334
336
339
340 for (; i != e; ++i, ++s) *i /= *s;
341 return *this;
342}
343
346 ARTS_ASSERT(nelem() == x.nelem());
347
348 ConstIterator1D s = x.begin();
349
351 const ComplexIterator1D e = end();
352
353 for (; i != e; ++i, ++s) *i /= *s;
354 return *this;
355}
356
359 const ConstComplexVectorView& x) {
360 ARTS_ASSERT(nelem() == x.nelem());
361
363
365 const ComplexIterator1D e = end();
366
367 for (; i != e; ++i, ++s) *i += *s;
368 return *this;
369}
370
373 ARTS_ASSERT(nelem() == x.nelem());
374
375 ConstIterator1D s = x.begin();
376
378 const ComplexIterator1D e = end();
379
380 for (; i != e; ++i, ++s) *i += *s;
381 return *this;
382}
383
386 const ConstComplexVectorView& x) {
387 ARTS_ASSERT(nelem() == x.nelem());
388
390
392 const ComplexIterator1D e = end();
393
394 for (; i != e; ++i, ++s) *i -= *s;
395 return *this;
396}
397
400 ARTS_ASSERT(nelem() == x.nelem());
401
402 ConstIterator1D s = x.begin();
403
405 const ComplexIterator1D e = end();
406
407 for (; i != e; ++i, ++s) *i -= *s;
408 return *this;
409}
410
412ComplexVectorView::operator ComplexMatrixView() {
413 // The old version (before 2013-01-18) of this was:
414 // return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
415 // Bus this was a bug! The problem is that the matrix index operator adds
416 // the mstart from both row and columm range object to mdata
417
418 return ComplexMatrixView(mdata, mrange, Range(0, 1));
419}
420
425 // Nothing to do here.
426}
427
431 : ConstComplexVectorView(data, range) {
432 // Nothing to do here.
433}
434
446 const Range& p,
447 const Range& n)
448 : ConstComplexVectorView(data, p, n) {
449 // Nothing to do here.
450}
451
457 const ConstComplexIterator1D& end,
458 ComplexIterator1D target) {
459 if (origin.mstride == 1 && target.mstride == 1)
460 memcpy((void*)target.mx,
461 (void*)origin.mx,
462 sizeof(Complex) * (end.mx - origin.mx));
463 else
464 for (; origin != end; ++origin, ++target) *target = *origin;
465}
466
468void copy(Complex x, ComplexIterator1D target, const ComplexIterator1D& end) {
469 for (; target != end; ++target) *target = x;
470}
471
472// Functions for Vector:
473// ---------------------
474
477 : ComplexVectorView(new Complex[n], Range(0, n)) {
478 // Nothing to do here.
479}
480
483 : ComplexVectorView(new Complex[n], Range(0, n)) {
484 // Here we can access the raw memory directly, for slightly
485 // increased efficiency:
486 const Complex* stop = mdata + n;
487 for (Complex* x = mdata; x < stop; ++x) *x = fill;
488}
489
492 : ComplexVectorView(new Complex[n], Range(0, n)) {
493 // Here we can access the raw memory directly, for slightly
494 // increased efficiency:
495 const Complex* stop = mdata + n;
496 for (Complex* x = mdata; x < stop; ++x) *x = fill;
497}
498
508 : ComplexVectorView(new Complex[extent], Range(0, extent)) {
509 // Fill with values:
510 Complex x = start;
512 const ComplexIterator1D e = end();
513 for (; i != e; ++i) {
514 *i = x;
515 x += stride;
516 }
517}
518
528 : ComplexVectorView(new Complex[extent], Range(0, extent)) {
529 // Fill with values:
530 Complex x = start;
532 const ComplexIterator1D e = end();
533 for (; i != e; ++i) {
534 *i = x;
535 x += stride;
536 }
537}
538
548 : ComplexVectorView(new Complex[extent], Range(0, extent)) {
549 // Fill with values:
550 Complex x = start;
552 const ComplexIterator1D e = end();
553 for (; i != e; ++i) {
554 *i = x;
555 x += stride;
556 }
557}
558
568 : ComplexVectorView(new Complex[extent], Range(0, extent)) {
569 // Fill with values:
570 Complex x = start;
572 const ComplexIterator1D e = end();
573 for (; i != e; ++i) {
574 *i = x;
575 x += stride;
576 }
577}
578
585 : ComplexVectorView(new Complex[v.nelem()], Range(0, v.nelem())) {
586 copy(v.begin(), v.end(), begin());
587}
588
592 : ComplexVectorView(new Complex[v.nelem()], Range(0, v.nelem())) {
593 copy(v.begin(), v.end(), begin());
594}
595
597 for (Index i=0; i<nelem(); i++) operator[](i) = Complex(v[i], 0);
598}
599
601ComplexVector::ComplexVector(const std::vector<Complex>& v)
602 : ComplexVectorView(new Complex[v.size()], Range(0, v.size())) {
603 auto vec_it_end = v.end();
604 ComplexIterator1D this_it = this->begin();
605 for (auto vec_it = v.begin();
606 vec_it != vec_it_end;
607 ++vec_it, ++this_it)
608 *this_it = *vec_it;
609}
610
612ComplexVector::ComplexVector(const std::vector<Numeric>& v)
613 : ComplexVectorView(new Complex[v.size()], Range(0, v.size())) {
614 auto vec_it_end = v.end();
615 ComplexIterator1D this_it = this->begin();
616 for (auto vec_it = v.begin();
617 vec_it != vec_it_end;
618 ++vec_it, ++this_it)
619 *this_it = *vec_it;
620}
621
623
648 if (this not_eq &v) {
649 resize(v.nelem());
651 }
652
653 return *this;
654}
655
657 return ConstVectorView(
658 reinterpret_cast<Numeric*>(mdata),
660}
661
663 return ConstVectorView(
664 reinterpret_cast<Numeric*>(mdata),
666}
667
669 return VectorView(
670 reinterpret_cast<Numeric*>(mdata),
672}
673
675 return VectorView(
676 reinterpret_cast<Numeric*>(mdata),
678}
679
681 return MatrixView(reinterpret_cast<Numeric*>(mdata),
682 Range(2 * mrr.mstart, mrr.mextent, mrr.mstride * 2),
683 Range(2 * mcr.mstart, mcr.mextent, mcr.mstride * 2));
684}
685
687 return MatrixView(reinterpret_cast<Numeric*>(mdata),
688 Range(2 * mrr.mstart, mrr.mextent, mrr.mstride * 2),
689 Range(2 * mcr.mstart + 1, mcr.mextent, mcr.mstride * 2));
690}
691
693 return ConstMatrixView(reinterpret_cast<Numeric*>(mdata),
694 Range(2 * mrr.mstart, mrr.mextent, mrr.mstride * 2),
695 Range(2 * mcr.mstart, mcr.mextent, mcr.mstride * 2));
696}
697
699 return ConstMatrixView(
700 reinterpret_cast<Numeric*>(mdata),
701 Range(2 * mrr.mstart, mrr.mextent, mrr.mstride * 2),
702 Range(2 * mcr.mstart + 1, mcr.mextent, mcr.mstride * 2));
703}
704
706
723 resize(x.nelem());
725 return *this;
726}
727
732 return *this;
733}
734
735// /** Assignment operator from VectorView. Assignment operators are not
736// inherited. */
737// inline Vector& Vector::operator=(const VectorView v)
738// {
739// cout << "Assigning Vector from Vector View.\n";
740// // Check that sizes are compatible:
741// ARTS_ASSERT(mrange.mextent==v.mrange.mextent);
742// VectorView::operator=(v);
743// return *this;
744// }
745
750 ARTS_ASSERT(0 <= n);
751 if (mrange.mextent != n) {
752 delete[] mdata;
753 mdata = new Complex[n];
755 mrange.mextent = n;
756 mrange.mstride = 1;
758}
759
761void swap(ComplexVector& v1, ComplexVector& v2) noexcept {
762 using std::swap;
763 swap(v1.mrange, v2.mrange);
764 swap(v1.mdata, v2.mdata);
765}
766
769ComplexVector::~ComplexVector() noexcept { delete[] mdata; }
770
771// Functions for ConstMatrixView:
772// ------------------------------
773
775// bool ConstComplexMatrixView::empty() const
776// {
777// return (nrows() == 0 || ncols() == 0);
778// }
779
781// Index ConstComplexMatrixView::nrows() const
782// {
783// return mrr.mextent;
784// }
785
787// Index ConstComplexMatrixView::ncols() const
788// {
789// return mcr.mextent;
790// }
791
796 const Range& r, const Range& c) const {
797 return ConstComplexMatrixView(mdata, mrr, mcr, r, c);
798}
799
806 Index c) const {
807 // Check that c is valid:
808 ARTS_ASSERT(0 <= c);
810
812}
813
820 Index r, const Range& c) const {
821 // Check that r is valid:
822 ARTS_ASSERT(0 <= r);
824
826}
827
831 mrr.mstride);
832}
833
838 mcr),
839 mrr.mstride);
840}
841
843
852 Index n = std::min(mrr.mextent, mcr.mextent);
854 Range(0, n, mrr.mstride + mcr.mstride));
855}
856
861 const Range& rr,
862 const Range& cr)
863 : mrr(rr), mcr(cr), mdata(data) {
864 // Nothing to do here.
865}
866
882 const Range& pr,
883 const Range& pc,
884 const Range& nr,
885 const Range& nc)
886 : mrr(pr, nr), mcr(pc, nc), mdata(data) {
887 // Nothing to do here.
888}
889
897std::ostream& operator<<(std::ostream& os, const ConstComplexMatrixView& v) {
898 // Row iterators:
899 ConstComplexIterator2D ir = v.begin();
900 const ConstComplexIterator2D end_row = v.end();
901
902 if (ir != end_row) {
903 ConstComplexIterator1D ic = ir->begin();
904 const ConstComplexIterator1D end_col = ir->end();
905
906 if (ic != end_col) {
907 os << std::setw(3) << *ic;
908 ++ic;
909 }
910 for (; ic != end_col; ++ic) {
911 os << " " << std::setw(3) << *ic;
912 }
913 ++ir;
914 }
915 for (; ir != end_row; ++ir) {
916 ConstComplexIterator1D ic = ir->begin();
917 const ConstComplexIterator1D end_col = ir->end();
918
919 os << "\n";
920 if (ic != end_col) {
921 os << std::setw(3) << *ic;
922 ++ic;
923 }
924 for (; ic != end_col; ++ic) {
925 os << " " << std::setw(3) << *ic;
926 }
927 }
928
929 return os;
930}
931
932// Functions for ComplexMatrixView:
933// -------------------------
934
939 const Range& c) {
940 return ComplexMatrixView(mdata, mrr, mcr, r, c);
941}
942
948 // Check that c is valid:
949 ARTS_ASSERT(0 <= c);
951
952 return ComplexVectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
953}
954
960 // Check that r is valid:
961 ARTS_ASSERT(0 <= r);
963
964 return ComplexVectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
965}
966
970 mrr.mstride);
971}
972
975 return ComplexIterator2D(
977 mrr.mstride);
978}
979
981
990 Index n = std::min(mrr.mextent, mcr.mextent);
992 Range(0, n, mrr.mstride + mcr.mstride));
993}
994
1000 const ConstComplexMatrixView& m) {
1001 // Check that sizes are compatible:
1004
1005 copy(m.begin(), m.end(), begin());
1006 return *this;
1007}
1008
1015 // Check that sizes are compatible:
1018
1019 copy(m.begin(), m.end(), begin());
1020 return *this;
1021}
1022
1027 // Check that sizes are compatible:
1030
1031 copy(m.begin(), m.end(), begin());
1032 return *this;
1033}
1034
1040 const ConstComplexVectorView& v) {
1041 // Check that sizes are compatible:
1042 ARTS_ASSERT(mrr.mextent == v.nelem());
1043 ARTS_ASSERT(mcr.mextent == 1);
1044 // dummy = ConstComplexMatrixView(v.mdata,v.mrange,Range(v.mrange.mstart,1));;
1046 copy(dummy.begin(), dummy.end(), begin());
1047 return *this;
1048}
1049
1053 copy(x, begin(), end());
1054 return *this;
1055}
1056
1059 const ComplexIterator2D er = end();
1060 for (ComplexIterator2D r = begin(); r != er; ++r) {
1061 const ComplexIterator1D ec = r->end();
1062 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c *= x;
1063 }
1064 return *this;
1065}
1066
1069 const ComplexIterator2D er = end();
1070 for (ComplexIterator2D r = begin(); r != er; ++r) {
1071 const ComplexIterator1D ec = r->end();
1072 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c *= x;
1073 }
1074 return *this;
1075}
1076
1079 const ComplexIterator2D er = end();
1080 for (ComplexIterator2D r = begin(); r != er; ++r) {
1081 const ComplexIterator1D ec = r->end();
1082 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c /= x;
1083 }
1084 return *this;
1085}
1086
1089 const ComplexIterator2D er = end();
1090 for (ComplexIterator2D r = begin(); r != er; ++r) {
1091 const ComplexIterator1D ec = r->end();
1092 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c /= x;
1093 }
1094 return *this;
1095}
1096
1099 const ComplexIterator2D er = end();
1100 for (ComplexIterator2D r = begin(); r != er; ++r) {
1101 const ComplexIterator1D ec = r->end();
1102 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c += x;
1103 }
1104 return *this;
1105}
1106
1109 const ComplexIterator2D er = end();
1110 for (ComplexIterator2D r = begin(); r != er; ++r) {
1111 const ComplexIterator1D ec = r->end();
1112 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c += x;
1113 }
1114 return *this;
1115}
1116
1119 const ComplexIterator2D er = end();
1120 for (ComplexIterator2D r = begin(); r != er; ++r) {
1121 const ComplexIterator1D ec = r->end();
1122 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c -= x;
1123 }
1124 return *this;
1125}
1126
1129 const ComplexIterator2D er = end();
1130 for (ComplexIterator2D r = begin(); r != er; ++r) {
1131 const ComplexIterator1D ec = r->end();
1132 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c -= x;
1133 }
1134 return *this;
1135}
1136
1139 const ConstComplexMatrixView& x) {
1140 ARTS_ASSERT(nrows() == x.nrows());
1141 ARTS_ASSERT(ncols() == x.ncols());
1144 const ComplexIterator2D er = end();
1145 for (; r != er; ++r, ++sr) {
1146 ConstComplexIterator1D sc = sr->begin();
1147 ComplexIterator1D c = r->begin();
1148 const ComplexIterator1D ec = r->end();
1149 for (; c != ec; ++c, ++sc) *c *= *sc;
1150 }
1151 return *this;
1152}
1153
1156 ARTS_ASSERT(nrows() == x.nrows());
1157 ARTS_ASSERT(ncols() == x.ncols());
1158 ConstIterator2D sr = x.begin();
1160 const ComplexIterator2D er = end();
1161 for (; r != er; ++r, ++sr) {
1162 ConstIterator1D sc = sr->begin();
1163 ComplexIterator1D c = r->begin();
1164 const ComplexIterator1D ec = r->end();
1165 for (; c != ec; ++c, ++sc) *c *= *sc;
1166 }
1167 return *this;
1168}
1169
1172 const ConstComplexMatrixView& x) {
1173 ARTS_ASSERT(nrows() == x.nrows());
1174 ARTS_ASSERT(ncols() == x.ncols());
1177 const ComplexIterator2D er = end();
1178 for (; r != er; ++r, ++sr) {
1179 ConstComplexIterator1D sc = sr->begin();
1180 ComplexIterator1D c = r->begin();
1181 const ComplexIterator1D ec = r->end();
1182 for (; c != ec; ++c, ++sc) *c /= *sc;
1183 }
1184 return *this;
1185}
1186
1189 ARTS_ASSERT(nrows() == x.nrows());
1190 ARTS_ASSERT(ncols() == x.ncols());
1191 ConstIterator2D sr = x.begin();
1193 const ComplexIterator2D er = end();
1194 for (; r != er; ++r, ++sr) {
1195 ConstIterator1D sc = sr->begin();
1196 ComplexIterator1D c = r->begin();
1197 const ComplexIterator1D ec = r->end();
1198 for (; c != ec; ++c, ++sc) *c /= *sc;
1199 }
1200 return *this;
1201}
1202
1205 const ConstComplexMatrixView& x) {
1206 ARTS_ASSERT(nrows() == x.nrows());
1207 ARTS_ASSERT(ncols() == x.ncols());
1210 const ComplexIterator2D er = end();
1211 for (; r != er; ++r, ++sr) {
1212 ConstComplexIterator1D sc = sr->begin();
1213 ComplexIterator1D c = r->begin();
1214 const ComplexIterator1D ec = r->end();
1215 for (; c != ec; ++c, ++sc) *c += *sc;
1216 }
1217 return *this;
1218}
1219
1222 ARTS_ASSERT(nrows() == x.nrows());
1223 ARTS_ASSERT(ncols() == x.ncols());
1224 ConstIterator2D sr = x.begin();
1226 const ComplexIterator2D er = end();
1227 for (; r != er; ++r, ++sr) {
1228 ConstIterator1D sc = sr->begin();
1229 ComplexIterator1D c = r->begin();
1230 const ComplexIterator1D ec = r->end();
1231 for (; c != ec; ++c, ++sc) *c += *sc;
1232 }
1233 return *this;
1234}
1235
1238 const ConstComplexMatrixView& x) {
1239 ARTS_ASSERT(nrows() == x.nrows());
1240 ARTS_ASSERT(ncols() == x.ncols());
1243 const ComplexIterator2D er = end();
1244 for (; r != er; ++r, ++sr) {
1245 ConstComplexIterator1D sc = sr->begin();
1246 ComplexIterator1D c = r->begin();
1247 const ComplexIterator1D ec = r->end();
1248 for (; c != ec; ++c, ++sc) *c -= *sc;
1249 }
1250 return *this;
1251}
1252
1255 ARTS_ASSERT(nrows() == x.nrows());
1256 ARTS_ASSERT(ncols() == x.ncols());
1257 ConstIterator2D sr = x.begin();
1259 const ComplexIterator2D er = end();
1260 for (; r != er; ++r, ++sr) {
1261 ConstIterator1D sc = sr->begin();
1262 ComplexIterator1D c = r->begin();
1263 const ComplexIterator1D ec = r->end();
1264 for (; c != ec; ++c, ++sc) *c -= *sc;
1265 }
1266 return *this;
1267}
1268
1272 // Nothing to do here.
1273}
1274
1279 const Range& rr,
1280 const Range& cr)
1281 : ConstComplexMatrixView(data, rr, cr) {
1282 // Nothing to do here.
1283}
1284
1300 const Range& pr,
1301 const Range& pc,
1302 const Range& nr,
1303 const Range& nc)
1304 : ConstComplexMatrixView(data, pr, pc, nr, nc) {
1305 // Nothing to do here.
1306}
1307
1318 const ConstComplexIterator2D& end,
1319 ComplexIterator2D target) {
1320 for (; origin != end; ++origin, ++target) {
1321 ConstComplexIterator1D o = origin->begin();
1322 const ConstComplexIterator1D e = origin->end();
1323 ComplexIterator1D t = target->begin();
1324 for (; o != e; ++o, ++t) *t = *o;
1325 }
1326}
1327
1330 for (; target != end; ++target) {
1331 ComplexIterator1D t = target->begin();
1332 const ComplexIterator1D e = target->end();
1333 for (; t != e; ++t) *t = x;
1334 }
1335}
1336
1339 for (; target != end; ++target) {
1340 ComplexIterator1D t = target->begin();
1341 const ComplexIterator1D e = target->end();
1342 for (; t != e; ++t) *t = x;
1343 }
1344}
1345
1346// Functions for ComplexMatrix:
1347// ---------------------
1348
1352 : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
1353 // Nothing to do here.
1354}
1355
1358 : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
1359 // Here we can access the raw memory directly, for slightly
1360 // increased efficiency:
1361 const Complex* stop = mdata + r * c;
1362 for (Complex* x = mdata; x < stop; ++x) *x = fill;
1363}
1364
1367 : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
1368 // Here we can access the raw memory directly, for slightly
1369 // increased efficiency:
1370 const Complex* stop = mdata + r * c;
1371 for (Complex* x = mdata; x < stop; ++x) *x = fill;
1372}
1373
1377 : ComplexMatrixView(new Complex[m.nrows() * m.ncols()],
1378 Range(0, m.nrows(), m.ncols()),
1379 Range(0, m.ncols())) {
1380 copy(m.begin(), m.end(), begin());
1381}
1382
1386 : ComplexMatrixView(new Complex[m.nrows() * m.ncols()],
1387 Range(0, m.nrows(), m.ncols()),
1388 Range(0, m.ncols())) {
1389 // There is a catch here: If m is an empty matrix, then it will have
1390 // 0 colunns. But this is used to initialize the stride of the row
1391 // Range! Thus, this method has to be consistent with the behaviour
1392 // of Range::Range. For now, Range::Range allows also stride 0.
1393 copy(m.begin(), m.end(), begin());
1394}
1395
1397
1422 swap(*this, m);
1423 return *this;
1424}
1425
1429 copy(x, begin(), end());
1430 return *this;
1431}
1432
1434
1447 resize(v.nelem(), 1);
1449 copy(dummy.begin(), dummy.end(), begin());
1450 return *this;
1451}
1452
1457 ARTS_ASSERT(0 <= r);
1458 ARTS_ASSERT(0 <= c);
1459
1460 if (mrr.mextent != r || mcr.mextent != c) {
1461 delete[] mdata;
1462 mdata = new Complex[r * c];
1463
1464 mrr.mstart = 0;
1465 mrr.mextent = r;
1466 mrr.mstride = c;
1467
1468 mcr.mstart = 0;
1469 mcr.mextent = c;
1470 mcr.mstride = 1;
1471 }
1472}
1473
1475void swap(ComplexMatrix& m1, ComplexMatrix& m2) noexcept {
1476 using std::swap;
1477 swap(m1.mrr, m2.mrr);
1478 swap(m1.mcr, m2.mcr);
1479 swap(m1.mdata, m2.mdata);
1480}
1481
1485 // cout << "Destroying a Matrix:\n"
1486 // << *this << "\n........................................\n";
1487 delete[] mdata;
1488}
1489
1491{
1492 ARTS_ASSERT(ncols() == nrows());
1493
1494 // help's internal variables are all mutable, so const is just for default parameter and to not use copies
1495 help.resize_if_smaller(int(ncols()));
1496
1497 // Compute LU decomposition using LAPACK dgetrf_.
1498 int info;
1499 lapack::zgetrf_(help.size(), help.size(), mdata, help.size(), help.ipivdata(), &info);
1500 lapack::zgetri_(help.size(), mdata, help.size(), help.ipivdata(), help.workdata(), help.lsize(), &info);
1501
1502 // Check for success.
1503 ARTS_USER_ERROR_IF (info not_eq 0,
1504 "Error inverting matrix: Matrix not of full rank.");
1505
1506 return *this;
1507}
1508
1511 return ConstComplexMatrixView(m.mdata, m.mcr, m.mrr);
1512}
1513
1517 return ComplexMatrixView(m.mdata, m.mcr, m.mrr);
1518}
1519
1523 return transpose((ComplexMatrixView)v);
1524}
1525
1536 // cout << "Assigning VectorView from Array<Numeric>.\n";
1537
1538 // Check that sizes are compatible:
1539 ARTS_ASSERT(mrange.mextent == v.nelem());
1540
1541 // Iterators for Array:
1542 auto i = v.begin();
1543 const auto e = v.end();
1544 // Iterator for Vector:
1545 ComplexIterator1D target = begin();
1546
1547 for (; i != e; ++i, ++target) *target = *i;
1548
1549 return *this;
1550}
1551
1552//Functions operating on the complex vectors
1553
1556 const ConstComplexVectorView& b) {
1557 // Check dimensions:
1558 ARTS_ASSERT(a.nelem() == b.nelem());
1559
1560 const ConstComplexIterator1D ae = a.end();
1561 ConstComplexIterator1D ai = a.begin();
1562 ConstComplexIterator1D bi = b.begin();
1563
1564 Complex res = 0;
1565 for (; ai != ae; ++ai, ++bi) res += (*ai) * (*bi);
1566
1567 return res;
1568}
1569
1571
1582 const ConstComplexVectorView& x) {
1583 ARTS_ASSERT(x.nelem() == M.nrows());
1584 ARTS_ASSERT(y.nelem() == M.ncols());
1585
1586 auto eigen_y = matpack::eigen::row_vec(y);
1587 if (y.mdata == x.mdata)
1589 else
1590 eigen_y.noalias() = matpack::eigen::mat(M) * matpack::eigen::row_vec(x);
1591}
1592
1594
1606 const ConstComplexMatrixView& B,
1607 const ConstComplexMatrixView& C) {
1608 ARTS_ASSERT(B.nrows() == A.nrows());
1609 ARTS_ASSERT(C.ncols() == A.ncols());
1610 ARTS_ASSERT(B.ncols() == C.nrows());
1611
1612 auto eigen_A = matpack::eigen::mat(A);
1613 if (A.mdata == B.mdata || A.mdata == C.mdata)
1614 eigen_A = matpack::eigen::mat(B) * matpack::eigen::mat(C);
1615 else
1616 eigen_A.noalias() = matpack::eigen::mat(B) * matpack::eigen::mat(C);
1617}
1618
1620 const ConstComplexMatrixView& B,
1621 const ConstMatrixView& C) {
1622 ARTS_ASSERT(B.nrows() == A.nrows());
1623 ARTS_ASSERT(C.ncols() == A.ncols());
1624 ARTS_ASSERT(B.ncols() == C.nrows());
1625
1626 auto eigen_A = matpack::eigen::mat(A);
1627 if (A.mdata == B.mdata)
1628 eigen_A = matpack::eigen::mat(B) * matpack::eigen::mat(C);
1629 else
1630 eigen_A.noalias() = matpack::eigen::mat(B) * matpack::eigen::mat(C);
1631}
1632
1634 const ConstMatrixView& B,
1635 const ConstComplexMatrixView& C) {
1636 ARTS_ASSERT(B.nrows() == A.nrows());
1637 ARTS_ASSERT(C.ncols() == A.ncols());
1638 ARTS_ASSERT(B.ncols() == C.nrows());
1639
1640 auto eigen_A = matpack::eigen::mat(A);
1641 if (A.mdata == C.mdata)
1642 eigen_A = matpack::eigen::mat(B) * matpack::eigen::mat(C);
1643 else
1644 eigen_A.noalias() = matpack::eigen::mat(B) * matpack::eigen::mat(C);
1645}
1646
1648 const ConstMatrixView& B,
1649 const ConstMatrixView& C) {
1650 ARTS_ASSERT(B.nrows() == A.nrows());
1651 ARTS_ASSERT(C.ncols() == A.ncols());
1652 ARTS_ASSERT(B.ncols() == C.nrows());
1653
1655}
1656
1658// Helper function for debugging
1659#ifndef NDEBUG
1660
1675 return mv(r, c);
1676}
1677
1678#endif
Interface for BLAS library.
This can be used to make arrays out of anything.
Definition: array.h:48
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
The iterator class for sub vectors.
Index mstride
Stride.
Complex * mx
Current position.
The row iterator class for sub matrices.
The ComplexMatrixView class.
ComplexMatrixView & operator/=(Complex x)
Division by scalar.
ComplexVectorView diagonal()
ComplexMatrix diagonal as vector.
Complex & operator()(Index r, Index c)
Plain index operator.
ComplexMatrixView & operator+=(Complex x)
Addition of scalar.
MatrixView imag()
Get a view of the imaginary parts of the matrix.
MatrixView real()
Get a view of the real part of the matrix.
ComplexIterator2D end()
Return iterator behind last row.
ComplexMatrixView & operator=(const ConstComplexMatrixView &v)
Assignment operator.
ComplexMatrixView & operator-=(Complex x)
Subtraction of scalar.
ComplexMatrixView()
Default constructor.
ComplexMatrixView & operator*=(Complex x)
Multiplication by scalar.
friend class ComplexVectorView
ComplexIterator2D begin()
Return iterator to first row.
The ComplexMatrix class.
void resize(Index r, Index c)
Resize function.
~ComplexMatrix() noexcept override
Destructor for Matrix.
ComplexMatrix & inv(const lapack_help::Inverse< Complex > &help=lapack_help::Inverse< Complex >{0})
friend void swap(ComplexMatrix &m1, ComplexMatrix &m2) noexcept
Swaps two objects.
ComplexMatrix & operator=(ComplexMatrix x)
Assignment operator from another matrix.
ComplexMatrix()=default
The ComplexVectorView class.
VectorView real()
Get a view of the real part of the vector.
ComplexVectorView operator-=(Complex x)
Subtraction of scalar.
ComplexVectorView operator+=(Complex x)
Addition of scalar.
VectorView imag()
Get a view of the imaginary part of the vector.
ComplexVectorView operator*=(Complex x)
Multiplication by scalar.
ComplexVectorView & operator=(const ConstComplexVectorView &v)
Assignment operator.
ComplexVectorView()=default
ComplexIterator1D end()
Return iterator behind last element.
Complex & operator[](Index n)
Plain Index operator.
ComplexVectorView operator/=(Complex x)
Division by scalar.
ComplexIterator1D begin()
Return iterator to first element.
The ComplexVector class.
ComplexVector()=default
~ComplexVector() noexcept override
Destructor for ComplexVector.
void resize(Index n)
Resize function.
ComplexVector & operator=(const ComplexVector &v)
Assignment from another Vector.
The constant iterator class for sub vectors.
const Complex * mx
Current position.
The const row iterator class for sub matrices.
A constant view of a ComplexMatrix.
friend class ConstComplexVectorView
ConstComplexIterator2D begin() const
Return const iterator to first row.
Range mcr
The column range of mdata that is actually used.
ConstComplexVectorView diagonal() const
ComplexMatrix diagonal as vector.
Complex operator()(Index r, Index c) const
Plain const index operator.
Complex * mdata
Pointer to the plain C array that holds the data.
Index ncols() const noexcept
ConstComplexMatrixView()=default
Range mrr
The row range of mdata that is actually used.
Index nrows() const noexcept
ConstMatrixView imag() const
Get a view of the imaginary part of the matrix.
ConstComplexIterator2D end() const
Return const iterator behind last row.
ConstMatrixView real() const
Get a view of the real part of the matrix.
A constant view of a ComplexVector.
Complex sum() const
Returns true if variable size is zero.
ConstVectorView imag() const
Get a view of the imaginary part of the vector.
ConstComplexIterator1D end() const
Return const iterator behind last element.
ConstVectorView real() const
Get a view of the real part of the vector.
Range mrange
The range of mdata that is actually used.
ConstComplexVectorView()=default
Complex * mdata
Pointer to the plain C array that holds the data.
ConstComplexIterator1D begin() const
Return const iterator to first element.
Index nelem() const noexcept
const Complex & operator[](Index n) const
Plain const index operator.
The constant iterator class for sub vectors.
Definition: matpackI.h:458
The const row iterator class for sub matrices.
Definition: matpackI.h:859
A constant view of a Matrix.
Definition: matpackI.h:1065
ConstIterator2D begin() const ARTS_NOEXCEPT
Return const iterator to first row.
Definition: matpackI.cc:449
Index nrows() const noexcept
Definition: matpackI.h:1079
Index ncols() const noexcept
Definition: matpackI.h:1080
A constant view of a Vector.
Definition: matpackI.h:521
ConstIterator1D begin() const ARTS_NOEXCEPT
Return const iterator to first element.
Definition: matpackI.cc:63
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
The MatrixView class.
Definition: matpackI.h:1188
The range class.
Definition: matpackI.h:160
Index mstart
The start index.
Definition: matpackI.h:377
Index mstride
The stride.
Definition: matpackI.h:382
Index mextent
The number of elements.
Definition: matpackI.h:380
The VectorView class.
Definition: matpackI.h:674
The Vector class.
Definition: matpackI.h:910
#define ARTS_ASSERT(condition,...)
Definition: debug.h:102
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:153
The declarations of all the exception classes.
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
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
void swap(ComplexVector &v1, ComplexVector &v2) noexcept
Swaps two objects.
Complex debug_matrixview_get_elem(ComplexMatrixView &mv, Index r, Index c)
Helper function to access matrix elements.
Complex operator*(const ConstComplexVectorView &a, const ConstComplexVectorView &b)
Scalar product.
std::ostream & operator<<(std::ostream &os, const ConstComplexVectorView &v)
Output operator.
void copy(ConstComplexIterator1D origin, const ConstComplexIterator1D &end, ComplexIterator1D target)
Copy data between begin and end to target.
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
std::complex< Numeric > Complex
void copy(ConstComplexIterator1D origin, const ConstComplexIterator1D &end, ComplexIterator1D target)
Copy data between begin and end to target.
void zgetri_(int *n, std::complex< double > *A, int *lda, int *ipiv, std::complex< double > *work, int *lwork, int *info)
void zgetrf_(int *m, int *n, std::complex< double > *A, int *lda, int *ipiv, int *info)
auto row_vec(matpack::vector auto &&x)
Map the input to a non-owning const-correct Eigen Map representing a column vector.
Definition: matpack_eigen.h:71
auto mat(matpack::matrix auto &&x)
Map the input to a non-owning const-correct Eigen Map representing a row matrix.
Definition: matpack_eigen.h:91
#define M
Definition: rng.cc:165
Struct cannot be const, but can be passed as const to allow defaults.
Definition: lapack.h:23
T * workdata() const noexcept
Definition: lapack.h:45
bool resize_if_smaller(Index N) const
Definition: lapack.h:31
int * lsize() const noexcept
Definition: lapack.h:44
int * size() const noexcept
Definition: lapack.h:43
int * ipivdata() const noexcept
Definition: lapack.h:46
#define v
#define a
#define c
#define b