complex.cc

Go to the documentation of this file.

/* Copyright (C) 2002-2012 Stefan Buehler <sbuehler@ltu.se>

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */

#include "complex.h"
#include <cmath>
#include <cstring>
#include "blas.h"
#include "exceptions.h"

using std::runtime_error;
using std::setw;

// Functions for ConstComplexVectorView:
// ------------------------------

// bool ConstComplexVectorView::empty() const
// {
//     return (nelem() == 0);
// }

// Index ConstComplexVectorView::nelem() const
// {
//   return mrange.mextent;
// }

Complex ConstComplexVectorView::sum() const {
  Complex s = 0;
  ConstComplexIterator1D i = begin();
  const ConstComplexIterator1D e = end();

  for (; i != e; ++i) s += *i;

  return s;
}

ConstComplexVectorView ConstComplexVectorView::operator[](
    const Range& r) const {
  return ConstComplexVectorView(mdata, mrange, r);
}

ConstComplexIterator1D ConstComplexVectorView::begin() const {
  return ConstComplexIterator1D(mdata + mrange.mstart, mrange.mstride);
}

ConstComplexIterator1D ConstComplexVectorView::end() const {
  return ConstComplexIterator1D(
      mdata + mrange.mstart + (mrange.mextent) * mrange.mstride,
      mrange.mstride);
}

ConstComplexVectorView::operator ConstComplexMatrixView() const {
  return ConstComplexMatrixView(mdata, mrange, Range(0, 1));
}

ConstComplexVectorView::ConstComplexVectorView(const Complex& a)
    : mrange(0, 1), mdata(&const_cast<Complex&>(a)) {
  // Nothing to do here.
}

ConstComplexVectorView::ConstComplexVectorView(Complex* data,
                                               const Range& range)
    : mrange(range), mdata(data) {
  // Nothing to do here.
}

ConstComplexVectorView::ConstComplexVectorView(Complex* data,
                                               const Range& p,
                                               const Range& n)
    : mrange(p, n), mdata(data) {
  // Nothing to do here.
}

std::ostream& operator<<(std::ostream& os, const ConstComplexVectorView& v) {
  ConstComplexIterator1D i = v.begin();
  const ConstComplexIterator1D end = v.end();

  if (i != end) {
    os << setw(3) << *i;
    ++i;
  }
  for (; i != end; ++i) {
    os << " " << setw(3) << *i;
  }

  return os;
}

// Functions for ComplexVectorView:
// ------------------------

ComplexVectorView::ComplexVectorView(const ComplexVector&) {
  throw runtime_error(
      "Creating a ComplexVectorView from a const ComplexVector is not allowed.");
  // This is not really a runtime error, but I don't want to start
  // producing direct output from inside matpack. And just exiting is
  // not so nice.
  // If you see this error, there is a bug in the code, not in the
  // ARTS input.
}

ComplexVectorView::ComplexVectorView(ComplexVector& v) {
  mdata = v.mdata;
  mrange = v.mrange;
}

ComplexVectorView ComplexVectorView::operator[](const Range& r) {
  return ComplexVectorView(mdata, mrange, r);
}

ComplexIterator1D ComplexVectorView::begin() {
  return ComplexIterator1D(mdata + mrange.mstart, mrange.mstride);
}

ComplexIterator1D ComplexVectorView::end() {
  return ComplexIterator1D(
      mdata + mrange.mstart + (mrange.mextent) * mrange.mstride,
      mrange.mstride);
}

ComplexVectorView& ComplexVectorView::operator=(
    const ConstComplexVectorView& v) {
  //  cout << "Assigning VectorView from ConstVectorView.\n";

  // Check that sizes are compatible:
  assert(mrange.mextent == v.mrange.mextent);

  copy(v.begin(), v.end(), begin());

  return *this;
}

ComplexVectorView& ComplexVectorView::operator=(const ComplexVectorView& v) {
  //  cout << "Assigning ComplexVectorView from ComplexVectorView.\n";

  // Check that sizes are compatible:
  assert(mrange.mextent == v.mrange.mextent);

  copy(v.begin(), v.end(), begin());

  return *this;
}

ComplexVectorView& ComplexVectorView::operator=(const ComplexVector& v) {
  //  cout << "Assigning ComplexVectorView from Vector.\n";

  // Check that sizes are compatible:
  assert(mrange.mextent == v.mrange.mextent);

  copy(v.begin(), v.end(), begin());

  return *this;
}

ComplexVectorView& ComplexVectorView::operator=(Complex x) {
  copy(x, begin(), end());
  return *this;
}

ComplexVectorView ComplexVectorView::operator*=(Complex x) {
  const ComplexIterator1D e = end();
  for (ComplexIterator1D i = begin(); i != e; ++i) *i *= x;
  return *this;
}

ComplexVectorView ComplexVectorView::operator*=(Numeric x) {
  const ComplexIterator1D e = end();
  for (ComplexIterator1D i = begin(); i != e; ++i) *i *= x;
  return *this;
}

