double zero_;

template<typename _Tp>
inline void
swap(_Tp& __a, _Tp& __b)
{
  _Tp __tmp = __a;
  __a = __b;
  __b = __tmp;
}

template<typename _Tp>
inline const _Tp&
min(const _Tp& __a, const _Tp& __b)
{
  if (__b < __a) return __b; return __a;
}

template<typename _Tp>
inline const _Tp&
max(const _Tp& __a, const _Tp& __b)
{
  if (__a < __b) return __b; return __a;
}

template<typename _InputIter, typename _OutputIter>
inline _OutputIter
copy(_InputIter __first, _InputIter __last,
     _OutputIter __result)
{
  for ( ; __first != __last; ++__result, ++__first)
    *__result = *__first;
  return __result;
}

template<typename _ForwardIter, typename _Tp>
void
fill(_ForwardIter __first, _ForwardIter __last, const _Tp& __value)
{
  for ( ; __first != __last; ++__first)
    *__first = __value;
}


struct dummy {
    void foo () {}
};

template<typename C>
struct container_reference {
    inline container_reference (): c_ (0) {}
    inline container_reference (C &c): c_ (&c) {}
    inline C &operator () () const { return *c_; }

    C *c_;
};

template <typename=void> struct unbounded_array {
    inline unbounded_array (): data_ (new double [9]) {}
    inline double& operator [] (unsigned int i) { return data_ [i]; }
    double* data_;
};

template<typename T1, class T2>
struct scalar_swap {
    inline void operator () (T1& t1, T2& t2) const {
      swap (t1, t2);
    }
};

struct row_major {
    static inline unsigned int element (unsigned int i, unsigned int size1, unsigned int j, unsigned int size2) {
      return i * size2 + j;
    }
};

template<typename E>
struct matrix_expression {
    inline E &operator () () { return *static_cast<E *> (this); }
};

template<typename E>
struct matrix_reference {
    inline matrix_reference (E &e): e_ (e) {}

    inline unsigned int size1 () const { return e_.size1 (); }
    inline unsigned int size2 () const { return e_.size2 (); }

    inline double& operator () (unsigned int i, unsigned int j) { return e_(i, j); }

    inline typename E::iterator1 find1 (int rank, unsigned int i, unsigned int j) {
      return e_.find1 (rank, i, j);
    }

    inline typename E::iterator2 find2 (int rank, unsigned int i, unsigned int j) {
      return e_.find2 (rank, i, j);
    }

    E& e_;
};

template<typename F1, class M, class E>
void matrix_swap (const F1 &f1, M &m, matrix_expression<E> &e) {
  typename M::iterator1 it1 (m.begin1 ());
  typename E::iterator1 it1e (e ().begin1 ());
  int size1 ((m.end1 () - it1));
  while (-- size1 >= 0) {

    typename M::iterator2 it2 (it1.begin ());
    typename E::iterator2 it2e (it1e.begin ());
    int size2 ((it1.end () - it2));

    F1 () (*it2, *it2e);
  }
}

template <typename=void>
struct banded_matrix {
    inline banded_matrix ():
		    size1_ (3), size2_ (3),
		    lower_ (1), upper_ (1) {}

    inline unsigned int size1 () const { return size1_; }
    inline unsigned int size2 () const { return size2_; }
    inline unsigned int lower () const { return lower_; }
    inline unsigned int upper () const { return upper_; }
    inline unbounded_array<> &data () { return data_; }

    inline double& operator () (unsigned int i, unsigned int j) {
      unsigned int k = j;
      unsigned int l = upper_ + i - j;
      if (k < size2_ &&
	  l < lower_ + 1 + upper_)
	return data () [row_major::element (k, size2_,
					    l, lower_ + 1 + upper_)];
      dummy ().foo ();
      return zero_;
    }

    struct iterator1;
    struct iterator2;

    inline iterator1 find1 (int rank, unsigned int i, unsigned int j) {
      if (rank == 1) {
	unsigned int upper_i = min (j + 1 + lower_, size1_);
	i = min (i, upper_i);
      }
      return iterator1 (*this, i, j);
    }

    inline iterator2 find2 (int rank, unsigned int i, unsigned int j) {
      if (rank == 1) {
	unsigned int lower_j = max (int (i - lower_), int (0));
	j = max (j, lower_j);
	unsigned int upper_j = min (i + 1 + upper_, size2_);
	j = min (j, upper_j);
      }
      return iterator2 (*this, i, j);
    }

    struct iterator1 : container_reference<banded_matrix> {
        inline iterator1 (banded_matrix &m, unsigned int it1, unsigned int it2):
			container_reference<banded_matrix> (m), it1_ (it1), it2_ (it2) {}

