ARTS  2.3.1285(git:92a29ea9-dirty)
matpackI.cc
Go to the documentation of this file.
1 /* Copyright (C) 2001-2019 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 
24 #include "matpackI.h"
25 #include <cmath>
26 #include <cstring>
27 #include "blas.h"
28 #include "exceptions.h"
29 
30 using std::cout;
31 using std::endl;
32 using std::runtime_error;
33 using std::setw;
34 
35 // Define the global joker object:
36 extern const Joker joker = Joker();
37 
38 static const Numeric RAD2DEG = 57.295779513082323;
39 
40 std::ostream& operator<<(std::ostream& os, const Range& r) {
41  os << "Range(" << r.get_start() << ", " << r.get_extent() << ", "
42  << r.get_stride() << ")";
43  return os;
44 }
45 
46 // Functions for ConstVectorView:
47 // ------------------------------
48 
49 bool ConstVectorView::empty() const { return (nelem() == 0); }
50 
52 
54  Numeric s = 0;
56  const ConstIterator1D e = end();
57 
58  for (; i != e; ++i) s += *i;
59 
60  return s;
61 }
62 
64  return ConstVectorView(mdata, mrange, r);
65 }
66 
69 }
70 
72  return ConstIterator1D(
74  mrange.mstride);
75 }
76 
77 ConstVectorView::operator ConstMatrixView() const {
78  // The old version (before 2013-01-18) of this was:
79  // return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
80  // Bus this was a bug! The problem is that the matrix index operator adds
81  // the mstart from both row and columm range object to mdata
82 
83  return ConstMatrixView(mdata, mrange, Range(0, 1));
84 }
85 
86 /* This one is a bit tricky: We have to cast away the arguments const
87  qualifier, because mdata is not const. This should be safe, since
88  there are no non-const methods for ConstVectorView. */
90  : mrange(0, 1), mdata(&const_cast<Numeric&>(a)) {
91  // Nothing to do here.
92 }
93 
95  : mrange(range), mdata(data) {
96  // Nothing to do here.
97 }
98 
100  : mrange(p, n), mdata(data) {
101  // Nothing to do here.
102 }
103 
104 /* Output operator. This demonstrates how iterators can be used to
105  traverse the vector. The iterators know which part of the vector
106  is `active', and also the stride. */
107 std::ostream& operator<<(std::ostream& os, const ConstVectorView& v) {
108  ConstIterator1D i = v.begin();
109  const ConstIterator1D end = v.end();
110 
111  if (i != end) {
112  os << setw(3) << *i;
113  ++i;
114  }
115  for (; i != end; ++i) {
116  os << " " << setw(3) << *i;
117  }
118 
119  return os;
120 }
121 
122 // Functions for VectorView:
123 // ------------------------
124 
126  throw runtime_error(
127  "Creating a VectorView from a const Vector is not allowed.");
128  // This is not really a runtime error, but I don't want to start
129  // producing direct output from inside matpack. And just exiting is
130  // not so nice.
131  // If you see this error, there is a bug in the code, not in the
132  // ARTS input.
133 }
134 
136  mdata = v.mdata;
137  mrange = v.mrange;
138 }
139 
141  return VectorView(mdata, mrange, r);
142 }
143 
146 }
147 
150  mrange.mstride);
151 }
152 
154  // cout << "Assigning VectorView from ConstVectorView.\n";
155 
156  // Check that sizes are compatible:
157  assert(mrange.mextent == v.mrange.mextent);
158 
159  copy(v.begin(), v.end(), begin());
160 
161  return *this;
162 }
163 
165  // cout << "Assigning VectorView from VectorView.\n";
166 
167  // Check that sizes are compatible:
168  assert(mrange.mextent == v.mrange.mextent);
169 
170  copy(v.begin(), v.end(), begin());
171 
172  return *this;
173 }
174 
176  // cout << "Assigning VectorView from Vector.\n";
177 
178  // Check that sizes are compatible:
179  assert(mrange.mextent == v.mrange.mextent);
180 
181  copy(v.begin(), v.end(), begin());
182 
183  return *this;
184 }
185 
187  copy(x, begin(), end());
188  return *this;
189 }
190 
192  const Iterator1D e = end();
193  for (Iterator1D i = begin(); i != e; ++i) *i *= x;
194  return *this;
195 }
196 
198  const Iterator1D e = end();
199  for (Iterator1D i = begin(); i != e; ++i) *i /= x;
200  return *this;
201 }
202 
204  const Iterator1D e = end();
205  for (Iterator1D i = begin(); i != e; ++i) *i += x;
206  return *this;
207 }
208 
210  const Iterator1D e = end();
211  for (Iterator1D i = begin(); i != e; ++i) *i -= x;
212  return *this;
213 }
214 
216  assert(nelem() == x.nelem());
217 
218  ConstIterator1D s = x.begin();
219 
220  Iterator1D i = begin();
221  const Iterator1D e = end();
222 
223  for (; i != e; ++i, ++s) *i *= *s;
224  return *this;
225 }
226 
228  assert(nelem() == x.nelem());
229 
230  ConstIterator1D s = x.begin();
231 
232  Iterator1D i = begin();
233  const Iterator1D e = end();
234 
235  for (; i != e; ++i, ++s) *i /= *s;
236  return *this;
237 }
238 
240  assert(nelem() == x.nelem());
241 
242  ConstIterator1D s = x.begin();
243 
244  Iterator1D i = begin();
245  const Iterator1D e = end();
246 
247  for (; i != e; ++i, ++s) *i += *s;
248  return *this;
249 }
250 
252  assert(nelem() == x.nelem());
253 
254  ConstIterator1D s = x.begin();
255 
256  Iterator1D i = begin();
257  const Iterator1D e = end();
258 
259  for (; i != e; ++i, ++s) *i -= *s;
260  return *this;
261 }
262 
263 VectorView::operator MatrixView() {
264  // The old version (before 2013-01-18) of this was:
265  // return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
266  // Bus this was a bug! The problem is that the matrix index operator adds
267  // the mstart from both row and columm range object to mdata
268 
269  return MatrixView(mdata, mrange, Range(0, 1));
270 }
271 
273  if (mrange.mstart != 0 || mrange.mstride != 1)
274  throw std::runtime_error(
275  "A VectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");
276 
277  return mdata;
278 }
279 
281  if (mrange.mstart != 0 || mrange.mstride != 1)
282  throw std::runtime_error(
283  "A VectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");
284 
285  return mdata;
286 }
287 
289  // Nothing to do here.
290 }
291 
293  : ConstVectorView(data, range) {
294  // Nothing to do here.
295 }
296 
298  : ConstVectorView(data, p, n) {
299  // Nothing to do here.
300 }
301 
302 void copy(ConstIterator1D origin,
303  const ConstIterator1D& end,
304  Iterator1D target) {
305  if (origin.mstride == 1 && target.mstride == 1)
306  memcpy((void*)target.mx,
307  (void*)origin.mx,
308  sizeof(Numeric) * (end.mx - origin.mx));
309  else
310  for (; origin != end; ++origin, ++target) *target = *origin;
311 }
312 
313 void copy(Numeric x, Iterator1D target, const Iterator1D& end) {
314  for (; target != end; ++target) *target = x;
315 }
316 
317 // Functions for Vector:
318 // ---------------------
319 
320 Vector::Vector(std::initializer_list<Numeric> init)
321  : VectorView(new Numeric[init.size()], Range(0, init.size())) {
322  std::copy(init.begin(), init.end(), begin());
323 }
324 
326  // Nothing to do here.
327 }
328 
330  : VectorView(new Numeric[n], Range(0, n)) {
331  // Here we can access the raw memory directly, for slightly
332  // increased efficiency:
333  std::fill_n(mdata, n, fill);
334 }
335 
337  : VectorView(new Numeric[extent], Range(0, extent)) {
338  // Fill with values:
339  Numeric x = start;
340  Iterator1D i = begin();
341  const Iterator1D e = end();
342  for (; i != e; ++i) {
343  *i = x;
344  x += stride;
345  }
346 }
347 
349  : VectorView(new Numeric[v.nelem()], Range(0, v.nelem())) {
350  copy(v.begin(), v.end(), begin());
351 }
352 
354  : VectorView(new Numeric[v.nelem()], Range(0, v.nelem())) {
355  std::memcpy(mdata, v.mdata, nelem() * sizeof(Numeric));
356 }
357 
358 Vector::Vector(const std::vector<Numeric>& v)
359  : VectorView(new Numeric[v.size()], Range(0, v.size())) {
360  std::vector<Numeric>::const_iterator vec_it_end = v.end();
361  Iterator1D this_it = this->begin();
362  for (std::vector<Numeric>::const_iterator vec_it = v.begin();
363  vec_it != vec_it_end;
364  ++vec_it, ++this_it)
365  *this_it = *vec_it;
366 }
367 
368 Vector& Vector::operator=(std::initializer_list<Numeric> v) {
369  resize(v.size());
370  std::copy(v.begin(), v.end(), begin());
371  return *this;
372 }
373 
375  if (this != &v) {
376  resize(v.nelem());
377  std::memcpy(mdata, v.mdata, nelem() * sizeof(Numeric));
378  }
379  return *this;
380 }
381 
383  if (this != &v) {
384  delete[] mdata;
385  mdata = v.mdata;
386  mrange = v.mrange;
387  v.mrange = Range(0, 0);
388  v.mdata = nullptr;
389  }
390  return *this;
391 }
392 
394  resize(x.nelem());
396  return *this;
397 }
398 
400  std::fill_n(mdata, nelem(), x);
401  return *this;
402 }
403 
405  assert(0 <= n);
406  if (mrange.mextent != n) {
407  delete[] mdata;
408  mdata = new Numeric[n];
409  mrange.mstart = 0;
410  mrange.mextent = n;
411  mrange.mstride = 1;
412  }
413 }
414 
415 void swap(Vector& v1, Vector& v2) {
416  std::swap(v1.mrange, v2.mrange);
417  std::swap(v1.mdata, v2.mdata);
418 }
419 
420 Vector::~Vector() { delete[] mdata; }
421 
422 // Functions for ConstMatrixView:
423 // ------------------------------
424 
426 bool ConstMatrixView::empty() const { return (nrows() == 0 || ncols() == 0); }
427 
429 Index ConstMatrixView::nrows() const { return mrr.mextent; }
430 
432 Index ConstMatrixView::ncols() const { return mcr.mextent; }
433 
438  const Range& c) const {
439  return ConstMatrixView(mdata, mrr, mcr, r, c);
440 }
441 
448  // Check that c is valid:
449  assert(0 <= c);
450  assert(c < mcr.mextent);
451 
452  return ConstVectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
453 }
454 
461  // Check that r is valid:
462  assert(0 <= r);
463  assert(r < mrr.mextent);
464 
465  return ConstVectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
466 }
467 
470  return ConstIterator2D(ConstVectorView(mdata + mrr.mstart, mcr), mrr.mstride);
471 }
472 
475  return ConstIterator2D(
476  ConstVectorView(mdata + mrr.mstart + (mrr.mextent) * mrr.mstride, mcr),
477  mrr.mstride);
478 }
479 
481 
490  Index n = std::min(mrr.mextent, mcr.mextent);
491  return ConstVectorView(mdata + mrr.mstart + mcr.mstart,
492  Range(0, n, mrr.mstride + mcr.mstride));
493 }
494 
499  const Range& rr,
500  const Range& cr)
501  : mrr(rr), mcr(cr), mdata(data) {
502  // Nothing to do here.
503 }
504 
520  const Range& pr,
521  const Range& pc,
522  const Range& nr,
523  const Range& nc)
524  : mrr(pr, nr), mcr(pc, nc), mdata(data) {
525  // Nothing to do here.
526 }
527 
535 std::ostream& operator<<(std::ostream& os, const ConstMatrixView& v) {
536  // Row iterators:
537  ConstIterator2D ir = v.begin();
538  const ConstIterator2D end_row = v.end();
539 
540  if (ir != end_row) {
541  ConstIterator1D ic = ir->begin();
542  const ConstIterator1D end_col = ir->end();
543 
544  if (ic != end_col) {
545  os << setw(3) << *ic;
546  ++ic;
547  }
548  for (; ic != end_col; ++ic) {
549  os << " " << setw(3) << *ic;
550  }
551  ++ir;
552  }
553  for (; ir != end_row; ++ir) {
554  ConstIterator1D ic = ir->begin();
555  const ConstIterator1D end_col = ir->end();
556 
557  os << "\n";
558  if (ic != end_col) {
559  os << setw(3) << *ic;
560  ++ic;
561  }
562  for (; ic != end_col; ++ic) {
563  os << " " << setw(3) << *ic;
564  }
565  }
566 
567  return os;
568 }
569 
570 // Functions for MatrixView:
571 // -------------------------
572 
577  return MatrixView(mdata, mrr, mcr, r, c);
578 }
579 
585  // Check that c is valid:
586  assert(0 <= c);
587  assert(c < mcr.mextent);
588 
589  return VectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
590 }
591 
597  // Check that r is valid:
598  assert(0 <= r);
599  assert(r < mrr.mextent);
600 
601  return VectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
602 }
603 
605 // hiden by the non-const operator of the derived class.*/
606 //ConstIterator2D MatrixView::begin() const
607 //{
608 // return ConstMatrixView::begin();
609 //}
610 //
612 //ConstIterator2D MatrixView::end() const
613 //{
614 // return ConstMatrixView::end();
615 //}
616 
620 }
621 
624  return Iterator2D(
626  mrr.mstride);
627 }
628 
634  // Check that sizes are compatible:
635  assert(mrr.mextent == m.mrr.mextent);
636  assert(mcr.mextent == m.mcr.mextent);
637 
638  copy(m.begin(), m.end(), begin());
639  return *this;
640 }
641 
648  // Check that sizes are compatible:
649  assert(mrr.mextent == m.mrr.mextent);
650  assert(mcr.mextent == m.mcr.mextent);
651 
652  copy(m.begin(), m.end(), begin());
653  return *this;
654 }
655 
660  // Check that sizes are compatible:
661  assert(mrr.mextent == m.mrr.mextent);
662  assert(mcr.mextent == m.mcr.mextent);
663 
664  copy(m.begin(), m.end(), begin());
665  return *this;
666 }
667 
673  // Check that sizes are compatible:
674  assert(mrr.mextent == v.nelem());
675  assert(mcr.mextent == 1);
676  // dummy = ConstMatrixView(v.mdata,v.mrange,Range(v.mrange.mstart,1));;
677  ConstMatrixView dummy(v);
678  copy(dummy.begin(), dummy.end(), begin());
679  return *this;
680 }
681 
685  copy(x, begin(), end());
686  return *this;
687 }
688 
691  const Iterator2D er = end();
692  for (Iterator2D r = begin(); r != er; ++r) {
693  const Iterator1D ec = r->end();
694  for (Iterator1D c = r->begin(); c != ec; ++c) *c *= x;
695  }
696  return *this;
697 }
698 
701  const Iterator2D er = end();
702  for (Iterator2D r = begin(); r != er; ++r) {
703  const Iterator1D ec = r->end();
704  for (Iterator1D c = r->begin(); c != ec; ++c) *c /= x;
705  }
706  return *this;
707 }
708 
711  const Iterator2D er = end();
712  for (Iterator2D r = begin(); r != er; ++r) {
713  const Iterator1D ec = r->end();
714  for (Iterator1D c = r->begin(); c != ec; ++c) *c += x;
715  }
716  return *this;
717 }
718 
721  const Iterator2D er = end();
722  for (Iterator2D r = begin(); r != er; ++r) {
723  const Iterator1D ec = r->end();
724  for (Iterator1D c = r->begin(); c != ec; ++c) *c -= x;
725  }
726  return *this;
727 }
728 
736  if (mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
737  mcr.mstride != 1)
738  throw std::runtime_error(
739  "A MatrixView can only be converted to a plain C-array if it's pointing to a continuous block of data");
740 
741  return mdata;
742 }
743 
751  if (mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
752  mcr.mstride != 1)
753  throw std::runtime_error(
754  "A MatrixView can only be converted to a plain C-array if it's pointing to a continuous block of data");
755 
756  return mdata;
757 }
758 
761  assert(nrows() == x.nrows());
762  assert(ncols() == x.ncols());
763  ConstIterator2D sr = x.begin();
764  Iterator2D r = begin();
765  const Iterator2D er = end();
766  for (; r != er; ++r, ++sr) {
767  ConstIterator1D sc = sr->begin();
768  Iterator1D c = r->begin();
769  const Iterator1D ec = r->end();
770  for (; c != ec; ++c, ++sc) *c *= *sc;
771  }
772  return *this;
773 }
774 
777  assert(nrows() == x.nrows());
778  assert(ncols() == x.ncols());
779  ConstIterator2D sr = x.begin();
780  Iterator2D r = begin();
781  const Iterator2D er = end();
782  for (; r != er; ++r, ++sr) {
783  ConstIterator1D sc = sr->begin();
784  Iterator1D c = r->begin();
785  const Iterator1D ec = r->end();
786  for (; c != ec; ++c, ++sc) *c /= *sc;
787  }
788  return *this;
789 }
790 
793  assert(nrows() == x.nrows());
794  assert(ncols() == x.ncols());
795  ConstIterator2D sr = x.begin();
796  Iterator2D r = begin();
797  const Iterator2D er = end();
798  for (; r != er; ++r, ++sr) {
799  ConstIterator1D sc = sr->begin();
800  Iterator1D c = r->begin();
801  const Iterator1D ec = r->end();
802  for (; c != ec; ++c, ++sc) *c += *sc;
803  }
804  return *this;
805 }
806 
809  assert(nrows() == x.nrows());
810  assert(ncols() == x.ncols());
811  ConstIterator2D sr = x.begin();
812  Iterator2D r = begin();
813  const Iterator2D er = end();
814  for (; r != er; ++r, ++sr) {
815  ConstIterator1D sc = sr->begin();
816  Iterator1D c = r->begin();
817  const Iterator1D ec = r->end();
818  for (; c != ec; ++c, ++sc) *c -= *sc;
819  }
820  return *this;
821 }
822 
825  assert(nrows() == x.nelem());
826  assert(ncols() == 1);
827  ConstIterator1D sc = x.begin();
828  Iterator2D r = begin();
829  const Iterator2D er = end();
830  for (; r != er; ++r, ++sc) {
831  Iterator1D c = r->begin();
832  *c *= *sc;
833  }
834  return *this;
835 }
836 
839  assert(nrows() == x.nelem());
840  assert(ncols() == 1);
841  ConstIterator1D sc = x.begin();
842  Iterator2D r = begin();
843  const Iterator2D er = end();
844  for (; r != er; ++r, ++sc) {
845  Iterator1D c = r->begin();
846  *c /= *sc;
847  }
848  return *this;
849 }
850 
853  assert(nrows() == x.nelem());
854  assert(ncols() == 1);
855  ConstIterator1D sc = x.begin();
856  Iterator2D r = begin();
857  const Iterator2D er = end();
858  for (; r != er; ++r, ++sc) {
859  Iterator1D c = r->begin();
860  *c += *sc;
861  }
862  return *this;
863 }
864 
867  assert(nrows() == x.nelem());
868  assert(ncols() == 1);
869  ConstIterator1D sc = x.begin();
870  Iterator2D r = begin();
871  const Iterator2D er = end();
872  for (; r != er; ++r, ++sc) {
873  Iterator1D c = r->begin();
874  *c -= *sc;
875  }
876  return *this;
877 }
878 
883  : ConstMatrixView(data, rr, cr) {
884  // Nothing to do here.
885 }
886 
902  const Range& pr,
903  const Range& pc,
904  const Range& nr,
905  const Range& nc)
906  : ConstMatrixView(data, pr, pc, nr, nc) {
907  // Nothing to do here.
908 }
909 
919 void copy(ConstIterator2D origin,
920  const ConstIterator2D& end,
921  Iterator2D target) {
922  for (; origin != end; ++origin, ++target) {
923  ConstIterator1D o = origin->begin();
924  const ConstIterator1D e = origin->end();
925  Iterator1D t = target->begin();
926  for (; o != e; ++o, ++t) *t = *o;
927  }
928 }
929 
931 void copy(Numeric x, Iterator2D target, const Iterator2D& end) {
932  for (; target != end; ++target) {
933  Iterator1D t = target->begin();
934  const Iterator1D e = target->end();
935  for (; t != e; ++t) *t = x;
936  }
937 }
938 
939 // Functions for Matrix:
940 // ---------------------
941 
945  : MatrixView(new Numeric[r * c], Range(0, r, c), Range(0, c)) {
946  // Nothing to do here.
947 }
948 
951  : MatrixView(new Numeric[r * c], Range(0, r, c), Range(0, c)) {
952  // Here we can access the raw memory directly, for slightly
953  // increased efficiency:
954  std::fill_n(mdata, r * c, fill);
955 }
956 
960  : MatrixView(new Numeric[m.nrows() * m.ncols()],
961  Range(0, m.nrows(), m.ncols()),
962  Range(0, m.ncols())) {
963  copy(m.begin(), m.end(), begin());
964 }
965 
969  : MatrixView(new Numeric[m.nrows() * m.ncols()],
970  Range(0, m.nrows(), m.ncols()),
971  Range(0, m.ncols())) {
972  // There is a catch here: If m is an empty matrix, then it will have
973  // 0 colunns. But this is used to initialize the stride of the row
974  // Range! Thus, this method has to be consistent with the behaviour
975  // of Range::Range. For now, Range::Range allows also stride 0.
976  std::memcpy(mdata, m.mdata, nrows() * ncols() * sizeof(Numeric));
977 }
978 
980 
1005  if (this != &m) {
1006  resize(m.nrows(), m.ncols());
1007  std::memcpy(mdata, m.mdata, nrows() * ncols() * sizeof(Numeric));
1008  }
1009  return *this;
1010 }
1011 
1014  if (this != &m) {
1015  delete[] mdata;
1016  mdata = m.mdata;
1017  mrr = m.mrr;
1018  mcr = m.mcr;
1019  m.mrr = Range(0, 0);
1020  m.mcr = Range(0, 0);
1021  m.mdata = nullptr;
1022  }
1023  return *this;
1024 }
1025 
1029  std::fill_n(mdata, nrows() * ncols(), x);
1030  return *this;
1031 }
1032 
1034 
1047  resize(v.nelem(), 1);
1048  ConstMatrixView dummy(v);
1049  copy(dummy.begin(), dummy.end(), begin());
1050  return *this;
1051 }
1052 
1057  assert(0 <= r);
1058  assert(0 <= c);
1059 
1060  if (mrr.mextent != r || mcr.mextent != c) {
1061  delete[] mdata;
1062  mdata = new Numeric[r * c];
1063 
1064  mrr.mstart = 0;
1065  mrr.mextent = r;
1066  mrr.mstride = c;
1067 
1068  mcr.mstart = 0;
1069  mcr.mextent = c;
1070  mcr.mstride = 1;
1071  }
1072 }
1073 
1075 void swap(Matrix& m1, Matrix& m2) {
1076  std::swap(m1.mrr, m2.mrr);
1077  std::swap(m1.mcr, m2.mcr);
1078  std::swap(m1.mdata, m2.mdata);
1079 }
1080 
1084  // cout << "Destroying a Matrix:\n"
1085  // << *this << "\n........................................\n";
1086  delete[] mdata;
1087 }
1088 
1089 // Some general Matrix Vector functions:
1090 
1093  // Check dimensions:
1094  assert(a.nelem() == b.nelem());
1095 
1096  const ConstIterator1D ae = a.end();
1097  ConstIterator1D ai = a.begin();
1098  ConstIterator1D bi = b.begin();
1099 
1100  Numeric res = 0;
1101  for (; ai != ae; ++ai, ++bi) res += (*ai) * (*bi);
1102 
1103  return res;
1104 }
1105 
1107 
1123 void mult(VectorView y, const ConstMatrixView& M, const ConstVectorView& x) {
1124  assert(y.mrange.get_extent() == M.mrr.get_extent());
1125  assert(M.mcr.get_extent() == x.mrange.get_extent());
1126  assert((M.mcr.get_extent() != 0) && (M.mrr.get_extent() != 0));
1127 
1128  if ((M.mcr.get_stride() == 1) || (M.mrr.get_stride() == 1)) {
1129  char trans;
1130  int m, n;
1131  double zero = 0.0;
1132  double one = 1.0;
1133  int LDA, incx, incy;
1134 
1135  if (M.mcr.get_stride() != 1) {
1136  trans = 'n';
1137  m = (int)M.mrr.get_extent();
1138  n = (int)M.mcr.get_extent();
1139  LDA = (int)M.mcr.get_stride();
1140  } else {
1141  trans = 't';
1142  m = (int)M.mcr.get_extent();
1143  n = (int)M.mrr.get_extent();
1144  LDA = (int)M.mrr.get_stride();
1145  if (M.mrr.get_stride() == 1) LDA = m;
1146  }
1147 
1148  incx = (int)x.mrange.get_stride();
1149  incy = (int)y.mrange.get_stride();
1150 
1151  double* mstart = M.mdata + M.mcr.get_start() + M.mrr.get_start();
1152  double* ystart = y.mdata + y.mrange.get_start();
1153  double* xstart = x.mdata + x.mrange.get_start();
1154 
1155  dgemv_(&trans,
1156  &m,
1157  &n,
1158  &one,
1159  mstart,
1160  &LDA,
1161  xstart,
1162  &incx,
1163  &zero,
1164  ystart,
1165  &incy);
1166 
1167  } else {
1168  mult_general(y, M, x);
1169  }
1170 }
1171 
1182  const ConstMatrixView& M,
1183  const ConstVectorView& x) {
1184  // Check dimensions:
1185  assert(y.mrange.mextent == M.mrr.mextent);
1186  assert(M.mcr.mextent == x.mrange.mextent);
1187  assert(M.mcr.mextent != 0 && M.mrr.mextent != 0);
1188 
1189  // Let's first find the pointers to the starting positions
1190  Numeric* mdata = M.mdata + M.mcr.mstart + M.mrr.mstart;
1191  Numeric* xdata = x.mdata + x.mrange.mstart;
1192  Numeric* yelem = y.mdata + y.mrange.mstart;
1193 
1194  Index i = M.mrr.mextent;
1195  while (i--) {
1196  Numeric* melem = mdata;
1197  Numeric* xelem = xdata; // Reset xelem to first element of source vector
1198 
1199  // Multiply first element of matrix row with first element of source
1200  // vector. This is done outside the loop to avoid initialization of the
1201  // target vector's element with zero (saves one assignment)
1202  *yelem = *melem * *xelem;
1203 
1204  Index j = M.mcr.mextent; // --j (instead of j-- like in the outer loop)
1205  while (--j) // is important here because we only want
1206  { // mextent-1 cycles
1207  melem += M.mcr.mstride;
1208  xelem += x.mrange.mstride;
1209  *yelem += *melem * *xelem;
1210  }
1211 
1212  mdata += M.mrr.mstride; // Jump to next matrix row
1213  yelem += y.mrange.mstride; // Jump to next element in target vector
1214  }
1215 }
1216 
1218 
1242 void mult(MatrixView A, const ConstMatrixView& B, const ConstMatrixView& C) {
1243  // Check dimensions:
1244  assert(A.nrows() == B.nrows());
1245  assert(A.ncols() == C.ncols());
1246  assert(B.ncols() == C.nrows());
1247 
1248  // Catch trivial case if one of the matrices is empty.
1249  if ((B.nrows() == 0) || (B.ncols() == 0) || (C.ncols() == 0)) return;
1250 
1251  // Matrices B and C must be continuous in at least on dimension, C
1252  // must be continuous along the second dimension.
1253  if (((B.mrr.get_stride() == 1) || (B.mcr.get_stride() == 1)) &&
1254  ((C.mrr.get_stride() == 1) || (C.mcr.get_stride() == 1)) &&
1255  (A.mcr.get_stride() == 1)) {
1256  // BLAS uses column-major order while arts uses row-major order.
1257  // Hence instead of C = A * B we compute C^T = A^T * B^T!
1258 
1259  int k, m, n;
1260 
1261  k = (int)B.ncols();
1262  m = (int)C.ncols();
1263  n = (int)B.nrows();
1264 
1265  // Note also the clash in nomenclature: BLAS uses C = A * B while
1266  // arts uses A = B * C. Taking into accout this and the difference in
1267  // memory layouts, we need to map the MatrixViews A, B and C to the BLAS
1268  // arguments as follows:
1269  // A (arts) -> C (BLAS)
1270  // B (arts) -> B (BLAS)
1271  // C (arts) -> A (BLAS)
1272 
1273  // Char indicating whether A (BLAS) or B (BLAS) should be transposed.
1274  char transa, transb;
1275  // Sizes of the matrices along the direction in which they are
1276  // traversed.
1277  int lda, ldb, ldc;
1278 
1279  // Check if C (arts) is transposed.
1280  if (C.mrr.get_stride() == 1) {
1281  transa = 'T';
1282  lda = (int)C.mcr.get_stride();
1283  } else {
1284  transa = 'N';
1285  lda = (int)C.mrr.get_stride();
1286  }
1287 
1288  // Check if B (arts) is transposed.
1289  if (B.mrr.get_stride() == 1) {
1290  transb = 'T';
1291  ldb = (int)B.mcr.get_stride();
1292  } else {
1293  transb = 'N';
1294  ldb = (int)B.mrr.get_stride();
1295  }
1296 
1297  // In case B (arts) has only one column, column and row stride are 1.
1298  // We therefore need to set ldb to k, since dgemm_ requires lda to be at
1299  // least k / m if A is non-transposed / transposed.
1300  if ((B.mcr.get_stride() == 1) && (B.mrr.get_stride() == 1)) {
1301  transb = 'N';
1302  ldb = k;
1303  }
1304 
1305  // The same holds for C (arts).
1306  if ((C.mcr.get_stride() == 1) && (C.mrr.get_stride() == 1)) {
1307  transa = 'N';
1308  lda = m;
1309  }
1310 
1311  ldc = (int)A.mrr.get_stride();
1312  // The same holds for A (arts).
1313  if ((A.mcr.get_stride() == 1) && (A.mrr.get_stride() == 1)) {
1314  ldc = m;
1315  }
1316  double alpha = 1.0, beta = 0.0;
1317 
1318  dgemm_(&transa,
1319  &transb,
1320  &m,
1321  &n,
1322  &k,
1323  &alpha,
1324  C.mdata + C.mrr.get_start() + C.mcr.get_start(),
1325  &lda,
1326  B.mdata + B.mrr.get_start() + B.mcr.get_start(),
1327  &ldb,
1328  &beta,
1329  A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1330  &ldc);
1331 
1332  } else {
1333  mult_general(A, B, C);
1334  }
1335 }
1336 
1338 
1347  const ConstMatrixView& B,
1348  const ConstMatrixView& C) {
1349  // Check dimensions:
1350  assert(A.nrows() == B.nrows());
1351  assert(A.ncols() == C.ncols());
1352  assert(B.ncols() == C.nrows());
1353 
1354  // Let's get the transpose of C, so that we can use 2D iterators to
1355  // access the columns (= rows of the transpose).
1356  ConstMatrixView CT = transpose(C);
1357 
1358  const Iterator2D ae = A.end();
1359  Iterator2D ai = A.begin();
1360  ConstIterator2D bi = B.begin();
1361 
1362  // This walks through the rows of A and B:
1363  for (; ai != ae; ++ai, ++bi) {
1364  const Iterator1D ace = ai->end();
1365  Iterator1D aci = ai->begin();
1366  ConstIterator2D cti = CT.begin();
1367 
1368  // This walks through the columns of A with a 1D iterator, and
1369  // at the same time through the rows of CT, which are the columns of
1370  // C, with a 2D iterator:
1371  for (; aci != ace; ++aci, ++cti) {
1372  // The operator * is used to compute the scalar product
1373  // between rows of B and rows of C.transpose().
1374  *aci = (*bi) * (*cti);
1375  }
1376  }
1377 }
1378 
1380 
1393 void cross3(VectorView c, const ConstVectorView& a, const ConstVectorView& b) {
1394  assert(a.nelem() == 3);
1395  assert(b.nelem() == 3);
1396  assert(c.nelem() == 3);
1397 
1398  c[0] = a[1] * b[2] - a[2] * b[1];
1399  c[1] = a[2] * b[0] - a[0] * b[2];
1400  c[2] = a[0] * b[1] - a[1] * b[0];
1401 }
1402 
1413  assert(a.nelem() == b.nelem());
1414  Numeric arg = (a * b) / sqrt(a * a) / sqrt(b * b);
1415 
1416  // Due to numerical inaccuracy, arg might be slightly larger than 1.
1417  // We catch those cases to avoid spurious returns of NaNs
1418  return fabs(arg) > 1. ? 0. : acos(arg) * RAD2DEG;
1419 }
1420 
1435  assert(a.nelem() == b.nelem());
1436  assert(a.nelem() == c.nelem());
1437 
1438  const Numeric C = (a * b) / (a * a);
1439  c = a;
1440  c *= C;
1441 };
1442 
1445  return ConstMatrixView(m.mdata, m.mcr, m.mrr);
1446 }
1447 
1451 
1455 
1476 void transform(VectorView y, double (&my_func)(double), ConstVectorView x) {
1477  // Check dimensions:
1478  assert(y.nelem() == x.nelem());
1479 
1480  const ConstIterator1D xe = x.end();
1481  ConstIterator1D xi = x.begin();
1482  Iterator1D yi = y.begin();
1483  for (; xi != xe; ++xi, ++yi) *yi = my_func(*xi);
1484 }
1485 
1504 void transform(MatrixView y, double (&my_func)(double), ConstMatrixView x) {
1505  // Check dimensions:
1506  assert(y.nrows() == x.nrows());
1507  assert(y.ncols() == x.ncols());
1508 
1509  const ConstIterator2D rxe = x.end();
1510  ConstIterator2D rx = x.begin();
1511  Iterator2D ry = y.begin();
1512  for (; rx != rxe; ++rx, ++ry) {
1513  const ConstIterator1D cxe = rx->end();
1514  ConstIterator1D cx = rx->begin();
1515  Iterator1D cy = ry->begin();
1516  for (; cx != cxe; ++cx, ++cy) *cy = my_func(*cx);
1517  }
1518 }
1519 
1522  // Initial value for max:
1523  Numeric max = x[0];
1524 
1525  const ConstIterator1D xe = x.end();
1526  ConstIterator1D xi = x.begin();
1527 
1528  for (; xi != xe; ++xi) {
1529  if (*xi > max) max = *xi;
1530  }
1531 
1532  return max;
1533 }
1534 
1537  // Initial value for max:
1538  Numeric max = x(0, 0);
1539 
1540  const ConstIterator2D rxe = x.end();
1541  ConstIterator2D rx = x.begin();
1542 
1543  for (; rx != rxe; ++rx) {
1544  const ConstIterator1D cxe = rx->end();
1545  ConstIterator1D cx = rx->begin();
1546 
1547  for (; cx != cxe; ++cx)
1548  if (*cx > max) max = *cx;
1549  }
1550 
1551  return max;
1552 }
1553 
1556  // Initial value for min:
1557  Numeric min = x[0];
1558 
1559  const ConstIterator1D xe = x.end();
1560  ConstIterator1D xi = x.begin();
1561 
1562  for (; xi != xe; ++xi) {
1563  if (*xi < min) min = *xi;
1564  }
1565 
1566  return min;
1567 }
1568 
1571  // Initial value for min:
1572  Numeric min = x(0, 0);
1573 
1574  const ConstIterator2D rxe = x.end();
1575  ConstIterator2D rx = x.begin();
1576 
1577  for (; rx != rxe; ++rx) {
1578  const ConstIterator1D cxe = rx->end();
1579  ConstIterator1D cx = rx->begin();
1580 
1581  for (; cx != cxe; ++cx)
1582  if (*cx < min) min = *cx;
1583  }
1584 
1585  return min;
1586 }
1587 
1590  // Initial value for mean:
1591  Numeric mean = 0;
1592 
1593  const ConstIterator1D xe = x.end();
1594  ConstIterator1D xi = x.begin();
1595 
1596  for (; xi != xe; ++xi) mean += *xi;
1597 
1598  mean /= (Numeric)x.nelem();
1599 
1600  return mean;
1601 }
1602 
1605  // Initial value for mean:
1606  Numeric mean = 0;
1607 
1608  const ConstIterator2D rxe = x.end();
1609  ConstIterator2D rx = x.begin();
1610 
1611  for (; rx != rxe; ++rx) {
1612  const ConstIterator1D cxe = rx->end();
1613  ConstIterator1D cx = rx->begin();
1614 
1615  for (; cx != cxe; ++cx) mean += *cx;
1616  }
1617 
1618  mean /= (Numeric)(x.nrows() * x.ncols());
1619 
1620  return mean;
1621 }
1622 
1633  // cout << "Assigning VectorView from Array<Numeric>.\n";
1634 
1635  // Check that sizes are compatible:
1636  assert(mrange.mextent == v.nelem());
1637 
1638  // Iterators for Array:
1639  Array<Numeric>::const_iterator i = v.begin();
1640  const Array<Numeric>::const_iterator e = v.end();
1641  // Iterator for Vector:
1642  Iterator1D target = begin();
1643 
1644  for (; i != e; ++i, ++target) *target = *i;
1645 
1646  return *this;
1647 }
1648 
1649 // Const
1650 
1651 // Converts constant matrix to constant eigen map
1653  return ConstMatrixViewMap(A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1654  A.nrows(),
1655  A.ncols(),
1656  StrideType(A.mrr.get_stride(), A.mcr.get_stride()));
1657 }
1658 
1659 // Converts constant vector to constant eigen row-view
1661  return ConstMatrixViewMap(A.mdata + A.mrange.get_start(),
1662  A.nelem(),
1663  1,
1664  StrideType(A.mrange.get_stride(), 1));
1665 }
1666 
1667 // Converts constant vector to constant eigen row-view
1669  return MapToEigen(A);
1670 }
1671 
1672 // Converts constant vector to constant eigen column-view
1674  return ConstMatrixViewMap(A.mdata + A.mrange.get_start(),
1675  1,
1676  A.nelem(),
1677  StrideType(1, A.mrange.get_stride()));
1678 }
1679 
1680 // Non- const
1681 
1682 // Converts matrix to eigen map
1684  return MatrixViewMap(A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1685  A.nrows(),
1686  A.ncols(),
1687  StrideType(A.mrr.get_stride(), A.mcr.get_stride()));
1688 }
1689 
1690 // Converts vector to eigen map row-view
1692  return MatrixViewMap(A.mdata + A.mrange.get_start(),
1693  A.nelem(),
1694  1,
1695  StrideType(A.mrange.get_stride(), 1));
1696 }
1697 
1698 // Converts vector to eigen map row-view
1700 
1701 // Converts vector to eigen map column-view
1703  return MatrixViewMap(A.mdata + A.mrange.get_start(),
1704  1,
1705  A.nelem(),
1706  StrideType(1, A.mrange.get_stride()));
1707 }
1708 
1709 // Special 4x4
1710 
1711 // Converts matrix to eigen map
1713  return Matrix4x4ViewMap(A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1714  4,
1715  4,
1716  StrideType(A.mrr.get_stride(), A.mcr.get_stride()));
1717 }
1718 
1719 // Converts constant matrix to constant eigen map
1721  return ConstMatrix4x4ViewMap(
1722  A.mdata + A.mrr.get_start() + A.mcr.get_start(),
1723  4,
1724  4,
1725  StrideType(A.mrr.get_stride(), A.mcr.get_stride()));
1726 }
1727 
1729 // Helper function for debugging
1730 #ifndef NDEBUG
1731 
1746  return mv(r, c);
1747 }
1748 
1749 #endif
1750 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
friend void swap(Matrix &m1, Matrix &m2)
Swaps two objects.
Definition: matpackI.cc:1075
Iterator2D begin()
Return const iterator to first row.
Definition: matpackI.cc:618
const Numeric * mx
Current position.
Definition: matpackI.h:460
void swap(Vector &v1, Vector &v2)
Definition: matpackI.cc:415
The VectorView class.
Definition: matpackI.h:610
Index mstride
The stride.
Definition: matpackI.h:355
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b)
cross3
Definition: matpackI.cc:1393
ConstIterator1D begin() const
Return const iterator to first element.
Definition: matpackI.cc:67
Index nelem() const
Number of elements.
Definition: array.h:195
bool empty() const
Returns true if variable size is zero.
Definition: matpackI.cc:426
Interface for BLAS library.
Index mstart
The start index.
Definition: matpackI.h:346
ConstVectorView diagonal() const
Matrix diagonal as vector.
Definition: matpackI.cc:489
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1081
Eigen::Map< MatrixType, 0, StrideType > MatrixViewMap
Definition: matpackI.h:110
MatrixView & operator*=(Numeric x)
Multiplication by scalar.
Definition: matpackI.cc:690
const Numeric * get_c_array() const
Conversion to plain C-array.
Definition: matpackI.cc:735
The Vector class.
Definition: matpackI.h:860
The MatrixView class.
Definition: matpackI.h:1093
void copy(ConstIterator1D origin, const ConstIterator1D &end, Iterator1D target)
Definition: matpackI.cc:302
VectorView operator*=(Numeric x)
Multiplication by scalar.
Definition: matpackI.cc:191
The range class.
Definition: matpackI.h:160
Numeric & operator[](Index n)
Plain Index operator.
Definition: matpackI.h:631
void dgemv_(char *trans, int *m, int *n, double *alpha, double *A, int *LDA, double *x, int *incx, double *beta, double *y, int *incy)
Matrix-Vector Multiplication.
Eigen::Map< Matrix4x4Type, 0, StrideType > Matrix4x4ViewMap
Definition: matpackI.h:113
friend class ConstIterator2D
Definition: matpackI.h:727
bool empty() const
Returns true if variable size is zero.
Definition: matpackI.cc:49
Index mstride
Stride.
Definition: matpackI.h:462
The row iterator class for sub matrices.
Definition: matpackI.h:765
Numeric operator*(const ConstVectorView &a, const ConstVectorView &b)
Scalar product.
Definition: matpackI.cc:1092
Iterator2D end()
Return iterator behind last row.
Definition: matpackI.cc:623
Numeric min(const ConstVectorView &x)
Min function, vector version.
Definition: matpackI.cc:1555
Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > StrideType
Definition: matpackI.h:104
VectorView operator+=(Numeric x)
Addition of scalar.
Definition: matpackI.cc:203
ConstIterator1D end() const
Return const iterator behind last element.
Definition: matpackI.cc:71
constexpr Index get_start() const
Returns the start index of the range.
Definition: matpackI.h:327
G0 G2 FVC Y DV Numeric Numeric Numeric Zeeman LowerQuantumNumbers void * data
Vector()=default
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
Matrix & operator=(const Matrix &x)
Assignment operator from another matrix.
Definition: matpackI.cc:1004
ConstIterator2D begin() const
Return const iterator to first row.
Definition: matpackI.cc:469
Iterator1D begin()
Return iterator to first element.
Definition: matpackI.cc:144
VectorView operator-=(Numeric x)
Subtraction of scalar.
Definition: matpackI.cc:209
The const row iterator class for sub matrices.
Definition: matpackI.h:809
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
Iterator1D end()
Return iterator behind last element.
Definition: matpackI.cc:148
const Numeric * get_c_array() const
Conversion to plain C-array, const-version.
Definition: matpackI.cc:272
friend class ConstMatrixView
Definition: matpackI.h:542
Range mrr
The row range of mdata that is actually used.
Definition: matpackI.h:1077
MatrixView()=default
The Joker class.
Definition: matpackI.h:126
friend class MatrixView
Definition: matpackI.h:1029
ConstMatrixViewMap MapToEigenCol(const ConstVectorView &A)
Definition: matpackI.cc:1673
friend void mult_general(VectorView, const ConstMatrixView &, const ConstVectorView &)
Matrix Vector multiplication.
Definition: matpackI.cc:1181
void proj(Vector &c, ConstVectorView a, ConstVectorView b)
Definition: matpackI.cc:1434
Numeric max(const ConstVectorView &x)
Max function, vector version.
Definition: matpackI.cc:1521
Numeric sum() const
The sum of all elements of a Vector.
Definition: matpackI.cc:53
constexpr Index get_extent() const
Returns the extent of the range.
Definition: matpackI.h:329
The declarations of all the exception classes.
Eigen::Map< const Matrix4x4Type, 0, StrideType > ConstMatrix4x4ViewMap
Definition: matpackI.h:114
Numeric operator[](Index n) const
Plain const index operator.
Definition: matpackI.h:507
ConstIterator2D end() const
Return const iterator behind last row.
Definition: matpackI.cc:474
virtual ~Vector()
Destructor for Vector.
Definition: matpackI.cc:420
#define beta
friend void swap(Vector &v1, Vector &v2)
Swaps two objects.
Definition: matpackI.cc:415
Numeric operator()(Index r, Index c) const
Plain const index operator.
Definition: matpackI.h:999
Eigen::Map< const MatrixType, 0, StrideType > ConstMatrixViewMap
Definition: matpackI.h:111
Matrix()=default
const Joker joker
friend class VectorView
Definition: matpackI.h:540
VectorView()=default
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
The Matrix class.
Definition: matpackI.h:1193
std::ostream & operator<<(std::ostream &os, const Range &r)
Definition: matpackI.cc:40
constexpr Index get_stride() const
Returns the stride of the range.
Definition: matpackI.h:331
Range mcr
The column range of mdata that is actually used.
Definition: matpackI.h:1079
Implementation of Matrix, Vector, and such stuff.
Index mextent
The number of elements.
Definition: matpackI.h:353
This can be used to make arrays out of anything.
Definition: array.h:40
friend ConstMatrix4x4ViewMap MapToEigen4x4(const ConstMatrixView &)
Definition: matpackI.cc:1720
MatrixView & operator=(const ConstMatrixView &v)
Assignment operator.
Definition: matpackI.cc:633
Workspace & init(Workspace &ws)
ConstMatrixView()=default
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
MatrixView & operator-=(Numeric x)
Subtraction of scalar.
Definition: matpackI.cc:720
friend ConstMatrixView transpose(ConstMatrixView m)
Const version of transpose.
Definition: matpackI.cc:1444
friend ConstMatrixViewMap MapToEigen(const ConstMatrixView &)
Definition: matpackI.cc:1652
A constant view of a Vector.
Definition: matpackI.h:476
Numeric vector_angle(ConstVectorView a, ConstVectorView b)
Definition: matpackI.cc:1412
VectorView & operator=(const ConstVectorView &v)
Assignment operator.
Definition: matpackI.cc:153
Vector & operator=(const Vector &v)
Assignment from another Vector.
Definition: matpackI.cc:374
Numeric mean(const ConstVectorView &x)
Mean function, vector version.
Definition: matpackI.cc:1589
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
Definition: zeemandata.h:77
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:596
A constant view of a Matrix.
Definition: matpackI.h:982
ConstVectorView()=default
Numeric debug_matrixview_get_elem(MatrixView &mv, Index r, Index c)
Helper function to access matrix elements.
Definition: matpackI.cc:1745
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:1476
MatrixView & operator/=(Numeric x)
Division by scalar.
Definition: matpackI.cc:700
Numeric * mx
Current position.
Definition: matpackI.h:413
Range mrange
The range of mdata that is actually used.
Definition: matpackI.h:594
The constant iterator class for sub vectors.
Definition: matpackI.h:420
Index mstride
Stride.
Definition: matpackI.h:415
MatrixView & operator+=(Numeric x)
Addition of scalar.
Definition: matpackI.cc:710
void dgemm_(char *transa, char *transb, int *m, int *n, int *k, double *alpha, double *A, int *lda, double *B, int *ldb, double *beta, double *C, int *ldc)
BLAS matrix multiplication.
The iterator class for sub vectors.
Definition: matpackI.h:360
friend void mult(MatrixView, const ConstMatrixView &, const ConstMatrixView &)
Matrix-Matrix Multiplication.
Definition: matpackI.cc:1242
Numeric & operator()(Index r, Index c)
Plain index operator.
Definition: matpackI.h:1108
ConstMatrixViewMap MapToEigenRow(const ConstVectorView &A)
Definition: matpackI.cc:1668
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
const Numeric RAD2DEG
Global constant, conversion from radians to degrees.
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
VectorView operator/=(Numeric x)
Division by scalar.
Definition: matpackI.cc:197
virtual ~Matrix()
Destructor for Matrix.
Definition: matpackI.cc:1083