ComplexVectorView ComplexVectorView::operator/=(Complex x) {
  const ComplexIterator1D e = end();
  for (ComplexIterator1D i = begin(); i != e; ++i) *i /= x;
  return *this;
}

ComplexVectorView ComplexVectorView::operator/=(Numeric x) {
  const ComplexIterator1D e = end();
  for (ComplexIterator1D i = begin(); i != e; ++i) *i /= x;
  return *this;
}

ComplexVectorView ComplexVectorView::operator+=(Complex x) {
  const ComplexIterator1D e = end();
  for (ComplexIterator1D i = begin(); i != e; ++i) *i += x;
  return *this;
}

ComplexVectorView ComplexVectorView::operator+=(Numeric x) {
  const ComplexIterator1D e = end();
  for (ComplexIterator1D i = begin(); i != e; ++i) *i += x;
  return *this;
}

ComplexVectorView ComplexVectorView::operator-=(Complex x) {
  const ComplexIterator1D e = end();
  for (ComplexIterator1D i = begin(); i != e; ++i) *i -= x;
  return *this;
}

ComplexVectorView ComplexVectorView::operator-=(Numeric x) {
  const ComplexIterator1D e = end();
  for (ComplexIterator1D i = begin(); i != e; ++i) *i -= x;
  return *this;
}

ComplexVectorView ComplexVectorView::operator*=(
    const ConstComplexVectorView& x) {
  assert(nelem() == x.nelem());

  ConstComplexIterator1D s = x.begin();

  ComplexIterator1D i = begin();
  const ComplexIterator1D e = end();

  for (; i != e; ++i, ++s) *i *= *s;
  return *this;
}

ComplexVectorView ComplexVectorView::operator*=(const ConstVectorView& x) {
  assert(nelem() == x.nelem());

  ConstIterator1D s = x.begin();

  ComplexIterator1D i = begin();
  const ComplexIterator1D e = end();

  for (; i != e; ++i, ++s) *i *= *s;
  return *this;
}

ComplexVectorView ComplexVectorView::operator/=(
    const ConstComplexVectorView& x) {
  assert(nelem() == x.nelem());

  ConstComplexIterator1D s = x.begin();

  ComplexIterator1D i = begin();
  const ComplexIterator1D e = end();

  for (; i != e; ++i, ++s) *i /= *s;
  return *this;
}

ComplexVectorView ComplexVectorView::operator/=(const ConstVectorView& x) {
  assert(nelem() == x.nelem());

  ConstIterator1D s = x.begin();

  ComplexIterator1D i = begin();
  const ComplexIterator1D e = end();

  for (; i != e; ++i, ++s) *i /= *s;
  return *this;
}

ComplexVectorView ComplexVectorView::operator+=(
    const ConstComplexVectorView& x) {
  assert(nelem() == x.nelem());

  ConstComplexIterator1D s = x.begin();

  ComplexIterator1D i = begin();
  const ComplexIterator1D e = end();

  for (; i != e; ++i, ++s) *i += *s;
  return *this;
}

ComplexVectorView ComplexVectorView::operator+=(const ConstVectorView& x) {
  assert(nelem() == x.nelem());

  ConstIterator1D s = x.begin();

  ComplexIterator1D i = begin();
  const ComplexIterator1D e = end();

  for (; i != e; ++i, ++s) *i += *s;
  return *this;
}

ComplexVectorView ComplexVectorView::operator-=(
    const ConstComplexVectorView& x) {
  assert(nelem() == x.nelem());

  ConstComplexIterator1D s = x.begin();

  ComplexIterator1D i = begin();
  const ComplexIterator1D e = end();

  for (; i != e; ++i, ++s) *i -= *s;
  return *this;
}

ComplexVectorView ComplexVectorView::operator-=(const ConstVectorView& x) {
  assert(nelem() == x.nelem());

  ConstIterator1D s = x.begin();

  ComplexIterator1D i = begin();
  const ComplexIterator1D e = end();

  for (; i != e; ++i, ++s) *i -= *s;
  return *this;
}

ComplexVectorView::operator ComplexMatrixView() {
  // The old version (before 2013-01-18) of this was:
  //    return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
  // Bus this was a bug! The problem is that the matrix index operator adds
  // the mstart from both row and columm range object to mdata

  return ComplexMatrixView(mdata, mrange, Range(0, 1));
}

const Complex* ComplexVectorView::get_c_array() const {
  if (mrange.mstart != 0 || mrange.mstride != 1)
    throw runtime_error(
        "A ComplexVectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");

  return mdata;
}

Complex* ComplexVectorView::get_c_array() {
  if (mrange.mstart != 0 || mrange.mstride != 1)
    throw runtime_error(
        "A VectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");

  return mdata;
}

ComplexVectorView::ComplexVectorView(Complex& a) : ConstComplexVectorView(a) {
  // Nothing to do here.
}

ComplexVectorView::ComplexVectorView(Complex* data, const Range& range)
    : ConstComplexVectorView(data, range) {
  // Nothing to do here.
}