        inline int operator - (const iterator1 &it) const {
	  return it1_ - it.it1_;
        }

        unsigned int it1_;
        unsigned int it2_;
    };

    struct iterator2 : container_reference<banded_matrix> {
        inline iterator2 (banded_matrix &m, unsigned int it1, unsigned int it2):
			container_reference<banded_matrix> (m), it1_ (it1), it2_ (it2) {}

        inline iterator2 &operator ++ () {
	  return *this;
        }

        inline int operator - (const iterator2 &it) const {
	  return it2_ - it.it2_;
        }

        inline double& operator * () const {
	  return (*this) () (it1_, it2_);
        }

        inline unsigned int index1 () const { return it1_; }
        inline unsigned int index2 () const { return it2_; }

        unsigned int it1_;
        unsigned int it2_;
    };

    unsigned int size1_;
    unsigned int size2_;
    unsigned int lower_;
    unsigned int upper_;
    unbounded_array<> data_;
};

template<typename M>
struct banded_adaptor : matrix_expression<banded_adaptor<M> > {
    typedef M matrix_type;
    typedef matrix_reference<banded_matrix<> > matrix_closure_type;
    typedef typename M::iterator1 iterator1_type;
    typedef typename M::iterator2 iterator2_type;

    inline banded_adaptor (matrix_type &data):
		    data_ (data), lower_ (1), upper_ (1) {}

    inline unsigned int size1 () const     { return data_.size1 (); }
    inline unsigned int size2 () const     { return data_.size2 (); }
    inline unsigned int lower () const     { return 1; }
    inline unsigned int upper () const     { return upper_; }
    inline matrix_closure_type &data () { return data_; }

    inline double& operator () (unsigned int i, unsigned int j) {
      unsigned int k = j;
      unsigned int l = upper_ + i - j;
      if (k < size2 () &&
	  l < lower_ + 1 + upper_)
	return data () (i, j);

      dummy ().foo ();
      return zero_;
    }

    inline void swap (banded_adaptor &m) {
      if (this != &m) {
	matrix_swap (scalar_swap<double, double> (), *this, m);
      }
    }

    struct iterator1;
    struct iterator2;

    inline iterator1 find1 (int rank, unsigned int i, unsigned int j) {
      if (rank == 1) {
	i = min (j + 1 + lower_, size1 ());
      }
      return iterator1 (*this, data ().find1 (rank, i, j));
    }

    inline iterator2 find2 (int rank, unsigned int i, unsigned int j) {
      if (rank == 1) {
	unsigned int upper_j = min (i + 1 + upper_, size2 ());
	j = min (j, upper_j);
      }
      return iterator2 (*this, data ().find2 (rank, i, j));
    }

    struct iterator1 : container_reference<banded_adaptor> {
        inline iterator1 (banded_adaptor &m, const iterator1_type &it1):
			container_reference<banded_adaptor> (m), it1_ (it1) {}

        inline int operator - (const iterator1 &it) const {
	  return it1_ - it.it1_;
        }

        inline iterator2 begin () const {
	  return (*this) ().find2 (1, index1 (), 0);
        }

        inline iterator2 end () const {
	  return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
        }

        inline unsigned int index1 () const { return it1_.it1_; }

        iterator1_type it1_;
    };

    inline iterator1 begin1 () { return find1 (0, 0, 0); }
    inline iterator1 end1 ()   { return find1 (0, size1 (), 0); }

    struct iterator2 : container_reference<banded_adaptor> {
        inline iterator2 (banded_adaptor &m, const iterator2_type &it2):
			container_reference<banded_adaptor> (m), it2_ (it2) {}

        inline iterator2 &operator ++ () { return *this; }

        inline int operator - (const iterator2 &it) const {
	  return it2_ - it.it2_;
        }

        inline double& operator * () const {
	  unsigned int i = it2_.index1 ();
	  unsigned int j = it2_.index2 ();
	  unsigned int l = (*this) ().upper () + i - j;
	  unsigned int q = (*this) ().size2 ();
	  if (j < q &&
	      l < (*this) ().lower () + 1 + (*this) ().upper ())
	    return *it2_;

	  return (*this) () (i, j);
        }

        iterator2_type it2_;
    };

    matrix_closure_type data_;
    unsigned int lower_;
    unsigned int upper_;
};

int main ()
{
  banded_matrix<> m1,m2;
  banded_adaptor<banded_matrix<> > bam1 (m1), bam2 (m2);
  bam1.swap (bam2);
}