double zero_;

template<typename T>
inline const T&
min(const T& a, const T& b)
{
  if (b < a) return b; return a;
}

struct barrier { barrier () {} };

template<typename C>
struct container_reference {
    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_;
};

inline unsigned int element (unsigned int i, unsigned int j) {
  return i + j;
}


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

    inline unsigned int size2 () const { return e_.size2 (); }
    inline double& foo (unsigned int i, unsigned int j) { return e_(i, j); }

    E& e_;
};

template <typename=void>
struct banded_matrix {
    inline banded_matrix () : size2_ (3) {}

    inline unsigned int size2 () const { return 3; }
    inline unbounded_array<> &data () { return data_; }

    inline double& operator () (unsigned int i, unsigned int j) {
      unsigned int dead1 = j;
      unsigned int dead2 = 1 + i - j;
      if (j < size2_ && i-j < 2)
	return data () [element (j,i-j+1)];
      barrier ();
      return zero_;
    }

    struct iterator1;
    struct iterator2;

    inline iterator1 find1 (unsigned int i) {
      return iterator1 (*this, i);
    }

    inline iterator2 find2 (unsigned int i) {
      return iterator2 (*this, i, min (i + 2, 3U));
    }

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

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

        unsigned int it1_;
        unsigned int fill_in;
    };

    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_;
    };

    const unsigned int size2_;
    unbounded_array<> data_;
};

template<typename M>
struct banded_adaptor {
    inline banded_adaptor (M &data) : data_ (data), 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_reference<banded_matrix<> > &data () { return data_; }

    inline double& operator () (unsigned int i, unsigned int j) {
      unsigned int dead1 = j;
      unsigned int dead2 = upper_ + i - j;
      if (j < size2 () && i-j < 1)
	return data ().foo (i, j);

      barrier ();
      return zero_;
    }

    struct iterator2;

    inline iterator2 find2 (unsigned int i) {
      return iterator2 (*this, data ().e_.find2 (i));
    }

    struct iterator1 : container_reference<banded_adaptor> {
        inline iterator1 (banded_adaptor &m, const typename M::iterator1 &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 (it1_.it1_);
        }

        inline iterator2 end () const {
	  return (*this) ().find2 (it1_.it1_);
        }

        typename M::iterator1 it1_;
    };

    inline iterator1 begin1 () { return iterator1 (*this, data ().e_.find1 (0)); }
    inline iterator1 end1 ()   { return iterator1 (*this, data ().e_.find1 (3)); }

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

        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);
        }

        typename M::iterator2 it2_;
    };

    matrix_reference<banded_matrix<> > data_;
    unsigned int upper_;
};

void matrix_swap (banded_adaptor<banded_matrix<> > &bam1,
		  banded_adaptor<banded_matrix<> > &bam2)
{
  banded_adaptor<banded_matrix<> >::iterator1 it1  = bam1.begin1 (),
					      it1e = bam2.begin1 ();
  int size1 = bam1.end1 () - it1;
  while (-- size1 >= 0) {
    banded_adaptor<banded_matrix<> >::iterator2 dummy (it1.begin());

    double x = *it1.begin();
    *it1.begin()  = *it1e.begin();
    *it1e.begin() = x;
  }
}

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