ComplexVectorView::ComplexVectorView(Complex* data,
                                     const Range& p,
                                     const Range& n)
    : ConstComplexVectorView(data, p, n) {
  // Nothing to do here.
}

void copy(ConstComplexIterator1D origin,
          const ConstComplexIterator1D& end,
          ComplexIterator1D target) {
  if (origin.mstride == 1 && target.mstride == 1)
    memcpy((void*)target.mx,
           (void*)origin.mx,
           sizeof(Complex) * (end.mx - origin.mx));
  else
    for (; origin != end; ++origin, ++target) *target = *origin;
}

void copy(Complex x, ComplexIterator1D target, const ComplexIterator1D& end) {
  for (; target != end; ++target) *target = x;
}

// Functions for Vector:
// ---------------------

ComplexVector::ComplexVector() {
  // Nothing to do here
}

ComplexVector::ComplexVector(Index n)
    : ComplexVectorView(new Complex[n], Range(0, n)) {
  // Nothing to do here.
}

ComplexVector::ComplexVector(Index n, Complex fill)
    : ComplexVectorView(new Complex[n], Range(0, n)) {
  // Here we can access the raw memory directly, for slightly
  // increased efficiency:
  const Complex* stop = mdata + n;
  for (Complex* x = mdata; x < stop; ++x) *x = fill;
}

ComplexVector::ComplexVector(Index n, Numeric fill)
    : ComplexVectorView(new Complex[n], Range(0, n)) {
  // Here we can access the raw memory directly, for slightly
  // increased efficiency:
  const Complex* stop = mdata + n;
  for (Complex* x = mdata; x < stop; ++x) *x = fill;
}

ComplexVector::ComplexVector(Complex start, Index extent, Complex stride)
    : ComplexVectorView(new Complex[extent], Range(0, extent)) {
  // Fill with values:
  Complex x = start;
  ComplexIterator1D i = begin();
  const ComplexIterator1D e = end();
  for (; i != e; ++i) {
    *i = x;
    x += stride;
  }
}

ComplexVector::ComplexVector(Numeric start, Index extent, Complex stride)
    : ComplexVectorView(new Complex[extent], Range(0, extent)) {
  // Fill with values:
  Complex x = start;
  ComplexIterator1D i = begin();
  const ComplexIterator1D e = end();
  for (; i != e; ++i) {
    *i = x;
    x += stride;
  }
}

ComplexVector::ComplexVector(Complex start, Index extent, Numeric stride)
    : ComplexVectorView(new Complex[extent], Range(0, extent)) {
  // Fill with values:
  Complex x = start;
  ComplexIterator1D i = begin();
  const ComplexIterator1D e = end();
  for (; i != e; ++i) {
    *i = x;
    x += stride;
  }
}

ComplexVector::ComplexVector(Numeric start, Index extent, Numeric stride)
    : ComplexVectorView(new Complex[extent], Range(0, extent)) {
  // Fill with values:
  Complex x = start;
  ComplexIterator1D i = begin();
  const ComplexIterator1D e = end();
  for (; i != e; ++i) {
    *i = x;
    x += stride;
  }
}

ComplexVector::ComplexVector(const ConstComplexVectorView& v)
    : ComplexVectorView(new Complex[v.nelem()], Range(0, v.nelem())) {
  copy(v.begin(), v.end(), begin());
}

ComplexVector::ComplexVector(const ComplexVector& v)
    : ComplexVectorView(new Complex[v.nelem()], Range(0, v.nelem())) {
  copy(v.begin(), v.end(), begin());
}

ComplexVector::ComplexVector(const std::vector<Complex>& v)
    : ComplexVectorView(new Complex[v.size()], Range(0, v.size())) {
  std::vector<Complex>::const_iterator vec_it_end = v.end();
  ComplexIterator1D this_it = this->begin();
  for (std::vector<Complex>::const_iterator vec_it = v.begin();
       vec_it != vec_it_end;
       ++vec_it, ++this_it)
    *this_it = *vec_it;
}

ComplexVector::ComplexVector(const std::vector<Numeric>& v)
    : ComplexVectorView(new Complex[v.size()], Range(0, v.size())) {
  std::vector<Numeric>::const_iterator vec_it_end = v.end();
  ComplexIterator1D this_it = this->begin();
  for (std::vector<Numeric>::const_iterator vec_it = v.begin();
       vec_it != vec_it_end;
       ++vec_it, ++this_it)
    *this_it = *vec_it;
}


ComplexVector& ComplexVector::operator=(ComplexVector v) {
  swap(*this, v);
  return *this;
}


ComplexVector& ComplexVector::operator=(const Array<Complex>& x) {
  resize(x.nelem());
  ComplexVectorView::operator=(x);
  return *this;
}

ComplexVector& ComplexVector::operator=(Complex x) {
  ComplexVectorView::operator=(x);
  return *this;
}

// /** Assignment operator from VectorView. Assignment operators are not
//     inherited. */
// inline Vector& Vector::operator=(const VectorView v)
// {
//   cout << "Assigning Vector from Vector View.\n";
//  // Check that sizes are compatible:
//   assert(mrange.mextent==v.mrange.mextent);
//   VectorView::operator=(v);
//   return *this;
// }

