ARTS 2.5.4 (git: 4c0d3b4d)
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 <algorithm>
30#include <cmath>
31#include <cstring>
32
33// Functions for ConstComplexVectorView:
34// ------------------------------
35
37// bool ConstComplexVectorView::empty() const
38// {
39// return (nelem() == 0);
40// }
41
51// Index ConstComplexVectorView::nelem() const
52// {
53// return mrange.mextent;
54// }
55
58 Complex s = 0;
60 const ConstComplexIterator1D e = end();
61
62 for (; i != e; ++i) s += *i;
63
64 return s;
65}
66
71 const Range& r) const {
73}
74
78}
79
85}
86
88ConstComplexVectorView::operator ConstComplexMatrixView() const {
89 return ConstComplexMatrixView(mdata, mrange, Range(0, 1));
90}
91
100 : mrange(0, 1), mdata(&const_cast<Complex&>(a)) {
101 // Nothing to do here.
102}
103
107 const Range& range)
108 : mrange(range), mdata(data) {
109 // Nothing to do here.
110}
111
123 const Range& p,
124 const Range& n)
125 : mrange(p, n), mdata(data) {
126 // Nothing to do here.
127}
128
132std::ostream& operator<<(std::ostream& os, const ConstComplexVectorView& v) {
133 ConstComplexIterator1D i = v.begin();
134 const ConstComplexIterator1D end = v.end();
135
136 if (i != end) {
137 os << std::setw(3) << *i;
138 ++i;
139 }
140 for (; i != end; ++i) {
141 os << " " << std::setw(3) << *i;
142 }
143
144 return os;
145}
146
147// Functions for ComplexVectorView:
148// ------------------------
149
153 ARTS_ASSERT (false,
154 "Creating a ComplexVectorView from a const ComplexVector is not allowed.\n"
155 "This is not really a runtime error, but I don't want to start\n"
156 "producing direct output from inside matpack. And just exiting is\n"
157 "not so nice.\n"
158 "If you see this error, there is a bug in the code, not in the\n"
159 "ARTS input.")
160}
161
164 mdata = v.mdata;
165 mrange = v.mrange;
166}
167
172 return ComplexVectorView(mdata, mrange, r);
173}
174
178}
179
182 return ComplexIterator1D(
185}
186
192 const ConstComplexVectorView& v) {
193 // cout << "Assigning VectorView from ConstVectorView.\n";
194
195 // Check that sizes are compatible:
196 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
197
198 copy(v.begin(), v.end(), begin());
199
200 return *this;
201}
202
209 // cout << "Assigning ComplexVectorView from ComplexVectorView.\n";
210
211 // Check that sizes are compatible:
212 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
213
214 copy(v.begin(), v.end(), begin());
215
216 return *this;
217}
218
222 // cout << "Assigning ComplexVectorView from Vector.\n";
223
224 // Check that sizes are compatible:
225 ARTS_ASSERT(mrange.mextent == v.mrange.mextent);
226
227 copy(v.begin(), v.end(), begin());
228
229 return *this;
230}
231
235 copy(x, begin(), end());
236 return *this;
237}
238
242 copy(x, begin(), end());
243 return *this;
244}
245
248 const ComplexIterator1D e = end();
249 for (ComplexIterator1D i = begin(); i != e; ++i) *i *= x;
250 return *this;
251}
252
255 const ComplexIterator1D e = end();
256 for (ComplexIterator1D i = begin(); i != e; ++i) *i *= x;
257 return *this;
258}
259
262 const ComplexIterator1D e = end();
263 for (ComplexIterator1D i = begin(); i != e; ++i) *i /= x;
264 return *this;
265}
266
269 const ComplexIterator1D e = end();
270 for (ComplexIterator1D i = begin(); i != e; ++i) *i /= x;
271 return *this;
272}
273
276 const ComplexIterator1D e = end();
277 for (ComplexIterator1D i = begin(); i != e; ++i) *i += x;
278 return *this;
279}
280
283 const ComplexIterator1D e = end();
284 for (ComplexIterator1D i = begin(); i != e; ++i) *i += x;
285 return *this;
286}
287
290 const ComplexIterator1D e = end();
291 for (ComplexIterator1D i = begin(); i != e; ++i) *i -= x;
292 return *this;
293}
294
297 const ComplexIterator1D e = end();
298 for (ComplexIterator1D i = begin(); i != e; ++i) *i -= x;
299 return *this;
300}
301
304 const ConstComplexVectorView& x) {
305 ARTS_ASSERT(nelem() == x.nelem());
306
308
310 const ComplexIterator1D e = end();
311
312 for (; i != e; ++i, ++s) *i *= *s;
313 return *this;
314}
315
318 ARTS_ASSERT(nelem() == x.nelem());
319
320 ConstIterator1D s = x.begin();
321
323 const ComplexIterator1D e = end();
324
325 for (; i != e; ++i, ++s) *i *= *s;
326 return *this;
327}
328
331 const ConstComplexVectorView& x) {
332 ARTS_ASSERT(nelem() == x.nelem());
333
335
337 const ComplexIterator1D e = end();
339 for (; i != e; ++i, ++s) *i /= *s;
340 return *this;
341}
342
345 ARTS_ASSERT(nelem() == x.nelem());
346
347 ConstIterator1D s = x.begin();
348
350 const ComplexIterator1D e = end();
352 for (; i != e; ++i, ++s) *i /= *s;
353 return *this;
354}
359 ARTS_ASSERT(nelem() == x.nelem());
360
362
364 const ComplexIterator1D e = end();
365
366 for (; i != e; ++i, ++s) *i += *s;
367 return *this;
368}
369
372 ARTS_ASSERT(nelem() == x.nelem());
373
374 ConstIterator1D s = x.begin();
375
377 const ComplexIterator1D e = end();
378
379 for (; i != e; ++i, ++s) *i += *s;
380 return *this;
381}
382
385 const ConstComplexVectorView& x) {
386 ARTS_ASSERT(nelem() == x.nelem());
387
389
391 const ComplexIterator1D e = end();
392
393 for (; i != e; ++i, ++s) *i -= *s;
394 return *this;
395}
396
399 ARTS_ASSERT(nelem() == x.nelem());
400
401 ConstIterator1D s = x.begin();
402
404 const ComplexIterator1D e = end();
405
406 for (; i != e; ++i, ++s) *i -= *s;
407 return *this;
408}
409
411ComplexVectorView::operator ComplexMatrixView() {
412 // The old version (before 2013-01-18) of this was:
413 // return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
414 // Bus this was a bug! The problem is that the matrix index operator adds
415 // the mstart from both row and columm range object to mdata
416
417 return ComplexMatrixView(mdata, mrange, Range(0, 1));
418}
419
427 ARTS_ASSERT (not (mrange.mstart != 0 || mrange.mstride != 1),
428 "A ComplexVectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");
429 return mdata;
430}
431
439 ARTS_ASSERT (not (mrange.mstart != 0 || mrange.mstride != 1),
440 "A ComplexVectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");
441 return mdata;
442}
443
448 // Nothing to do here.
449}
450
454 : ConstComplexVectorView(data, range) {
455 // Nothing to do here.
456}
457
469 const Range& p,
470 const Range& n)
471 : ConstComplexVectorView(data, p, n) {
472 // Nothing to do here.
473}
474
481 ComplexIterator1D target) {
482 if (origin.mstride == 1 && target.mstride == 1)
483 memcpy((void*)target.mx,
484 (void*)origin.mx,
485 sizeof(Complex) * (end.mx - origin.mx));
486 else
487 for (; origin != end; ++origin, ++target) *target = *origin;
488}
489
492 for (; target != end; ++target) *target = x;
493}
494
495// Functions for Vector:
496// ---------------------
497
500 : ComplexVectorView(new Complex[n], Range(0, n)) {
501 // Nothing to do here.
502}
503
506 : ComplexVectorView(new Complex[n], Range(0, n)) {
507 // Here we can access the raw memory directly, for slightly
508 // increased efficiency:
509 const Complex* stop = mdata + n;
510 for (Complex* x = mdata; x < stop; ++x) *x = fill;
511}
512
515 : ComplexVectorView(new Complex[n], Range(0, n)) {
516 // Here we can access the raw memory directly, for slightly
517 // increased efficiency:
518 const Complex* stop = mdata + n;
519 for (Complex* x = mdata; x < stop; ++x) *x = fill;
520}
521
531 : ComplexVectorView(new Complex[extent], Range(0, extent)) {
532 // Fill with values:
533 Complex x = start;
535 const ComplexIterator1D e = end();
536 for (; i != e; ++i) {
537 *i = x;
538 x += stride;
539 }
540}
541
551 : ComplexVectorView(new Complex[extent], Range(0, extent)) {
552 // Fill with values:
553 Complex x = start;
555 const ComplexIterator1D e = end();
556 for (; i != e; ++i) {
557 *i = x;
558 x += stride;
559 }
560}
561
571 : ComplexVectorView(new Complex[extent], Range(0, extent)) {
572 // Fill with values:
573 Complex x = start;
575 const ComplexIterator1D e = end();
576 for (; i != e; ++i) {
577 *i = x;
578 x += stride;
579 }
580}
581
591 : ComplexVectorView(new Complex[extent], Range(0, extent)) {
592 // Fill with values:
593 Complex x = start;
595 const ComplexIterator1D e = end();
596 for (; i != e; ++i) {
597 *i = x;
598 x += stride;
599 }
600}
601
608 : ComplexVectorView(new Complex[v.nelem()], Range(0, v.nelem())) {
609 copy(v.begin(), v.end(), begin());
610}
611
615 : ComplexVectorView(new Complex[v.nelem()], Range(0, v.nelem())) {
616 copy(v.begin(), v.end(), begin());
617}
618
620 for (Index i=0; i<nelem(); i++) operator[](i) = Complex(v[i], 0);
621}
622
624ComplexVector::ComplexVector(const std::vector<Complex>& v)
625 : ComplexVectorView(new Complex[v.size()], Range(0, v.size())) {
626 auto vec_it_end = v.end();
627 ComplexIterator1D this_it = this->begin();
628 for (auto vec_it = v.begin();
629 vec_it != vec_it_end;
630 ++vec_it, ++this_it)
631 *this_it = *vec_it;
632}
633
635ComplexVector::ComplexVector(const std::vector<Numeric>& v)
636 : ComplexVectorView(new Complex[v.size()], Range(0, v.size())) {
637 auto vec_it_end = v.end();
638 ComplexIterator1D this_it = this->begin();
639 for (auto vec_it = v.begin();
640 vec_it != vec_it_end;
641 ++vec_it, ++this_it)
642 *this_it = *vec_it;
643}
644
646
671 const auto n = v.nelem();
672 resize(n);
673
674 auto* in = v.mdata + v.mrange.mstart;
675 auto* end = v.mdata + n;
676 auto* out = mdata;
677
678 while (in < end) {
679 out = in;
680 out++;
681 in += v.mrange.mstride;
682 }
683
684 return *this;
685}
686
688
705 resize(x.nelem());
707 return *this;
708}
709
714 return *this;
715}
717// /** Assignment operator from VectorView. Assignment operators are not
718// inherited. */
719// inline Vector& Vector::operator=(const VectorView v)
720// {
721// cout << "Assigning Vector from Vector View.\n";
722// // Check that sizes are compatible:
723// ARTS_ASSERT(mrange.mextent==v.mrange.mextent);
724// VectorView::operator=(v);
725// return *this;
726// }
727
732 ARTS_ASSERT(0 <= n);
733 if (mrange.mextent != n) {
734 delete[] mdata;
735 mdata = new Complex[n];
736 mrange.mstart = 0;
737 mrange.mextent = n;
738 mrange.mstride = 1;
741
744 std::swap(v1.mrange, v2.mrange);
745 std::swap(v1.mdata, v2.mdata);
746}
747
751
752// Functions for ConstMatrixView:
753// ------------------------------
754
756// bool ConstComplexMatrixView::empty() const
757// {
758// return (nrows() == 0 || ncols() == 0);
759// }
760
762// Index ConstComplexMatrixView::nrows() const
763// {
764// return mrr.mextent;
765// }
766
768// Index ConstComplexMatrixView::ncols() const
769// {
770// return mcr.mextent;
771// }
772
777 const Range& r, const Range& c) const {
778 return ConstComplexMatrixView(mdata, mrr, mcr, r, c);
779}
780
787 Index c) const {
788 // Check that c is valid:
789 ARTS_ASSERT(0 <= c);
791
793}
794
801 Index r, const Range& c) const {
802 // Check that r is valid:
803 ARTS_ASSERT(0 <= r);
805
807}
808
812 mrr.mstride);
813}
814
819 mcr),
820 mrr.mstride);
821}
822
824
835 Range(0, n, mrr.mstride + mcr.mstride));
836}
837
842 const Range& rr,
843 const Range& cr)
844 : mrr(rr), mcr(cr), mdata(data) {
845 // Nothing to do here.
846}
847
863 const Range& pr,
864 const Range& pc,
865 const Range& nr,
866 const Range& nc)
867 : mrr(pr, nr), mcr(pc, nc), mdata(data) {
868 // Nothing to do here.
869}
870
878std::ostream& operator<<(std::ostream& os, const ConstComplexMatrixView& v) {
879 // Row iterators:
880 ConstComplexIterator2D ir = v.begin();
881 const ConstComplexIterator2D end_row = v.end();
882
883 if (ir != end_row) {
884 ConstComplexIterator1D ic = ir->begin();
885 const ConstComplexIterator1D end_col = ir->end();
886
887 if (ic != end_col) {
888 os << std::setw(3) << *ic;
889 ++ic;
890 }
891 for (; ic != end_col; ++ic) {
892 os << " " << std::setw(3) << *ic;
893 }
894 ++ir;
895 }
896 for (; ir != end_row; ++ir) {
897 ConstComplexIterator1D ic = ir->begin();
898 const ConstComplexIterator1D end_col = ir->end();
899
900 os << "\n";
901 if (ic != end_col) {
902 os << std::setw(3) << *ic;
903 ++ic;
904 }
905 for (; ic != end_col; ++ic) {
906 os << " " << std::setw(3) << *ic;
907 }
908 }
909
910 return os;
911}
912
913// Functions for ComplexMatrixView:
914// -------------------------
915
920 const Range& c) {
921 return ComplexMatrixView(mdata, mrr, mcr, r, c);
922}
923
929 // Check that c is valid:
930 ARTS_ASSERT(0 <= c);
932
933 return ComplexVectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
934}
935
941 // Check that r is valid:
942 ARTS_ASSERT(0 <= r);
944
945 return ComplexVectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
946}
947
951 mrr.mstride);
952}
953
956 return ComplexIterator2D(
958 mrr.mstride);
959}
960
962
973 Range(0, n, mrr.mstride + mcr.mstride));
974}
975
981 const ConstComplexMatrixView& m) {
982 // Check that sizes are compatible:
985
986 copy(m.begin(), m.end(), begin());
987 return *this;
988}
989
996 // Check that sizes are compatible:
999
1000 copy(m.begin(), m.end(), begin());
1001 return *this;
1002}
1003
1008 // Check that sizes are compatible:
1011
1012 copy(m.begin(), m.end(), begin());
1013 return *this;
1014}
1015
1021 const ConstComplexVectorView& v) {
1022 // Check that sizes are compatible:
1023 ARTS_ASSERT(mrr.mextent == v.nelem());
1024 ARTS_ASSERT(mcr.mextent == 1);
1025 // dummy = ConstComplexMatrixView(v.mdata,v.mrange,Range(v.mrange.mstart,1));;
1027 copy(dummy.begin(), dummy.end(), begin());
1028 return *this;
1029}
1030
1034 copy(x, begin(), end());
1035 return *this;
1036}
1037
1040 const ComplexIterator2D er = end();
1041 for (ComplexIterator2D r = begin(); r != er; ++r) {
1042 const ComplexIterator1D ec = r->end();
1043 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c *= x;
1044 }
1045 return *this;
1046}
1047
1050 const ComplexIterator2D er = end();
1051 for (ComplexIterator2D r = begin(); r != er; ++r) {
1052 const ComplexIterator1D ec = r->end();
1053 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c *= x;
1054 }
1055 return *this;
1056}
1057
1060 const ComplexIterator2D er = end();
1061 for (ComplexIterator2D r = begin(); r != er; ++r) {
1062 const ComplexIterator1D ec = r->end();
1063 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c /= x;
1064 }
1065 return *this;
1066}
1067
1070 const ComplexIterator2D er = end();
1071 for (ComplexIterator2D r = begin(); r != er; ++r) {
1072 const ComplexIterator1D ec = r->end();
1073 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c /= x;
1074 }
1075 return *this;
1076}
1077
1080 const ComplexIterator2D er = end();
1081 for (ComplexIterator2D r = begin(); r != er; ++r) {
1082 const ComplexIterator1D ec = r->end();
1083 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c += x;
1084 }
1085 return *this;
1086}
1087
1090 const ComplexIterator2D er = end();
1091 for (ComplexIterator2D r = begin(); r != er; ++r) {
1092 const ComplexIterator1D ec = r->end();
1093 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c += x;
1094 }
1095 return *this;
1096}
1097
1100 const ComplexIterator2D er = end();
1101 for (ComplexIterator2D r = begin(); r != er; ++r) {
1102 const ComplexIterator1D ec = r->end();
1103 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c -= x;
1104 }
1105 return *this;
1106}
1107
1110 const ComplexIterator2D er = end();
1111 for (ComplexIterator2D r = begin(); r != er; ++r) {
1112 const ComplexIterator1D ec = r->end();
1113 for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c -= x;
1114 }
1115 return *this;
1116}
1117
1125 ARTS_ASSERT (not (mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
1126 mcr.mstride != 1),
1127 "A MatrixView can only be converted to a plain C-array if it's pointing to a continuous block of data");
1128 return mdata;
1129}
1130
1138 ARTS_ASSERT (not (mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
1139 mcr.mstride != 1),
1140 "A MatrixView can only be converted to a plain C-array if it's pointing to a continuous block of data");
1141 return mdata;
1142}
1143
1146 const ConstComplexMatrixView& x) {
1147 ARTS_ASSERT(nrows() == x.nrows());
1148 ARTS_ASSERT(ncols() == x.ncols());
1151 const ComplexIterator2D er = end();
1152 for (; r != er; ++r, ++sr) {
1153 ConstComplexIterator1D 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 ARTS_ASSERT(nrows() == x.nrows());
1164 ARTS_ASSERT(ncols() == x.ncols());
1165 ConstIterator2D sr = x.begin();
1167 const ComplexIterator2D er = end();
1168 for (; r != er; ++r, ++sr) {
1169 ConstIterator1D 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 const ConstComplexMatrixView& x) {
1180 ARTS_ASSERT(nrows() == x.nrows());
1181 ARTS_ASSERT(ncols() == x.ncols());
1184 const ComplexIterator2D er = end();
1185 for (; r != er; ++r, ++sr) {
1186 ConstComplexIterator1D 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 ARTS_ASSERT(nrows() == x.nrows());
1197 ARTS_ASSERT(ncols() == x.ncols());
1198 ConstIterator2D sr = x.begin();
1200 const ComplexIterator2D er = end();
1201 for (; r != er; ++r, ++sr) {
1202 ConstIterator1D 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 const ConstComplexMatrixView& x) {
1213 ARTS_ASSERT(nrows() == x.nrows());
1214 ARTS_ASSERT(ncols() == x.ncols());
1217 const ComplexIterator2D er = end();
1218 for (; r != er; ++r, ++sr) {
1219 ConstComplexIterator1D 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
1229 ARTS_ASSERT(nrows() == x.nrows());
1230 ARTS_ASSERT(ncols() == x.ncols());
1231 ConstIterator2D sr = x.begin();
1233 const ComplexIterator2D er = end();
1234 for (; r != er; ++r, ++sr) {
1235 ConstIterator1D 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 const ConstComplexMatrixView& x) {
1246 ARTS_ASSERT(nrows() == x.nrows());
1247 ARTS_ASSERT(ncols() == x.ncols());
1250 const ComplexIterator2D er = end();
1251 for (; r != er; ++r, ++sr) {
1252 ConstComplexIterator1D sc = sr->begin();
1253 ComplexIterator1D c = r->begin();
1254 const ComplexIterator1D ec = r->end();
1255 for (; c != ec; ++c, ++sc) *c -= *sc;
1256 }
1257 return *this;
1258}
1259
1262 ARTS_ASSERT(nrows() == x.nrows());
1263 ARTS_ASSERT(ncols() == x.ncols());
1264 ConstIterator2D sr = x.begin();
1266 const ComplexIterator2D er = end();
1267 for (; r != er; ++r, ++sr) {
1268 ConstIterator1D sc = sr->begin();
1269 ComplexIterator1D c = r->begin();
1270 const ComplexIterator1D ec = r->end();
1271 for (; c != ec; ++c, ++sc) *c -= *sc;
1272 }
1273 return *this;
1274}
1275
1279 // Nothing to do here.
1280}
1281
1286 const Range& rr,
1287 const Range& cr)
1288 : ConstComplexMatrixView(data, rr, cr) {
1289 // Nothing to do here.
1290}
1291
1307 const Range& pr,
1308 const Range& pc,
1309 const Range& nr,
1310 const Range& nc)
1311 : ConstComplexMatrixView(data, pr, pc, nr, nc) {
1312 // Nothing to do here.
1313}
1314
1326 ComplexIterator2D target) {
1327 for (; origin != end; ++origin, ++target) {
1328 ConstComplexIterator1D o = origin->begin();
1329 const ConstComplexIterator1D e = origin->end();
1330 ComplexIterator1D t = target->begin();
1331 for (; o != e; ++o, ++t) *t = *o;
1332 }
1333}
1334
1337 for (; target != end; ++target) {
1338 ComplexIterator1D t = target->begin();
1339 const ComplexIterator1D e = target->end();
1340 for (; t != e; ++t) *t = x;
1341 }
1342}
1343
1346 for (; target != end; ++target) {
1347 ComplexIterator1D t = target->begin();
1348 const ComplexIterator1D e = target->end();
1349 for (; t != e; ++t) *t = x;
1350 }
1351}
1352
1353// Functions for ComplexMatrix:
1354// ---------------------
1355
1359 : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
1360 // Nothing to do here.
1361}
1362
1365 : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
1366 // Here we can access the raw memory directly, for slightly
1367 // increased efficiency:
1368 const Complex* stop = mdata + r * c;
1369 for (Complex* x = mdata; x < stop; ++x) *x = fill;
1370}
1371
1374 : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
1375 // Here we can access the raw memory directly, for slightly
1376 // increased efficiency:
1377 const Complex* stop = mdata + r * c;
1378 for (Complex* x = mdata; x < stop; ++x) *x = fill;
1379}
1380
1384 : ComplexMatrixView(new Complex[m.nrows() * m.ncols()],
1385 Range(0, m.nrows(), m.ncols()),
1386 Range(0, m.ncols())) {
1387 copy(m.begin(), m.end(), begin());
1388}
1389
1393 : ComplexMatrixView(new Complex[m.nrows() * m.ncols()],
1394 Range(0, m.nrows(), m.ncols()),
1395 Range(0, m.ncols())) {
1396 // There is a catch here: If m is an empty matrix, then it will have
1397 // 0 colunns. But this is used to initialize the stride of the row
1398 // Range! Thus, this method has to be consistent with the behaviour
1399 // of Range::Range. For now, Range::Range allows also stride 0.
1400 copy(m.begin(), m.end(), begin());
1401}
1402
1404
1429 swap(*this, m);
1430 return *this;
1431}
1432
1436 copy(x, begin(), end());
1437 return *this;
1438}
1439
1441
1454 resize(v.nelem(), 1);
1456 copy(dummy.begin(), dummy.end(), begin());
1457 return *this;
1458}
1459
1464 ARTS_ASSERT(0 <= r);
1465 ARTS_ASSERT(0 <= c);
1466
1467 if (mrr.mextent != r || mcr.mextent != c) {
1468 delete[] mdata;
1469 mdata = new Complex[r * c];
1470
1471 mrr.mstart = 0;
1472 mrr.mextent = r;
1473 mrr.mstride = c;
1474
1475 mcr.mstart = 0;
1476 mcr.mextent = c;
1477 mcr.mstride = 1;
1478 }
1479}
1480
1483 std::swap(m1.mrr, m2.mrr);
1484 std::swap(m1.mcr, m2.mcr);
1485 std::swap(m1.mdata, m2.mdata);
1486}
1487
1491 // cout << "Destroying a Matrix:\n"
1492 // << *this << "\n........................................\n";
1493 delete[] mdata;
1494}
1495
1497{
1498 ARTS_ASSERT(ncols() == nrows());
1499
1500 // help's internal variables are all mutable, so const is just for default parameter and to not use copies
1501 help.resize_if_smaller(int(ncols()));
1502
1503 // Compute LU decomposition using LAPACK dgetrf_.
1504 int info;
1505 lapack::zgetrf_(help.size(), help.size(), mdata, help.size(), help.ipivdata(), &info);
1506 lapack::zgetri_(help.size(), mdata, help.size(), help.ipivdata(), help.workdata(), help.lsize(), &info);
1507
1508 // Check for success.
1509 ARTS_USER_ERROR_IF (info not_eq 0,
1510 "Error inverting matrix: Matrix not of full rank.");
1511
1512 return *this;
1513}
1514
1517 return ConstComplexMatrixView(m.mdata, m.mcr, m.mrr);
1518}
1519
1523 return ComplexMatrixView(m.mdata, m.mcr, m.mrr);
1524}
1525
1529 return transpose((ComplexMatrixView)v);
1530}
1531
1542 // cout << "Assigning VectorView from Array<Numeric>.\n";
1543
1544 // Check that sizes are compatible:
1545 ARTS_ASSERT(mrange.mextent == v.nelem());
1546
1547 // Iterators for Array:
1548 auto i = v.begin();
1549 const auto e = v.end();
1550 // Iterator for Vector:
1551 ComplexIterator1D target = begin();
1552
1553 for (; i != e; ++i, ++target) *target = *i;
1554
1555 return *this;
1556}
1557
1558//Functions operating on the complex vectors
1559
1562 const ConstComplexVectorView& b) {
1563 // Check dimensions:
1564 ARTS_ASSERT(a.nelem() == b.nelem());
1565
1566 const ConstComplexIterator1D ae = a.end();
1567 ConstComplexIterator1D ai = a.begin();
1568 ConstComplexIterator1D bi = b.begin();
1569
1570 Complex res = 0;
1571 for (; ai != ae; ++ai, ++bi) res += (*ai) * (*bi);
1572
1573 return res;
1574}
1575
1577
1588 const ConstComplexVectorView& x) {
1589 ARTS_ASSERT(x.nelem() == M.nrows());
1590 ARTS_ASSERT(y.nelem() == M.ncols());
1591
1593 if (y.mdata == x.mdata)
1594 eigen_y = MapToEigen(M) * MapToEigenRow(x);
1595 else
1596 eigen_y.noalias() = MapToEigen(M) * MapToEigenRow(x);
1597}
1598
1600
1612 const ConstComplexMatrixView& B,
1613 const ConstComplexMatrixView& C) {
1614 ARTS_ASSERT(B.nrows() == A.nrows());
1615 ARTS_ASSERT(C.ncols() == A.ncols());
1616 ARTS_ASSERT(B.ncols() == C.nrows());
1617
1618 ComplexMatrixViewMap eigen_A = MapToEigen(A);
1619 if (A.mdata == B.mdata || A.mdata == C.mdata)
1620 eigen_A = MapToEigen(B) * MapToEigen(C);
1621 else
1622 eigen_A.noalias() = MapToEigen(B) * MapToEigen(C);
1623}
1625 const ConstComplexMatrixView& B,
1626 const ConstMatrixView& C) {
1627 ARTS_ASSERT(B.nrows() == A.nrows());
1628 ARTS_ASSERT(C.ncols() == A.ncols());
1629 ARTS_ASSERT(B.ncols() == C.nrows());
1630
1631 ComplexMatrixViewMap eigen_A = MapToEigen(A);
1632 if (A.mdata == B.mdata)
1633 eigen_A = MapToEigen(B) * MapToEigen(C);
1634 else
1635 eigen_A.noalias() = MapToEigen(B) * MapToEigen(C);
1636}
1638 const ConstMatrixView& B,
1639 const ConstComplexMatrixView& C) {
1640 ARTS_ASSERT(B.nrows() == A.nrows());
1641 ARTS_ASSERT(C.ncols() == A.ncols());
1642 ARTS_ASSERT(B.ncols() == C.nrows());
1643
1644 ComplexMatrixViewMap eigen_A = MapToEigen(A);
1645 if (A.mdata == C.mdata)
1646 eigen_A = MapToEigen(B) * MapToEigen(C);
1647 else
1648 eigen_A.noalias() = MapToEigen(B) * MapToEigen(C);
1649}
1651 const ConstMatrixView& B,
1652 const ConstMatrixView& C) {
1653 ARTS_ASSERT(B.nrows() == A.nrows());
1654 ARTS_ASSERT(C.ncols() == A.ncols());
1655 ARTS_ASSERT(B.ncols() == C.nrows());
1656
1657 ComplexMatrixViewMap eigen_A = MapToEigen(A);
1658 eigen_A.noalias() = MapToEigen(B) * MapToEigen(C);
1659}
1660
1661// Converts constant matrix to constant eigen map
1664 A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1665 A.nrows(),
1666 A.ncols(),
1668}
1669
1670// Converts constant vector to constant eigen row-view
1673 A.nelem(),
1674 1,
1675 StrideType(A.mrange.get_stride(), 1));
1676}
1677
1678// Converts constant vector to constant eigen row-view
1680 return MapToEigen(A);
1681}
1682
1683// Converts constant vector to constant eigen column-view
1686 1,
1687 A.nelem(),
1688 StrideType(1, A.mrange.get_stride()));
1689}
1690
1691// Converts matrix to eigen map
1693 return ComplexMatrixViewMap(
1694 A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1695 A.nrows(),
1696 A.ncols(),
1698}
1699
1700// Converts vector to eigen map row-view
1703 A.nelem(),
1704 1,
1705 StrideType(A.mrange.get_stride(), 1));
1706}
1707
1708// Converts vector to eigen map row-view
1710 return MapToEigen(A);
1711}
1712
1713// Converts vector to eigen map column-view
1716 1,
1717 A.nelem(),
1718 StrideType(1, A.mrange.get_stride()));
1719}
1720
1722// Helper function for debugging
1723#ifndef NDEBUG
1724
1739 return mv(r, c);
1740}
1741
1742#endif
Interface for BLAS library.
This can be used to make arrays out of anything.
Definition: array.h:108
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:197
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()=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:454
The const row iterator class for sub matrices.
Definition: matpackI.h:857
A constant view of a Matrix.
Definition: matpackI.h:1050
ConstIterator2D begin() const ARTS_NOEXCEPT
Return const iterator to first row.
Definition: matpackI.cc:469
Index nrows() const noexcept
Definition: matpackI.h:1061
Index ncols() const noexcept
Definition: matpackI.h:1062
A constant view of a Vector.
Definition: matpackI.h:517
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:541
The range class.
Definition: matpackI.h:166
Index mstart
The start index.
Definition: matpackI.h:368
Index mstride
The stride.
Definition: matpackI.h:377
constexpr Index get_start() const noexcept
Returns the start index of the range.
Definition: matpackI.h:347
Index mextent
The number of elements.
Definition: matpackI.h:375
constexpr Index get_stride() const noexcept
Returns the stride of the range.
Definition: matpackI.h:351
The Vector class.
Definition: matpackI.h:908
#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:113
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.
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: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)
#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