ARTS 2.5.0 (git: 9ee3ac6c)
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 "matpack_complex.h"
27#include "blas.h"
28#include "exceptions.h"
29#include <cmath>
30#include <cstring>
31
32// Functions for ConstComplexVectorView:
33// ------------------------------
34
36// bool ConstComplexVectorView::empty() const
37// {
38// return (nelem() == 0);
39// }
40
50// Index ConstComplexVectorView::nelem() const
51// {
52// return mrange.mextent;
53// }
54
57 Complex s = 0;
59 const ConstComplexIterator1D e = end();
60
61 for (; i != e; ++i) s += *i;
62
63 return s;
64}
65
70 const Range& r) const {
72}
73
77}
78
84}
85
87ConstComplexVectorView::operator ConstComplexMatrixView() const {
88 return ConstComplexMatrixView(mdata, mrange, Range(0, 1));
89}
90
99 : mrange(0, 1), mdata(&const_cast<Complex&>(a)) {
100 // Nothing to do here.
101}
102
106 const Range& range)
107 : mrange(range), mdata(data) {
108 // Nothing to do here.
109}
110
122 const Range& p,
123 const Range& n)
124 : mrange(p, n), mdata(data) {
125 // Nothing to do here.
126}
127
131std::ostream& operator<<(std::ostream& os, const ConstComplexVectorView& v) {
132 ConstComplexIterator1D i = v.begin();
133 const ConstComplexIterator1D end = v.end();
134
135 if (i != end) {
136 os << std::setw(3) << *i;
137 ++i;
138 }
139 for (; i != end; ++i) {
140 os << " " << std::setw(3) << *i;
141 }
142
143 return os;
144}
145
146// Functions for ComplexVectorView:
147// ------------------------
148
152 ARTS_ASSERT (false,
153 "Creating a ComplexVectorView from a const ComplexVector is not allowed.\n"
154 "This is not really a runtime error, but I don't want to start\n"
155 "producing direct output from inside matpack. And just exiting is\n"
156 "not so nice.\n"
157 "If you see this error, there is a bug in the code, not in the\n"
158 "ARTS input.")
159}
160
163 mdata = v.mdata;
164 mrange = v.mrange;
165}
166
171 return ComplexVectorView(mdata, mrange, r);
172}
173
177}
178
181 return ComplexIterator1D(
184}
185
191 const ConstComplexVectorView& v) {
192 // cout << "Assigning VectorView from ConstVectorView.\n";
193
194 // Check that sizes are compatible:
195 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
196
197 copy(v.begin(), v.end(), begin());
198
199 return *this;
200}
201
208 // cout << "Assigning ComplexVectorView from ComplexVectorView.\n";
209
210 // Check that sizes are compatible:
211 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
212
213 copy(v.begin(), v.end(), begin());
214
215 return *this;
216}
217
221 // cout << "Assigning ComplexVectorView from Vector.\n";
222
223 // Check that sizes are compatible:
224 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
225
226 copy(v.begin(), v.end(), begin());
228 return *this;
229}
230
234 copy(x, begin(), end());
235 return *this;
236}
237
240 const ComplexIterator1D e = end();
241 for (ComplexIterator1D i = begin(); i != e; ++i) *i *= x;
242 return *this;
243}
244
247 const ComplexIterator1D e = end();
248 for (ComplexIterator1D i = begin(); i != e; ++i) *i *= x;
249 return *this;
250}
251
254 const ComplexIterator1D e = end();
255 for (ComplexIterator1D i = begin(); i != e; ++i) *i /= x;
256 return *this;
257}
258
261 const ComplexIterator1D e = end();
262 for (ComplexIterator1D i = begin(); i != e; ++i) *i /= x;
263 return *this;
264}
265
268 const ComplexIterator1D e = end();
269 for (ComplexIterator1D i = begin(); i != e; ++i) *i += x;
270 return *this;
271}
272
275 const ComplexIterator1D e = end();
276 for (ComplexIterator1D i = begin(); i != e; ++i) *i += x;
277 return *this;
278}
279
282 const ComplexIterator1D e = end();
283 for (ComplexIterator1D i = begin(); i != e; ++i) *i -= x;
284 return *this;
285}
286
289 const ComplexIterator1D e = end();
290 for (ComplexIterator1D i = begin(); i != e; ++i) *i -= x;
291 return *this;
292}
293
296 const ConstComplexVectorView& x) {
297 ARTS_ASSERT(nelem() == x.nelem());
298
300
302 const ComplexIterator1D e = end();
303
304 for (; i != e; ++i, ++s) *i *= *s;
305 return *this;
306}
307
310 ARTS_ASSERT(nelem() == x.nelem());
311
312 ConstIterator1D s = x.begin();
313
315 const ComplexIterator1D e = end();
316
317 for (; i != e; ++i, ++s) *i *= *s;
318 return *this;
319}
320
323 const ConstComplexVectorView& x) {
324 ARTS_ASSERT(nelem() == x.nelem());
325
327
329 const ComplexIterator1D e = end();
330
331 for (; i != e; ++i, ++s) *i /= *s;
332 return *this;
333}
337 ARTS_ASSERT(nelem() == x.nelem());
340
342 const ComplexIterator1D e = end();
343
344 for (; i != e; ++i, ++s) *i /= *s;
345 return *this;
346}
347
350 const ConstComplexVectorView& x) {
352
354
358 for (; i != e; ++i, ++s) *i += *s;
359 return *this;
360}
361
364 ARTS_ASSERT(nelem() == x.nelem());
365
366 ConstIterator1D s = x.begin();
367
369 const ComplexIterator1D e = end();
370
371 for (; i != e; ++i, ++s) *i += *s;
372 return *this;
373}
374
377 const ConstComplexVectorView& x) {
378 ARTS_ASSERT(nelem() == x.nelem());
379
381
383 const ComplexIterator1D e = end();
384
385 for (; i != e; ++i, ++s) *i -= *s;
386 return *this;
387}
388
391 ARTS_ASSERT(nelem() == x.nelem());
392
393 ConstIterator1D s = x.begin();
394
396 const ComplexIterator1D e = end();
397
398 for (; i != e; ++i, ++s) *i -= *s;
399 return *this;
400}
401
403ComplexVectorView::operator ComplexMatrixView() {
404 // The old version (before 2013-01-18) of this was:
405 // return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
406 // Bus this was a bug! The problem is that the matrix index operator adds
407 // the mstart from both row and columm range object to mdata
408
409 return ComplexMatrixView(mdata, mrange, Range(0, 1));
410}
411
419 ARTS_ASSERT (not (mrange.mstart != 0 || mrange.mstride != 1),
420 "A ComplexVectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");
421 return mdata;
422}
423
431 ARTS_ASSERT (not (mrange.mstart != 0 || mrange.mstride != 1),
432 "A ComplexVectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");
433 return mdata;
434}
435
440 // Nothing to do here.
441}
442
446 : ConstComplexVectorView(data, range) {
447 // Nothing to do here.
448}
449
461 const Range& p,
462 const Range& n)
464 // Nothing to do here.
465}
466
473 ComplexIterator1D target) {
474 if (origin.mstride == 1 && target.mstride == 1)
475 memcpy((void*)target.mx,
476 (void*)origin.mx,
477 sizeof(Complex) * (end.mx - origin.mx));
478 else
479 for (; origin != end; ++origin, ++target) *target = *origin;
480}
481
484 for (; target != end; ++target) *target = x;
485}
486
487// Functions for Vector:
488// ---------------------
489
492
495 : ComplexVectorView(new Complex[n], Range(0, n)) {
496 // Nothing to do here.
497}
498
501 : ComplexVectorView(new Complex[n], Range(0, n)) {
502 // Here we can access the raw memory directly, for slightly
503 // increased efficiency:
504 const Complex* stop = mdata + n;
505 for (Complex* x = mdata; x < stop; ++x) *x = fill;
506}
507
510 : ComplexVectorView(new Complex[n], Range(0, n)) {
511 // Here we can access the raw memory directly, for slightly
512 // increased efficiency:
513 const Complex* stop = mdata + n;
514 for (Complex* x = mdata; x < stop; ++x) *x = fill;
515}
516
526 : ComplexVectorView(new Complex[extent], Range(0, extent)) {
527 // Fill with values:
528 Complex x = start;
530 const ComplexIterator1D e = end();
531 for (; i != e; ++i) {
532 *i = x;
533 x += stride;
534 }
535}
536
546 : ComplexVectorView(new Complex[extent], Range(0, extent)) {
547 // Fill with values:
548 Complex x = start;
550 const ComplexIterator1D e = end();
551 for (; i != e; ++i) {
552 *i = x;
553 x += stride;
554 }
555}
556
566 : ComplexVectorView(new Complex[extent], Range(0, extent)) {
567 // Fill with values:
568 Complex x = start;
570 const ComplexIterator1D e = end();
571 for (; i != e; ++i) {
572 *i = x;
573 x += stride;
574 }
575}
576
586 : ComplexVectorView(new Complex[extent], Range(0, extent)) {
587 // Fill with values:
588 Complex x = start;
590 const ComplexIterator1D e = end();
591 for (; i != e; ++i) {
592 *i = x;
593 x += stride;
594 }
595}
596
603 : ComplexVectorView(new Complex[v.nelem()], Range(0, v.nelem())) {
604 copy(v.begin(), v.end(), begin());
605}
606
610 : ComplexVectorView(new Complex[v.nelem()], Range(0, v.nelem())) {
611 copy(v.begin(), v.end(), begin());
612}
613
615 for (Index i=0; i<nelem(); i++) operator[](i) = Complex(v[i], 0);
616}
617
619ComplexVector::ComplexVector(const std::vector<Complex>& v)
620 : ComplexVectorView(new Complex[v.size()], Range(0, v.size())) {
621 auto vec_it_end = v.end();
622 ComplexIterator1D this_it = this->begin();
623 for (auto vec_it = v.begin();
624 vec_it != vec_it_end;
625 ++vec_it, ++this_it)
626 *this_it = *vec_it;
627}
628
630ComplexVector::ComplexVector(const std::vector<Numeric>& v)
631 : ComplexVectorView(new Complex[v.size()], Range(0, v.size())) {
632 auto vec_it_end = v.end();
633 ComplexIterator1D this_it = this->begin();
634 for (auto vec_it = v.begin();
635 vec_it != vec_it_end;
636 ++vec_it, ++this_it)
637 *this_it = *vec_it;
638}
639
641
666 swap(*this, v);
667 return *this;
668}
669
671
688 resize(x.nelem());
690 return *this;
691}
692
697 return *this;
700// /** Assignment operator from VectorView. Assignment operators are not
701// inherited. */
702// inline Vector& Vector::operator=(const VectorView v)
703// {
704// cout << "Assigning Vector from Vector View.\n";
705// // Check that sizes are compatible:
706// ARTS_ASSERT(mrange.mextent==v.mrange.mextent);
707// VectorView::operator=(v);
708// return *this;
709// }
715 ARTS_ASSERT(0 <= n);
716 if (mrange.mextent != n) {
717 delete[] mdata;
718 mdata = new Complex[n];
719 mrange.mstart = 0;
721 mrange.mstride = 1;
722 }
723}
724
727 std::swap(v1.mrange, v2.mrange);
728 std::swap(v1.mdata, v2.mdata);
729}
730
734
735// Functions for ConstMatrixView:
736// ------------------------------
737
739// bool ConstComplexMatrixView::empty() const
740// {
741// return (nrows() == 0 || ncols() == 0);
742// }
743
745// Index ConstComplexMatrixView::nrows() const
746// {
747// return mrr.mextent;
748// }
749
751// Index ConstComplexMatrixView::ncols() const
752// {
753// return mcr.mextent;
754// }
755
760 const Range& r, const Range& c) const {
761 return ConstComplexMatrixView(mdata, mrr, mcr, r, c);
762}
763
770 Index c) const {
771 // Check that c is valid:
772 ARTS_ASSERT(0 <= c);
774
776}
777
784 Index r, const Range& c) const {
785 // Check that r is valid:
786 ARTS_ASSERT(0 <= r);
788
790}
791
795 mrr.mstride);
796}
797
802 mcr),
803 mrr.mstride);
804}
805
807
818 Range(0, n, mrr.mstride + mcr.mstride));
819}
820
825 const Range& rr,
826 const Range& cr)
827 : mrr(rr), mcr(cr), mdata(data) {
828 // Nothing to do here.
829}
830
846 const Range& pr,
847 const Range& pc,
848 const Range& nr,
849 const Range& nc)
850 : mrr(pr, nr), mcr(pc, nc), mdata(data) {
851 // Nothing to do here.
852}
853
861std::ostream& operator<<(std::ostream& os, const ConstComplexMatrixView& v) {
862 // Row iterators:
864 const ConstComplexIterator2D end_row = v.end();
865
866 if (ir != end_row) {
867 ConstComplexIterator1D ic = ir->begin();
868 const ConstComplexIterator1D end_col = ir->end();
869
870 if (ic != end_col) {
871 os << std::setw(3) << *ic;
872 ++ic;
873 }
874 for (; ic != end_col; ++ic) {
875 os << " " << std::setw(3) << *ic;
876 }
877 ++ir;
878 }
879 for (; ir != end_row; ++ir) {
880 ConstComplexIterator1D ic = ir->begin();
881 const ConstComplexIterator1D end_col = ir->end();
882
883 os << "\n";
884 if (ic != end_col) {
885 os << std::setw(3) << *ic;
886 ++ic;
887 }
888 for (; ic != end_col; ++ic) {
889 os << " " << std::setw(3) << *ic;
890 }
891 }
892
893 return os;
894}
895
896// Functions for ComplexMatrixView:
897// -------------------------
898
903 const Range& c) {
904 return ComplexMatrixView(mdata, mrr, mcr, r, c);
905}
906
912 // Check that c is valid:
913 ARTS_ASSERT(0 <= c);
915
916 return ComplexVectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
917}
918
924 // Check that r is valid:
925 ARTS_ASSERT(0 <= r);
927
928 return ComplexVectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
929}
930
934 mrr.mstride);
935}
936
939 return ComplexIterator2D(
941 mrr.mstride);
942}
943
945
956 Range(0, n, mrr.mstride + mcr.mstride));
957}
958
964 const ConstComplexMatrixView& m) {
965 // Check that sizes are compatible:
968
969 copy(m.begin(), m.end(), begin());
970 return *this;
971}
972
979 // Check that sizes are compatible:
982
983 copy(m.begin(), m.end(), begin());
984 return *this;
985}
986
991 // Check that sizes are compatible:
994
995 copy(m.begin(), m.end(), begin());
996 return *this;
997}
998
1004 const ConstComplexVectorView& v) {
1005 // Check that sizes are compatible:
1006 ARTS_ASSERT(mrr.mextent == v.nelem());
1007 ARTS_ASSERT(mcr.mextent == 1);
1008 // dummy = ConstComplexMatrixView(v.mdata,v.mrange,Range(v.mrange.mstart,1));;
1010 copy(dummy.begin(), dummy.end(), begin());
1011 return *this;
1012}
1013
1017 copy(x, begin(), end());
1018 return *this;
1019}
1020
1023 const ComplexIterator2D er = end();
1024 for (ComplexIterator2D r = begin(); r != er; ++r) {
1025 const ComplexIterator1D ec = r->end();
1026 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c *= x;
1027 }
1028 return *this;
1029}
1030
1033 const ComplexIterator2D er = end();
1034 for (ComplexIterator2D r = begin(); r != er; ++r) {
1035 const ComplexIterator1D ec = r->end();
1036 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c *= x;
1037 }
1038 return *this;
1039}
1040
1043 const ComplexIterator2D er = end();
1044 for (ComplexIterator2D r = begin(); r != er; ++r) {
1045 const ComplexIterator1D ec = r->end();
1046 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c /= x;
1047 }
1048 return *this;
1049}
1050
1053 const ComplexIterator2D er = end();
1054 for (ComplexIterator2D r = begin(); r != er; ++r) {
1055 const ComplexIterator1D ec = r->end();
1056 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c /= x;
1057 }
1058 return *this;
1059}
1060
1063 const ComplexIterator2D er = end();
1064 for (ComplexIterator2D r = begin(); r != er; ++r) {
1065 const ComplexIterator1D ec = r->end();
1066 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c += x;
1067 }
1068 return *this;
1069}
1070
1073 const ComplexIterator2D er = end();
1074 for (ComplexIterator2D r = begin(); r != er; ++r) {
1075 const ComplexIterator1D ec = r->end();
1076 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c += x;
1077 }
1078 return *this;
1079}
1080
1083 const ComplexIterator2D er = end();
1084 for (ComplexIterator2D r = begin(); r != er; ++r) {
1085 const ComplexIterator1D ec = r->end();
1086 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c -= x;
1087 }
1088 return *this;
1089}
1090
1093 const ComplexIterator2D er = end();
1094 for (ComplexIterator2D r = begin(); r != er; ++r) {
1095 const ComplexIterator1D ec = r->end();
1096 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c -= x;
1097 }
1098 return *this;
1099}
1100
1108 ARTS_ASSERT (not (mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
1109 mcr.mstride != 1),
1110 "A MatrixView can only be converted to a plain C-array if it's pointing to a continuous block of data");
1111 return mdata;
1112}
1113
1121 ARTS_ASSERT (not (mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
1122 mcr.mstride != 1),
1123 "A MatrixView can only be converted to a plain C-array if it's pointing to a continuous block of data");
1124 return mdata;
1125}
1126
1129 const ConstComplexMatrixView& x) {
1130 ARTS_ASSERT(nrows() == x.nrows());
1131 ARTS_ASSERT(ncols() == x.ncols());
1134 const ComplexIterator2D er = end();
1135 for (; r != er; ++r, ++sr) {
1136 ConstComplexIterator1D sc = sr->begin();
1137 ComplexIterator1D c = r->begin();
1138 const ComplexIterator1D ec = r->end();
1139 for (; c != ec; ++c, ++sc) *c *= *sc;
1140 }
1141 return *this;
1142}
1143
1146 ARTS_ASSERT(nrows() == x.nrows());
1147 ARTS_ASSERT(ncols() == x.ncols());
1148 ConstIterator2D sr = x.begin();
1150 const ComplexIterator2D er = end();
1151 for (; r != er; ++r, ++sr) {
1152 ConstIterator1D sc = sr->begin();
1153 ComplexIterator1D c = r->begin();
1154 const ComplexIterator1D ec = r->end();
1155 for (; c != ec; ++c, ++sc) *c *= *sc;
1156 }
1157 return *this;
1158}
1159
1162 const ConstComplexMatrixView& x) {
1163 ARTS_ASSERT(nrows() == x.nrows());
1164 ARTS_ASSERT(ncols() == x.ncols());
1167 const ComplexIterator2D er = end();
1168 for (; r != er; ++r, ++sr) {
1169 ConstComplexIterator1D sc = sr->begin();
1170 ComplexIterator1D c = r->begin();
1171 const ComplexIterator1D ec = r->end();
1172 for (; c != ec; ++c, ++sc) *c /= *sc;
1173 }
1174 return *this;
1175}
1176
1179 ARTS_ASSERT(nrows() == x.nrows());
1180 ARTS_ASSERT(ncols() == x.ncols());
1181 ConstIterator2D sr = x.begin();
1183 const ComplexIterator2D er = end();
1184 for (; r != er; ++r, ++sr) {
1185 ConstIterator1D sc = sr->begin();
1186 ComplexIterator1D c = r->begin();
1187 const ComplexIterator1D ec = r->end();
1188 for (; c != ec; ++c, ++sc) *c /= *sc;
1189 }
1190 return *this;
1191}
1192
1195 const ConstComplexMatrixView& x) {
1196 ARTS_ASSERT(nrows() == x.nrows());
1197 ARTS_ASSERT(ncols() == x.ncols());
1200 const ComplexIterator2D er = end();
1201 for (; r != er; ++r, ++sr) {
1202 ConstComplexIterator1D sc = sr->begin();
1203 ComplexIterator1D c = r->begin();
1204 const ComplexIterator1D ec = r->end();
1205 for (; c != ec; ++c, ++sc) *c += *sc;
1206 }
1207 return *this;
1208}
1209
1212 ARTS_ASSERT(nrows() == x.nrows());
1213 ARTS_ASSERT(ncols() == x.ncols());
1214 ConstIterator2D sr = x.begin();
1216 const ComplexIterator2D er = end();
1217 for (; r != er; ++r, ++sr) {
1218 ConstIterator1D sc = sr->begin();
1219 ComplexIterator1D c = r->begin();
1220 const ComplexIterator1D ec = r->end();
1221 for (; c != ec; ++c, ++sc) *c += *sc;
1222 }
1223 return *this;
1224}
1225
1228 const ConstComplexMatrixView& x) {
1229 ARTS_ASSERT(nrows() == x.nrows());
1230 ARTS_ASSERT(ncols() == x.ncols());
1233 const ComplexIterator2D er = end();
1234 for (; r != er; ++r, ++sr) {
1235 ConstComplexIterator1D sc = sr->begin();
1236 ComplexIterator1D c = r->begin();
1237 const ComplexIterator1D ec = r->end();
1238 for (; c != ec; ++c, ++sc) *c -= *sc;
1239 }
1240 return *this;
1241}
1242
1245 ARTS_ASSERT(nrows() == x.nrows());
1246 ARTS_ASSERT(ncols() == x.ncols());
1247 ConstIterator2D sr = x.begin();
1249 const ComplexIterator2D er = end();
1250 for (; r != er; ++r, ++sr) {
1251 ConstIterator1D sc = sr->begin();
1252 ComplexIterator1D c = r->begin();
1253 const ComplexIterator1D ec = r->end();
1254 for (; c != ec; ++c, ++sc) *c -= *sc;
1255 }
1256 return *this;
1257}
1258
1262 // Nothing to do here.
1263}
1264
1269 const Range& rr,
1270 const Range& cr)
1271 : ConstComplexMatrixView(data, rr, cr) {
1272 // Nothing to do here.
1273}
1274
1290 const Range& pr,
1291 const Range& pc,
1292 const Range& nr,
1293 const Range& nc)
1294 : ConstComplexMatrixView(data, pr, pc, nr, nc) {
1295 // Nothing to do here.
1296}
1297
1309 ComplexIterator2D target) {
1310 for (; origin != end; ++origin, ++target) {
1311 ConstComplexIterator1D o = origin->begin();
1312 const ConstComplexIterator1D e = origin->end();
1313 ComplexIterator1D t = target->begin();
1314 for (; o != e; ++o, ++t) *t = *o;
1315 }
1316}
1317
1320 for (; target != end; ++target) {
1321 ComplexIterator1D t = target->begin();
1322 const ComplexIterator1D e = target->end();
1323 for (; t != e; ++t) *t = x;
1324 }
1325}
1326
1329 for (; target != end; ++target) {
1330 ComplexIterator1D t = target->begin();
1331 const ComplexIterator1D e = target->end();
1332 for (; t != e; ++t) *t = x;
1333 }
1334}
1335
1336// Functions for ComplexMatrix:
1337// ---------------------
1338
1342 : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
1343 // Nothing to do here.
1344}
1345
1348 : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
1349 // Here we can access the raw memory directly, for slightly
1350 // increased efficiency:
1351 const Complex* stop = mdata + r * c;
1352 for (Complex* x = mdata; x < stop; ++x) *x = fill;
1353}
1354
1357 : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
1358 // Here we can access the raw memory directly, for slightly
1359 // increased efficiency:
1360 const Complex* stop = mdata + r * c;
1361 for (Complex* x = mdata; x < stop; ++x) *x = fill;
1362}
1363
1367 : ComplexMatrixView(new Complex[m.nrows() * m.ncols()],
1368 Range(0, m.nrows(), m.ncols()),
1369 Range(0, m.ncols())) {
1370 copy(m.begin(), m.end(), begin());
1371}
1372
1376 : ComplexMatrixView(new Complex[m.nrows() * m.ncols()],
1377 Range(0, m.nrows(), m.ncols()),
1378 Range(0, m.ncols())) {
1379 // There is a catch here: If m is an empty matrix, then it will have
1380 // 0 colunns. But this is used to initialize the stride of the row
1381 // Range! Thus, this method has to be consistent with the behaviour
1382 // of Range::Range. For now, Range::Range allows also stride 0.
1383 copy(m.begin(), m.end(), begin());
1384}
1385
1387
1412 swap(*this, m);
1413 return *this;
1414}
1415
1419 copy(x, begin(), end());
1420 return *this;
1421}
1422
1424
1437 resize(v.nelem(), 1);
1439 copy(dummy.begin(), dummy.end(), begin());
1440 return *this;
1441}
1442
1447 ARTS_ASSERT(0 <= r);
1448 ARTS_ASSERT(0 <= c);
1449
1450 if (mrr.mextent != r || mcr.mextent != c) {
1451 delete[] mdata;
1452 mdata = new Complex[r * c];
1453
1454 mrr.mstart = 0;
1455 mrr.mextent = r;
1456 mrr.mstride = c;
1457
1458 mcr.mstart = 0;
1459 mcr.mextent = c;
1460 mcr.mstride = 1;
1461 }
1462}
1463
1466 std::swap(m1.mrr, m2.mrr);
1467 std::swap(m1.mcr, m2.mcr);
1468 std::swap(m1.mdata, m2.mdata);
1469}
1470
1474 // cout << "Destroying a Matrix:\n"
1475 // << *this << "\n........................................\n";
1476 delete[] mdata;
1477}
1478
1480{
1481 ARTS_ASSERT(ncols() == nrows());
1482
1483 // help's internal variables are all mutable, so const is just for default parameter and to not use copies
1484 help.resize_if_smaller(int(ncols()));
1485
1486 // Compute LU decomposition using LAPACK dgetrf_.
1487 int info;
1488 lapack::zgetrf_(help.size(), help.size(), mdata, help.size(), help.ipivdata(), &info);
1489 lapack::zgetri_(help.size(), mdata, help.size(), help.ipivdata(), help.workdata(), help.lsize(), &info);
1490
1491 // Check for success.
1492 ARTS_USER_ERROR_IF (info not_eq 0,
1493 "Error inverting matrix: Matrix not of full rank.");
1494
1495 return *this;
1496}
1497
1500 return ConstComplexMatrixView(m.mdata, m.mcr, m.mrr);
1501}
1502
1506 return ComplexMatrixView(m.mdata, m.mcr, m.mrr);
1507}
1508
1512 return transpose((ComplexMatrixView)v);
1513}
1514
1525 // cout << "Assigning VectorView from Array<Numeric>.\n";
1526
1527 // Check that sizes are compatible:
1528 ARTS_ASSERT(mrange.mextent == v.nelem());
1529
1530 // Iterators for Array:
1531 auto i = v.begin();
1532 const auto e = v.end();
1533 // Iterator for Vector:
1534 ComplexIterator1D target = begin();
1535
1536 for (; i != e; ++i, ++target) *target = *i;
1537
1538 return *this;
1539}
1540
1541//Functions operating on the complex vectors
1542
1545 const ConstComplexVectorView& b) {
1546 // Check dimensions:
1547 ARTS_ASSERT(a.nelem() == b.nelem());
1548
1549 const ConstComplexIterator1D ae = a.end();
1550 ConstComplexIterator1D ai = a.begin();
1551 ConstComplexIterator1D bi = b.begin();
1552
1553 Complex res = 0;
1554 for (; ai != ae; ++ai, ++bi) res += (*ai) * (*bi);
1555
1556 return res;
1557}
1558
1560
1571 const ConstComplexVectorView& x) {
1572 ARTS_ASSERT(x.nelem() == M.nrows());
1573 ARTS_ASSERT(y.nelem() == M.ncols());
1574
1576 if (y.mdata == x.mdata)
1577 eigen_y = MapToEigen(M) * MapToEigenRow(x);
1578 else
1579 eigen_y.noalias() = MapToEigen(M) * MapToEigenRow(x);
1580}
1581
1583
1595 const ConstComplexMatrixView& B,
1596 const ConstComplexMatrixView& C) {
1597 ARTS_ASSERT(B.nrows() == A.nrows());
1598 ARTS_ASSERT(C.ncols() == A.ncols());
1599 ARTS_ASSERT(B.ncols() == C.nrows());
1600
1601 ComplexMatrixViewMap eigen_A = MapToEigen(A);
1602 if (A.mdata == B.mdata || A.mdata == C.mdata)
1603 eigen_A = MapToEigen(B) * MapToEigen(C);
1604 else
1605 eigen_A.noalias() = MapToEigen(B) * MapToEigen(C);
1606}
1608 const ConstComplexMatrixView& B,
1609 const ConstMatrixView& C) {
1610 ARTS_ASSERT(B.nrows() == A.nrows());
1611 ARTS_ASSERT(C.ncols() == A.ncols());
1612 ARTS_ASSERT(B.ncols() == C.nrows());
1613
1614 ComplexMatrixViewMap eigen_A = MapToEigen(A);
1615 if (A.mdata == B.mdata)
1616 eigen_A = MapToEigen(B) * MapToEigen(C);
1617 else
1618 eigen_A.noalias() = MapToEigen(B) * MapToEigen(C);
1619}
1621 const ConstMatrixView& B,
1622 const ConstComplexMatrixView& C) {
1623 ARTS_ASSERT(B.nrows() == A.nrows());
1624 ARTS_ASSERT(C.ncols() == A.ncols());
1625 ARTS_ASSERT(B.ncols() == C.nrows());
1626
1627 ComplexMatrixViewMap eigen_A = MapToEigen(A);
1628 if (A.mdata == C.mdata)
1629 eigen_A = MapToEigen(B) * MapToEigen(C);
1630 else
1631 eigen_A.noalias() = MapToEigen(B) * MapToEigen(C);
1632}
1634 const ConstMatrixView& B,
1635 const ConstMatrixView& C) {
1636 ARTS_ASSERT(B.nrows() == A.nrows());
1637 ARTS_ASSERT(C.ncols() == A.ncols());
1638 ARTS_ASSERT(B.ncols() == C.nrows());
1639
1640 ComplexMatrixViewMap eigen_A = MapToEigen(A);
1641 eigen_A.noalias() = MapToEigen(B) * MapToEigen(C);
1642}
1643
1644// Converts constant matrix to constant eigen map
1647 A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1648 A.nrows(),
1649 A.ncols(),
1651}
1652
1653// Converts constant vector to constant eigen row-view
1656 A.nelem(),
1657 1,
1658 StrideType(A.mrange.get_stride(), 1));
1659}
1660
1661// Converts constant vector to constant eigen row-view
1663 return MapToEigen(A);
1664}
1665
1666// Converts constant vector to constant eigen column-view
1669 1,
1670 A.nelem(),
1671 StrideType(1, A.mrange.get_stride()));
1672}
1673
1674// Converts matrix to eigen map
1676 return ComplexMatrixViewMap(
1677 A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1678 A.nrows(),
1679 A.ncols(),
1681}
1682
1683// Converts vector to eigen map row-view
1686 A.nelem(),
1687 1,
1688 StrideType(A.mrange.get_stride(), 1));
1689}
1690
1691// Converts vector to eigen map row-view
1693 return MapToEigen(A);
1694}
1695
1696// Converts vector to eigen map column-view
1699 1,
1700 A.nelem(),
1701 StrideType(1, A.mrange.get_stride()));
1702}
1703
1705// Helper function for debugging
1706#ifndef NDEBUG
1707
1722 return mv(r, c);
1723}
1724
1725#endif
Index nrows
void * data
Index ncols
Interface for BLAS library.
This can be used to make arrays out of anything.
Definition: array.h:107
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
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.
const Complex * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array.
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.
const Complex * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array.
ComplexIterator1D begin()
Return iterator to first element.
The ComplexVector class.
ComplexVector & operator=(ComplexVector v)
Assignment from another Vector.
ComplexVector()
Default constructor.
void resize(Index n)
Resize function.
~ComplexVector() override
Destructor for ComplexVector.
friend void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
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:425
The const row iterator class for sub matrices.
Definition: matpackI.h:825
A constant view of a Matrix.
Definition: matpackI.h:1014
ConstIterator2D begin() const ARTS_NOEXCEPT
Return const iterator to first row.
Definition: matpackI.cc:473
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
A constant view of a Vector.
Definition: matpackI.h:489
ConstIterator1D begin() const ARTS_NOEXCEPT
Return const iterator to first element.
Definition: matpackI.cc:66
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The range class.
Definition: matpackI.h:165
constexpr Index get_stride() const ARTS_NOEXCEPT
Returns the stride of the range.
Definition: matpackI.h:336
Index mstart
The start index.
Definition: matpackI.h:351
constexpr Index get_start() const ARTS_NOEXCEPT
Returns the start index of the range.
Definition: matpackI.h:332
Index mstride
The stride.
Definition: matpackI.h:360
Index mextent
The number of elements.
Definition: matpackI.h:358
The Vector class.
Definition: matpackI.h:876
#define ARTS_NOEXCEPT
Definition: debug.h:80
#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)
Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > StrideType
Definition: matpackI.h:109
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
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.
ConstComplexMatrixViewMap MapToEigenRow(const ConstComplexVectorView &A)
ConstComplexMatrixViewMap MapToEigen(const ConstComplexMatrixView &A)
void copy(ConstComplexIterator1D origin, const ConstComplexIterator1D &end, ComplexIterator1D target)
Copy data between begin and end to target.
ConstComplexMatrixViewMap MapToEigenCol(const ConstComplexVectorView &A)
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
Eigen::Map< const ComplexMatrixType, 0, StrideType > ConstComplexMatrixViewMap
Eigen::Map< ComplexMatrixType, 0, StrideType > ComplexMatrixViewMap
Index nelem(const Lines &l)
Number of lines.
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
Definition: zeemandata.h:78
constexpr Rational end(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the largest M for a polarization type of this transition.
Definition: zeemandata.h:109
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)
#define v
#define a
#define c
#define b
#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