void ComplexVector::resize(Index n) {
  assert(0 <= n);
  if (mrange.mextent != n) {
    delete[] mdata;
    mdata = new Complex[n];
    mrange.mstart = 0;
    mrange.mextent = n;
    mrange.mstride = 1;
  }
}

void swap(ComplexVector& v1, ComplexVector& v2) {
  std::swap(v1.mrange, v2.mrange);
  std::swap(v1.mdata, v2.mdata);
}

ComplexVector::~ComplexVector() { delete[] mdata; }

// Functions for ConstMatrixView:
// ------------------------------

// bool ConstComplexMatrixView::empty() const
// {
//     return (nrows() == 0 || ncols() == 0);
// }

// Index ConstComplexMatrixView::nrows() const
// {
//   return mrr.mextent;
// }

// Index ConstComplexMatrixView::ncols() const
// {
//   return mcr.mextent;
// }

ConstComplexMatrixView ConstComplexMatrixView::operator()(
    const Range& r, const Range& c) const {
  return ConstComplexMatrixView(mdata, mrr, mcr, r, c);
}

ConstComplexVectorView ConstComplexMatrixView::operator()(const Range& r,
                                                          Index c) const {
  // Check that c is valid:
  assert(0 <= c);
  assert(c < mcr.mextent);

  return ConstComplexVectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
}

ConstComplexVectorView ConstComplexMatrixView::operator()(
    Index r, const Range& c) const {
  // Check that r is valid:
  assert(0 <= r);
  assert(r < mrr.mextent);

  return ConstComplexVectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
}

ConstComplexIterator2D ConstComplexMatrixView::begin() const {
  return ConstComplexIterator2D(ConstComplexVectorView(mdata + mrr.mstart, mcr),
                                mrr.mstride);
}

ConstComplexIterator2D ConstComplexMatrixView::end() const {
  return ConstComplexIterator2D(
      ConstComplexVectorView(mdata + mrr.mstart + (mrr.mextent) * mrr.mstride,
                             mcr),
      mrr.mstride);
}


ConstComplexVectorView ConstComplexMatrixView::diagonal() const {
  Index n = std::min(mrr.mextent, mcr.mextent);
  return ConstComplexVectorView(mdata + mrr.mstart + mcr.mstart,
                                Range(0, n, mrr.mstride + mcr.mstride));
}

ConstComplexMatrixView::ConstComplexMatrixView(Complex* data,
                                               const Range& rr,
                                               const Range& cr)
    : mrr(rr), mcr(cr), mdata(data) {
  // Nothing to do here.
}

ConstComplexMatrixView::ConstComplexMatrixView(Complex* data,
                                               const Range& pr,
                                               const Range& pc,
                                               const Range& nr,
                                               const Range& nc)
    : mrr(pr, nr), mcr(pc, nc), mdata(data) {
  // Nothing to do here.
}

std::ostream& operator<<(std::ostream& os, const ConstComplexMatrixView& v) {
  // Row iterators:
  ConstComplexIterator2D ir = v.begin();
  const ConstComplexIterator2D end_row = v.end();

  if (ir != end_row) {
    ConstComplexIterator1D ic = ir->begin();
    const ConstComplexIterator1D end_col = ir->end();

    if (ic != end_col) {
      os << setw(3) << *ic;
      ++ic;
    }
    for (; ic != end_col; ++ic) {
      os << " " << setw(3) << *ic;
    }
    ++ir;
  }
  for (; ir != end_row; ++ir) {
    ConstComplexIterator1D ic = ir->begin();
    const ConstComplexIterator1D end_col = ir->end();

    os << "\n";
    if (ic != end_col) {
      os << setw(3) << *ic;
      ++ic;
    }
    for (; ic != end_col; ++ic) {
      os << " " << setw(3) << *ic;
    }
  }

  return os;
}

// Functions for ComplexMatrixView:
// -------------------------

ComplexMatrixView ComplexMatrixView::operator()(const Range& r,
                                                const Range& c) {
  return ComplexMatrixView(mdata, mrr, mcr, r, c);
}

ComplexVectorView ComplexMatrixView::operator()(const Range& r, Index c) {
  // Check that c is valid:
  assert(0 <= c);
  assert(c < mcr.mextent);

  return ComplexVectorView(mdata + mcr.mstart + c * mcr.mstride, mrr, r);
}

ComplexVectorView ComplexMatrixView::operator()(Index r, const Range& c) {
  // Check that r is valid:
  assert(0 <= r);
  assert(r < mrr.mextent);

  return ComplexVectorView(mdata + mrr.mstart + r * mrr.mstride, mcr, c);
}

ComplexIterator2D ComplexMatrixView::begin() {
  return ComplexIterator2D(ComplexVectorView(mdata + mrr.mstart, mcr),
                           mrr.mstride);
}

