ARTS  2.0.49
matpackIV.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2008 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 "matpackIV.h"
26 #include "exceptions.h"
27 
30 Tensor3View* Iterator4D::operator->()
31 {
32  return &msv;
33 }
34 
36 Tensor3View& Iterator4D::operator*()
37 {
38  return msv;
39 }
40 
43 const ConstTensor3View* ConstIterator4D::operator->() const
44 {
45  return &msv;
46 }
47 
49 const ConstTensor3View& ConstIterator4D::operator*() const
50 {
51  return msv;
52 }
53 
54 
55 
56 // Functions for ConstTensor4View:
57 // ------------------------------
58 
60 Index ConstTensor4View::nbooks() const
61 {
62  return mbr.mextent;
63 }
64 
66 Index ConstTensor4View::npages() const
67 {
68  return mpr.mextent;
69 }
70 
72 Index ConstTensor4View::nrows() const
73 {
74  return mrr.mextent;
75 }
76 
78 Index ConstTensor4View::ncols() const
79 {
80  return mcr.mextent;
81 }
82 
86 ConstTensor4View ConstTensor4View::operator()(const Range& b,
87  const Range& p,
88  const Range& r,
89  const Range& c) const
90 {
91  return ConstTensor4View(mdata,
92  mbr, mpr, mrr, mcr,
93  b, p, r, c);
94 }
95 
98 ConstTensor3View ConstTensor4View::operator()(const Range& b,
99  const Range& p,
100  const Range& r,
101  Index c ) const
102 {
103  // Check that c is valid:
104  assert( 0 <= c );
105  assert( c < mcr.mextent );
106 
107  return ConstTensor3View(mdata +
108  mcr.mstart + c * mcr.mstride,
109  mbr, mpr, mrr,
110  b, p, r);
111 }
112 
115 ConstTensor3View ConstTensor4View::operator()(const Range& b,
116  const Range& p,
117  Index r,
118  const Range& c) const
119 {
120  // Check that r is valid:
121  assert( 0 <= r );
122  assert( r < mrr.mextent );
123 
124  return ConstTensor3View(mdata +
125  mrr.mstart + r * mrr.mstride,
126  mbr, mpr, mcr,
127  b, p, c);
128 }
129 
132 ConstTensor3View ConstTensor4View::operator()(const Range& b,
133  Index p,
134  const Range& r,
135  const Range& c) const
136 {
137  // Check that p is valid:
138  assert( 0 <= p );
139  assert( p < mpr.mextent );
140 
141  return ConstTensor3View(mdata +
142  mpr.mstart + p * mpr.mstride,
143  mbr, mrr, mcr,
144  b, r, c);
145 }
146 
149 ConstTensor3View ConstTensor4View::operator()(Index b,
150  const Range& p,
151  const Range& r,
152  const Range& c) const
153 {
154  // Check that b is valid:
155  assert( 0 <= b );
156  assert( b < mbr.mextent );
157 
158  return ConstTensor3View(mdata +
159  mbr.mstart + b * mbr.mstride,
160  mpr, mrr, mcr,
161  p, r, c);
162 }
163 
166 ConstMatrixView ConstTensor4View::operator()(const Range& b,
167  const Range& p,
168  Index r,
169  Index c ) const
170 {
171  // Check that r and c are valid:
172  assert( 0 <= r );
173  assert( 0 <= c );
174  assert( r < mrr.mextent );
175  assert( c < mcr.mextent );
176 
177  return ConstMatrixView(mdata +
178  mrr.mstart + r * mrr.mstride +
179  mcr.mstart + c * mcr.mstride,
180  mbr, mpr,
181  b, p);
182 }
183 
186 ConstMatrixView ConstTensor4View::operator()(const Range& b,
187  Index p,
188  const Range& r,
189  Index c ) const
190 {
191  // Check that p and c are valid:
192  assert( 0 <= p );
193  assert( 0 <= c );
194  assert( p < mpr.mextent );
195  assert( c < mcr.mextent );
196 
197  return ConstMatrixView(mdata +
198  mpr.mstart + p * mpr.mstride +
199  mcr.mstart + c * mcr.mstride,
200  mbr, mrr,
201  b, r);
202 }
203 
206 ConstMatrixView ConstTensor4View::operator()(const Range& b,
207  Index p,
208  Index r,
209  const Range& c) const
210 {
211  // Check that p and r are valid:
212  assert( 0 <= p );
213  assert( 0 <= r );
214  assert( p < mpr.mextent );
215  assert( r < mrr.mextent );
216 
217  return ConstMatrixView(mdata +
218  mpr.mstart + p * mpr.mstride +
219  mrr.mstart + r * mrr.mstride,
220  mbr, mcr,
221  b, c);
222 }
223 
226 ConstMatrixView ConstTensor4View::operator()(Index b,
227  const Range& p,
228  Index r,
229  const Range& c) const
230 {
231  // Check that b and r are valid:
232  assert( 0 <= b );
233  assert( 0 <= r );
234  assert( b < mbr.mextent );
235  assert( r < mrr.mextent );
236 
237  return ConstMatrixView(mdata +
238  mbr.mstart + b * mbr.mstride +
239  mrr.mstart + r * mrr.mstride,
240  mpr, mcr,
241  p, c);
242 }
243 
246 ConstMatrixView ConstTensor4View::operator()(Index b,
247  const Range& p,
248  const Range& r,
249  Index c ) const
250 {
251  // Check that b and c are valid:
252  assert( 0 <= b );
253  assert( 0 <= c );
254  assert( b < mbr.mextent );
255  assert( c < mcr.mextent );
256 
257  return ConstMatrixView(mdata +
258  mbr.mstart + b * mbr.mstride +
259  mcr.mstart + c * mcr.mstride,
260  mpr, mrr,
261  p, r);
262 }
263 
266 ConstMatrixView ConstTensor4View::operator()(Index b,
267  Index p,
268  const Range& r,
269  const Range& c) const
270 {
271  // Check that b and p are valid:
272  assert( 0 <= b );
273  assert( 0 <= p );
274  assert( b < mbr.mextent );
275  assert( p < mpr.mextent );
276 
277  return ConstMatrixView(mdata +
278  mbr.mstart + b * mbr.mstride +
279  mpr.mstart + p * mpr.mstride,
280  mrr, mcr,
281  r, c);
282 }
283 
286 ConstVectorView ConstTensor4View::operator()(const Range& b,
287  Index p,
288  Index r,
289  Index c ) const
290 {
291  // Check that p, r and c are valid:
292  assert( 0 <= p );
293  assert( 0 <= r );
294  assert( 0 <= c );
295  assert( p < mpr.mextent );
296  assert( r < mrr.mextent );
297  assert( c < mcr.mextent );
298 
299  return ConstVectorView(mdata +
300  mpr.mstart + p * mpr.mstride +
301  mrr.mstart + r * mrr.mstride +
302  mcr.mstart + c * mcr.mstride,
303  mbr,
304  b);
305 }
306 
309 ConstVectorView ConstTensor4View::operator()(Index b,
310  const Range& p,
311  Index r,
312  Index c ) const
313 {
314  // Check that b, r and c are valid:
315  assert( 0 <= b );
316  assert( 0 <= r );
317  assert( 0 <= c );
318  assert( b < mbr.mextent );
319  assert( r < mrr.mextent );
320  assert( c < mcr.mextent );
321 
322  return ConstVectorView(mdata +
323  mbr.mstart + b * mbr.mstride +
324  mrr.mstart + r * mrr.mstride +
325  mcr.mstart + c * mcr.mstride,
326  mpr,
327  p);
328 }
329 
332 ConstVectorView ConstTensor4View::operator()(Index b,
333  Index p,
334  const Range& r,
335  Index c ) const
336 {
337  // Check that b, p and c are valid:
338  assert( 0 <= b );
339  assert( 0 <= p );
340  assert( 0 <= c );
341  assert( b < mbr.mextent );
342  assert( p < mpr.mextent );
343  assert( c < mcr.mextent );
344 
345  return ConstVectorView(mdata +
346  mbr.mstart + b * mbr.mstride +
347  mpr.mstart + p * mpr.mstride +
348  mcr.mstart + c * mcr.mstride,
349  mrr,
350  r);
351 }
352 
355 ConstVectorView ConstTensor4View::operator()(Index b,
356  Index p,
357  Index r,
358  const Range& c) const
359 {
360  // Check that b, p and r are valid:
361  assert( 0 <= b );
362  assert( 0 <= p );
363  assert( 0 <= r );
364  assert( b < mbr.mextent );
365  assert( p < mpr.mextent );
366  assert( r < mrr.mextent );
367 
368  return ConstVectorView(mdata +
369  mbr.mstart + b * mbr.mstride +
370  mpr.mstart + p * mpr.mstride +
371  mrr.mstart + r * mrr.mstride,
372  mcr,
373  c);
374 }
375 
382 Numeric *Tensor4View::get_c_array()
383 {
384  if (mbr.mstart != 0 || mbr.mstride != mpr.mextent * mrr.mextent * mcr.mextent
385  || mpr.mstart != 0 || mpr.mstride != mrr.mextent * mcr.mextent
386  || mrr.mstart != 0 || mrr.mstride != mcr.mextent
387  || mcr.mstart != 0 || mcr.mstride != 1)
388  throw runtime_error("A Tensor4View can only be converted to a plain C-array if mbr.mstart == 0 and mbr.mstride == mrp.extent*mrr.extent*mcr.extent and mpr.mstart == 0 and mpr.mstride == mrr.extent*mcr.extent and mrr.mstart == 0 and mrr.mstride == mcr.extent and mcr.mstart == 0 and mcr.mstride == 1");
389 
390  return mdata;
391 }
392 
399 const Numeric *Tensor4View::get_c_array() const
400 {
401  if (mbr.mstart != 0 || mbr.mstride != mpr.mextent * mrr.mextent * mcr.mextent
402  || mpr.mstart != 0 || mpr.mstride != mrr.mextent * mcr.mextent
403  || mrr.mstart != 0 || mrr.mstride != mcr.mextent
404  || mcr.mstart != 0 || mcr.mstride != 1)
405  throw runtime_error("A Tensor4View can only be converted to a plain C-array if mbr.mstart == 0 and mbr.mstride == mrp.extent*mrr.extent*mcr.extent and mpr.mstart == 0 and mpr.mstride == mrr.extent*mcr.extent and mrr.mstart == 0 and mrr.mstride == mcr.extent and mcr.mstart == 0 and mcr.mstride == 1");
406 
407  return mdata;
408 }
409 
411 ConstIterator4D ConstTensor4View::begin() const
412 {
413  return ConstIterator4D( ConstTensor3View(mdata + mbr.mstart,
414  mpr, mrr, mcr),
415  mbr.mstride );
416 }
417 
419 ConstIterator4D ConstTensor4View::end() const
420 {
421  return ConstIterator4D( ConstTensor3View(mdata + mbr.mstart +
422  (mbr.mextent) * mbr.mstride,
423  mpr, mrr, mcr),
424  mbr.mstride );
425 }
426 
428 ConstTensor4View::ConstTensor4View(const ConstTensor3View& a) :
429  mbr(0,1,a.mpr.mextent*a.mrr.mextent*a.mcr.mextent),
430  mpr(a.mpr),
431  mrr(a.mrr),
432  mcr(a.mcr),
433  mdata(a.mdata)
434 {
435  // Nothing to do here.
436 }
437 
440 ConstTensor4View::ConstTensor4View() :
441  mbr(0,0,1),
442  mpr(0,0,1),
443  mrr(0,0,1),
444  mcr(0,0,1),
445  mdata(NULL)
446 {
447  // Nothing to do here.
448 }
449 
454 ConstTensor4View::ConstTensor4View(Numeric *data,
455  const Range& br,
456  const Range& pr,
457  const Range& rr,
458  const Range& cr) :
459  mbr(br),
460  mpr(pr),
461  mrr(rr),
462  mcr(cr),
463  mdata(data)
464 {
465  // Nothing to do here.
466 }
467 
475 ConstTensor4View::ConstTensor4View(Numeric *data,
476  const Range& pb,
477  const Range& pp,
478  const Range& pr,
479  const Range& pc,
480  const Range& nb,
481  const Range& np,
482  const Range& nr,
483  const Range& nc) :
484  mbr(pb,nb),
485  mpr(pp,np),
486  mrr(pr,nr),
487  mcr(pc,nc),
488  mdata(data)
489 {
490  // Nothing to do here.
491 }
492 
496 ostream& operator<<(ostream& os, const ConstTensor4View& v)
497 {
498  // Page iterators:
499  ConstIterator4D ib = v.begin();
500  const ConstIterator4D end_book = v.end();
501 
502  if ( ib != end_book ) {
503  os << *ib;
504  ++ib;
505  }
506 
507  for ( ; ib != end_book; ++ib ) {
508  os << "\n\n";
509  os << *ib;
510  }
511 
512  return os;
513 }
514 
515 
516 // Functions for Tensor4View:
517 // -------------------------
518 
523 ConstTensor4View Tensor4View::operator()(const Range& b,
524  const Range& p,
525  const Range& r,
526  const Range& c) const
527 {
528  return ConstTensor4View::operator()(b,p,r,c);
529 }
530 
535 ConstTensor3View Tensor4View::operator()(const Range& b,
536  const Range& p,
537  const Range& r,
538  Index c ) const
539 {
540  return ConstTensor4View::operator()(b,p,r,c);
541 }
542 
547 ConstTensor3View Tensor4View::operator()(const Range& b,
548  const Range& p,
549  Index r,
550  const Range& c) const
551 {
552  return ConstTensor4View::operator()(b,p,r,c);
553 }
554 
559 ConstTensor3View Tensor4View::operator()(const Range& b,
560  Index p,
561  const Range& r,
562  const Range& c) const
563 {
564  return ConstTensor4View::operator()(b,p,r,c);
565 }
566 
571 ConstTensor3View Tensor4View::operator()(Index b,
572  const Range& p,
573  const Range& r,
574  const Range& c) const
575 {
576  return ConstTensor4View::operator()(b,p,r,c);
577 }
578 
583 ConstMatrixView Tensor4View::operator()(const Range& b,
584  const Range& p,
585  Index r,
586  Index c ) const
587 {
588  return ConstTensor4View::operator()(b,p,r,c);
589 }
590 
595 ConstMatrixView Tensor4View::operator()(const Range& b,
596  Index p,
597  const Range& r,
598  Index c ) const
599 {
600  return ConstTensor4View::operator()(b,p,r,c);
601 }
602 
607 ConstMatrixView Tensor4View::operator()(const Range& b,
608  Index p,
609  Index r,
610  const Range& c) const
611 {
612  return ConstTensor4View::operator()(b,p,r,c);
613 }
614 
619 ConstMatrixView Tensor4View::operator()(Index b,
620  const Range& p,
621  Index r,
622  const Range& c) const
623 {
624  return ConstTensor4View::operator()(b,p,r,c);
625 }
626 
631 ConstMatrixView Tensor4View::operator()(Index b,
632  const Range& p,
633  const Range& r,
634  Index c ) const
635 {
636  return ConstTensor4View::operator()(b,p,r,c);
637 }
638 
643 ConstMatrixView Tensor4View::operator()(Index b,
644  Index p,
645  const Range& r,
646  const Range& c) const
647 {
648  return ConstTensor4View::operator()(b,p,r,c);
649 }
650 
655 ConstVectorView Tensor4View::operator()(const Range& b,
656  Index p,
657  Index r,
658  Index c ) const
659 {
660  return ConstTensor4View::operator()(b,p,r,c);
661 }
662 
667 ConstVectorView Tensor4View::operator()(Index b,
668  const Range& p,
669  Index r,
670  Index c ) const
671 {
672  return ConstTensor4View::operator()(b,p,r,c);
673 }
674 
679 ConstVectorView Tensor4View::operator()(Index b,
680  Index p,
681  const Range& r,
682  Index c ) const
683 {
684  return ConstTensor4View::operator()(b,p,r,c);
685 }
686 
691 ConstVectorView Tensor4View::operator()(Index b,
692  Index p,
693  Index r,
694  const Range& c) const
695 {
696  return ConstTensor4View::operator()(b,p,r,c);
697 }
698 
702 Tensor4View Tensor4View::operator()(const Range& b,
703  const Range& p,
704  const Range& r,
705  const Range& c)
706 {
707  return Tensor4View(mdata,
708  mbr, mpr, mrr, mcr,
709  b, p, r, c);
710 }
711 
714 Tensor3View Tensor4View::operator()(const Range& b,
715  const Range& p,
716  const Range& r,
717  Index c)
718 {
719  // Check that c is valid:
720  assert( 0 <= c );
721  assert( c < mcr.mextent );
722 
723  return Tensor3View(mdata +
724  mcr.mstart + c * mcr.mstride,
725  mbr, mpr, mrr,
726  b, p, r);
727 }
728 
731 Tensor3View Tensor4View::operator()(const Range& b,
732  const Range& p,
733  Index r,
734  const Range& c)
735 {
736  // Check that r is valid:
737  assert( 0 <= r );
738  assert( r < mrr.mextent );
739 
740  return Tensor3View(mdata +
741  mrr.mstart + r * mrr.mstride,
742  mbr, mpr, mcr,
743  b, p, c);
744 }
745 
748 Tensor3View Tensor4View::operator()(const Range& b,
749  Index p,
750  const Range& r,
751  const Range& c)
752 {
753  // Check that p is valid:
754  assert( 0 <= p );
755  assert( p < mpr.mextent );
756 
757  return Tensor3View(mdata +
758  mpr.mstart + p * mpr.mstride,
759  mbr, mrr, mcr,
760  b, r, c);
761 }
762 
765 Tensor3View Tensor4View::operator()(Index b,
766  const Range& p,
767  const Range& r,
768  const Range& c)
769 {
770  // Check that b is valid:
771  assert( 0 <= b );
772  assert( b < mbr.mextent );
773 
774  return Tensor3View(mdata +
775  mbr.mstart + b * mbr.mstride,
776  mpr, mrr, mcr,
777  p, r, c);
778 }
779 
782 MatrixView Tensor4View::operator()(const Range& b,
783  const Range& p,
784  Index r,
785  Index c)
786 {
787  // Check that r and c are valid:
788  assert( 0 <= r );
789  assert( 0 <= c );
790  assert( r < mrr.mextent );
791  assert( c < mcr.mextent );
792 
793  return MatrixView(mdata +
794  mrr.mstart + r * mrr.mstride +
795  mcr.mstart + c * mcr.mstride,
796  mbr, mpr,
797  b, p);
798 }
799 
802 MatrixView Tensor4View::operator()(const Range& b,
803  Index p,
804  const Range& r,
805  Index c)
806 {
807  // Check that p and c are valid:
808  assert( 0 <= p );
809  assert( 0 <= c );
810  assert( p < mpr.mextent );
811  assert( c < mcr.mextent );
812 
813  return MatrixView(mdata +
814  mpr.mstart + p * mpr.mstride +
815  mcr.mstart + c * mcr.mstride,
816  mbr, mrr,
817  b, r);
818 }
819 
822 MatrixView Tensor4View::operator()(const Range& b,
823  Index p,
824  Index r,
825  const Range& c)
826 {
827  // Check that p and r are valid:
828  assert( 0 <= p );
829  assert( 0 <= r );
830  assert( p < mpr.mextent );
831  assert( r < mrr.mextent );
832 
833  return MatrixView(mdata +
834  mpr.mstart + p * mpr.mstride +
835  mrr.mstart + r * mrr.mstride,
836  mbr, mcr,
837  b, c);
838 }
839 
842 MatrixView Tensor4View::operator()(Index b,
843  const Range& p,
844  Index r,
845  const Range& c)
846 {
847  // Check that b and r are valid:
848  assert( 0 <= b );
849  assert( 0 <= r );
850  assert( b < mbr.mextent );
851  assert( r < mrr.mextent );
852 
853  return MatrixView(mdata +
854  mbr.mstart + b * mbr.mstride +
855  mrr.mstart + r * mrr.mstride,
856  mpr, mcr,
857  p, c);
858 }
859 
862 MatrixView Tensor4View::operator()(Index b,
863  const Range& p,
864  const Range& r,
865  Index c)
866 {
867  // Check that b and c are valid:
868  assert( 0 <= b );
869  assert( 0 <= c );
870  assert( b < mbr.mextent );
871  assert( c < mcr.mextent );
872 
873  return MatrixView(mdata +
874  mbr.mstart + b * mbr.mstride +
875  mcr.mstart + c * mcr.mstride,
876  mpr, mrr,
877  p, r);
878 }
879 
882 MatrixView Tensor4View::operator()(Index b,
883  Index p,
884  const Range& r,
885  const Range& c)
886 {
887  // Check that b and p are valid:
888  assert( 0 <= b );
889  assert( 0 <= p );
890  assert( b < mbr.mextent );
891  assert( p < mpr.mextent );
892 
893  return MatrixView(mdata +
894  mbr.mstart + b * mbr.mstride +
895  mpr.mstart + p * mpr.mstride,
896  mrr, mcr,
897  r, c);
898 }
899 
902 VectorView Tensor4View::operator()(const Range& b,
903  Index p,
904  Index r,
905  Index c)
906 {
907  // Check that p, r and c are valid:
908  assert( 0 <= p );
909  assert( 0 <= r );
910  assert( 0 <= c );
911  assert( p < mpr.mextent );
912  assert( r < mrr.mextent );
913  assert( c < mcr.mextent );
914 
915  return VectorView(mdata +
916  mpr.mstart + p * mpr.mstride +
917  mrr.mstart + r * mrr.mstride +
918  mcr.mstart + c * mcr.mstride,
919  mbr,
920  b);
921 }
922 
925 VectorView Tensor4View::operator()(Index b,
926  const Range& p,
927  Index r,
928  Index c)
929 {
930  // Check that b, r and c are valid:
931  assert( 0 <= b );
932  assert( 0 <= r );
933  assert( 0 <= c );
934  assert( b < mbr.mextent );
935  assert( r < mrr.mextent );
936  assert( c < mcr.mextent );
937 
938  return VectorView(mdata +
939  mbr.mstart + b * mbr.mstride +
940  mrr.mstart + r * mrr.mstride +
941  mcr.mstart + c * mcr.mstride,
942  mpr,
943  p);
944 }
945 
948 VectorView Tensor4View::operator()(Index b,
949  Index p,
950  const Range& r,
951  Index c)
952 {
953  // Check that b, p and c are valid:
954  assert( 0 <= b );
955  assert( 0 <= p );
956  assert( 0 <= c );
957  assert( b < mbr.mextent );
958  assert( p < mpr.mextent );
959  assert( c < mcr.mextent );
960 
961  return VectorView(mdata +
962  mbr.mstart + b * mbr.mstride +
963  mpr.mstart + p * mpr.mstride +
964  mcr.mstart + c * mcr.mstride,
965  mrr,
966  r);
967 }
968 
971 VectorView Tensor4View::operator()(Index b,
972  Index p,
973  Index r,
974  const Range& c)
975 {
976  // Check that b, p and r are valid:
977  assert( 0 <= b );
978  assert( 0 <= p );
979  assert( 0 <= r );
980  assert( b < mbr.mextent );
981  assert( p < mpr.mextent );
982  assert( r < mrr.mextent );
983 
984  return VectorView(mdata +
985  mbr.mstart + b * mbr.mstride +
986  mpr.mstart + p * mpr.mstride +
987  mrr.mstart + r * mrr.mstride,
988  mcr,
989  c);
990 }
991 
994 ConstIterator4D Tensor4View::begin() const
995 {
996  return ConstTensor4View::begin();
997 }
998 
1000 ConstIterator4D Tensor4View::end() const
1001 {
1002  return ConstTensor4View::end();
1003 }
1004 
1006 Iterator4D Tensor4View::begin()
1007 {
1008  return Iterator4D( Tensor3View(mdata + mbr.mstart,
1009  mpr, mrr, mcr),
1010  mbr.mstride );
1011 }
1012 
1014 Iterator4D Tensor4View::end()
1015 {
1016  return Iterator4D( Tensor3View(mdata + mbr.mstart +
1017  (mbr.mextent) * mbr.mstride,
1018  mpr, mrr, mcr),
1019  mbr.mstride );
1020 }
1021 
1026 Tensor4View& Tensor4View::operator=(const ConstTensor4View& m)
1027 {
1028  // Check that sizes are compatible:
1029  assert( mbr.mextent == m.mbr.mextent );
1030  assert( mpr.mextent == m.mpr.mextent );
1031  assert( mrr.mextent == m.mrr.mextent );
1032  assert( mcr.mextent == m.mcr.mextent );
1033 
1034  copy( m.begin(), m.end(), begin() );
1035  return *this;
1036 }
1037 
1043 Tensor4View& Tensor4View::operator=(const Tensor4View& m)
1044 {
1045  // Check that sizes are compatible:
1046  assert( mbr.mextent == m.mbr.mextent );
1047  assert( mpr.mextent == m.mpr.mextent );
1048  assert( mrr.mextent == m.mrr.mextent );
1049  assert( mcr.mextent == m.mcr.mextent );
1050 
1051  copy( m.begin(), m.end(), begin() );
1052  return *this;
1053 }
1054 
1058 Tensor4View& Tensor4View::operator=(const Tensor4& m)
1059 {
1060  // Check that sizes are compatible:
1061  assert( mbr.mextent == m.mbr.mextent );
1062  assert( mpr.mextent == m.mpr.mextent );
1063  assert( mrr.mextent == m.mrr.mextent );
1064  assert( mcr.mextent == m.mcr.mextent );
1065 
1066  copy( m.begin(), m.end(), begin() );
1067  return *this;
1068 }
1069 
1072 Tensor4View& Tensor4View::operator=(Numeric x)
1073 {
1074  copy( x, begin(), end() );
1075  return *this;
1076 }
1077 
1078 // Some little helper functions:
1079 //------------------------------
1080 
1082 Tensor4View& Tensor4View::operator*=(Numeric x)
1083 {
1084  const Iterator4D eb = end();
1085  for ( Iterator4D b = begin(); b != eb ; ++b )
1086  {
1087  *b *= x;
1088  }
1089  return *this;
1090 }
1091 
1093 Tensor4View& Tensor4View::operator/=(Numeric x)
1094 {
1095  const Iterator4D eb = end();
1096  for ( Iterator4D b = begin(); b != eb ; ++b )
1097  {
1098  *b /= x;
1099  }
1100  return *this;
1101 }
1102 
1104 Tensor4View& Tensor4View::operator+=(Numeric x)
1105 {
1106  const Iterator4D eb = end();
1107  for ( Iterator4D b = begin(); b != eb ; ++b )
1108  {
1109  *b += x;
1110  }
1111  return *this;
1112 }
1113 
1115 Tensor4View& Tensor4View::operator-=(Numeric x)
1116 {
1117  const Iterator4D eb = end();
1118  for ( Iterator4D b = begin(); b != eb ; ++b )
1119  {
1120  *b -= x;
1121  }
1122  return *this;
1123 }
1124 
1126 Tensor4View& Tensor4View::operator*=(const ConstTensor4View& x)
1127 {
1128  assert( nbooks() == x.nbooks() );
1129  assert( npages() == x.npages() );
1130  assert( nrows() == x.nrows() );
1131  assert( ncols() == x.ncols() );
1132  ConstIterator4D xb = x.begin();
1133  Iterator4D b = begin();
1134  const Iterator4D eb = end();
1135  for ( ; b != eb ; ++b, ++xb )
1136  {
1137  *b *= *xb;
1138  }
1139  return *this;
1140 }
1141 
1143 Tensor4View& Tensor4View::operator/=(const ConstTensor4View& x)
1144 {
1145  assert( nbooks() == x.nbooks() );
1146  assert( npages() == x.npages() );
1147  assert( nrows() == x.nrows() );
1148  assert( ncols() == x.ncols() );
1149  ConstIterator4D xb = x.begin();
1150  Iterator4D b = begin();
1151  const Iterator4D eb = end();
1152  for ( ; b != eb ; ++b, ++xb )
1153  {
1154  *b /= *xb;
1155  }
1156  return *this;
1157 }
1158 
1160 Tensor4View& Tensor4View::operator+=(const ConstTensor4View& x)
1161 {
1162  assert( nbooks() == x.nbooks() );
1163  assert( npages() == x.npages() );
1164  assert( nrows() == x.nrows() );
1165  assert( ncols() == x.ncols() );
1166  ConstIterator4D xb = x.begin();
1167  Iterator4D b = begin();
1168  const Iterator4D eb = end();
1169  for ( ; b != eb ; ++b, ++xb )
1170  {
1171  *b += *xb;
1172  }
1173  return *this;
1174 }
1175 
1177 Tensor4View& Tensor4View::operator-=(const ConstTensor4View& x)
1178 {
1179  assert( nbooks() == x.nbooks() );
1180  assert( npages() == x.npages() );
1181  assert( nrows() == x.nrows() );
1182  assert( ncols() == x.ncols() );
1183  ConstIterator4D xb = x.begin();
1184  Iterator4D b = begin();
1185  const Iterator4D eb = end();
1186  for ( ; b != eb ; ++b, ++xb )
1187  {
1188  *b -= *xb;
1189  }
1190  return *this;
1191 }
1192 
1194 Tensor4View::Tensor4View(const Tensor3View& a) :
1195  ConstTensor4View( a.mdata,
1196  Range(0,1,
1197  a.mpr.mextent*
1198  a.mrr.mextent*
1199  a.mcr.mextent),
1200  a.mpr,
1201  a.mrr,
1202  a.mcr )
1203 {
1204  // Nothing to do here.
1205 }
1206 
1209 Tensor4View::Tensor4View() :
1210  ConstTensor4View()
1211 {
1212  // Nothing to do here.
1213 }
1214 
1218 Tensor4View::Tensor4View(Numeric *data,
1219  const Range& br,
1220  const Range& pr,
1221  const Range& rr,
1222  const Range& cr) :
1223  ConstTensor4View(data, br, pr, rr, cr)
1224 {
1225  // Nothing to do here.
1226 }
1227 
1246 Tensor4View::Tensor4View(Numeric *data,
1247  const Range& pb,
1248  const Range& pp,
1249  const Range& pr,
1250  const Range& pc,
1251  const Range& nb,
1252  const Range& np,
1253  const Range& nr,
1254  const Range& nc) :
1255  ConstTensor4View(data, pb, pp, pr, pc, nb, np, nr, nc)
1256 {
1257  // Nothing to do here.
1258 }
1259 
1264 void copy(ConstIterator4D origin,
1265  const ConstIterator4D& end,
1266  Iterator4D target)
1267 {
1268  for ( ; origin != end ; ++origin, ++target )
1269  {
1270  // We use the copy function for the next smaller rank of tensor
1271  // recursively:
1272  copy( origin->begin(), origin->end(), target->begin() );
1273  }
1274 }
1275 
1277 void copy(Numeric x,
1278  Iterator4D target,
1279  const Iterator4D& end)
1280 {
1281  for ( ; target != end ; ++target )
1282  {
1283  // We use the copy function for the next smaller rank of tensor
1284  // recursively:
1285  copy( x, target->begin(), target->end() );
1286  }
1287 }
1288 
1289 
1290 // Functions for Tensor4:
1291 // ---------------------
1292 
1294 Tensor4::Tensor4() :
1295  Tensor4View::Tensor4View()
1296 {
1297  // Nothing to do here. However, note that the default constructor
1298  // for Tensor4View has been called in the initializer list. That is
1299  // crucial, otherwise internal range objects will not be properly
1300  // initialized.
1301 }
1302 
1305 Tensor4::Tensor4(Index b, Index p, Index r, Index c) :
1306  Tensor4View( new Numeric[b*p*r*c],
1307  Range( 0, b, p*r*c ),
1308  Range( 0, p, r*c ),
1309  Range( 0, r, c ),
1310  Range( 0, c) )
1311 {
1312  // Nothing to do here.
1313 }
1314 
1316 Tensor4::Tensor4(Index b, Index p, Index r, Index c, Numeric fill) :
1317  Tensor4View( new Numeric[b*p*r*c],
1318  Range( 0, b, p*r*c ),
1319  Range( 0, p, r*c ),
1320  Range( 0, r, c ),
1321  Range( 0, c) )
1322 {
1323  // Here we can access the raw memory directly, for slightly
1324  // increased efficiency:
1325  const Numeric *stop = mdata + b*p*r*c;
1326  for ( Numeric *x = mdata; x < stop; ++x )
1327  *x = fill;
1328 }
1329 
1332 Tensor4::Tensor4(const ConstTensor4View& m) :
1333  Tensor4View( new Numeric[m.nbooks()*m.npages()*m.nrows()*m.ncols()],
1334  Range( 0, m.nbooks(), m.npages()*m.nrows()*m.ncols() ),
1335  Range( 0, m.npages(), m.nrows()*m.ncols() ),
1336  Range( 0, m.nrows(), m.ncols() ),
1337  Range( 0, m.ncols() ) )
1338 {
1339  copy( m.begin(), m.end(), begin() );
1340 }
1341 
1344 Tensor4::Tensor4(const Tensor4& m) :
1345  Tensor4View( new Numeric[m.nbooks()*m.npages()*m.nrows()*m.ncols()],
1346  Range( 0, m.nbooks(), m.npages()*m.nrows()*m.ncols() ),
1347  Range( 0, m.npages(), m.nrows()*m.ncols() ),
1348  Range( 0, m.nrows(), m.ncols() ),
1349  Range( 0, m.ncols() ) )
1350 {
1351  // There is a catch here: If m is an empty tensor, then it will have
1352  // dimensions of size 0. But these are used to initialize the stride
1353  // for higher dimensions! Thus, this method has to be consistent
1354  // with the behaviour of Range::Range. For now, Range::Range allows
1355  // also stride 0.
1356  copy( m.begin(), m.end(), begin());
1357 }
1358 
1360 
1383 Tensor4& Tensor4::operator=(const Tensor4& m)
1384 {
1385  // cout << "Tensor4 copy: m = " << m.nrows() << " " << m.ncols() << "\n";
1386  // cout << " n = " << nrows() << " " << ncols() << "\n";
1387 
1388  resize( m.mbr.mextent, m.mpr.mextent, m.mrr.mextent, m.mcr.mextent );
1389  copy( m.begin(), m.end(), begin() );
1390  return *this;
1391 }
1392 
1395 Tensor4& Tensor4::operator=(Numeric x)
1396 {
1397  copy( x, begin(), end() );
1398  return *this;
1399 }
1400 
1404 void Tensor4::resize(Index b, Index p, Index r, Index c)
1405 {
1406  assert( 0 <= b );
1407  assert( 0 <= p );
1408  assert( 0 <= r );
1409  assert( 0 <= c );
1410 
1411  if ( mbr.mextent != b ||
1412  mpr.mextent != p ||
1413  mrr.mextent != r ||
1414  mcr.mextent != c )
1415  {
1416  delete[] mdata;
1417  mdata = new Numeric[b*p*r*c];
1418 
1419  mbr.mstart = 0;
1420  mbr.mextent = b;
1421  mbr.mstride = p*r*c;
1422 
1423  mpr.mstart = 0;
1424  mpr.mextent = p;
1425  mpr.mstride = r*c;
1426 
1427  mrr.mstart = 0;
1428  mrr.mextent = r;
1429  mrr.mstride = c;
1430 
1431  mcr.mstart = 0;
1432  mcr.mextent = c;
1433  mcr.mstride = 1;
1434  }
1435 }
1436 
1439 Tensor4::~Tensor4()
1440 {
1441 // cout << "Destroying a Tensor4:\n"
1442 // << *this << "\n........................................\n";
1443  delete[] mdata;
1444 }
1445 
1446 
1462 void transform( Tensor4View y,
1463  double (&my_func)(double),
1464  ConstTensor4View x )
1465 {
1466  // Check dimensions:
1467  assert( y.nbooks() == x.nbooks() );
1468  assert( y.npages() == x.npages() );
1469  assert( y.nrows() == x.nrows() );
1470  assert( y.ncols() == x.ncols() );
1471 
1472  const ConstIterator4D xe = x.end();
1473  ConstIterator4D xi = x.begin();
1474  Iterator4D yi = y.begin();
1475  for ( ; xi != xe; ++xi, ++yi )
1476  {
1477  // Use the transform function of lower dimensional tensors
1478  // recursively:
1479  transform( *yi, my_func, *xi );
1480  }
1481 }
1482 
1484 Numeric max(const ConstTensor4View& x)
1485 {
1486  const ConstIterator4D xe = x.end();
1487  ConstIterator4D xi = x.begin();
1488 
1489  // Initial value for max:
1490  Numeric themax = max(*xi);
1491  ++xi;
1492 
1493  for ( ; xi != xe ; ++xi )
1494  {
1495  // Use the max function of lower dimensional tensors
1496  // recursively:
1497  Numeric maxi = max(*xi);
1498  if ( maxi > themax )
1499  themax = maxi;
1500  }
1501 
1502  return themax;
1503 }
1504 
1506 Numeric min(const ConstTensor4View& x)
1507 {
1508  const ConstIterator4D xe = x.end();
1509  ConstIterator4D xi = x.begin();
1510 
1511  // Initial value for min:
1512  Numeric themin = min(*xi);
1513  ++xi;
1514 
1515  for ( ; xi != xe ; ++xi )
1516  {
1517  // Use the min function of lower dimensional tensors
1518  // recursively:
1519  Numeric mini = min(*xi);
1520  if ( mini < themin )
1521  themin = mini;
1522  }
1523 
1524  return themin;
1525 }
1526 
1527 
1529 // Helper function for debugging
1530 #ifndef NDEBUG
1531 
1547 Numeric debug_tensor4view_get_elem (Tensor4View& tv,
1548  Index b, Index p, Index r, Index c)
1549 {
1550  return tv(b, p, r, c);
1551 }
1552 
1553 #endif
1554 
Tensor3View::end
ConstIterator3D end() const
Return const iterator behind last row.
Definition: matpackIII.cc:470
Tensor4::resize
void resize(Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackIV.cc:1404
ConstTensor4View::mpr
Range mpr
The page range of mdata that is actually used.
Definition: matpackIV.h:227
MatrixView
The MatrixView class.
Definition: matpackI.h:668
exceptions.h
The declarations of all the exception classes.
ConstTensor4View::end
ConstIterator4D end() const
Return const iterator behind last book.
Definition: matpackIV.cc:419
Tensor4View::Tensor4View
Tensor4View()
Default constructor.
Definition: matpackIV.cc:1209
min
Numeric min(const ConstTensor4View &x)
Min function, tensor version.
Definition: matpackIV.cc:1506
transform
void transform(Tensor4View y, double(&my_func)(double), ConstTensor4View x)
A generic transform function for tensors, which can be used to implement mathematical functions opera...
Definition: matpackIV.cc:1462
Tensor4View::operator*=
Tensor4View & operator*=(Numeric x)
Multiplication by scalar.
Definition: matpackIV.cc:1082
ConstTensor4View::mbr
Range mbr
The book range of mdata that is actually used.
Definition: matpackIV.h:225
ConstTensor4View::operator()
ConstTensor4View operator()(const Range &b, const Range &p, const Range &r, const Range &c) const
Const index operator for subrange.
Definition: matpackIV.cc:86
ConstTensor4View::mcr
Range mcr
The column range of mdata that is actually used.
Definition: matpackIV.h:231
Tensor4::~Tensor4
virtual ~Tensor4()
Destructor for Tensor4.
Definition: matpackIV.cc:1439
Tensor4::operator=
Tensor4 & operator=(const Tensor4 &x)
Assignment operator from another tensor.
Definition: matpackIV.cc:1383
Range::mstart
Index mstart
The start index.
Definition: matpackI.h:204
Tensor4
The Tensor4 class.
Definition: matpackIV.h:375
Tensor4::Tensor4
Tensor4()
Default constructor.
Definition: matpackIV.cc:1294
matpackIV.h
Iterator4D::operator->
Tensor3View * operator->()
The -> operator is needed, so that we can write i->begin() to get the 3D iterators.
Definition: matpackIV.cc:30
ConstIterator4D
Const version of Iterator4D.
Definition: matpackIV.h:84
ConstTensor4View
A constant view of a Tensor4.
Definition: matpackIV.h:149
Tensor3View
The Tensor3View class.
Definition: matpackIII.h:234
Tensor4View::operator+=
Tensor4View & operator+=(Numeric x)
Addition of scalar.
Definition: matpackIV.cc:1104
ConstIterator4D::operator->
const ConstTensor3View * operator->() const
The -> operator is needed, so that we can write i->begin() to get the 3D iterators.
Definition: matpackIV.cc:43
ConstTensor3View::begin
ConstIterator3D begin() const
Return const iterator to first page.
Definition: matpackIII.cc:141
eb
#define eb
Definition: continua.cc:14566
debug_tensor4view_get_elem
Numeric debug_tensor4view_get_elem(Tensor4View &tv, Index b, Index p, Index r, Index c)
Helper function to access tensor elements.
Definition: matpackIV.cc:1547
VectorView
The VectorView class.
Definition: matpackI.h:373
ConstTensor4View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:78
ConstTensor4View::begin
ConstIterator4D begin() const
Return const iterator to first book.
Definition: matpackIV.cc:411
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
ConstTensor4View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:66
ConstTensor4View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:60
Tensor4View::operator=
Tensor4View & operator=(const ConstTensor4View &v)
Assignment operator.
Definition: matpackIV.cc:1026
Tensor4View::end
ConstIterator4D end() const
Return const iterator behind last book.
Definition: matpackIV.cc:1000
Iterator4D::msv
Tensor3View msv
Current position.
Definition: matpackIV.h:78
Tensor4View::operator/=
Tensor4View & operator/=(Numeric x)
Division by scalar.
Definition: matpackIV.cc:1093
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:591
Tensor4View::begin
ConstIterator4D begin() const
Return const iterator to first book.
Definition: matpackIV.cc:994
ConstTensor4View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackIV.cc:72
ConstTensor3View::end
ConstIterator3D end() const
Return const iterator behind last page.
Definition: matpackIII.cc:150
ConstIterator4D::operator*
const ConstTensor3View & operator*() const
Dereferencing.
Definition: matpackIV.cc:49
Range
The range class.
Definition: matpackI.h:148
Range::mextent
Index mextent
The number of elements.
Definition: matpackI.h:207
Iterator4D::operator*
Tensor3View & operator*()
Dereferencing.
Definition: matpackIV.cc:36
Iterator4D
Implementation of Tensors of Rank 4.
Definition: matpackIV.h:38
Tensor3View::begin
ConstIterator3D begin() const
Return const iterator to first row.
Definition: matpackIII.cc:464
copy
void copy(ConstIterator4D origin, const ConstIterator4D &end, Iterator4D target)
Copy data between begin and end to target.
Definition: matpackIV.cc:1264
Tensor4View::operator()
ConstTensor4View operator()(const Range &b, const Range &p, const Range &r, const Range &c) const
Const index operator for subrange.
Definition: matpackIV.cc:523
ConstTensor4View::mrr
Range mrr
The row range of mdata that is actually used.
Definition: matpackIV.h:229
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:147
Tensor4View
The Tensor4View class.
Definition: matpackIV.h:245
ConstIterator4D::msv
ConstTensor3View msv
Current position.
Definition: matpackIV.h:127
Tensor4View::operator-=
Tensor4View & operator-=(Numeric x)
Subtraction of scalar.
Definition: matpackIV.cc:1115
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
operator<<
ostream & operator<<(ostream &os, const ConstTensor4View &v)
Output operator.
Definition: matpackIV.cc:496
max
Numeric max(const ConstTensor4View &x)
Max function, tensor version.
Definition: matpackIV.cc:1484
Range::mstride
Index mstride
The stride.
Definition: matpackI.h:209
ConstTensor4View::ConstTensor4View
ConstTensor4View()
Default constructor.
Definition: matpackIV.cc:440
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:300
Tensor4View::get_c_array
const Numeric * get_c_array() const
Conversion to plain C-array.
Definition: matpackIV.cc:399
ConstTensor4View::mdata
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackIV.h:233