ARTS 2.5.9 (git: 825fa5f2)
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
27#include "debug.h"
28#include "exceptions.h"
29
30// Functions for ConstTensor3View:
31// ------------------------------
32
37 const Range& r,
38 const Range& c) const {
39 return ConstTensor3View(mdata, mpr, mrr, mcr, p, r, c);
40}
41
45 const Range& r,
46 Index c) const {
47 // Check that c is valid:
48 ARTS_ASSERT(0 <= c);
50
51 return ConstMatrixView(mdata + mcr.mstart + c * mcr.mstride, mpr, mrr, p, r);
52}
53
57 Index r,
58 const Range& c) const {
59 // Check that r is valid:
60 ARTS_ASSERT(0 <= r);
62
63 return ConstMatrixView(mdata + mrr.mstart + r * mrr.mstride, mpr, mcr, p, c);
64}
65
69 const Range& r,
70 const Range& c) const {
71 // Check that p is valid:
72 ARTS_ASSERT(0 <= p);
74
75 return ConstMatrixView(mdata + mpr.mstart + p * mpr.mstride, mrr, mcr, r, c);
76}
77
81 Index r,
82 const Range& c) const {
83 // Check that p and r are valid:
84 ARTS_ASSERT(0 <= p);
85 ARTS_ASSERT(0 <= r);
88
89 return ConstVectorView(
90 mdata + mpr.mstart + p * mpr.mstride + mrr.mstart + r * mrr.mstride,
91 mcr,
92 c);
93}
94
98 const Range& r,
99 Index c) const {
100 // Check that p and c are valid:
101 ARTS_ASSERT(0 <= p);
102 ARTS_ASSERT(0 <= c);
105
106 return ConstVectorView(
108 mrr,
109 r);
110}
111
115 Index r,
116 Index c) const {
117 // Check that r and c are valid:
118 ARTS_ASSERT(0 <= r);
119 ARTS_ASSERT(0 <= c);
122
123 return ConstVectorView(
125 mpr,
126 p);
127}
128
132 mpr.mstride);
133}
134
137 return ConstIterator3D(
140 mpr.mstride);
141}
142
145 : mpr(0, 1, a.mrr.mextent * a.mcr.mextent),
146 mrr(a.mrr),
147 mcr(a.mcr),
148 mdata(a.mdata) {
149 // Nothing to do here.
150}
151
157 const Range& pr,
158 const Range& rr,
159 const Range& cr)
160 : mpr(pr), mrr(rr), mcr(cr), mdata(data) {
161 // Nothing to do here.
162}
163
172 const Range& pp,
173 const Range& pr,
174 const Range& pc,
175 const Range& np,
176 const Range& nr,
177 const Range& nc)
178 : mpr(pp, np), mrr(pr, nr), mcr(pc, nc), mdata(data) {
179 // Nothing to do here.
180}
181
185std::ostream& operator<<(std::ostream& os, const ConstTensor3View& v) {
186 // Page iterators:
187 ConstIterator3D ip = v.begin();
188 const ConstIterator3D end_page = v.end();
189
190 if (ip != end_page) {
191 os << *ip;
192 ++ip;
193 }
194
195 for (; ip != end_page; ++ip) {
196 os << "\n\n";
197 os << *ip;
198 }
199
200 return os;
201}
202
203// Functions for Tensor3View:
204// -------------------------
205
210 const Range& r,
211 const Range& c) {
212 return Tensor3View(mdata, mpr, mrr, mcr, p, r, c);
213}
214
218 // Check that c is valid:
219 ARTS_ASSERT(0 <= c);
221
222 return MatrixView(mdata + mcr.mstart + c * mcr.mstride, mpr, mrr, p, r);
223}
224
228 // Check that r is valid:
229 ARTS_ASSERT(0 <= r);
231
232 return MatrixView(mdata + mrr.mstart + r * mrr.mstride, mpr, mcr, p, c);
233}
234
238 // Check that p is valid:
239 ARTS_ASSERT(0 <= p);
241
242 return MatrixView(mdata + mpr.mstart + p * mpr.mstride, mrr, mcr, r, c);
243}
244
248 // Check that p and r are valid:
249 ARTS_ASSERT(0 <= p);
250 ARTS_ASSERT(0 <= r);
253
254 return VectorView(
255 mdata + mpr.mstart + p * mpr.mstride + mrr.mstart + r * mrr.mstride,
256 mcr,
257 c);
258}
259
263 // Check that p and c are valid:
264 ARTS_ASSERT(0 <= p);
265 ARTS_ASSERT(0 <= c);
268
269 return VectorView(
271 mrr,
272 r);
273}
274
278 // Check that r and r are valid:
279 ARTS_ASSERT(0 <= r);
280 ARTS_ASSERT(0 <= c);
283
284 return VectorView(
286 mpr,
287 p);
288}
289
297 ARTS_ASSERT(mpr.mstart == 0 and
298 (mpr.mstride == mcr.mextent * mrr.mextent or size() == 0),
299 "Page ",
300 mpr)
301 ARTS_ASSERT(mrr.mstart == 0 and (mrr.mstride == mcr.mextent or size() == 0),
302 "Row ",
303 mrr)
305 mcr.mstart == 0 and (mcr.mstride == 1 or size() == 0), "Column ", mcr)
306 return mdata;
307}
308
316 ARTS_ASSERT(mpr.mstart == 0 and
317 (mpr.mstride == mcr.mextent * mrr.mextent or size() == 0),
318 "Page ",
319 mpr)
320 ARTS_ASSERT(mrr.mstart == 0 and (mrr.mstride == mcr.mextent or size() == 0),
321 "Row ",
322 mrr)
324 mcr.mstart == 0 and (mcr.mstride == 1 or size() == 0), "Column ", mcr)
325 return mdata;
326}
327
331}
332
335 return Iterator3D(
337 mpr.mstride);
338}
339
345 // Check that sizes are compatible:
349
350 copy(m.begin(), m.end(), begin());
351 return *this;
352}
353
360 // Check that sizes are compatible:
364
365 copy(m.begin(), m.end(), begin());
366 return *this;
367}
368
373 // Check that sizes are compatible:
377
378 copy(m.begin(), m.end(), begin());
379 return *this;
380}
381
385 copy(x, begin(), end());
386 return *this;
387}
388
389// Some little helper functions:
390//------------------------------
391
392Numeric add(Numeric x, Numeric y) { return x + y; }
393
396 const Iterator3D ep = end();
397 for (Iterator3D p = begin(); p != ep; ++p) {
398 *p *= x;
399 }
400 return *this;
401}
402
405 const Iterator3D ep = end();
406 for (Iterator3D p = begin(); p != ep; ++p) {
407 *p /= x;
408 }
409 return *this;
410}
411
414 const Iterator3D ep = end();
415 for (Iterator3D p = begin(); p != ep; ++p) {
416 *p += x;
417 }
418 return *this;
419}
420
423 const Iterator3D ep = end();
424 for (Iterator3D p = begin(); p != ep; ++p) {
425 *p -= x;
426 }
427 return *this;
428}
429
432 ARTS_ASSERT(npages() == x.npages());
433 ARTS_ASSERT(nrows() == x.nrows());
434 ARTS_ASSERT(ncols() == x.ncols());
435 ConstIterator3D xp = x.begin();
436 Iterator3D p = begin();
437 const Iterator3D ep = end();
438 for (; p != ep; ++p, ++xp) {
439 *p *= *xp;
440 }
441 return *this;
442}
443
446 ARTS_ASSERT(npages() == x.npages());
447 ARTS_ASSERT(nrows() == x.nrows());
448 ARTS_ASSERT(ncols() == x.ncols());
449 ConstIterator3D xp = x.begin();
450 Iterator3D p = begin();
451 const Iterator3D ep = end();
452 for (; p != ep; ++p, ++xp) {
453 *p /= *xp;
454 }
455 return *this;
456}
457
460 ARTS_ASSERT(npages() == x.npages());
461 ARTS_ASSERT(nrows() == x.nrows());
462 ARTS_ASSERT(ncols() == x.ncols());
463 ConstIterator3D xp = x.begin();
464 Iterator3D p = begin();
465 const Iterator3D ep = end();
466 for (; p != ep; ++p, ++xp) {
467 *p += *xp;
468 }
469 return *this;
470}
471
474 ARTS_ASSERT(npages() == x.npages());
475 ARTS_ASSERT(nrows() == x.nrows());
476 ARTS_ASSERT(ncols() == x.ncols());
477 ConstIterator3D xp = x.begin();
478 Iterator3D p = begin();
479 const Iterator3D ep = end();
480 for (; p != ep; ++p, ++xp) {
481 *p -= *xp;
482 }
483 return *this;
484}
485
489 a.mdata, Range(0, 1, a.mrr.mextent * a.mcr.mextent), a.mrr, a.mcr) {
490 // Nothing to do here.
491}
492
497 const Range& pr,
498 const Range& rr,
499 const Range& cr)
500 : ConstTensor3View(data, pr, rr, cr) {
501 // Nothing to do here.
502}
503
521 const Range& pp,
522 const Range& pr,
523 const Range& pc,
524 const Range& np,
525 const Range& nr,
526 const Range& nc)
527 : ConstTensor3View(data, pp, pr, pc, np, nr, nc) {
528 // Nothing to do here.
529}
530
536 const ConstIterator3D& end,
537 Iterator3D target) {
538 for (; origin != end; ++origin, ++target) {
539 // We use the copy function for the next smaller rank of tensor
540 // recursively:
541 copy(origin->begin(), origin->end(), target->begin());
542 }
543}
544
546void copy(Numeric x, Iterator3D target, const Iterator3D& end) {
547 for (; target != end; ++target) {
548 // We use the copy function for the next smaller rank of tensor
549 // recursively:
550 copy(x, target->begin(), target->end());
551 }
552}
553
554// Functions for Tensor3:
555// ---------------------
556
560 : Tensor3View(new Numeric[p * r * c],
561 Range(0, p, r * c),
562 Range(0, r, c),
563 Range(0, c)) {
564 // Nothing to do here.
565}
566
569 : Tensor3View(new Numeric[p * r * c],
570 Range(0, p, r * c),
571 Range(0, r, c),
572 Range(0, c)) {
573 // Here we can access the raw memory directly, for slightly
574 // increased efficiency:
575 std::fill_n(mdata, p * r * c, fill);
576}
577
581 : Tensor3View(new Numeric[m.npages() * m.nrows() * m.ncols()],
582 Range(0, m.npages(), m.nrows() * m.ncols()),
583 Range(0, m.nrows(), m.ncols()),
584 Range(0, m.ncols())) {
585 copy(m.begin(), m.end(), begin());
586}
587
591 : Tensor3View(new Numeric[m.npages() * m.nrows() * m.ncols()],
592 Range(0, m.npages(), m.nrows() * m.ncols()),
593 Range(0, m.nrows(), m.ncols()),
594 Range(0, m.ncols())) {
595 // There is a catch here: If m is an empty tensor, then it will have
596 // dimensions of size 0. But these are used to initialize the stride
597 // for higher dimensions! Thus, this method has to be consistent
598 // with the behaviour of Range::Range. For now, Range::Range allows
599 // also stride 0.
600 std::memcpy(mdata, m.mdata, npages() * nrows() * ncols() * sizeof(Numeric));
601}
602
604
628 if (this != &x) {
629 resize(x.npages(), x.nrows(), x.ncols());
630 std::memcpy(mdata, x.mdata, npages() * nrows() * ncols() * sizeof(Numeric));
631 }
632 return *this;
633}
634
637 if (this != &x) {
638 delete[] mdata;
639 mdata = x.mdata;
640 mpr = x.mpr;
641 mrr = x.mrr;
642 mcr = x.mcr;
643 x.mpr = Range(0, 0);
644 x.mrr = Range(0, 0);
645 x.mcr = Range(0, 0);
646 x.mdata = nullptr;
647 }
648 return *this;
649}
650
654 std::fill_n(mdata, npages() * nrows() * ncols(), x);
655 return *this;
656}
657
662 ARTS_ASSERT(0 <= p);
663 ARTS_ASSERT(0 <= r);
664 ARTS_ASSERT(0 <= c);
665
666 if (mpr.mextent != p || mrr.mextent != r || mcr.mextent != c) {
667 delete[] mdata;
668 mdata = new Numeric[p * r * c];
669
670 mpr.mstart = 0;
671 mpr.mextent = p;
672 mpr.mstride = r * c;
673
674 mrr.mstart = 0;
675 mrr.mextent = r;
676 mrr.mstride = c;
677
678 mcr.mstart = 0;
679 mcr.mextent = c;
680 mcr.mstride = 1;
681 }
682}
683
685void swap(Tensor3& t1, Tensor3& t2) noexcept {
686 using std::swap;
687 swap(t1.mpr, t2.mpr);
688 swap(t1.mrr, t2.mrr);
689 swap(t1.mcr, t2.mcr);
690 swap(t1.mdata, t2.mdata);
691}
692
696 // cout << "Destroying a Tensor3:\n"
697 // << *this << "\n........................................\n";
698 delete[] mdata;
699}
700
716void transform(Tensor3View y, double (&my_func)(double), ConstTensor3View x) {
717 // Check dimensions:
718 ARTS_ASSERT(y.npages() == x.npages());
719 ARTS_ASSERT(y.nrows() == x.nrows());
720 ARTS_ASSERT(y.ncols() == x.ncols());
721
722 const ConstIterator3D xe = x.end();
723 ConstIterator3D xi = x.begin();
724 Iterator3D yi = y.begin();
725 for (; xi != xe; ++xi, ++yi) {
726 // Use the transform function of lower dimensional tensors
727 // recursively:
728 transform(*yi, my_func, *xi);
729 }
730}
731
734 const ConstIterator3D xe = x.end();
735 ConstIterator3D xi = x.begin();
736
737 // Initial value for max:
738 Numeric themax = max(*xi);
739 ++xi;
740
741 for (; xi != xe; ++xi) {
742 // Use the max function of lower dimensional tensors
743 // recursively:
744 Numeric maxi = max(*xi);
745 if (maxi > themax) themax = maxi;
746 }
747
748 return themax;
749}
750
753 const ConstIterator3D xe = x.end();
754 ConstIterator3D xi = x.begin();
755
756 // Initial value for min:
757 Numeric themin = min(*xi);
758 ++xi;
759
760 for (; xi != xe; ++xi) {
761 // Use the min function of lower dimensional tensors
762 // recursively:
763 Numeric mini = min(*xi);
764 if (mini < themin) themin = mini;
765 }
766
767 return themin;
768}
769
771// Helper function for debugging
772#ifndef NDEBUG
773
789 return tv(p, r, c);
790}
791
792#endif
794
796
812 ARTS_ASSERT(A.npages() == B.nelem());
813 ARTS_ASSERT(A.nrows() == C.nrows());
814 ARTS_ASSERT(A.ncols() == C.ncols());
815
816 for (Index ii = 0; ii < B.nelem(); ii++) {
817 A(ii, joker, joker) = C;
818 A(ii, joker, joker) *= B[ii];
819 }
820}
base max(const Array< base > &x)
Max function.
Definition: array.h:145
base min(const Array< base > &x)
Min function.
Definition: array.h:161
Const version of Iterator3D.
Definition: matpackIII.h:78
A constant view of a Matrix.
Definition: matpackI.h:1065
ConstIterator2D begin() const ARTS_NOEXCEPT
Return const iterator to first row.
Definition: matpackI.cc:449
Index nrows() const noexcept
Definition: matpackI.h:1079
ConstIterator2D end() const ARTS_NOEXCEPT
Return const iterator behind last row.
Definition: matpackI.cc:454
Index ncols() const noexcept
Definition: matpackI.h:1080
A constant view of a Tensor3.
Definition: matpackIII.h:133
ConstTensor3View operator()(const Range &p, const Range &r, const Range &c) const
Const index operator for subrange.
Definition: matpackIII.cc:36
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:145
ConstIterator3D begin() const
Return const iterator to first page.
Definition: matpackIII.cc:130
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackIII.h:239
ConstIterator3D end() const
Return const iterator behind last page.
Definition: matpackIII.cc:136
ConstTensor3View()=default
Range mrr
The row range of mdata that is actually used.
Definition: matpackIII.h:235
Range mcr
The column range of mdata that is actually used.
Definition: matpackIII.h:237
Index size() const noexcept
Definition: matpackIII.h:154
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:148
Range mpr
The page range of mdata that is actually used.
Definition: matpackIII.h:233
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:151
A constant view of a Vector.
Definition: matpackI.h:521
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
Implementation of Tensors of Rank 3.
Definition: matpackIII.h:37
The MatrixView class.
Definition: matpackI.h:1188
Iterator2D end() ARTS_NOEXCEPT
Return iterator behind last row.
Definition: matpackI.cc:606
Iterator2D begin() ARTS_NOEXCEPT
‍** Return const iterator to first row. Has to be redefined here, since it is
Definition: matpackI.cc:601
The range class.
Definition: matpackI.h:160
Index mstart
The start index.
Definition: matpackI.h:377
Index mstride
The stride.
Definition: matpackI.h:382
Index mextent
The number of elements.
Definition: matpackI.h:380
The Tensor3View class.
Definition: matpackIII.h:251
Iterator3D begin()
Return iterator to first page.
Definition: matpackIII.cc:329
Tensor3View & operator/=(Numeric x)
Division by scalar.
Definition: matpackIII.cc:404
Tensor3View & operator-=(Numeric x)
Subtraction of scalar.
Definition: matpackIII.cc:422
Tensor3View & operator=(const ConstTensor3View &v)
Assignment operator.
Definition: matpackIII.cc:344
Tensor3View & operator*=(Numeric x)
Multiplication by scalar.
Definition: matpackIII.cc:395
Tensor3View & operator+=(Numeric x)
Addition of scalar.
Definition: matpackIII.cc:413
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array.
Definition: matpackIII.cc:315
Tensor3View operator()(const Range &p, const Range &r, const Range &c)
Index operator for subrange.
Definition: matpackIII.cc:209
Tensor3View()=default
Iterator3D end()
Return iterator behind last page.
Definition: matpackIII.cc:334
The Tensor3 class.
Definition: matpackIII.h:352
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:661
~Tensor3() noexcept override
Destructor for Tensor3.
Definition: matpackIII.cc:695
Tensor3 & operator=(const matpack::tensor3_like_not_tensor3 auto &init)
Set from a tensor type.
Definition: matpackIII.h:372
Tensor3()=default
The VectorView class.
Definition: matpackI.h:674
Helper macros for debugging.
#define ARTS_NOEXCEPT
Definition: debug.h:79
#define ARTS_ASSERT(condition,...)
Definition: debug.h:82
The declarations of all the exception classes.
void mult(Tensor3View A, const ConstVectorView B, const ConstMatrixView C)
mult Tensor3
Definition: matpackIII.cc:811
Numeric debug_tensor3view_get_elem(Tensor3View &tv, Index p, Index r, Index c)
Helper function to access tensor elements.
Definition: matpackIII.cc:788
Numeric add(Numeric x, Numeric y)
Definition: matpackIII.cc:392
void swap(Tensor3 &t1, Tensor3 &t2) noexcept
Swaps two objects.
Definition: matpackIII.cc:685
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:716
void copy(ConstIterator3D origin, const ConstIterator3D &end, Iterator3D target)
Copy data between begin and end to target.
Definition: matpackIII.cc:535
std::ostream & operator<<(std::ostream &os, const ConstTensor3View &v)
Output operator.
Definition: matpackIII.cc:185
void copy(ConstIterator3D origin, const ConstIterator3D &end, Iterator3D target)
Copy data between begin and end to target.
Definition: matpackIII.cc:535
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
const Joker joker
#define v
#define a
#define c