ARTS 2.5.0 (git: 9ee3ac6c)
matpackIII.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
25#include "matpackIII.h"
26#include "exceptions.h"
27
28// Functions for ConstTensor3View:
29// ------------------------------
30
32
37 return (npages() == 0 || nrows() == 0 || ncols() == 0);
38}
39
44 const Range& r,
45 const Range& c) const {
46 return ConstTensor3View(mdata, mpr, mrr, mcr, p, r, c);
47}
48
52 const Range& r,
53 Index c) const {
54 // Check that c is valid:
55 ARTS_ASSERT(0 <= c);
57
58 return ConstMatrixView(mdata + mcr.mstart + c * mcr.mstride, mpr, mrr, p, r);
59}
60
64 Index r,
65 const Range& c) const {
66 // Check that r is valid:
67 ARTS_ASSERT(0 <= r);
69
70 return ConstMatrixView(mdata + mrr.mstart + r * mrr.mstride, mpr, mcr, p, c);
71}
72
76 const Range& r,
77 const Range& c) const {
78 // Check that p is valid:
79 ARTS_ASSERT(0 <= p);
81
82 return ConstMatrixView(mdata + mpr.mstart + p * mpr.mstride, mrr, mcr, r, c);
83}
84
88 Index r,
89 const Range& c) const {
90 // Check that p and r are valid:
91 ARTS_ASSERT(0 <= p);
92 ARTS_ASSERT(0 <= r);
95
96 return ConstVectorView(
97 mdata + mpr.mstart + p * mpr.mstride + mrr.mstart + r * mrr.mstride,
98 mcr,
99 c);
100}
101
105 const Range& r,
106 Index c) const {
107 // Check that p and c are valid:
108 ARTS_ASSERT(0 <= p);
109 ARTS_ASSERT(0 <= c);
112
113 return ConstVectorView(
115 mrr,
116 r);
117}
118
122 Index r,
123 Index c) const {
124 // Check that r and c are valid:
125 ARTS_ASSERT(0 <= r);
126 ARTS_ASSERT(0 <= c);
129
130 return ConstVectorView(
132 mpr,
133 p);
134}
135
139 mpr.mstride);
140}
141
144 return ConstIterator3D(
147 mpr.mstride);
148}
149
152 : mpr(0, 1, a.mrr.mextent * a.mcr.mextent),
153 mrr(a.mrr),
154 mcr(a.mcr),
155 mdata(a.mdata) {
156 // Nothing to do here.
157}
158
164 const Range& pr,
165 const Range& rr,
166 const Range& cr)
167 : mpr(pr), mrr(rr), mcr(cr), mdata(data) {
168 // Nothing to do here.
169}
170
179 const Range& pp,
180 const Range& pr,
181 const Range& pc,
182 const Range& np,
183 const Range& nr,
184 const Range& nc)
185 : mpr(pp, np), mrr(pr, nr), mcr(pc, nc), mdata(data) {
186 // Nothing to do here.
188
192std::ostream& operator<<(std::ostream& os, const ConstTensor3View& v) {
193 // Page iterators:
194 ConstIterator3D ip = v.begin();
195 const ConstIterator3D end_page = v.end();
196
197 if (ip != end_page) {
198 os << *ip;
199 ++ip;
200 }
201
202 for (; ip != end_page; ++ip) {
203 os << "\n\n";
204 os << *ip;
205 }
206
207 return os;
208}
209
210// Functions for Tensor3View:
211// -------------------------
212
217 const Range& r,
218 const Range& c) {
219 return Tensor3View(mdata, mpr, mrr, mcr, p, r, c);
220}
221
225 // Check that c is valid:
226 ARTS_ASSERT(0 <= c);
228
229 return MatrixView(mdata + mcr.mstart + c * mcr.mstride, mpr, mrr, p, r);
230}
231
235 // Check that r is valid:
236 ARTS_ASSERT(0 <= r);
238
239 return MatrixView(mdata + mrr.mstart + r * mrr.mstride, mpr, mcr, p, c);
240}
241
245 // Check that p is valid:
246 ARTS_ASSERT(0 <= p);
248
249 return MatrixView(mdata + mpr.mstart + p * mpr.mstride, mrr, mcr, r, c);
250}
251
255 // Check that p and r are valid:
256 ARTS_ASSERT(0 <= p);
257 ARTS_ASSERT(0 <= r);
260
261 return VectorView(
262 mdata + mpr.mstart + p * mpr.mstride + mrr.mstart + r * mrr.mstride,
263 mcr,
264 c);
265}
266
270 // Check that p and c are valid:
271 ARTS_ASSERT(0 <= p);
272 ARTS_ASSERT(0 <= c);
275
276 return VectorView(
278 mrr,
279 r);
280}
281
285 // Check that r and r are valid:
286 ARTS_ASSERT(0 <= r);
287 ARTS_ASSERT(0 <= c);
290
291 return VectorView(
293 mpr,
294 p);
295}
296
304 ARTS_ASSERT (not (mpr.mstart != 0 || mpr.mstride != mrr.mextent * mcr.mextent ||
305 mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
306 mcr.mstride != 1),
307 "A Tensor3View can only be converted to a plain C-array if it's pointing to a continuous block of data");
308 return mdata;
309}
310
318 ARTS_ASSERT (not (mpr.mstart != 0 || mpr.mstride != mrr.mextent * mcr.mextent ||
319 mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
320 mcr.mstride != 1),
321 "A Tensor3View can only be converted to a plain C-array if it's pointing to a continuous block of data");
322 return mdata;
323}
324
328}
329
332 return Iterator3D(
334 mpr.mstride);
335}
336
342 // Check that sizes are compatible:
346
347 copy(m.begin(), m.end(), begin());
348 return *this;
349}
350
357 // Check that sizes are compatible:
361
362 copy(m.begin(), m.end(), begin());
363 return *this;
364}
365
370 // Check that sizes are compatible:
374
375 copy(m.begin(), m.end(), begin());
376 return *this;
377}
378
382 copy(x, begin(), end());
383 return *this;
384}
385
386// Some little helper functions:
387//------------------------------
388
389Numeric add(Numeric x, Numeric y) { return x + y; }
390
393 const Iterator3D ep = end();
394 for (Iterator3D p = begin(); p != ep; ++p) {
395 *p *= x;
396 }
397 return *this;
398}
399
402 const Iterator3D ep = end();
403 for (Iterator3D p = begin(); p != ep; ++p) {
404 *p /= x;
405 }
406 return *this;
407}
408
411 const Iterator3D ep = end();
412 for (Iterator3D p = begin(); p != ep; ++p) {
413 *p += x;
414 }
415 return *this;
416}
417
420 const Iterator3D ep = end();
421 for (Iterator3D p = begin(); p != ep; ++p) {
422 *p -= x;
423 }
424 return *this;
425}
426
429 ARTS_ASSERT(npages() == x.npages());
430 ARTS_ASSERT(nrows() == x.nrows());
431 ARTS_ASSERT(ncols() == x.ncols());
432 ConstIterator3D xp = x.begin();
433 Iterator3D p = begin();
434 const Iterator3D ep = end();
435 for (; p != ep; ++p, ++xp) {
436 *p *= *xp;
437 }
438 return *this;
439}
440
443 ARTS_ASSERT(npages() == x.npages());
444 ARTS_ASSERT(nrows() == x.nrows());
445 ARTS_ASSERT(ncols() == x.ncols());
446 ConstIterator3D xp = x.begin();
447 Iterator3D p = begin();
448 const Iterator3D ep = end();
449 for (; p != ep; ++p, ++xp) {
450 *p /= *xp;
451 }
452 return *this;
453}
454
457 ARTS_ASSERT(npages() == x.npages());
458 ARTS_ASSERT(nrows() == x.nrows());
459 ARTS_ASSERT(ncols() == x.ncols());
460 ConstIterator3D xp = x.begin();
461 Iterator3D p = begin();
462 const Iterator3D ep = end();
463 for (; p != ep; ++p, ++xp) {
464 *p += *xp;
465 }
466 return *this;
467}
468
471 ARTS_ASSERT(npages() == x.npages());
472 ARTS_ASSERT(nrows() == x.nrows());
473 ARTS_ASSERT(ncols() == x.ncols());
474 ConstIterator3D xp = x.begin();
475 Iterator3D p = begin();
476 const Iterator3D ep = end();
477 for (; p != ep; ++p, ++xp) {
478 *p -= *xp;
479 }
480 return *this;
481}
482
486 a.mdata, Range(0, 1, a.mrr.mextent * a.mcr.mextent), a.mrr, a.mcr) {
487 // Nothing to do here.
488}
489
494 const Range& pr,
495 const Range& rr,
496 const Range& cr)
497 : ConstTensor3View(data, pr, rr, cr) {
498 // Nothing to do here.
499}
500
518 const Range& pp,
519 const Range& pr,
520 const Range& pc,
521 const Range& np,
522 const Range& nr,
523 const Range& nc)
524 : ConstTensor3View(data, pp, pr, pc, np, nr, nc) {
525 // Nothing to do here.
526}
527
533 const ConstIterator3D& end,
534 Iterator3D target) {
535 for (; origin != end; ++origin, ++target) {
536 // We use the copy function for the next smaller rank of tensor
537 // recursively:
538 copy(origin->begin(), origin->end(), target->begin());
539 }
540}
541
543void copy(Numeric x, Iterator3D target, const Iterator3D& end) {
544 for (; target != end; ++target) {
545 // We use the copy function for the next smaller rank of tensor
546 // recursively:
547 copy(x, target->begin(), target->end());
548 }
549}
550
551// Functions for Tensor3:
552// ---------------------
553
557 : Tensor3View(new Numeric[p * r * c],
558 Range(0, p, r * c),
559 Range(0, r, c),
560 Range(0, c)) {
561 // Nothing to do here.
562}
563
566 : Tensor3View(new Numeric[p * r * c],
567 Range(0, p, r * c),
568 Range(0, r, c),
569 Range(0, c)) {
570 // Here we can access the raw memory directly, for slightly
571 // increased efficiency:
572 std::fill_n(mdata, p * r * c, fill);
573}
574
578 : Tensor3View(new Numeric[m.npages() * m.nrows() * m.ncols()],
579 Range(0, m.npages(), m.nrows() * m.ncols()),
580 Range(0, m.nrows(), m.ncols()),
581 Range(0, m.ncols())) {
582 copy(m.begin(), m.end(), begin());
583}
584
588 : Tensor3View(new Numeric[m.npages() * m.nrows() * m.ncols()],
589 Range(0, m.npages(), m.nrows() * m.ncols()),
590 Range(0, m.nrows(), m.ncols()),
591 Range(0, m.ncols())) {
592 // There is a catch here: If m is an empty tensor, then it will have
593 // dimensions of size 0. But these are used to initialize the stride
594 // for higher dimensions! Thus, this method has to be consistent
595 // with the behaviour of Range::Range. For now, Range::Range allows
596 // also stride 0.
597 std::memcpy(mdata, m.mdata, npages() * nrows() * ncols() * sizeof(Numeric));
598}
599
601
625 if (this != &x) {
626 resize(x.npages(), x.nrows(), x.ncols());
627 std::memcpy(mdata, x.mdata, npages() * nrows() * ncols() * sizeof(Numeric));
628 }
629 return *this;
630}
631
634 if (this != &x) {
635 delete[] mdata;
636 mdata = x.mdata;
637 mpr = x.mpr;
638 mrr = x.mrr;
639 mcr = x.mcr;
640 x.mpr = Range(0, 0);
641 x.mrr = Range(0, 0);
642 x.mcr = Range(0, 0);
643 x.mdata = nullptr;
644 }
645 return *this;
646}
647
651 std::fill_n(mdata, npages() * nrows() * ncols(), x);
652 return *this;
653}
654
659 ARTS_ASSERT(0 <= p);
660 ARTS_ASSERT(0 <= r);
661 ARTS_ASSERT(0 <= c);
662
663 if (mpr.mextent != p || mrr.mextent != r || mcr.mextent != c) {
664 delete[] mdata;
665 mdata = new Numeric[p * r * c];
666
667 mpr.mstart = 0;
668 mpr.mextent = p;
669 mpr.mstride = r * c;
670
671 mrr.mstart = 0;
672 mrr.mextent = r;
673 mrr.mstride = c;
674
675 mcr.mstart = 0;
676 mcr.mextent = c;
677 mcr.mstride = 1;
678 }
679}
680
682void swap(Tensor3& t1, Tensor3& t2) {
683 std::swap(t1.mpr, t2.mpr);
684 std::swap(t1.mrr, t2.mrr);
685 std::swap(t1.mcr, t2.mcr);
686 std::swap(t1.mdata, t2.mdata);
687}
688
692 // cout << "Destroying a Tensor3:\n"
693 // << *this << "\n........................................\n";
694 delete[] mdata;
695}
696
712void transform(Tensor3View y, double (&my_func)(double), ConstTensor3View x) {
713 // Check dimensions:
714 ARTS_ASSERT(y.npages() == x.npages());
715 ARTS_ASSERT(y.nrows() == x.nrows());
716 ARTS_ASSERT(y.ncols() == x.ncols());
717
718 const ConstIterator3D xe = x.end();
719 ConstIterator3D xi = x.begin();
720 Iterator3D yi = y.begin();
721 for (; xi != xe; ++xi, ++yi) {
722 // Use the transform function of lower dimensional tensors
723 // recursively:
724 transform(*yi, my_func, *xi);
725 }
726}
727
730 const ConstIterator3D xe = x.end();
731 ConstIterator3D xi = x.begin();
732
733 // Initial value for max:
734 Numeric themax = max(*xi);
735 ++xi;
736
737 for (; xi != xe; ++xi) {
738 // Use the max function of lower dimensional tensors
739 // recursively:
740 Numeric maxi = max(*xi);
741 if (maxi > themax) themax = maxi;
742 }
743
744 return themax;
745}
746
749 const ConstIterator3D xe = x.end();
750 ConstIterator3D xi = x.begin();
751
752 // Initial value for min:
753 Numeric themin = min(*xi);
754 ++xi;
755
756 for (; xi != xe; ++xi) {
757 // Use the min function of lower dimensional tensors
758 // recursively:
759 Numeric mini = min(*xi);
760 if (mini < themin) themin = mini;
761 }
762
763 return themin;
764}
765
767// Helper function for debugging
768#ifndef NDEBUG
769
785 return tv(p, r, c);
786}
787
788#endif
790
792
808 ARTS_ASSERT(A.npages() == B.nelem());
809 ARTS_ASSERT(A.nrows() == C.nrows());
810 ARTS_ASSERT(A.ncols() == C.ncols());
811
812 for (Index ii = 0; ii < B.nelem(); ii++) {
813 A(ii, joker, joker) = C;
814 A(ii, joker, joker) *= B[ii];
815 }
816}
Index nrows
Index npages
void * data
Index ncols
Const version of Iterator3D.
Definition: matpackIII.h:76
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
ConstIterator2D end() const ARTS_NOEXCEPT
Return const iterator behind last row.
Definition: matpackI.cc:478
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 Tensor3.
Definition: matpackIII.h:132
ConstTensor3View operator()(const Range &p, const Range &r, const Range &c) const
Const index operator for subrange.
Definition: matpackIII.cc:43
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
ConstIterator3D begin() const
Return const iterator to first page.
Definition: matpackIII.cc:137
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackIII.h:227
ConstIterator3D end() const
Return const iterator behind last page.
Definition: matpackIII.cc:143
ConstTensor3View()=default
Range mrr
The row range of mdata that is actually used.
Definition: matpackIII.h:223
Range mcr
The column range of mdata that is actually used.
Definition: matpackIII.h:225
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:147
Range mpr
The page range of mdata that is actually used.
Definition: matpackIII.h:221
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
bool empty() const
Check if variable is empty.
Definition: matpackIII.cc:36
A constant view of a Vector.
Definition: matpackI.h:489
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
Implementation of Tensors of Rank 3.
Definition: matpackIII.h:34
The MatrixView class.
Definition: matpackI.h:1125
Iterator2D end() ARTS_NOEXCEPT
Return iterator behind last row.
Definition: matpackI.cc:627
Iterator2D begin() ARTS_NOEXCEPT
‍** Return const iterator to first row. Has to be redefined here, since it is
Definition: matpackI.cc:622
The range class.
Definition: matpackI.h:165
Index mstart
The start index.
Definition: matpackI.h:351
Index mstride
The stride.
Definition: matpackI.h:360
Index mextent
The number of elements.
Definition: matpackI.h:358
The Tensor3View class.
Definition: matpackIII.h:239
Iterator3D begin()
Return iterator to first page.
Definition: matpackIII.cc:326
Tensor3View & operator/=(Numeric x)
Division by scalar.
Definition: matpackIII.cc:401
Tensor3View & operator-=(Numeric x)
Subtraction of scalar.
Definition: matpackIII.cc:419
Tensor3View & operator=(const ConstTensor3View &v)
Assignment operator.
Definition: matpackIII.cc:341
Tensor3View & operator*=(Numeric x)
Multiplication by scalar.
Definition: matpackIII.cc:392
Tensor3View & operator+=(Numeric x)
Addition of scalar.
Definition: matpackIII.cc:410
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array.
Definition: matpackIII.cc:317
Tensor3View operator()(const Range &p, const Range &r, const Range &c)
Index operator for subrange.
Definition: matpackIII.cc:216
Tensor3View()=default
Iterator3D end()
Return iterator behind last page.
Definition: matpackIII.cc:331
The Tensor3 class.
Definition: matpackIII.h:339
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:658
virtual ~Tensor3()
Destructor for Tensor3.
Definition: matpackIII.cc:691
Tensor3()=default
Tensor3 & operator=(const Tensor3 &x)
Assignment operator from another tensor.
Definition: matpackIII.cc:624
The VectorView class.
Definition: matpackI.h:626
#define ARTS_NOEXCEPT
Definition: debug.h:80
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
The declarations of all the exception classes.
void mult(Tensor3View A, const ConstVectorView B, const ConstMatrixView C)
mult Tensor3
Definition: matpackIII.cc:807
Numeric debug_tensor3view_get_elem(Tensor3View &tv, Index p, Index r, Index c)
Helper function to access tensor elements.
Definition: matpackIII.cc:784
Numeric add(Numeric x, Numeric y)
Definition: matpackIII.cc:389
Numeric min(const ConstTensor3View &x)
Min function, tensor version.
Definition: matpackIII.cc:748
void transform(Tensor3View y, double(&my_func)(double), ConstTensor3View x)
A generic transform function for tensors, which can be used to implement mathematical functions opera...
Definition: matpackIII.cc:712
void swap(Tensor3 &t1, Tensor3 &t2)
Swaps two objects.
Definition: matpackIII.cc:682
void copy(ConstIterator3D origin, const ConstIterator3D &end, Iterator3D target)
Copy data between begin and end to target.
Definition: matpackIII.cc:532
std::ostream & operator<<(std::ostream &os, const ConstTensor3View &v)
Output operator.
Definition: matpackIII.cc:192
Numeric max(const ConstTensor3View &x)
Max function, tensor version.
Definition: matpackIII.cc:729
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
const Joker joker
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
#define v
#define a
#define c