ARTS  2.2.66
matpackI.cc
Go to the documentation of this file.
1 /* Copyright (C) 2001-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 <cstring>
26 #include <cmath>
27 #include "matpackI.h"
28 #include "exceptions.h"
29 
30 using std::setw;
31 using std::runtime_error;
32 
33 
34 // Define the global joker object:
35 extern const Joker joker = Joker();
36 
37 static const Numeric RAD2DEG = 57.295779513082323;
38 
39 // Functions for Range:
40 // --------------------
41 
42 /* Explicit constructor.
43 
44  \param Start must be >= 0.
45 
46  \param Extent also. Although internally negative extent means "to the end",
47  this can not be created this way, only with the joker. Zero
48  extent is allowed, though, which corresponds to an empty range.
49 
50  \param Stride can be anything. It can be omitted, in which case the
51  default value is 1. */
52 Range::Range(Index start, Index extent, Index stride) :
53  mstart(start), mextent(extent), mstride(stride)
54 {
55  // Start must be >= 0:
56  assert( 0<=mstart );
57  // Extent also. Although internally negative extent means "to the end",
58  // this can not be created this way, only with the joker. Zero
59  // extent is allowed, though, which corresponds to an empty range.
60  assert( 0<=mextent );
61  // Stride can be anything except 0.
62  // SAB 2001-09-21: Allow 0 stride.
63  // assert( 0!=mstride);
64 }
65 
68 Range::Range(Index start, Joker, Index stride) :
69  mstart(start), mextent(-1), mstride(stride)
70 {
71  // Start must be >= 0:
72  assert( 0<=mstart );
73 }
74 
79  mstart(0), mextent(-1), mstride(stride)
80 {
81  // Nothing to do here.
82 }
83 
89 Range::Range(Index max_size, const Range& r) :
90  mstart(r.mstart),
91  mextent(r.mextent),
92  mstride(r.mstride)
93 {
94  // Start must be >= 0:
95  assert( 0<=mstart );
96  // ... and < max_size:
97  assert( mstart<max_size );
98 
99  // Stride must be != 0:
100  assert( 0!=mstride);
101 
102  // Convert negative extent (joker) to explicit extent
103  if ( mextent<0 )
104  {
105  if ( 0<mstride )
106  mextent = 1 + (max_size-1-mstart)/mstride;
107  else
108  mextent = 1 + (0-mstart)/mstride;
109  }
110  else
111  {
112 #ifndef NDEBUG
113  // Check that extent is ok:
114  Index fin = mstart+(mextent-1)*mstride;
115  assert( 0 <= fin );
116  assert( fin < max_size );
117 #endif
118  }
119 }
120 
126 Range::Range(const Range& p, const Range& n) :
127  mstart(p.mstart + n.mstart*p.mstride),
128  mextent(n.mextent),
129  mstride(p.mstride*n.mstride)
130 {
131  // We have to juggle here a bit with previous, new, and resulting
132  // quantities. I.e.;
133  // p.mstride: Previous stride
134  // n.mstride: New stride (as specified)
135  // mstride: Resulting stride (old*new)
136 
137  // Get the previous final element:
138  Index prev_fin = p.mstart+(p.mextent-1)*p.mstride;
139 
140  // Resulting start must be >= previous start:
141  assert( p.mstart<=mstart );
142  // and <= prev_fin:
143  assert( mstart<=prev_fin );
144 
145  // Resulting stride must be != 0:
146  assert( 0!=mstride);
147 
148  // Convert negative extent (joker) to explicit extent
149  if ( mextent<0 )
150  {
151  if ( 0<mstride )
152  mextent = 1 + (prev_fin-mstart)/mstride;
153  else
154  mextent = 1 + (p.mstart-mstart)/mstride;
155  }
156  else
157  {
158 #ifndef NDEBUG
159  // Check that extent is ok:
160  Index fin = mstart+(mextent-1)*mstride;
161  assert( p.mstart <= fin );
162  assert( fin <= prev_fin );
163 #endif
164  }
165 
166 }
167 
168 // Functions for ConstVectorView:
169 // ------------------------------
170 
181 {
182  return mrange.mextent;
183 }
184 
187 {
188  Numeric s=0;
189  ConstIterator1D i = begin();
190  const ConstIterator1D e = end();
191 
192  for ( ; i!=e; ++i )
193  s += *i;
194 
195  return s;
196 }
197 
202 {
203  return ConstVectorView(mdata, mrange, r);
204 }
205 
208 {
209  return ConstIterator1D(mdata+mrange.mstart, mrange.mstride);
210 }
211 
214 {
215  return ConstIterator1D( mdata +
216  mrange.mstart +
217  (mrange.mextent)*mrange.mstride,
218  mrange.mstride );
219 }
220 
222 ConstVectorView::operator ConstMatrixView() const
223 {
224  // The old version (before 2013-01-18) of this was:
225  // return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
226  // Bus this was a bug! The problem is that the matrix index operator adds
227  // the mstart from both row and columm range object to mdata
228 
229  return ConstMatrixView(mdata,mrange,Range(0,1));
230 }
231 
240  mrange(0,1), mdata(&const_cast<Numeric&>(a))
241 {
242  // Nothing to do here.
243 }
244 
248  mrange(0,0), mdata(NULL)
249 {
250  // Nothing to do here.
251 }
252 
256  const Range& range) :
257  mrange(range),
258  mdata(data)
259 {
260  // Nothing to do here.
261 }
262 
274  const Range& p,
275  const Range& n) :
276  mrange(p,n),
277  mdata(data)
278 {
279  // Nothing to do here.
280 }
281 
285 std::ostream& operator<<(std::ostream& os, const ConstVectorView& v)
286 {
287  ConstIterator1D i=v.begin();
288  const ConstIterator1D end=v.end();
289 
290  if ( i!=end )
291  {
292  os << setw(3) << *i;
293  ++i;
294  }
295  for ( ; i!=end ; ++i )
296  {
297  os << " " << setw(3) << *i;
298  }
299 
300  return os;
301 }
302 
303 
304 // Functions for VectorView:
305 // ------------------------
306 
310 {
311  throw runtime_error("Creating a VectorView from a const Vector is not allowed.");
312  // This is not really a runtime error, but I don't want to start
313  // producing direct output from inside matpack. And just exiting is
314  // not so nice.
315  // If you see this error, there is a bug in the code, not in the
316  // ARTS input.
317 }
318 
321 {
322  mdata = v.mdata;
323  mrange = v.mrange;
324 }
325 
331 {
332  return ConstVectorView::operator[](r);
333 }
334 
339 {
340  return VectorView(mdata, mrange, r);
341 }
342 
347 {
348  return ConstVectorView::begin();
349 }
350 
355 {
356  return ConstVectorView::end();
357 }
358 
361 {
363 }
364 
367 {
368  return Iterator1D( mdata +
369  mrange.mstart +
371  mrange.mstride );
372 }
373 
379 {
380  // cout << "Assigning VectorView from ConstVectorView.\n";
381 
382  // Check that sizes are compatible:
383  assert(mrange.mextent==v.mrange.mextent);
384 
385  copy( v.begin(), v.end(), begin() );
386 
387  return *this;
388 }
389 
396 {
397  // cout << "Assigning VectorView from VectorView.\n";
398 
399  // Check that sizes are compatible:
400  assert(mrange.mextent==v.mrange.mextent);
401 
402  copy( v.begin(), v.end(), begin() );
403 
404  return *this;
405 }
406 
410 {
411  // cout << "Assigning VectorView from Vector.\n";
412 
413  // Check that sizes are compatible:
414  assert(mrange.mextent==v.mrange.mextent);
415 
416  copy( v.begin(), v.end(), begin() );
417 
418  return *this;
419 }
420 
424 {
425  copy( x, begin(), end() );
426  return *this;
427 }
428 
431 {
432  const Iterator1D e=end();
433  for ( Iterator1D i=begin(); i!=e ; ++i )
434  *i *= x;
435  return *this;
436 }
437 
440 {
441  const Iterator1D e=end();
442  for ( Iterator1D i=begin(); i!=e ; ++i )
443  *i /= x;
444  return *this;
445 }
446 
449 {
450  const Iterator1D e=end();
451  for ( Iterator1D i=begin(); i!=e ; ++i )
452  *i += x;
453  return *this;
454 }
455 
458 {
459  const Iterator1D e=end();
460  for ( Iterator1D i=begin(); i!=e ; ++i )
461  *i -= x;
462  return *this;
463 }
464 
467 {
468  assert( nelem()==x.nelem() );
469 
470  ConstIterator1D s=x.begin();
471 
472  Iterator1D i=begin();
473  const Iterator1D e=end();
474 
475  for ( ; i!=e ; ++i,++s )
476  *i *= *s;
477  return *this;
478 }
479 
482 {
483  assert( nelem()==x.nelem() );
484 
485  ConstIterator1D s=x.begin();
486 
487  Iterator1D i=begin();
488  const Iterator1D e=end();
489 
490  for ( ; i!=e ; ++i,++s )
491  *i /= *s;
492  return *this;
493 }
494 
497 {
498  assert( nelem()==x.nelem() );
499 
500  ConstIterator1D s=x.begin();
501 
502  Iterator1D i=begin();
503  const Iterator1D e=end();
504 
505  for ( ; i!=e ; ++i,++s )
506  *i += *s;
507  return *this;
508 }
509 
512 {
513  assert( nelem()==x.nelem() );
514 
515  ConstIterator1D s=x.begin();
516 
517  Iterator1D i=begin();
518  const Iterator1D e=end();
519 
520  for ( ; i!=e ; ++i,++s )
521  *i -= *s;
522  return *this;
523 }
524 
526 VectorView::operator MatrixView()
527 {
528  // The old version (before 2013-01-18) of this was:
529  // return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
530  // Bus this was a bug! The problem is that the matrix index operator adds
531  // the mstart from both row and columm range object to mdata
532 
533  return MatrixView(mdata,mrange,Range(0,1));
534 }
535 
543 {
544  if (mrange.mstart != 0 || mrange.mstride != 1)
545  throw std::runtime_error("A VectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");
546 
547  return mdata;
548 }
549 
557 {
558  if (mrange.mstart != 0 || mrange.mstride != 1)
559  throw std::runtime_error("A VectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");
560 
561  return mdata;
562 }
563 
568  ConstVectorView(a)
569 {
570  // Nothing to do here.
571 }
572 
577 {
578  // Nothing to do here.
579 }
580 
584  const Range& range) :
585  ConstVectorView(data,range)
586 {
587  // Nothing to do here.
588 }
589 
601  const Range& p,
602  const Range& n) :
603  ConstVectorView(data,p,n)
604 {
605  // Nothing to do here.
606 }
607 
612 void copy(ConstIterator1D origin,
613  const ConstIterator1D& end,
614  Iterator1D target)
615 {
616  if (origin.mstride == 1 && target.mstride == 1)
617  memcpy((void *)target.mx,
618  (void *)origin.mx,
619  sizeof(Numeric) * (end.mx - origin.mx));
620  else
621  for ( ; origin!=end ; ++origin,++target )
622  *target = *origin;
623 }
624 
626 void copy(Numeric x,
627  Iterator1D target,
628  const Iterator1D& end)
629 {
630  for ( ; target!=end ; ++target )
631  *target = x;
632 }
633 
634 
635 // Functions for Vector:
636 // ---------------------
637 
640 {
641  // Nothing to do here
642 }
643 
646  VectorView( new Numeric[n],
647  Range(0,n))
648 {
649  // Nothing to do here.
650 }
651 
654  VectorView( new Numeric[n],
655  Range(0,n))
656 {
657  // Here we can access the raw memory directly, for slightly
658  // increased efficiency:
659  const Numeric *stop = mdata+n;
660  for ( Numeric *x=mdata; x<stop; ++x )
661  *x = fill;
662 }
663 
672 Vector::Vector(Numeric start, Index extent, Numeric stride) :
673  VectorView( new Numeric[extent],
674  Range(0,extent))
675 {
676  // Fill with values:
677  Numeric x = start;
678  Iterator1D i=begin();
679  const Iterator1D e=end();
680  for ( ; i!=e; ++i )
681  {
682  *i = x;
683  x += stride;
684  }
685 }
686 
693  VectorView( new Numeric[v.nelem()],
694  Range(0,v.nelem()))
695 {
696  copy(v.begin(),v.end(),begin());
697 }
698 
702  VectorView( new Numeric[v.nelem()],
703  Range(0,v.nelem()))
704 {
705  copy(v.begin(),v.end(),begin());
706 }
707 
709 Vector::Vector(const std::vector<Numeric>& v) :
710  VectorView( new Numeric[v.size()],
711  Range(0,v.size()))
712 {
713  std::vector<Numeric>::const_iterator vec_it_end = v.end();
714  Iterator1D this_it = this->begin();
715  for (std::vector<Numeric>::const_iterator vec_it = v.begin();
716  vec_it != vec_it_end;
717  ++vec_it, ++this_it)
718  *this_it = *vec_it;
719 }
720 
722 
747 {
748  swap(*this, v);
749  return *this;
750 }
751 
753 
770 {
771  resize( x.nelem() );
773  return *this;
774 }
775 
779 {
781  return *this;
782 }
783 
784 // /** Assignment operator from VectorView. Assignment operators are not
785 // inherited. */
786 // inline Vector& Vector::operator=(const VectorView v)
787 // {
788 // cout << "Assigning Vector from Vector View.\n";
789 // // Check that sizes are compatible:
790 // assert(mrange.mextent==v.mrange.mextent);
791 // VectorView::operator=(v);
792 // return *this;
793 // }
794 
799 {
800  assert( 0<=n );
801  if ( mrange.mextent != n )
802  {
803  delete[] mdata;
804  mdata = new Numeric[n];
805  mrange.mstart = 0;
806  mrange.mextent = n;
807  mrange.mstride = 1;
808  }
809 }
810 
811 
813 void swap(Vector& v1, Vector& v2)
814 {
815  std::swap(v1.mrange, v2.mrange);
816  std::swap(v1.mdata, v2.mdata);
817 }
818 
819 
823 {
824  delete[] mdata;
825 }
826 
827 
828 // Functions for ConstMatrixView:
829 // ------------------------------
830 
833 {
834  return mrr.mextent;
835 }
836 
839 {
840  return mcr.mextent;
841 }
842 
847  const Range& c) const
848 {
849  return ConstMatrixView(mdata, mrr, mcr, r, c);
850 }
851 
858 {
859  // Check that c is valid:
860  assert( 0 <= c );
861  assert( c < mcr.mextent );
862 
863  return ConstVectorView(mdata + mcr.mstart + c*mcr.mstride,
864  mrr, r);
865 }
866 
873 {
874  // Check that r is valid:
875  assert( 0 <= r );
876  assert( r < mrr.mextent );
877 
878  return ConstVectorView(mdata + mrr.mstart + r*mrr.mstride,
879  mcr, c);
880 }
881 
884 {
885  return ConstIterator2D(ConstVectorView(mdata+mrr.mstart,
886  mcr),
887  mrr.mstride);
888 }
889 
892 {
893  return ConstIterator2D( ConstVectorView(mdata + mrr.mstart + (mrr.mextent)*mrr.mstride,
894  mcr),
895  mrr.mstride );
896 }
897 
901  mrr(0,0,1), mcr(0,0,1), mdata(NULL)
902 {
903  // Nothing to do here.
904 }
905 
910  const Range& rr,
911  const Range& cr) :
912  mrr(rr),
913  mcr(cr),
914  mdata(data)
915 {
916  // Nothing to do here.
917 }
918 
934  const Range& pr, const Range& pc,
935  const Range& nr, const Range& nc) :
936  mrr(pr,nr),
937  mcr(pc,nc),
938  mdata(data)
939 {
940  // Nothing to do here.
941 }
942 
950 std::ostream& operator<<(std::ostream& os, const ConstMatrixView& v)
951 {
952  // Row iterators:
953  ConstIterator2D ir=v.begin();
954  const ConstIterator2D end_row=v.end();
955 
956  if ( ir!=end_row )
957  {
958  ConstIterator1D ic = ir->begin();
959  const ConstIterator1D end_col = ir->end();
960 
961  if ( ic!=end_col )
962  {
963  os << setw(3) << *ic;
964  ++ic;
965  }
966  for ( ; ic!=end_col ; ++ic )
967  {
968  os << " " << setw(3) << *ic;
969  }
970  ++ir;
971  }
972  for ( ; ir!=end_row ; ++ir )
973  {
974  ConstIterator1D ic = ir->begin();
975  const ConstIterator1D end_col = ir->end();
976 
977  os << "\n";
978  if ( ic!=end_col )
979  {
980  os << setw(3) << *ic;
981  ++ic;
982  }
983  for ( ; ic!=end_col ; ++ic )
984  {
985  os << " " << setw(3) << *ic;
986  }
987  }
988 
989  return os;
990 }
991 
992 
993 // Functions for MatrixView:
994 // -------------------------
995 
1001 {
1002  return ConstMatrixView::operator()(r,c);
1003 }
1004 
1012 {
1013  return ConstMatrixView::operator()(r,c);
1014 }
1015 
1023 {
1024  return ConstMatrixView::operator()(r,c);
1025 }
1026 
1031 {
1032  return MatrixView(mdata, mrr, mcr, r, c);
1033 }
1034 
1040 {
1041  // Check that c is valid:
1042  assert( 0 <= c );
1043  assert( c < mcr.mextent );
1044 
1045  return VectorView(mdata + mcr.mstart + c*mcr.mstride,
1046  mrr, r);
1047 }
1048 
1054 {
1055  // Check that r is valid:
1056  assert( 0 <= r );
1057  assert( r < mrr.mextent );
1058 
1059  return VectorView(mdata + mrr.mstart + r*mrr.mstride,
1060  mcr, c);
1061 }
1062 
1066 {
1067  return ConstMatrixView::begin();
1068 }
1069 
1072 {
1073  return ConstMatrixView::end();
1074 }
1075 
1078 {
1080  mrr.mstride);
1081 }
1082 
1085 {
1087  mcr),
1088  mrr.mstride );
1089 }
1090 
1096 {
1097  // Check that sizes are compatible:
1098  assert(mrr.mextent==m.mrr.mextent);
1099  assert(mcr.mextent==m.mcr.mextent);
1100 
1101  copy( m.begin(), m.end(), begin() );
1102  return *this;
1103 }
1104 
1111 {
1112  // Check that sizes are compatible:
1113  assert(mrr.mextent==m.mrr.mextent);
1114  assert(mcr.mextent==m.mcr.mextent);
1115 
1116  copy( m.begin(), m.end(), begin() );
1117  return *this;
1118 }
1119 
1124 {
1125  // Check that sizes are compatible:
1126  assert(mrr.mextent==m.mrr.mextent);
1127  assert(mcr.mextent==m.mcr.mextent);
1128 
1129  copy( m.begin(), m.end(), begin() );
1130  return *this;
1131 }
1132 
1138 {
1139  // Check that sizes are compatible:
1140  assert( mrr.mextent==v.nelem() );
1141  assert( mcr.mextent==1 );
1142  // dummy = ConstMatrixView(v.mdata,v.mrange,Range(v.mrange.mstart,1));;
1143  ConstMatrixView dummy(v);
1144  copy( dummy.begin(), dummy.end(), begin() );
1145  return *this;
1146 }
1147 
1151 {
1152  copy( x, begin(), end() );
1153  return *this;
1154 }
1155 
1158 {
1159  const Iterator2D er=end();
1160  for ( Iterator2D r=begin(); r!=er ; ++r )
1161  {
1162  const Iterator1D ec = r->end();
1163  for ( Iterator1D c = r->begin(); c!=ec ; ++c )
1164  *c *= x;
1165  }
1166  return *this;
1167 }
1168 
1171 {
1172  const Iterator2D er=end();
1173  for ( Iterator2D r=begin(); r!=er ; ++r )
1174  {
1175  const Iterator1D ec = r->end();
1176  for ( Iterator1D c = r->begin(); c!=ec ; ++c )
1177  *c /= x;
1178  }
1179  return *this;
1180 }
1181 
1184 {
1185  const Iterator2D er=end();
1186  for ( Iterator2D r=begin(); r!=er ; ++r )
1187  {
1188  const Iterator1D ec = r->end();
1189  for ( Iterator1D c = r->begin(); c!=ec ; ++c )
1190  *c += x;
1191  }
1192  return *this;
1193 }
1194 
1197 {
1198  const Iterator2D er=end();
1199  for ( Iterator2D r=begin(); r!=er ; ++r )
1200  {
1201  const Iterator1D ec = r->end();
1202  for ( Iterator1D c = r->begin(); c!=ec ; ++c )
1203  *c -= x;
1204  }
1205  return *this;
1206 }
1207 
1215 {
1216  if (mrr.mstart != 0 || mrr.mstride != mcr.mextent
1217  || mcr.mstart != 0 || mcr.mstride != 1)
1218  throw std::runtime_error("A MatrixView can only be converted to a plain C-array if it's pointing to a continuous block of data");
1219 
1220  return mdata;
1221 }
1222 
1230 {
1231  if (mrr.mstart != 0 || mrr.mstride != mcr.mextent
1232  || mcr.mstart != 0 || mcr.mstride != 1)
1233  throw std::runtime_error("A MatrixView can only be converted to a plain C-array if it's pointing to a continuous block of data");
1234 
1235  return mdata;
1236 }
1237 
1240 {
1241  assert(nrows()==x.nrows());
1242  assert(ncols()==x.ncols());
1243  ConstIterator2D sr = x.begin();
1244  Iterator2D r = begin();
1245  const Iterator2D er = end();
1246  for ( ; r!=er ; ++r,++sr )
1247  {
1248  ConstIterator1D sc = sr->begin();
1249  Iterator1D c = r->begin();
1250  const Iterator1D ec = r->end();
1251  for ( ; c!=ec ; ++c,++sc )
1252  *c *= *sc;
1253  }
1254  return *this;
1255 }
1256 
1259 {
1260  assert(nrows()==x.nrows());
1261  assert(ncols()==x.ncols());
1262  ConstIterator2D sr = x.begin();
1263  Iterator2D r = begin();
1264  const Iterator2D er = end();
1265  for ( ; r!=er ; ++r,++sr )
1266  {
1267  ConstIterator1D sc = sr->begin();
1268  Iterator1D c = r->begin();
1269  const Iterator1D ec = r->end();
1270  for ( ; c!=ec ; ++c,++sc )
1271  *c /= *sc;
1272  }
1273  return *this;
1274 }
1275 
1278 {
1279  assert(nrows()==x.nrows());
1280  assert(ncols()==x.ncols());
1281  ConstIterator2D sr = x.begin();
1282  Iterator2D r = begin();
1283  const Iterator2D er = end();
1284  for ( ; r!=er ; ++r,++sr )
1285  {
1286  ConstIterator1D sc = sr->begin();
1287  Iterator1D c = r->begin();
1288  const Iterator1D ec = r->end();
1289  for ( ; c!=ec ; ++c,++sc )
1290  *c += *sc;
1291  }
1292  return *this;
1293 }
1294 
1297 {
1298  assert(nrows()==x.nrows());
1299  assert(ncols()==x.ncols());
1300  ConstIterator2D sr = x.begin();
1301  Iterator2D r = begin();
1302  const Iterator2D er = end();
1303  for ( ; r!=er ; ++r,++sr )
1304  {
1305  ConstIterator1D sc = sr->begin();
1306  Iterator1D c = r->begin();
1307  const Iterator1D ec = r->end();
1308  for ( ; c!=ec ; ++c,++sc )
1309  *c -= *sc;
1310  }
1311  return *this;
1312 }
1313 
1316 {
1317  assert(nrows()==x.nelem());
1318  assert(ncols()==1);
1319  ConstIterator1D sc = x.begin();
1320  Iterator2D r = begin();
1321  const Iterator2D er = end();
1322  for ( ; r!=er ; ++r,++sc )
1323  {
1324  Iterator1D c = r->begin();
1325  *c *= *sc;
1326  }
1327  return *this;
1328 }
1329 
1332 {
1333  assert(nrows()==x.nelem());
1334  assert(ncols()==1);
1335  ConstIterator1D sc = x.begin();
1336  Iterator2D r = begin();
1337  const Iterator2D er = end();
1338  for ( ; r!=er ; ++r,++sc )
1339  {
1340  Iterator1D c = r->begin();
1341  *c /= *sc;
1342  }
1343  return *this;
1344 }
1345 
1348 {
1349  assert(nrows()==x.nelem());
1350  assert(ncols()==1);
1351  ConstIterator1D sc = x.begin();
1352  Iterator2D r = begin();
1353  const Iterator2D er = end();
1354  for ( ; r!=er ; ++r,++sc )
1355  {
1356  Iterator1D c = r->begin();
1357  *c += *sc;
1358  }
1359  return *this;
1360 }
1361 
1364 {
1365  assert(nrows()==x.nelem());
1366  assert(ncols()==1);
1367  ConstIterator1D sc = x.begin();
1368  Iterator2D r = begin();
1369  const Iterator2D er = end();
1370  for ( ; r!=er ; ++r,++sc )
1371  {
1372  Iterator1D c = r->begin();
1373  *c -= *sc;
1374  }
1375  return *this;
1376 }
1377 
1381  ConstMatrixView()
1382 {
1383  // Nothing to do here.
1384 }
1385 
1390  const Range& rr,
1391  const Range& cr) :
1392  ConstMatrixView(data, rr, cr)
1393 {
1394  // Nothing to do here.
1395 }
1396 
1412  const Range& pr, const Range& pc,
1413  const Range& nr, const Range& nc) :
1414  ConstMatrixView(data,pr,pc,nr,nc)
1415 {
1416  // Nothing to do here.
1417 }
1418 
1428 void copy(ConstIterator2D origin,
1429  const ConstIterator2D& end,
1430  Iterator2D target)
1431 {
1432  for ( ; origin!=end ; ++origin,++target )
1433  {
1434  ConstIterator1D o = origin->begin();
1435  const ConstIterator1D e = origin->end();
1436  Iterator1D t = target->begin();
1437  for ( ; o!=e ; ++o,++t )
1438  *t = *o;
1439  }
1440 }
1441 
1443 void copy(Numeric x,
1444  Iterator2D target,
1445  const Iterator2D& end)
1446 {
1447  for ( ; target!=end ; ++target )
1448  {
1449  Iterator1D t = target->begin();
1450  const Iterator1D e = target->end();
1451  for ( ; t!=e ; ++t )
1452  *t = x;
1453  }
1454 }
1455 
1456 
1457 // Functions for Matrix:
1458 // ---------------------
1459 
1463 {
1464  // Nothing to do here. However, note that the default constructor
1465  // for MatrixView has been called in the initializer list. That is
1466  // crucial, otherwise internal range objects will not be properly
1467  // initialized.
1468 }
1469 
1473  MatrixView( new Numeric[r*c],
1474  Range(0,r,c),
1475  Range(0,c))
1476 {
1477  // Nothing to do here.
1478 }
1479 
1482  MatrixView( new Numeric[r*c],
1483  Range(0,r,c),
1484  Range(0,c))
1485 {
1486  // Here we can access the raw memory directly, for slightly
1487  // increased efficiency:
1488  const Numeric *stop = mdata+r*c;
1489  for ( Numeric *x=mdata; x<stop; ++x )
1490  *x = fill;
1491 }
1492 
1496  MatrixView( new Numeric[m.nrows()*m.ncols()],
1497  Range( 0, m.nrows(), m.ncols() ),
1498  Range( 0, m.ncols() ) )
1499 {
1500  copy(m.begin(),m.end(),begin());
1501 }
1502 
1506  MatrixView( new Numeric[m.nrows()*m.ncols()],
1507  Range( 0, m.nrows(), m.ncols() ),
1508  Range( 0, m.ncols() ) )
1509 {
1510  // There is a catch here: If m is an empty matrix, then it will have
1511  // 0 colunns. But this is used to initialize the stride of the row
1512  // Range! Thus, this method has to be consistent with the behaviour
1513  // of Range::Range. For now, Range::Range allows also stride 0.
1514  copy(m.begin(),m.end(),begin());
1515 }
1516 
1518 
1543 {
1544  swap(*this, m);
1545  return *this;
1546 }
1547 
1551 {
1552  copy( x, begin(), end() );
1553  return *this;
1554 }
1555 
1557 
1570 {
1571  resize( v.nelem(), 1 );
1572  ConstMatrixView dummy(v);
1573  copy( dummy.begin(), dummy.end(), begin() );
1574  return *this;
1575 }
1576 
1581 {
1582  assert( 0<=r );
1583  assert( 0<=c );
1584 
1585  if ( mrr.mextent!=r || mcr.mextent!=c )
1586  {
1587  delete[] mdata;
1588  mdata = new Numeric[r*c];
1589 
1590  mrr.mstart = 0;
1591  mrr.mextent = r;
1592  mrr.mstride = c;
1593 
1594  mcr.mstart = 0;
1595  mcr.mextent = c;
1596  mcr.mstride = 1;
1597  }
1598 }
1599 
1600 
1602 void swap(Matrix& m1, Matrix& m2)
1603 {
1604  std::swap(m1.mrr, m2.mrr);
1605  std::swap(m1.mcr, m2.mcr);
1606  std::swap(m1.mdata, m2.mdata);
1607 }
1608 
1609 
1613 {
1614 // cout << "Destroying a Matrix:\n"
1615 // << *this << "\n........................................\n";
1616  delete[] mdata;
1617 }
1618 
1619 
1620 // Some general Matrix Vector functions:
1621 
1624 {
1625  // Check dimensions:
1626  assert( a.nelem() == b.nelem() );
1627 
1628  const ConstIterator1D ae = a.end();
1629  ConstIterator1D ai = a.begin();
1630  ConstIterator1D bi = b.begin();
1631 
1632  Numeric res = 0;
1633  for ( ; ai!=ae ; ++ai, ++bi )
1634  res += (*ai) * (*bi);
1635 
1636  return res;
1637 }
1638 
1649  const ConstMatrixView& M,
1650  const ConstVectorView& x )
1651 {
1652  // Check dimensions:
1653  assert( y.mrange.mextent == M.mrr.mextent );
1654  assert( M.mcr.mextent == x.mrange.mextent );
1655  assert( M.mcr.mextent != 0 && M.mrr.mextent != 0);
1656 
1657  // Let's first find the pointers to the starting positions
1658  Numeric *mdata = M.mdata + M.mcr.mstart + M.mrr.mstart;
1659  Numeric *xdata = x.mdata + x.mrange.mstart;
1660  Numeric *yelem = y.mdata + y.mrange.mstart;
1661 
1662  Index i = M.mrr.mextent;
1663  while (i--)
1664  {
1665  Numeric *melem = mdata;
1666  Numeric *xelem = xdata; // Reset xelem to first element of source vector
1667 
1668  // Multiply first element of matrix row with first element of source
1669  // vector. This is done outside the loop to avoid initialization of the
1670  // target vector's element with zero (saves one assignment)
1671  *yelem = *melem * *xelem;
1672 
1673  Index j = M.mcr.mextent; // --j (instead of j-- like in the outer loop)
1674  while (--j) // is important here because we only want
1675  { // mextent-1 cycles
1676  melem += M.mcr.mstride;
1677  xelem += x.mrange.mstride;
1678  *yelem += *melem * *xelem;
1679  }
1680 
1681  mdata += M.mrr.mstride; // Jump to next matrix row
1682  yelem += y.mrange.mstride; // Jump to next element in target vector
1683  }
1684 }
1685 
1686 
1694  const ConstMatrixView& B,
1695  const ConstMatrixView& C )
1696 {
1697  // Check dimensions:
1698  assert( A.nrows() == B.nrows() );
1699  assert( A.ncols() == C.ncols() );
1700  assert( B.ncols() == C.nrows() );
1701 
1702  // Let's get the transpose of C, so that we can use 2D iterators to
1703  // access the columns (= rows of the transpose).
1704  ConstMatrixView CT = transpose(C);
1705 
1706  const Iterator2D ae = A.end();
1707  Iterator2D ai = A.begin();
1708  ConstIterator2D bi = B.begin();
1709 
1710  // This walks through the rows of A and B:
1711  for ( ; ai!=ae ; ++ai, ++bi )
1712  {
1713  const Iterator1D ace = ai->end();
1714  Iterator1D aci = ai->begin();
1715  ConstIterator2D cti = CT.begin();
1716 
1717  // This walks through the columns of A with a 1D iterator, and
1718  // at the same time through the rows of CT, which are the columns of
1719  // C, with a 2D iterator:
1720  for ( ; aci!=ace ; ++aci, ++cti )
1721  {
1722  // The operator * is used to compute the scalar product
1723  // between rows of B and rows of C.transpose().
1724  *aci = (*bi) * (*cti);
1725  }
1726  }
1727 }
1728 
1730 
1744 {
1745  assert( a.nelem() == 3 );
1746  assert( b.nelem() == 3 );
1747  assert( c.nelem() == 3 );
1748 
1749  c[0] = a[1]*b[2] - a[2]*b[1];
1750  c[1] = a[2]*b[0] - a[0]*b[2];
1751  c[2] = a[0]*b[1] - a[1]*b[0];
1752 }
1753 
1764 {
1765  assert( a.nelem() == b.nelem() );
1766  Numeric arg = (a*b) / sqrt(a*a) / sqrt(b*b);
1767 
1768  // Due to numerical inaccuracy, arg might be slightly larger than 1.
1769  // We catch those cases to avoid spurious returns of NaNs
1770  return fabs(arg) > 1. ? 0. : acos(arg) * RAD2DEG;
1771 }
1772 
1773 
1788 {
1789  assert( a.nelem() == b.nelem() );
1790  assert( a.nelem() == c.nelem() );
1791 
1792  const Numeric C = (a*b) / (a*a);
1793  c = a;
1794  c *= C;
1795 };
1796 
1797 
1800 {
1801  return ConstMatrixView(m.mdata, m.mcr, m.mrr);
1802 }
1803 
1807 {
1808  return MatrixView(m.mdata, m.mcr, m.mrr);
1809 }
1810 
1814 {
1815  return transpose((MatrixView)v);
1816 }
1817 
1839  double (&my_func)(double),
1840  ConstVectorView x )
1841 {
1842  // Check dimensions:
1843  assert( y.nelem()==x.nelem() );
1844 
1845  const ConstIterator1D xe = x.end();
1846  ConstIterator1D xi = x.begin();
1847  Iterator1D yi = y.begin();
1848  for ( ; xi!=xe; ++xi, ++yi )
1849  *yi = my_func(*xi);
1850 }
1851 
1871  double (&my_func)(double),
1872  ConstMatrixView x )
1873 {
1874  // Check dimensions:
1875  assert( y.nrows()==x.nrows() );
1876  assert( y.ncols()==x.ncols() );
1877 
1878  const ConstIterator2D rxe = x.end();
1879  ConstIterator2D rx = x.begin();
1880  Iterator2D ry = y.begin();
1881  for ( ; rx!=rxe; ++rx, ++ry )
1882  {
1883  const ConstIterator1D cxe = rx->end();
1884  ConstIterator1D cx = rx->begin();
1885  Iterator1D cy = ry->begin();
1886  for ( ; cx!=cxe; ++cx, ++cy )
1887  *cy = my_func(*cx);
1888  }
1889 }
1890 
1893 {
1894  // Initial value for max:
1895  Numeric max = x[0];
1896 
1897  const ConstIterator1D xe = x.end();
1898  ConstIterator1D xi = x.begin();
1899 
1900  for ( ; xi!=xe ; ++xi )
1901  {
1902  if ( *xi > max )
1903  max = *xi;
1904  }
1905 
1906  return max;
1907 }
1908 
1911 {
1912  // Initial value for max:
1913  Numeric max = x(0,0);
1914 
1915  const ConstIterator2D rxe = x.end();
1916  ConstIterator2D rx = x.begin();
1917 
1918  for ( ; rx!=rxe ; ++rx )
1919  {
1920  const ConstIterator1D cxe = rx->end();
1921  ConstIterator1D cx = rx->begin();
1922 
1923  for ( ; cx!=cxe ; ++cx )
1924  if ( *cx > max )
1925  max = *cx;
1926  }
1927 
1928  return max;
1929 }
1930 
1933 {
1934  // Initial value for min:
1935  Numeric min = x[0];
1936 
1937  const ConstIterator1D xe = x.end();
1938  ConstIterator1D xi = x.begin();
1939 
1940  for ( ; xi!=xe ; ++xi )
1941  {
1942  if ( *xi < min )
1943  min = *xi;
1944  }
1945 
1946  return min;
1947 }
1948 
1951 {
1952  // Initial value for min:
1953  Numeric min = x(0,0);
1954 
1955  const ConstIterator2D rxe = x.end();
1956  ConstIterator2D rx = x.begin();
1957 
1958  for ( ; rx!=rxe ; ++rx )
1959  {
1960  const ConstIterator1D cxe = rx->end();
1961  ConstIterator1D cx = rx->begin();
1962 
1963  for ( ; cx!=cxe ; ++cx )
1964  if ( *cx < min )
1965  min = *cx;
1966  }
1967 
1968  return min;
1969 }
1970 
1973 {
1974  // Initial value for mean:
1975  Numeric mean = 0;
1976 
1977  const ConstIterator1D xe = x.end();
1978  ConstIterator1D xi = x.begin();
1979 
1980  for ( ; xi!=xe ; ++xi ) mean += *xi;
1981 
1982  mean /= (Numeric)x.nelem();
1983 
1984  return mean;
1985 }
1986 
1989 {
1990  // Initial value for mean:
1991  Numeric mean = 0;
1992 
1993  const ConstIterator2D rxe = x.end();
1994  ConstIterator2D rx = x.begin();
1995 
1996  for ( ; rx!=rxe ; ++rx )
1997  {
1998  const ConstIterator1D cxe = rx->end();
1999  ConstIterator1D cx = rx->begin();
2000 
2001  for ( ; cx!=cxe ; ++cx ) mean += *cx;
2002  }
2003 
2004  mean /= (Numeric)(x.nrows()*x.ncols());
2005 
2006  return mean;
2007 }
2008 
2019 {
2020  // cout << "Assigning VectorView from Array<Numeric>.\n";
2021 
2022  // Check that sizes are compatible:
2023  assert(mrange.mextent==v.nelem());
2024 
2025  // Iterators for Array:
2026  Array<Numeric>::const_iterator i=v.begin();
2027  const Array<Numeric>::const_iterator e=v.end();
2028  // Iterator for Vector:
2029  Iterator1D target = begin();
2030 
2031  for ( ; i!=e ; ++i,++target )
2032  *target = *i;
2033 
2034  return *this;
2035 }
2036 
2038 // Helper function for debugging
2039 #ifndef NDEBUG
2040 
2055 {
2056  return mv(r, c);
2057 }
2058 
2059 #endif
2060 
2062 
friend class ConstVectorView
Definition: matpackI.h:158
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
friend void swap(Matrix &m1, Matrix &m2)
Swaps two objects.
Definition: matpackI.cc:1602
const Numeric * mx
Current position.
Definition: matpackI.h:276
void swap(Vector &v1, Vector &v2)
Swaps two objects.
Definition: matpackI.cc:813
The VectorView class.
Definition: matpackI.h:372
Matrix & operator=(Matrix x)
Assignment operator from another matrix.
Definition: matpackI.cc:1542
Index mstride
The stride.
Definition: matpackI.h:209
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b)
cross3
Definition: matpackI.cc:1743
ConstIterator1D begin() const
Return const iterator to first element.
Definition: matpackI.cc:207
#define C(a, b)
Definition: Faddeeva.cc:254
#define M
Definition: rng.cc:196
Index nelem() const
Number of elements.
Definition: array.h:176
friend void mult(VectorView, const ConstMatrixView &, const ConstVectorView &)
Matrix Vector multiplication.
Definition: matpackI.cc:1648
Index mstart
The start index.
Definition: matpackI.h:204
ConstIterator2D begin() const
Return const iterator to first row.
Definition: matpackI.cc:1065
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:667
Matrix()
Default constructor.
Definition: matpackI.cc:1461
MatrixView & operator*=(Numeric x)
Multiplication by scalar.
Definition: matpackI.cc:1157
const Numeric * get_c_array() const
Conversion to plain C-array.
Definition: matpackI.cc:1214
The Vector class.
Definition: matpackI.h:556
The MatrixView class.
Definition: matpackI.h:679
void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Copy data between begin and end to target.
Definition: matpackI.cc:612
Numeric operator()(Index r, Index c) const
Plain const index operator.
Definition: matpackI.h:687
ConstIterator2D end() const
Return const iterator behind last row.
Definition: matpackI.cc:1071
VectorView operator*=(Numeric x)
Multiplication by scalar.
Definition: matpackI.cc:430
The range class.
Definition: matpackI.h:148
Vector()
Default constructor.
Definition: matpackI.cc:639
friend class ConstIterator2D
Definition: matpackI.h:445
ConstVectorView()
Default constructor.
Definition: matpackI.cc:247
ConstMatrixView()
Default constructor.
Definition: matpackI.cc:900
Index mstride
Stride.
Definition: matpackI.h:278
friend class ConstMatrixView
Definition: matpackI.h:161
The row iterator class for sub matrices.
Definition: matpackI.h:468
Numeric operator*(const ConstVectorView &a, const ConstVectorView &b)
Scalar product.
Definition: matpackI.cc:1623
Numeric min(const ConstVectorView &x)
Min function, vector version.
Definition: matpackI.cc:1932
ConstIterator1D begin() const
Return const iterator to first element.
Definition: matpackI.cc:346
ConstIterator1D end() const
Return const iterator behind last element.
Definition: matpackI.cc:354
VectorView operator+=(Numeric x)
Addition of scalar.
Definition: matpackI.cc:448
ConstIterator1D end() const
Return const iterator behind last element.
Definition: matpackI.cc:213
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
Numeric operator[](Index n) const
Plain const index operator.
Definition: matpackI.h:383
ConstIterator2D begin() const
Return const iterator to first row.
Definition: matpackI.cc:883
Range(Index start, Index extent, Index stride=1)
Definition: matpackI.cc:52
VectorView operator-=(Numeric x)
Subtraction of scalar.
Definition: matpackI.cc:457
The const row iterator class for sub matrices.
Definition: matpackI.h:508
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
VectorView()
Default constructor.
Definition: matpackI.cc:575
const Numeric * get_c_array() const
Conversion to plain C-array.
Definition: matpackI.cc:542
friend class ConstMatrixView
Definition: matpackI.h:334
Range mrr
The row range of mdata that is actually used.
Definition: matpackI.h:663
The Joker class.
Definition: matpackI.h:114
friend class MatrixView
Definition: matpackI.h:639
void proj(Vector &c, ConstVectorView a, ConstVectorView b)
Definition: matpackI.cc:1787
Numeric max(const ConstVectorView &x)
Max function, vector version.
Definition: matpackI.cc:1892
Vector & operator=(Vector v)
Assignment from another Vector.
Definition: matpackI.cc:746
Numeric sum() const
The sum of all elements of a Vector.
Definition: matpackI.cc:186
The declarations of all the exception classes.
Numeric operator[](Index n) const
Plain const index operator.
Definition: matpackI.h:303
ConstIterator2D end() const
Return const iterator behind last row.
Definition: matpackI.cc:891
virtual ~Vector()
Destructor for Vector.
Definition: matpackI.cc:822
friend void swap(Vector &v1, Vector &v2)
Swaps two objects.
Definition: matpackI.cc:813
Numeric operator()(Index r, Index c) const
Plain const index operator.
Definition: matpackI.h:607
std::ostream & operator<<(std::ostream &os, const ConstVectorView &v)
Output operator.
Definition: matpackI.cc:285
const Joker joker
MatrixView()
Default constructor.
Definition: matpackI.cc:1380
friend class VectorView
Definition: matpackI.h:332
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
The Matrix class.
Definition: matpackI.h:788
Range mcr
The column range of mdata that is actually used.
Definition: matpackI.h:665
Index mextent
The number of elements.
Definition: matpackI.h:207
MatrixView & operator=(const ConstMatrixView &v)
Assignment operator.
Definition: matpackI.cc:1095
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
MatrixView & operator-=(Numeric x)
Subtraction of scalar.
Definition: matpackI.cc:1196
friend ConstMatrixView transpose(ConstMatrixView m)
Const version of transpose.
Definition: matpackI.cc:1799
A constant view of a Vector.
Definition: matpackI.h:292
Numeric vector_angle(ConstVectorView a, ConstVectorView b)
Definition: matpackI.cc:1763
VectorView & operator=(const ConstVectorView &v)
Assignment operator.
Definition: matpackI.cc:378
Numeric mean(const ConstVectorView &x)
Mean function, vector version.
Definition: matpackI.cc:1972
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:358
A constant view of a Matrix.
Definition: matpackI.h:596
Numeric debug_matrixview_get_elem(MatrixView &mv, Index r, Index c)
Helper function to access matrix elements.
Definition: matpackI.cc:2054
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1838
MatrixView & operator/=(Numeric x)
Division by scalar.
Definition: matpackI.cc:1170
Numeric * mx
Current position.
Definition: matpackI.h:242
Range mrange
The range of mdata that is actually used.
Definition: matpackI.h:356
The constant iterator class for sub vectors.
Definition: matpackI.h:249
Index mstride
Stride.
Definition: matpackI.h:244
MatrixView & operator+=(Numeric x)
Addition of scalar.
Definition: matpackI.cc:1183
The iterator class for sub vectors.
Definition: matpackI.h:214
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
const Numeric RAD2DEG
Global constant, conversion from radians to degrees.
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580
VectorView operator/=(Numeric x)
Division by scalar.
Definition: matpackI.cc:439
virtual ~Matrix()
Destructor for Matrix.
Definition: matpackI.cc:1612