ARTS 2.5.4 (git: 31ce4f0e)
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;
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) {
333 ARTS_ASSERT(nelem() == x.nelem());
336
338 const ComplexIterator1D e = end();
339
340 for (; i != e; ++i, ++s) *i /= *s;
341 return *this;
342}
343
346 ARTS_ASSERT(nelem() == x.nelem());
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
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
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 const auto n = v.nelem();
649 resize(n);
650
651 auto* in = v.mdata + v.mrange.mstart;
652 auto* end = v.mdata + n;
653 auto* out = mdata;
654
655 while (in < end) {
656 out = in;
657 out++;
658 in += v.mrange.mstride;
659 }
660
661 return *this;
662}
663
665
682 resize(x.nelem());
684 return *this;
685}
686
691 return *this;
692}
693
694// /** Assignment operator from VectorView. Assignment operators are not
695// inherited. */
696// inline Vector& Vector::operator=(const VectorView v)
697// {
698// cout << "Assigning Vector from Vector View.\n";
699// // Check that sizes are compatible:
700// ARTS_ASSERT(mrange.mextent==v.mrange.mextent);
701// VectorView::operator=(v);
702// return *this;
703// }
704
709 ARTS_ASSERT(0 <= n);
710 if (mrange.mextent != n) {
711 delete[] mdata;
712 mdata = new Complex[n];
713 mrange.mstart = 0;
714 mrange.mextent = n;
715 mrange.mstride = 1;
716 }
717}
721 std::swap(v1.mrange, v2.mrange);
723}
724
728
729// Functions for ConstMatrixView:
730// ------------------------------
731
733// bool ConstComplexMatrixView::empty() const
734// {
735// return (nrows() == 0 || ncols() == 0);
736// }
739// Index ConstComplexMatrixView::nrows() const
740// {
741// return mrr.mextent;
742// }
743
745// Index ConstComplexMatrixView::ncols() const
746// {
747// return mcr.mextent;
748// }
749
754 const Range& r, const Range& c) const {
755 return ConstComplexMatrixView(mdata, mrr, mcr, r, c);
756}
757
764 Index c) const {
765 // Check that c is valid:
766 ARTS_ASSERT(0 <= c);
768
770}
771
778 Index r, const Range& c) const {
779 // Check that r is valid:
780 ARTS_ASSERT(0 <= r);
782
784}
785
789 mrr.mstride);
790}
791
796 mcr),
797 mrr.mstride);
798}
799
801
812 Range(0, n, mrr.mstride + mcr.mstride));
813}
814
819 const Range& rr,
820 const Range& cr)
821 : mrr(rr), mcr(cr), mdata(data) {
822 // Nothing to do here.
823}
824
840 const Range& pr,
841 const Range& pc,
842 const Range& nr,
843 const Range& nc)
844 : mrr(pr, nr), mcr(pc, nc), mdata(data) {
845 // Nothing to do here.
846}
847
855std::ostream& operator<<(std::ostream& os, const ConstComplexMatrixView& v) {
856 // Row iterators:
857 ConstComplexIterator2D ir = v.begin();
858 const ConstComplexIterator2D end_row = v.end();
859
860 if (ir != end_row) {
861 ConstComplexIterator1D ic = ir->begin();
862 const ConstComplexIterator1D end_col = ir->end();
863
864 if (ic != end_col) {
865 os << std::setw(3) << *ic;
866 ++ic;
867 }
868 for (; ic != end_col; ++ic) {
869 os << " " << std::setw(3) << *ic;
870 }
871 ++ir;
872 }
873 for (; ir != end_row; ++ir) {
874 ConstComplexIterator1D ic = ir->begin();
875 const ConstComplexIterator1D end_col = ir->end();
876
877 os << "\n";
878 if (ic != end_col) {
879 os << std::setw(3) << *ic;
880 ++ic;
881 }
882 for (; ic != end_col; ++ic) {
883 os << " " << std::setw(3) << *ic;
884 }
885 }
886
887 return os;
888}
889
890// Functions for ComplexMatrixView:
891// -------------------------
892
897 const Range& c) {
898 return ComplexMatrixView(mdata, mrr, mcr, r, c);
899}
900
906 // Check that c is valid:
907 ARTS_ASSERT(0 <= c);
909
910 return ComplexVectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
911}
912
918 // Check that r is valid:
919 ARTS_ASSERT(0 <= r);
921
922 return ComplexVectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
923}
924
928 mrr.mstride);
929}
930
933 return ComplexIterator2D(
935 mrr.mstride);
936}
937
939
950 Range(0, n, mrr.mstride + mcr.mstride));
951}
952
958 const ConstComplexMatrixView& m) {
959 // Check that sizes are compatible:
962
963 copy(m.begin(), m.end(), begin());
964 return *this;
965}
966
973 // Check that sizes are compatible:
976
977 copy(m.begin(), m.end(), begin());
978 return *this;
979}
980
985 // Check that sizes are compatible:
988
989 copy(m.begin(), m.end(), begin());
990 return *this;
991}
992
998 const ConstComplexVectorView& v) {
999 // Check that sizes are compatible:
1000 ARTS_ASSERT(mrr.mextent == v.nelem());
1001 ARTS_ASSERT(mcr.mextent == 1);
1002 // dummy = ConstComplexMatrixView(v.mdata,v.mrange,Range(v.mrange.mstart,1));;
1004 copy(dummy.begin(), dummy.end(), begin());
1005 return *this;
1006}
1007
1011 copy(x, begin(), end());
1012 return *this;
1013}
1014
1017 const ComplexIterator2D er = end();
1018 for (ComplexIterator2D r = begin(); r != er; ++r) {
1019 const ComplexIterator1D ec = r->end();
1020 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c *= x;
1021 }
1022 return *this;
1023}
1024
1027 const ComplexIterator2D er = end();
1028 for (ComplexIterator2D r = begin(); r != er; ++r) {
1029 const ComplexIterator1D ec = r->end();
1030 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c *= x;
1031 }
1032 return *this;
1033}
1034
1037 const ComplexIterator2D er = end();
1038 for (ComplexIterator2D r = begin(); r != er; ++r) {
1039 const ComplexIterator1D ec = r->end();
1040 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c /= x;
1041 }
1042 return *this;
1043}
1044
1047 const ComplexIterator2D er = end();
1048 for (ComplexIterator2D r = begin(); r != er; ++r) {
1049 const ComplexIterator1D ec = r->end();
1050 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c /= x;
1051 }
1052 return *this;
1053}
1054
1057 const ComplexIterator2D er = end();
1058 for (ComplexIterator2D r = begin(); r != er; ++r) {
1059 const ComplexIterator1D ec = r->end();
1060 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c += x;
1061 }
1062 return *this;
1063}
1064
1067 const ComplexIterator2D er = end();
1068 for (ComplexIterator2D r = begin(); r != er; ++r) {
1069 const ComplexIterator1D ec = r->end();
1070 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c += x;
1071 }
1072 return *this;
1073}
1074
1077 const ComplexIterator2D er = end();
1078 for (ComplexIterator2D r = begin(); r != er; ++r) {
1079 const ComplexIterator1D ec = r->end();
1080 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c -= x;
1081 }
1082 return *this;
1083}
1084
1087 const ComplexIterator2D er = end();
1088 for (ComplexIterator2D r = begin(); r != er; ++r) {
1089 const ComplexIterator1D ec = r->end();
1090 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c -= x;
1091 }
1092 return *this;
1093}
1094
1097 const ConstComplexMatrixView& x) {
1098 ARTS_ASSERT(nrows() == x.nrows());
1099 ARTS_ASSERT(ncols() == x.ncols());
1102 const ComplexIterator2D er = end();
1103 for (; r != er; ++r, ++sr) {
1104 ConstComplexIterator1D sc = sr->begin();
1105 ComplexIterator1D c = r->begin();
1106 const ComplexIterator1D ec = r->end();
1107 for (; c != ec; ++c, ++sc) *c *= *sc;
1108 }
1109 return *this;
1110}
1111
1114 ARTS_ASSERT(nrows() == x.nrows());
1115 ARTS_ASSERT(ncols() == x.ncols());
1116 ConstIterator2D sr = x.begin();
1118 const ComplexIterator2D er = end();
1119 for (; r != er; ++r, ++sr) {
1120 ConstIterator1D sc = sr->begin();
1121 ComplexIterator1D c = r->begin();
1122 const ComplexIterator1D ec = r->end();
1123 for (; c != ec; ++c, ++sc) *c *= *sc;
1124 }
1125 return *this;
1126}
1127
1130 const ConstComplexMatrixView& x) {
1131 ARTS_ASSERT(nrows() == x.nrows());
1132 ARTS_ASSERT(ncols() == x.ncols());
1135 const ComplexIterator2D er = end();
1136 for (; r != er; ++r, ++sr) {
1137 ConstComplexIterator1D sc = sr->begin();
1138 ComplexIterator1D c = r->begin();
1139 const ComplexIterator1D ec = r->end();
1140 for (; c != ec; ++c, ++sc) *c /= *sc;
1141 }
1142 return *this;
1143}
1144
1147 ARTS_ASSERT(nrows() == x.nrows());
1148 ARTS_ASSERT(ncols() == x.ncols());
1149 ConstIterator2D sr = x.begin();
1151 const ComplexIterator2D er = end();
1152 for (; r != er; ++r, ++sr) {
1153 ConstIterator1D sc = sr->begin();
1154 ComplexIterator1D c = r->begin();
1155 const ComplexIterator1D ec = r->end();
1156 for (; c != ec; ++c, ++sc) *c /= *sc;
1157 }
1158 return *this;
1159}
1160
1163 const ConstComplexMatrixView& x) {
1164 ARTS_ASSERT(nrows() == x.nrows());
1165 ARTS_ASSERT(ncols() == x.ncols());
1168 const ComplexIterator2D er = end();
1169 for (; r != er; ++r, ++sr) {
1170 ConstComplexIterator1D sc = sr->begin();
1171 ComplexIterator1D c = r->begin();
1172 const ComplexIterator1D ec = r->end();
1173 for (; c != ec; ++c, ++sc) *c += *sc;
1174 }
1175 return *this;
1176}
1177
1180 ARTS_ASSERT(nrows() == x.nrows());
1181 ARTS_ASSERT(ncols() == x.ncols());
1182 ConstIterator2D sr = x.begin();
1184 const ComplexIterator2D er = end();
1185 for (; r != er; ++r, ++sr) {
1186 ConstIterator1D sc = sr->begin();
1187 ComplexIterator1D c = r->begin();
1188 const ComplexIterator1D ec = r->end();
1189 for (; c != ec; ++c, ++sc) *c += *sc;
1190 }
1191 return *this;
1192}
1193
1196 const ConstComplexMatrixView& x) {
1197 ARTS_ASSERT(nrows() == x.nrows());
1198 ARTS_ASSERT(ncols() == x.ncols());
1201 const ComplexIterator2D er = end();
1202 for (; r != er; ++r, ++sr) {
1203 ConstComplexIterator1D sc = sr->begin();
1204 ComplexIterator1D c = r->begin();
1205 const ComplexIterator1D ec = r->end();
1206 for (; c != ec; ++c, ++sc) *c -= *sc;
1207 }
1208 return *this;
1209}
1210
1213 ARTS_ASSERT(nrows() == x.nrows());
1214 ARTS_ASSERT(ncols() == x.ncols());
1215 ConstIterator2D sr = x.begin();
1217 const ComplexIterator2D er = end();
1218 for (; r != er; ++r, ++sr) {
1219 ConstIterator1D sc = sr->begin();
1220 ComplexIterator1D c = r->begin();
1221 const ComplexIterator1D ec = r->end();
1222 for (; c != ec; ++c, ++sc) *c -= *sc;
1223 }
1224 return *this;
1225}
1226
1230 // Nothing to do here.
1231}
1232
1237 const Range& rr,
1238 const Range& cr)
1239 : ConstComplexMatrixView(data, rr, cr) {
1240 // Nothing to do here.
1241}
1242
1258 const Range& pr,
1259 const Range& pc,
1260 const Range& nr,
1261 const Range& nc)
1262 : ConstComplexMatrixView(data, pr, pc, nr, nc) {
1263 // Nothing to do here.
1264}
1265
1277 ComplexIterator2D target) {
1278 for (; origin != end; ++origin, ++target) {
1279 ConstComplexIterator1D o = origin->begin();
1280 const ConstComplexIterator1D e = origin->end();
1281 ComplexIterator1D t = target->begin();
1282 for (; o != e; ++o, ++t) *t = *o;
1283 }
1284}
1285
1288 for (; target != end; ++target) {
1289 ComplexIterator1D t = target->begin();
1290 const ComplexIterator1D e = target->end();
1291 for (; t != e; ++t) *t = x;
1292 }
1293}
1294
1297 for (; target != end; ++target) {
1298 ComplexIterator1D t = target->begin();
1299 const ComplexIterator1D e = target->end();
1300 for (; t != e; ++t) *t = x;
1301 }
1302}
1303
1304// Functions for ComplexMatrix:
1305// ---------------------
1306
1310 : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
1311 // Nothing to do here.
1312}
1313
1316 : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
1317 // Here we can access the raw memory directly, for slightly
1318 // increased efficiency:
1319 const Complex* stop = mdata + r * c;
1320 for (Complex* x = mdata; x < stop; ++x) *x = fill;
1321}
1322
1325 : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
1326 // Here we can access the raw memory directly, for slightly
1327 // increased efficiency:
1328 const Complex* stop = mdata + r * c;
1329 for (Complex* x = mdata; x < stop; ++x) *x = fill;
1330}
1331
1335 : ComplexMatrixView(new Complex[m.nrows() * m.ncols()],
1336 Range(0, m.nrows(), m.ncols()),
1337 Range(0, m.ncols())) {
1338 copy(m.begin(), m.end(), begin());
1339}
1340
1344 : ComplexMatrixView(new Complex[m.nrows() * m.ncols()],
1345 Range(0, m.nrows(), m.ncols()),
1346 Range(0, m.ncols())) {
1347 // There is a catch here: If m is an empty matrix, then it will have
1348 // 0 colunns. But this is used to initialize the stride of the row
1349 // Range! Thus, this method has to be consistent with the behaviour
1350 // of Range::Range. For now, Range::Range allows also stride 0.
1351 copy(m.begin(), m.end(), begin());
1352}
1353
1355
1380 swap(*this, m);
1381 return *this;
1382}
1383
1387 copy(x, begin(), end());
1388 return *this;
1389}
1390
1392
1405 resize(v.nelem(), 1);
1407 copy(dummy.begin(), dummy.end(), begin());
1408 return *this;
1409}
1410
1415 ARTS_ASSERT(0 <= r);
1416 ARTS_ASSERT(0 <= c);
1417
1418 if (mrr.mextent != r || mcr.mextent != c) {
1419 delete[] mdata;
1420 mdata = new Complex[r * c];
1421
1422 mrr.mstart = 0;
1423 mrr.mextent = r;
1424 mrr.mstride = c;
1425
1426 mcr.mstart = 0;
1427 mcr.mextent = c;
1428 mcr.mstride = 1;
1429 }
1430}
1431
1434 std::swap(m1.mrr, m2.mrr);
1435 std::swap(m1.mcr, m2.mcr);
1436 std::swap(m1.mdata, m2.mdata);
1437}
1438
1442 // cout << "Destroying a Matrix:\n"
1443 // << *this << "\n........................................\n";
1444 delete[] mdata;
1445}
1446
1448{
1449 ARTS_ASSERT(ncols() == nrows());
1450
1451 // help's internal variables are all mutable, so const is just for default parameter and to not use copies
1452 help.resize_if_smaller(int(ncols()));
1453
1454 // Compute LU decomposition using LAPACK dgetrf_.
1455 int info;
1456 lapack::zgetrf_(help.size(), help.size(), mdata, help.size(), help.ipivdata(), &info);
1457 lapack::zgetri_(help.size(), mdata, help.size(), help.ipivdata(), help.workdata(), help.lsize(), &info);
1458
1459 // Check for success.
1460 ARTS_USER_ERROR_IF (info not_eq 0,
1461 "Error inverting matrix: Matrix not of full rank.");
1462
1463 return *this;
1464}
1465
1468 return ConstComplexMatrixView(m.mdata, m.mcr, m.mrr);
1469}
1470
1474 return ComplexMatrixView(m.mdata, m.mcr, m.mrr);
1475}
1476
1480 return transpose((ComplexMatrixView)v);
1481}
1482
1493 // cout << "Assigning VectorView from Array<Numeric>.\n";
1494
1495 // Check that sizes are compatible:
1496 ARTS_ASSERT(mrange.mextent == v.nelem());
1497
1498 // Iterators for Array:
1499 auto i = v.begin();
1500 const auto e = v.end();
1501 // Iterator for Vector:
1502 ComplexIterator1D target = begin();
1503
1504 for (; i != e; ++i, ++target) *target = *i;
1505
1506 return *this;
1507}
1508
1509//Functions operating on the complex vectors
1510
1513 const ConstComplexVectorView& b) {
1514 // Check dimensions:
1515 ARTS_ASSERT(a.nelem() == b.nelem());
1516
1517 const ConstComplexIterator1D ae = a.end();
1518 ConstComplexIterator1D ai = a.begin();
1519 ConstComplexIterator1D bi = b.begin();
1520
1521 Complex res = 0;
1522 for (; ai != ae; ++ai, ++bi) res += (*ai) * (*bi);
1523
1524 return res;
1525}
1526
1528
1539 const ConstComplexVectorView& x) {
1540 ARTS_ASSERT(x.nelem() == M.nrows());
1541 ARTS_ASSERT(y.nelem() == M.ncols());
1542
1543 auto eigen_y = matpack::eigen::row_vec(y);
1544 if (y.mdata == x.mdata)
1546 else
1547 eigen_y.noalias() = matpack::eigen::mat(M) * matpack::eigen::row_vec(x);
1548}
1549
1551
1563 const ConstComplexMatrixView& B,
1564 const ConstComplexMatrixView& C) {
1565 ARTS_ASSERT(B.nrows() == A.nrows());
1566 ARTS_ASSERT(C.ncols() == A.ncols());
1567 ARTS_ASSERT(B.ncols() == C.nrows());
1568
1569 auto eigen_A = matpack::eigen::mat(A);
1570 if (A.mdata == B.mdata || A.mdata == C.mdata)
1571 eigen_A = matpack::eigen::mat(B) * matpack::eigen::mat(C);
1572 else
1573 eigen_A.noalias() = matpack::eigen::mat(B) * matpack::eigen::mat(C);
1574}
1575
1577 const ConstComplexMatrixView& B,
1578 const ConstMatrixView& C) {
1579 ARTS_ASSERT(B.nrows() == A.nrows());
1580 ARTS_ASSERT(C.ncols() == A.ncols());
1581 ARTS_ASSERT(B.ncols() == C.nrows());
1582
1583 auto eigen_A = matpack::eigen::mat(A);
1584 if (A.mdata == B.mdata)
1585 eigen_A = matpack::eigen::mat(B) * matpack::eigen::mat(C);
1586 else
1587 eigen_A.noalias() = matpack::eigen::mat(B) * matpack::eigen::mat(C);
1588}
1589
1591 const ConstMatrixView& B,
1592 const ConstComplexMatrixView& C) {
1593 ARTS_ASSERT(B.nrows() == A.nrows());
1594 ARTS_ASSERT(C.ncols() == A.ncols());
1595 ARTS_ASSERT(B.ncols() == C.nrows());
1596
1597 auto eigen_A = matpack::eigen::mat(A);
1598 if (A.mdata == C.mdata)
1599 eigen_A = matpack::eigen::mat(B) * matpack::eigen::mat(C);
1600 else
1601 eigen_A.noalias() = matpack::eigen::mat(B) * matpack::eigen::mat(C);
1602}
1603
1605 const ConstMatrixView& B,
1606 const ConstMatrixView& C) {
1607 ARTS_ASSERT(B.nrows() == A.nrows());
1608 ARTS_ASSERT(C.ncols() == A.ncols());
1609 ARTS_ASSERT(B.ncols() == C.nrows());
1610
1612}
1613
1615// Helper function for debugging
1616#ifndef NDEBUG
1617
1632 return mv(r, c);
1633}
1634
1635#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.
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 & inv(const lapack_help::Inverse< Complex > &help=lapack_help::Inverse< Complex >{0})
~ComplexMatrix() override
Destructor for Matrix.
ComplexMatrix & operator=(ComplexMatrix x)
Assignment operator from another matrix.
ComplexMatrix()=default
friend void swap(ComplexMatrix &m1, ComplexMatrix &m2)
Swaps two objects.
The ComplexVectorView class.
ComplexVectorView operator-=(Complex x)
Subtraction of scalar.
ComplexVectorView operator+=(Complex x)
Addition of scalar.
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
void resize(Index n)
Resize function.
ComplexVector & operator=(const ComplexVector &v)
Assignment from another Vector.
~ComplexVector() override
Destructor for ComplexVector.
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
ConstComplexIterator2D end() const
Return const iterator behind last row.
A constant view of a ComplexVector.
Complex sum() const
Returns true if variable size is zero.
ConstComplexIterator1D end() const
Return const iterator behind last element.
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:449
The const row iterator class for sub matrices.
Definition: matpackI.h:848
A constant view of a Matrix.
Definition: matpackI.h:1043
ConstIterator2D begin() const ARTS_NOEXCEPT
Return const iterator to first row.
Definition: matpackI.cc:448
Index nrows() const noexcept
Definition: matpackI.h:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
A constant view of a Vector.
Definition: matpackI.h:512
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:536
The range class.
Definition: matpackI.h:159
Index mstart
The start index.
Definition: matpackI.h:367
Index mstride
The stride.
Definition: matpackI.h:372
Index mextent
The number of elements.
Definition: matpackI.h:370
The Vector class.
Definition: matpackI.h:899
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
The declarations of all the exception classes.
#define min(a, b)
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.
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.
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
std::complex< Numeric > Complex
Index nelem(const Lines &l)
Number of lines.
constexpr Numeric e
Elementary charge convenience name [C].
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
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 row vector.
Definition: matpack_eigen.h:63
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:76
#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