ComplexIterator2D ComplexMatrixView::end() {
  return ComplexIterator2D(
      ComplexVectorView(mdata + mrr.mstart + (mrr.mextent) * mrr.mstride, mcr),
      mrr.mstride);
}


ComplexVectorView ComplexMatrixView::diagonal() {
  Index n = std::min(mrr.mextent, mcr.mextent);
  return ComplexVectorView(mdata + mrr.mstart + mcr.mstart,
                           Range(0, n, mrr.mstride + mcr.mstride));
}

ComplexMatrixView& ComplexMatrixView::operator=(
    const ConstComplexMatrixView& m) {
  // Check that sizes are compatible:
  assert(mrr.mextent == m.mrr.mextent);
  assert(mcr.mextent == m.mcr.mextent);

  copy(m.begin(), m.end(), begin());
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator=(const ComplexMatrixView& m) {
  // Check that sizes are compatible:
  assert(mrr.mextent == m.mrr.mextent);
  assert(mcr.mextent == m.mcr.mextent);

  copy(m.begin(), m.end(), begin());
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator=(const ComplexMatrix& m) {
  // Check that sizes are compatible:
  assert(mrr.mextent == m.mrr.mextent);
  assert(mcr.mextent == m.mcr.mextent);

  copy(m.begin(), m.end(), begin());
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator=(
    const ConstComplexVectorView& v) {
  // Check that sizes are compatible:
  assert(mrr.mextent == v.nelem());
  assert(mcr.mextent == 1);
  //  dummy = ConstComplexMatrixView(v.mdata,v.mrange,Range(v.mrange.mstart,1));;
  ConstComplexMatrixView dummy(v);
  copy(dummy.begin(), dummy.end(), begin());
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator=(Complex x) {
  copy(x, begin(), end());
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator*=(Complex x) {
  const ComplexIterator2D er = end();
  for (ComplexIterator2D r = begin(); r != er; ++r) {
    const ComplexIterator1D ec = r->end();
    for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c *= x;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator*=(Numeric x) {
  const ComplexIterator2D er = end();
  for (ComplexIterator2D r = begin(); r != er; ++r) {
    const ComplexIterator1D ec = r->end();
    for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c *= x;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator/=(Complex x) {
  const ComplexIterator2D er = end();
  for (ComplexIterator2D r = begin(); r != er; ++r) {
    const ComplexIterator1D ec = r->end();
    for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c /= x;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator/=(Numeric x) {
  const ComplexIterator2D er = end();
  for (ComplexIterator2D r = begin(); r != er; ++r) {
    const ComplexIterator1D ec = r->end();
    for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c /= x;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator+=(Complex x) {
  const ComplexIterator2D er = end();
  for (ComplexIterator2D r = begin(); r != er; ++r) {
    const ComplexIterator1D ec = r->end();
    for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c += x;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator+=(Numeric x) {
  const ComplexIterator2D er = end();
  for (ComplexIterator2D r = begin(); r != er; ++r) {
    const ComplexIterator1D ec = r->end();
    for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c += x;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator-=(Complex x) {
  const ComplexIterator2D er = end();
  for (ComplexIterator2D r = begin(); r != er; ++r) {
    const ComplexIterator1D ec = r->end();
    for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c -= x;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator-=(Numeric x) {
  const ComplexIterator2D er = end();
  for (ComplexIterator2D r = begin(); r != er; ++r) {
    const ComplexIterator1D ec = r->end();
    for (ComplexIterator1D c = r->begin(); c != ec; ++c) *c -= x;
  }
  return *this;
}

const Complex* ComplexMatrixView::get_c_array() const {
  if (mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
      mcr.mstride != 1)
    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");

  return mdata;
}

Complex* ComplexMatrixView::get_c_array() {
  if (mrr.mstart != 0 || mrr.mstride != mcr.mextent || mcr.mstart != 0 ||
      mcr.mstride != 1)
    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");

  return mdata;
}

ComplexMatrixView& ComplexMatrixView::operator*=(
    const ConstComplexMatrixView& x) {
  assert(nrows() == x.nrows());
  assert(ncols() == x.ncols());
  ConstComplexIterator2D sr = x.begin();
  ComplexIterator2D r = begin();
  const ComplexIterator2D er = end();
  for (; r != er; ++r, ++sr) {
    ConstComplexIterator1D sc = sr->begin();
    ComplexIterator1D c = r->begin();
    const ComplexIterator1D ec = r->end();
    for (; c != ec; ++c, ++sc) *c *= *sc;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator*=(const ConstMatrixView& x) {
  assert(nrows() == x.nrows());
  assert(ncols() == x.ncols());
  ConstIterator2D sr = x.begin();
  ComplexIterator2D r = begin();
  const ComplexIterator2D er = end();
  for (; r != er; ++r, ++sr) {
    ConstIterator1D sc = sr->begin();
    ComplexIterator1D c = r->begin();
    const ComplexIterator1D ec = r->end();
    for (; c != ec; ++c, ++sc) *c *= *sc;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator/=(
    const ConstComplexMatrixView& x) {
  assert(nrows() == x.nrows());
  assert(ncols() == x.ncols());
  ConstComplexIterator2D sr = x.begin();
  ComplexIterator2D r = begin();
  const ComplexIterator2D er = end();
  for (; r != er; ++r, ++sr) {
    ConstComplexIterator1D sc = sr->begin();
    ComplexIterator1D c = r->begin();
    const ComplexIterator1D ec = r->end();
    for (; c != ec; ++c, ++sc) *c /= *sc;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator/=(const ConstMatrixView& x) {
  assert(nrows() == x.nrows());
  assert(ncols() == x.ncols());
  ConstIterator2D sr = x.begin();
  ComplexIterator2D r = begin();
  const ComplexIterator2D er = end();
  for (; r != er; ++r, ++sr) {
    ConstIterator1D sc = sr->begin();
    ComplexIterator1D c = r->begin();
    const ComplexIterator1D ec = r->end();
    for (; c != ec; ++c, ++sc) *c /= *sc;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator+=(
    const ConstComplexMatrixView& x) {
  assert(nrows() == x.nrows());
  assert(ncols() == x.ncols());
  ConstComplexIterator2D sr = x.begin();
  ComplexIterator2D r = begin();
  const ComplexIterator2D er = end();
  for (; r != er; ++r, ++sr) {
    ConstComplexIterator1D sc = sr->begin();
    ComplexIterator1D c = r->begin();
    const ComplexIterator1D ec = r->end();
    for (; c != ec; ++c, ++sc) *c += *sc;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator+=(const ConstMatrixView& x) {
  assert(nrows() == x.nrows());
  assert(ncols() == x.ncols());
  ConstIterator2D sr = x.begin();
  ComplexIterator2D r = begin();
  const ComplexIterator2D er = end();
  for (; r != er; ++r, ++sr) {
    ConstIterator1D sc = sr->begin();
    ComplexIterator1D c = r->begin();
    const ComplexIterator1D ec = r->end();
    for (; c != ec; ++c, ++sc) *c += *sc;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator-=(
    const ConstComplexMatrixView& x) {
  assert(nrows() == x.nrows());
  assert(ncols() == x.ncols());
  ConstComplexIterator2D sr = x.begin();
  ComplexIterator2D r = begin();
  const ComplexIterator2D er = end();
  for (; r != er; ++r, ++sr) {
    ConstComplexIterator1D sc = sr->begin();
    ComplexIterator1D c = r->begin();
    const ComplexIterator1D ec = r->end();
    for (; c != ec; ++c, ++sc) *c -= *sc;
  }
  return *this;
}

ComplexMatrixView& ComplexMatrixView::operator-=(const ConstMatrixView& x) {
  assert(nrows() == x.nrows());
  assert(ncols() == x.ncols());
  ConstIterator2D sr = x.begin();
  ComplexIterator2D r = begin();
  const ComplexIterator2D er = end();
  for (; r != er; ++r, ++sr) {
    ConstIterator1D sc = sr->begin();
    ComplexIterator1D c = r->begin();
    const ComplexIterator1D ec = r->end();
    for (; c != ec; ++c, ++sc) *c -= *sc;
  }
  return *this;
}

ComplexMatrixView::ComplexMatrixView() : ConstComplexMatrixView() {
  // Nothing to do here.
}

ComplexMatrixView::ComplexMatrixView(Complex* data,
                                     const Range& rr,
                                     const Range& cr)
    : ConstComplexMatrixView(data, rr, cr) {
  // Nothing to do here.
}

ComplexMatrixView::ComplexMatrixView(Complex* data,
                                     const Range& pr,
                                     const Range& pc,
                                     const Range& nr,
                                     const Range& nc)
    : ConstComplexMatrixView(data, pr, pc, nr, nc) {
  // Nothing to do here.
}

void copy(ConstComplexIterator2D origin,
          const ConstComplexIterator2D& end,
          ComplexIterator2D target) {
  for (; origin != end; ++origin, ++target) {
    ConstComplexIterator1D o = origin->begin();
    const ConstComplexIterator1D e = origin->end();
    ComplexIterator1D t = target->begin();
    for (; o != e; ++o, ++t) *t = *o;
  }
}

void copy(Complex x, ComplexIterator2D target, const ComplexIterator2D& end) {
  for (; target != end; ++target) {
    ComplexIterator1D t = target->begin();
    const ComplexIterator1D e = target->end();
    for (; t != e; ++t) *t = x;
  }
}

void copy(Numeric x, ComplexIterator2D target, const ComplexIterator2D& end) {
  for (; target != end; ++target) {
    ComplexIterator1D t = target->begin();
    const ComplexIterator1D e = target->end();
    for (; t != e; ++t) *t = x;
  }
}

// Functions for ComplexMatrix:
// ---------------------

ComplexMatrix::ComplexMatrix(Index r, Index c)
    : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
  // Nothing to do here.
}

ComplexMatrix::ComplexMatrix(Index r, Index c, Complex fill)
    : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
  // Here we can access the raw memory directly, for slightly
  // increased efficiency:
  const Complex* stop = mdata + r * c;
  for (Complex* x = mdata; x < stop; ++x) *x = fill;
}

ComplexMatrix::ComplexMatrix(Index r, Index c, Numeric fill)
    : ComplexMatrixView(new Complex[r * c], Range(0, r, c), Range(0, c)) {
  // Here we can access the raw memory directly, for slightly
  // increased efficiency:
  const Complex* stop = mdata + r * c;
  for (Complex* x = mdata; x < stop; ++x) *x = fill;
}

ComplexMatrix::ComplexMatrix(const ConstComplexMatrixView& m)
    : ComplexMatrixView(new Complex[m.nrows() * m.ncols()],
                        Range(0, m.nrows(), m.ncols()),
                        Range(0, m.ncols())) {
  copy(m.begin(), m.end(), begin());
}

ComplexMatrix::ComplexMatrix(const ComplexMatrix& m)
    : ComplexMatrixView(new Complex[m.nrows() * m.ncols()],
                        Range(0, m.nrows(), m.ncols()),
                        Range(0, m.ncols())) {
  // There is a catch here: If m is an empty matrix, then it will have
  // 0 colunns. But this is used to initialize the stride of the row
  // Range! Thus, this method has to be consistent with the behaviour
  // of Range::Range. For now, Range::Range allows also stride 0.
  copy(m.begin(), m.end(), begin());
}


ComplexMatrix& ComplexMatrix::operator=(ComplexMatrix m) {
  swap(*this, m);
  return *this;
}

ComplexMatrix& ComplexMatrix::operator=(Complex x) {
  copy(x, begin(), end());
  return *this;
}


ComplexMatrix& ComplexMatrix::operator=(const ConstComplexVectorView& v) {
  resize(v.nelem(), 1);
  ConstComplexMatrixView dummy(v);
  copy(dummy.begin(), dummy.end(), begin());
  return *this;
}

void ComplexMatrix::resize(Index r, Index c) {
  assert(0 <= r);
  assert(0 <= c);

  if (mrr.mextent != r || mcr.mextent != c) {
    delete[] mdata;
    mdata = new Complex[r * c];

    mrr.mstart = 0;
    mrr.mextent = r;
    mrr.mstride = c;

    mcr.mstart = 0;
    mcr.mextent = c;
    mcr.mstride = 1;
  }
}

void swap(ComplexMatrix& m1, ComplexMatrix& m2) {
  std::swap(m1.mrr, m2.mrr);
  std::swap(m1.mcr, m2.mcr);
  std::swap(m1.mdata, m2.mdata);
}

ComplexMatrix::~ComplexMatrix() {
  //   cout << "Destroying a Matrix:\n"
  //        << *this << "\n........................................\n";
  delete[] mdata;
}

ComplexMatrix& ComplexMatrix::inv(const lapack_help::Inverse<Complex>& help)
{
  assert(ncols() == nrows());
  
  // help's internal variables are all mutable, so const is just for default parameter and to not use copies
  help.resize_if_smaller(int(ncols()));
  
  // Compute LU decomposition using LAPACK dgetrf_.
  int info;
  lapack::zgetrf_(help.size(), help.size(), mdata, help.size(), help.ipivdata(), &info);
  lapack::zgetri_(help.size(), mdata, help.size(), help.ipivdata(), help.workdata(), help.lsize(), &info);
  
  // Check for success.
  if (info not_eq 0) {
    throw runtime_error("Error inverting matrix: Matrix not of full rank.");
  }
  
  return *this;
}

ConstComplexMatrixView transpose(ConstComplexMatrixView m) {
  return ConstComplexMatrixView(m.mdata, m.mcr, m.mrr);
}

ComplexMatrixView transpose(ComplexMatrixView m) {
  return ComplexMatrixView(m.mdata, m.mcr, m.mrr);
}

ComplexMatrixView transpose(ComplexVector v) {
  return transpose((ComplexMatrixView)v);
}

ComplexVectorView& ComplexVectorView::operator=(const Array<Complex>& v) {
  //  cout << "Assigning VectorView from Array<Numeric>.\n";

  // Check that sizes are compatible:
  assert(mrange.mextent == v.nelem());

  // Iterators for Array:
  Array<Complex>::const_iterator i = v.begin();
  const Array<Complex>::const_iterator e = v.end();
  // Iterator for Vector:
  ComplexIterator1D target = begin();

  for (; i != e; ++i, ++target) *target = *i;

  return *this;
}

//Functions operating on the complex vectors

Complex operator*(const ConstComplexVectorView& a,
                  const ConstComplexVectorView& b) {
  // Check dimensions:
  assert(a.nelem() == b.nelem());

  const ConstComplexIterator1D ae = a.end();
  ConstComplexIterator1D ai = a.begin();
  ConstComplexIterator1D bi = b.begin();

  Complex res = 0;
  for (; ai != ae; ++ai, ++bi) res += (*ai) * (*bi);

  return res;
}


void mult(ComplexVectorView y,
          const ConstComplexMatrixView& M,
          const ConstComplexVectorView& x) {
  assert(x.nelem() == M.nrows());
  assert(y.nelem() == M.ncols());

  ComplexMatrixViewMap eigen_y = MapToEigenRow(y);
  if (y.mdata == x.mdata)
    eigen_y = MapToEigen(M) * MapToEigenRow(x);
  else
    eigen_y.noalias() = MapToEigen(M) * MapToEigenRow(x);
}


void mult(ComplexMatrixView A,
          const ConstComplexMatrixView& B,
          const ConstComplexMatrixView& C) {
  assert(B.nrows() == A.nrows());
  assert(C.ncols() == A.ncols());
  assert(B.ncols() == C.nrows());

  ComplexMatrixViewMap eigen_A = MapToEigen(A);
  if (A.mdata == B.mdata || A.mdata == C.mdata)
    eigen_A = MapToEigen(B) * MapToEigen(C);
  else
    eigen_A.noalias() = MapToEigen(B) * MapToEigen(C);
}
void mult(ComplexMatrixView A,
          const ConstComplexMatrixView& B,
          const ConstMatrixView& C) {
  assert(B.nrows() == A.nrows());
  assert(C.ncols() == A.ncols());
  assert(B.ncols() == C.nrows());

  ComplexMatrixViewMap eigen_A = MapToEigen(A);
  if (A.mdata == B.mdata)
    eigen_A = MapToEigen(B) * MapToEigen(C);
  else
    eigen_A.noalias() = MapToEigen(B) * MapToEigen(C);
}
void mult(ComplexMatrixView A,
          const ConstMatrixView& B,
          const ConstComplexMatrixView& C) {
  assert(B.nrows() == A.nrows());
  assert(C.ncols() == A.ncols());
  assert(B.ncols() == C.nrows());

  ComplexMatrixViewMap eigen_A = MapToEigen(A);
  if (A.mdata == C.mdata)
    eigen_A = MapToEigen(B) * MapToEigen(C);
  else
    eigen_A.noalias() = MapToEigen(B) * MapToEigen(C);
}
void mult(ComplexMatrixView A,
          const ConstMatrixView& B,
          const ConstMatrixView& C) {
  assert(B.nrows() == A.nrows());
  assert(C.ncols() == A.ncols());
  assert(B.ncols() == C.nrows());

  ComplexMatrixViewMap eigen_A = MapToEigen(A);
  eigen_A.noalias() = MapToEigen(B) * MapToEigen(C);
}

// Converts constant matrix to constant eigen map
ComplexConstMatrixViewMap MapToEigen(const ConstComplexMatrixView& A) {
  return ComplexConstMatrixViewMap(
      A.mdata + A.mrr.get_start() + A.mcr.get_start(),
      A.nrows(),
      A.ncols(),
      ComplexStrideType(A.mrr.get_stride(), A.mcr.get_stride()));
}

// Converts constant vector to constant eigen row-view
ComplexConstMatrixViewMap MapToEigen(const ConstComplexVectorView& A) {
  return ComplexConstMatrixViewMap(A.mdata + A.mrange.get_start(),
                                   A.nelem(),
                                   1,
                                   ComplexStrideType(A.mrange.get_stride(), 1));
}

// Converts constant vector to constant eigen row-view
ComplexConstMatrixViewMap MapToEigenRow(const ConstComplexVectorView& A) {
  return MapToEigen(A);
}

// Converts constant vector to constant eigen column-view
ComplexConstMatrixViewMap MapToEigenCol(const ConstComplexVectorView& A) {
  return ComplexConstMatrixViewMap(A.mdata + A.mrange.get_start(),
                                   1,
                                   A.nelem(),
                                   ComplexStrideType(1, A.mrange.get_stride()));
}

// Converts matrix to eigen map
ComplexMatrixViewMap MapToEigen(ComplexMatrixView& A) {
  return ComplexMatrixViewMap(
      A.mdata + A.mrr.get_start() + A.mcr.get_start(),
      A.nrows(),
      A.ncols(),
      ComplexStrideType(A.mrr.get_stride(), A.mcr.get_stride()));
}

// Converts vector to eigen map row-view
ComplexMatrixViewMap MapToEigen(ComplexVectorView& A) {
  return ComplexMatrixViewMap(A.mdata + A.mrange.get_start(),
                              A.nelem(),
                              1,
                              ComplexStrideType(A.mrange.get_stride(), 1));
}

// Converts vector to eigen map row-view
ComplexMatrixViewMap MapToEigenRow(ComplexVectorView& A) {
  return MapToEigen(A);
}

// Converts vector to eigen map column-view
ComplexMatrixViewMap MapToEigenCol(ComplexVectorView& A) {
  return ComplexMatrixViewMap(A.mdata + A.mrange.get_start(),
                              1,
                              A.nelem(),
                              ComplexStrideType(1, A.mrange.get_stride()));
}

// Helper function for debugging
#ifndef NDEBUG

Complex debug_matrixview_get_elem(ComplexMatrixView& mv, Index r, Index c) {
  return mv(r, c);
}

#endif