code:
//  Copyright (c) 2000-2011 Joerg Walter, Mathias Koch, David Bellot
//
//  Distributed under the Boost Software License, Version 1.0. (See
//  accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt)
//
//  The authors gratefully acknowledge the support of
//  GeNeSys mbH & Co. KG in producing this work.

#ifndef _BOOST_UBLAS_BLAS_
#define _BOOST_UBLAS_BLAS_

#include <boost/numeric/ublas/traits.hpp>

namespace boost { namespace numeric { namespace ublas {

    /** Interface and implementation of BLAS level 1
     * This includes functions which perform \b vector-vector operations.
     * More information about BLAS can be found at
     * http://en.wikipedia.org/wiki/BLAS
     */
    namespace blas_1 {

        /** 1-Norm: \f$\sum_i |x_i|\f$ (also called \f$\mathcal{L}_1\f$ or Manhattan norm)
     *
     * \param v a vector or vector expression
     * \return the 1-Norm with type of the vector's type
     *
     * \tparam V type of the vector (not needed by default)
     */
        template<class V>
        typename type_traits<typename V::value_type>::real_type
        asum (const V &v) {
            return norm_1 (v);
        }

        /** 2-Norm: \f$\sum_i |x_i|^2\f$ (also called \f$\mathcal{L}_2\f$ or Euclidean norm)
     *
     * \param v a vector or vector expression
     * \return the 2-Norm with type of the vector's type
     *
     * \tparam V type of the vector (not needed by default)
     */
        template<class V>
        typename type_traits<typename V::value_type>::real_type
        nrm2 (const V &v) {
            return norm_2 (v);
        }

        /** Infinite-norm: \f$\max_i |x_i|\f$ (also called \f$\mathcal{L}_\infty\f$ norm)
     *
     * \param v a vector or vector expression
     * \return the Infinite-Norm with type of the vector's type
     *
     * \tparam V type of the vector (not needed by default)
     */
        template<class V>
        typename type_traits<typename V::value_type>::real_type
        amax (const V &v) {
            return norm_inf (v);
        }

        /** Inner product of vectors \f$v_1\f$ and \f$v_2\f$
     *
     * \param v1 first vector of the inner product
     * \param v2 second vector of the inner product
     * \return the inner product of the type of the most generic type of the 2 vectors
     *
     * \tparam V1 type of first vector (not needed by default)
     * \tparam V2 type of second vector (not needed by default)
     */
        template<class V1, class V2>
        typename promote_traits<typename V1::value_type, typename V2::value_type>::promote_type
        dot (const V1 &v1, const V2 &v2) {
            return inner_prod (v1, v2);
        }

        /** Copy vector \f$v_2\f$ to \f$v_1\f$
     *
     * \param v1 target vector
     * \param v2 source vector
     * \return a reference to the target vector
     *
     * \tparam V1 type of first vector (not needed by default)
     * \tparam V2 type of second vector (not needed by default)
     */
        template<class V1, class V2>
        V1 & copy (V1 &v1, const V2 &v2)
    {
            return v1.assign (v2);
        }

        /** Swap vectors \f$v_1\f$ and \f$v_2\f$
     *
     * \param v1 first vector
     * \param v2 second vector
     *
         * \tparam V1 type of first vector (not needed by default)
     * \tparam V2 type of second vector (not needed by default)
     */
    template<class V1, class V2>
        void swap (V1 &v1, V2 &v2)
    {
            v1.swap (v2);
        }

        /** scale vector \f$v\f$ with scalar \f$t\f$
     *
     * \param v vector to be scaled
     * \param t the scalar
     * \return \c t*v
     *
     * \tparam V type of the vector (not needed by default)
     * \tparam T type of the scalar (not needed by default)
     */
        template<class V, class T>
        V & scal (V &v, const T &t)
    {
            return v *= t;
        }

        /** Compute \f$v_1= v_1 +  t.v_2\f$
     *
     * \param v1 target and first vector
     * \param t the scalar
     * \param v2 second vector
     * \return a reference to the first and target vector
     *
     * \tparam V1 type of the first vector (not needed by default)
     * \tparam T type of the scalar (not needed by default)
     * \tparam V2 type of the second vector (not needed by default)
     */
        template<class V1, class T, class V2>
        V1 & axpy (V1 &v1, const T &t, const V2 &v2)
    {
            return v1.plus_assign (t * v2);
        }

    /** Performs rotation of points in the plane and assign the result to the first vector
     *
     * Each point is defined as a pair \c v1(i) and \c v2(i), being respectively
     * the \f$x\f$ and \f$y\f$ coordinates. The parameters \c t1 and \t2 are respectively
     * the cosine and sine of the angle of the rotation.
     * Results are not returned but directly written into \c v1.
     *
     * \param t1 cosine of the rotation
     * \param v1 vector of \f$x\f$ values
     * \param t2 sine of the rotation
     * \param v2 vector of \f$y\f$ values
     *
     * \tparam T1 type of the cosine value (not needed by default)
     * \tparam V1 type of the \f$x\f$ vector (not needed by default)
     * \tparam T2 type of the sine value (not needed by default)
     * \tparam V2 type of the \f$y\f$ vector (not needed by default)
     */
        template<class T1, class V1, class T2, class V2>
        void rot (const T1 &t1, V1 &v1, const T2 &t2, V2 &v2)
    {
            typedef typename promote_traits<typename V1::value_type, typename V2::value_type>::promote_type promote_type;
            vector<promote_type> vt (t1 * v1 + t2 * v2);
            v2.assign (- t2 * v1 + t1 * v2);
            v1.assign (vt);
        }

    }

    /** \brief Interface and implementation of BLAS level 2
     * This includes functions which perform \b matrix-vector operations.
     * More information about BLAS can be found at
     * http://en.wikipedia.org/wiki/BLAS
     */
    namespace blas_2 {

       /** \brief multiply vector \c v with triangular matrix \c m
    *
    * \param v a vector
    * \param m a triangular matrix
    * \return the result of the product
    *
    * \tparam V type of the vector (not needed by default)
    * \tparam M type of the matrix (not needed by default)
        */
        template<class V, class M>
        V & tmv (V &v, const M &m)
    {
            return v = prod (m, v);
        }

        /** \brief solve \f$m.x = v\f$ in place, where \c m is a triangular matrix
     *
     * \param v a vector
     * \param m a matrix
     * \param C (this parameter is not needed)
     * \return a result vector from the above operation
     *
     * \tparam V type of the vector (not needed by default)
     * \tparam M type of the matrix (not needed by default)
     * \tparam C n/a
         */
        template<class V, class M, class C>
        V & tsv (V &v, const M &m, C)
    {
            return v = solve (m, v, C ());
        }

        /** \brief compute \f$ v_1 = t_1.v_1 + t_2.(m.v_2)\f$, a general matrix-vector product
     *
     * \param v1 a vector
     * \param t1 a scalar
     * \param t2 another scalar
     * \param m a matrix
     * \param v2 another vector
     * \return the vector \c v1 with the result from the above operation
     *
     * \tparam V1 type of first vector (not needed by default)
     * \tparam T1 type of first scalar (not needed by default)
     * \tparam T2 type of second scalar (not needed by default)
     * \tparam M type of matrix (not needed by default)
     * \tparam V2 type of second vector (not needed by default)
         */
        template<class V1, class T1, class T2, class M, class V2>
        V1 & gmv (V1 &v1, const T1 &t1, const T2 &t2, const M &m, const V2 &v2)
    {
            return v1 = t1 * v1 + t2 * prod (m, v2);
        }

        /** \brief Rank 1 update: \f$ m = m + t.(v_1.v_2^T)\f$
     *
     * \param m a matrix
     * \param t a scalar
     * \param v1 a vector
     * \param v2 another vector
     * \return a matrix with the result from the above operation
     *
     * \tparam M type of matrix (not needed by default)
     * \tparam T type of scalar (not needed by default)
     * \tparam V1 type of first vector (not needed by default)
     * \tparam V2type of second vector (not needed by default)
     */
        template<class M, class T, class V1, class V2>
        M & gr (M &m, const T &t, const V1 &v1, const V2 &v2)
    {
#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
            return m += t * outer_prod (v1, v2);
#else
            return m = m + t * outer_prod (v1, v2);
#endif
        }

        /** \brief symmetric rank 1 update: \f$m = m + t.(v.v^T)\f$
     *
     * \param m a matrix
     * \param t a scalar
     * \param v a vector
     * \return a matrix with the result from the above operation
     *
     * \tparam M type of matrix (not needed by default)
     * \tparam T type of scalar (not needed by default)
     * \tparam V type of vector (not needed by default)
     */
        template<class M, class T, class V>
        M & sr (M &m, const T &t, const V &v)
    {
#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
            return m += t * outer_prod (v, v);
#else
            return m = m + t * outer_prod (v, v);
#endif
        }

        /** \brief hermitian rank 1 update: \f$m = m + t.(v.v^H)\f$
     *
     * \param m a matrix
     * \param t a scalar
     * \param v a vector
     * \return a matrix with the result from the above operation
     *
     * \tparam M type of matrix (not needed by default)
     * \tparam T type of scalar (not needed by default)
     * \tparam V type of vector (not needed by default)
     */
        template<class M, class T, class V>
        M & hr (M &m, const T &t, const V &v)
    {
#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
            return m += t * outer_prod (v, conj (v));
#else
            return m = m + t * outer_prod (v, conj (v));
#endif
        }

         /** \brief symmetric rank 2 update: \f$ m=m+ t.(v_1.v_2^T + v_2.v_1^T)\f$
      *
      * \param m a matrix
      * \param t a scalar
      * \param v1 a vector
      * \param v2 another vector
      * \return a matrix with the result from the above operation
      *
      * \tparam M type of matrix (not needed by default)
      * \tparam T type of scalar (not needed by default)
      * \tparam V1 type of first vector (not needed by default)
      * \tparam V2type of second vector (not needed by default)
          */
        template<class M, class T, class V1, class V2>
        M & sr2 (M &m, const T &t, const V1 &v1, const V2 &v2)
    {
#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
            return m += t * (outer_prod (v1, v2) + outer_prod (v2, v1));
#else
            return m = m + t * (outer_prod (v1, v2) + outer_prod (v2, v1));
#endif
        }

        /** \brief hermitian rank 2 update: \f$m=m+t.(v_1.v_2^H) + v_2.(t.v_1)^H)\f$
     *
     * \param m a matrix
     * \param t a scalar
     * \param v1 a vector
     * \param v2 another vector
     * \return a matrix with the result from the above operation
     *
     * \tparam M type of matrix (not needed by default)
     * \tparam T type of scalar (not needed by default)
     * \tparam V1 type of first vector (not needed by default)
     * \tparam V2type of second vector (not needed by default)
         */
        template<class M, class T, class V1, class V2>
        M & hr2 (M &m, const T &t, const V1 &v1, const V2 &v2)
    {
#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
            return m += t * outer_prod (v1, conj (v2)) + type_traits<T>::conj (t) * outer_prod (v2, conj (v1));
#else
            return m = m + t * outer_prod (v1, conj (v2)) + type_traits<T>::conj (t) * outer_prod (v2, conj (v1));
#endif
        }

    }

    /** \brief Interface and implementation of BLAS level 3
     * This includes functions which perform \b matrix-matrix operations.
     * More information about BLAS can be found at
     * http://en.wikipedia.org/wiki/BLAS
     */
    namespace blas_3 {

        /** \brief triangular matrix multiplication \f$m_1=t.m_2.m_3\f$ where \f$m_2\f$ and \f$m_3\f$ are triangular
     *
     * \param m1 a matrix for storing result
     * \param t a scalar
     * \param m2 a triangular matrix
     * \param m3 a triangular matrix
     * \return the matrix \c m1
     *
     * \tparam M1 type of the result matrix (not needed by default)
     * \tparam T type of the scalar (not needed by default)
     * \tparam M2 type of the first triangular matrix (not needed by default)
     * \tparam M3 type of the second triangular matrix (not needed by default)
     *
        */
        template<class M1, class T, class M2, class M3>
        M1 & tmm (M1 &m1, const T &t, const M2 &m2, const M3 &m3)
    {
            return m1 = t * prod (m2, m3);
        }

        /** \brief triangular solve \f$ m_2.x = t.m_1\f$ in place, \f$m_2\f$ is a triangular matrix
     *
     * \param m1 a matrix
     * \param t a scalar
     * \param m2 a triangular matrix
     * \param C (not used)
     * \return the \f$m_1\f$ matrix
     *
     * \tparam M1 type of the first matrix (not needed by default)
     * \tparam T type of the scalar (not needed by default)
     * \tparam M2 type of the triangular matrix (not needed by default)
     * \tparam C (n/a)
         */
        template<class M1, class T, class M2, class C>
        M1 & tsm (M1 &m1, const T &t, const M2 &m2, C)
    {
            return m1 = solve (m2, t * m1, C ());
        }

        /** \brief general matrix multiplication \f$m_1=t_1.m_1 + t_2.m_2.m_3\f$
     *
     * \param m1 first matrix
     * \param t1 first scalar
     * \param t2 second scalar
     * \param m2 second matrix
     * \param m3 third matrix
     * \return the matrix \c m1
     *
     * \tparam M1 type of the first matrix (not needed by default)
     * \tparam T1 type of the first scalar (not needed by default)
     * \tparam T2 type of the second scalar (not needed by default)
     * \tparam M2 type of the second matrix (not needed by default)
     * \tparam M3 type of the third matrix (not needed by default)
         */
        template<class M1, class T1, class T2, class M2, class M3>
        M1 & gmm (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3)
    {
            return m1 = t1 * m1 + t2 * prod (m2, m3);
        }

        /** \brief symmetric rank \a k update: \f$m_1=t.m_1+t_2.(m_2.m_2^T)\f$
     *
     * \param m1 first matrix
     * \param t1 first scalar
     * \param t2 second scalar
     * \param m2 second matrix
     * \return matrix \c m1
     *
     * \tparam M1 type of the first matrix (not needed by default)
     * \tparam T1 type of the first scalar (not needed by default)
     * \tparam T2 type of the second scalar (not needed by default)
     * \tparam M2 type of the second matrix (not needed by default)
     * \todo use opb_prod()
         */
        template<class M1, class T1, class T2, class M2>
        M1 & srk (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2)
    {
            return m1 = t1 * m1 + t2 * prod (m2, trans (m2));
        }

        /** \brief hermitian rank \a k update: \f$m_1=t.m_1+t_2.(m_2.m2^H)\f$
     *
     * \param m1 first matrix
     * \param t1 first scalar
     * \param t2 second scalar
     * \param m2 second matrix
     * \return matrix \c m1
     *
     * \tparam M1 type of the first matrix (not needed by default)
     * \tparam T1 type of the first scalar (not needed by default)
     * \tparam T2 type of the second scalar (not needed by default)
     * \tparam M2 type of the second matrix (not needed by default)
     * \todo use opb_prod()
         */
        template<class M1, class T1, class T2, class M2>
        M1 & hrk (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2)
    {
            return m1 = t1 * m1 + t2 * prod (m2, herm (m2));
        }

        /** \brief generalized symmetric rank \a k update: \f$m_1=t_1.m_1+t_2.(m_2.m3^T)+t_2.(m_3.m2^T)\f$
     *
     * \param m1 first matrix
     * \param t1 first scalar
     * \param t2 second scalar
     * \param m2 second matrix
     * \param m3 third matrix
     * \return matrix \c m1
     *
     * \tparam M1 type of the first matrix (not needed by default)
     * \tparam T1 type of the first scalar (not needed by default)
     * \tparam T2 type of the second scalar (not needed by default)
     * \tparam M2 type of the second matrix (not needed by default)
     * \tparam M3 type of the third matrix (not needed by default)
     * \todo use opb_prod()
         */
        template<class M1, class T1, class T2, class M2, class M3>
        M1 & sr2k (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3)
    {
            return m1 = t1 * m1 + t2 * (prod (m2, trans (m3)) + prod (m3, trans (m2)));
        }

        /** \brief generalized hermitian rank \a k update: * \f$m_1=t_1.m_1+t_2.(m_2.m_3^H)+(m_3.(t_2.m_2)^H)\f$
     *
     * \param m1 first matrix
     * \param t1 first scalar
     * \param t2 second scalar
     * \param m2 second matrix
     * \param m3 third matrix
     * \return matrix \c m1
     *
     * \tparam M1 type of the first matrix (not needed by default)
     * \tparam T1 type of the first scalar (not needed by default)
     * \tparam T2 type of the second scalar (not needed by default)
     * \tparam M2 type of the second matrix (not needed by default)
     * \tparam M3 type of the third matrix (not needed by default)
     * \todo use opb_prod()
         */
        template<class M1, class T1, class T2, class M2, class M3>
        M1 & hr2k (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3)
    {
            return m1 =
              t1 * m1
            + t2 * prod (m2, herm (m3))
            + type_traits<T2>::conj (t2) * prod (m3, herm (m2));
        }

    }

}}}

#endif

license: mit


code:
//
//  Copyright (c) 2000-2002
//  Joerg Walter, Mathias Koch
//
//  Distributed under the Boost Software License, Version 1.0. (See
//  accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt)
//
//  The authors gratefully acknowledge the support of
//  GeNeSys mbH & Co. KG in producing this work.
//

#ifndef _BOOST_UBLAS_CONCEPTS_
#define _BOOST_UBLAS_CONCEPTS_

#include <boost/concept_check.hpp>

// Concept checks based on ideas of Jeremy Siek

namespace boost { namespace numeric { namespace ublas {

    template<class I>
    struct Indexed1DIteratorConcept {
        typedef I iterator_type;

        void constraints () {
            iterator_type it = iterator_type ();
            // Index
            it.index ();
        }
    };

    template<class I>
    struct IndexedBidirectional1DIteratorConcept {
        typedef I iterator_type;

        void constraints () {
            function_requires< BidirectionalIteratorConcept<iterator_type> >();
            function_requires< Indexed1DIteratorConcept<iterator_type> >();
        }
    };

    template<class I>
    struct Mutable_IndexedBidirectional1DIteratorConcept {
        typedef I iterator_type;

        void constraints () {
            function_requires< Mutable_BidirectionalIteratorConcept<iterator_type> >();
            function_requires< Indexed1DIteratorConcept<iterator_type> >();
        }
    };

    template<class I>
    struct IndexedRandomAccess1DIteratorConcept {
        typedef I iterator_type;

        void constraints () {
            function_requires< RandomAccessIteratorConcept<iterator_type> >();
            function_requires< Indexed1DIteratorConcept<iterator_type> >();
        }
    };

    template<class I>
    struct Mutable_IndexedRandomAccess1DIteratorConcept {
        typedef I iterator_type;

        void constraints () {
            function_requires< Mutable_RandomAccessIteratorConcept<iterator_type> >();
            function_requires< Indexed1DIteratorConcept<iterator_type> >();
        }
    };

    template<class I>
    struct Indexed2DIteratorConcept {
        typedef I iterator_type;
        typedef typename I::dual_iterator_type dual_iterator_type;
        typedef typename I::dual_reverse_iterator_type dual_reverse_iterator_type;

        void constraints () {
            iterator_type it = iterator_type ();
            // Indices
            it.index1 ();
            it.index2 ();
            // Iterator begin/end
            dual_iterator_type it_begin (it.begin ());
            dual_iterator_type it_end (it.end ());
            // Reverse iterator begin/end
            dual_reverse_iterator_type it_rbegin (it.rbegin ());
            dual_reverse_iterator_type it_rend (it.rend ());
            ignore_unused_variable_warning (it_begin);
            ignore_unused_variable_warning (it_end);
            ignore_unused_variable_warning (it_rbegin);
            ignore_unused_variable_warning (it_rend);
        }
    };

    template<class I1, class I2>
    struct IndexedBidirectional2DIteratorConcept {
        typedef I1 subiterator1_type;
        typedef I2 subiterator2_type;

        void constraints () {
            function_requires< BidirectionalIteratorConcept<subiterator1_type> >();
            function_requires< BidirectionalIteratorConcept<subiterator2_type> >();
            function_requires< Indexed2DIteratorConcept<subiterator1_type> >();
            function_requires< Indexed2DIteratorConcept<subiterator2_type> >();
        }
    };

    template<class I1, class I2>
    struct Mutable_IndexedBidirectional2DIteratorConcept {
        typedef I1 subiterator1_type;
        typedef I2 subiterator2_type;

        void constraints () {
            function_requires< Mutable_BidirectionalIteratorConcept<subiterator1_type> >();
            function_requires< Mutable_BidirectionalIteratorConcept<subiterator2_type> >();
            function_requires< Indexed2DIteratorConcept<subiterator1_type> >();
            function_requires< Indexed2DIteratorConcept<subiterator2_type> >();
        }
    };

    template<class I1, class I2>
    struct IndexedRandomAccess2DIteratorConcept {
        typedef I1 subiterator1_type;
        typedef I2 subiterator2_type;

        void constraints () {
            function_requires< RandomAccessIteratorConcept<subiterator1_type> >();
            function_requires< RandomAccessIteratorConcept<subiterator2_type> >();
            function_requires< Indexed2DIteratorConcept<subiterator1_type> >();
            function_requires< Indexed2DIteratorConcept<subiterator2_type> >();
        }
    };

    template<class I1, class I2>
    struct Mutable_IndexedRandomAccess2DIteratorConcept {
        typedef I1 subiterator1_type;
        typedef I2 subiterator2_type;

        void constraints () {
            function_requires< Mutable_RandomAccessIteratorConcept<subiterator1_type> >();
            function_requires< Mutable_RandomAccessIteratorConcept<subiterator2_type> >();
            function_requires< Indexed2DIteratorConcept<subiterator1_type> >();
            function_requires< Indexed2DIteratorConcept<subiterator2_type> >();
        }
    };

    template<class C>
    struct StorageArrayConcept {
        typedef C container_type;
        typedef typename C::size_type size_type;
        typedef typename C::value_type value_type;

        void constraints () {
            function_requires< RandomAccessContainerConcept<container_type> >();
            size_type n (0);
            // Sizing constructor
            container_type c = container_type (n);
            // Initialised sizing constructor
            container_type (n, value_type (5));
            ignore_unused_variable_warning (c);
        }
    };

    template<class C>
    struct Mutable_StorageArrayConcept {
        typedef C container_type;
        typedef typename C::size_type size_type;
        typedef typename C::value_type value_type;
        typedef typename C::iterator iterator_type;

        void constraints () {
            function_requires< Mutable_RandomAccessContainerConcept<container_type> > ();
            size_type n (0);
            // Sizing constructor
            container_type c = container_type (n);
            // Initialised sizing constructor
            c = container_type (n, value_type (3));
            // Resize
            c.resize (n, value_type (5));
            // Resize - none preserving
            c.resize (n);
        }
    };

    template<class C>
    struct StorageSparseConcept {
        typedef C container_type;
        typedef typename C::size_type size_type;

        void constraints () {
            function_requires< ReversibleContainerConcept<container_type> > ();
        }
    };

    template<class C>
    struct Mutable_StorageSparseConcept {
        typedef C container_type;
        typedef typename C::size_type size_type;
        typedef typename C::value_type value_type;
        typedef typename C::iterator iterator_type;

        void constraints () {
            // NOTE - Not Mutable_ReversibleContainerConcept
            function_requires< ReversibleContainerConcept<container_type> >();
            container_type c = container_type ();
            value_type t = value_type ();
            iterator_type it = iterator_type (), it1 = iterator_type (), it2 = iterator_type ();
            // Insert
            c.insert (it, t);
            // Erase
            c.erase (it);
            // Range erase
            c.erase (it1, it2);
            // Clear
            c.clear ();
        }
    };

    template<class G>
    struct IndexSetConcept {
        typedef G generator_type;
        typedef typename G::size_type size_type;
        typedef typename G::value_type value_type;

        void constraints () {
            function_requires< AssignableConcept<generator_type> >();
            function_requires< ReversibleContainerConcept<generator_type> >();
            generator_type g = generator_type ();
            size_type n (0);
            value_type t;
            // Element access
            t = g (n);
            ignore_unused_variable_warning (t);
        }
    };

    /** \brief Scalar expression concept.
     *
     * requirements
     * \li \c SE::value_type is the type of the scalar expression
     * \li \c SE must be convertable to \c SE::value_type
     * \li the constant \c SE::complexity must exist
     *
     * \param SE the type of the scalar expression
     */
    template<class SE>
    struct ScalarExpressionConcept {
        typedef SE scalar_expression_type;
        typedef typename SE::value_type value_type;

        static const unsigned complexity = SE::complexity;

        void constraints () {
            scalar_expression_type *sp;
            scalar_expression_type s = *sp;
            value_type t;
            // Conversion
            t = s;
            ignore_unused_variable_warning (t);
        }
    };

    /** \brief Vector expression concept.
     *
     * requirements
     * \li \c VE::value_type is the type of the elements
     * \li \c VE::const_reference The return type when accessing an element of a constant vector
     * expression. Must be convertable to a \c value_type.
     * \li \c VE::size_type is the (unsigned) type of the indices
     * \li \c VE::difference_type is the (signed) type of distances between indices
     * \li \c VE::category
     *
     * \li the constant \c SE::complexity must exist
     *
     * \param SE the type of the scalar expression
     */
    template<class VE>
    struct VectorExpressionConcept {
        typedef VE vector_expression_type;
        typedef typename VE::type_category type_category;
        typedef typename VE::size_type size_type;
        typedef typename VE::difference_type difference_type;
        typedef typename VE::value_type value_type;
        typedef typename VE::const_reference const_reference;
        typedef typename VE::const_iterator const_iterator_type;
        typedef typename VE::const_reverse_iterator const_reverse_iterator_type;

        void constraints () {
            vector_expression_type *vp;
            const vector_expression_type *cvp;
            vector_expression_type v = *vp;
            const vector_expression_type cv = *cvp;
            size_type n (0), i (0);
            value_type t;
            // Find (internal?)
            const_iterator_type cit (v.find (i));
            // Beginning of range
            const_iterator_type cit_begin (v.begin ());
            // End of range
            const_iterator_type cit_end (v.end ());
            // Size
            n = v.size ();
            // Beginning of reverse range
            const_reverse_iterator_type crit_begin (cv.rbegin ());
            // End of reverse range
            const_reverse_iterator_type crit_end (cv.rend ());
            // Element access
            t = v (i);
            ignore_unused_variable_warning (n);
            ignore_unused_variable_warning (cit);
            ignore_unused_variable_warning (cit_begin);
            ignore_unused_variable_warning (cit_end);
            ignore_unused_variable_warning (crit_begin);
            ignore_unused_variable_warning (crit_end);
            ignore_unused_variable_warning (t);
        }
    };

    template<class VE>
    struct Mutable_VectorExpressionConcept {
        typedef VE vector_expression_type;
        typedef typename VE::size_type size_type;
        typedef typename VE::value_type value_type;
        typedef typename VE::iterator iterator_type;
        typedef typename VE::reverse_iterator reverse_iterator_type;

        void constraints () {
            function_requires< AssignableConcept<vector_expression_type> >();
            function_requires< VectorExpressionConcept<vector_expression_type> >();
            vector_expression_type *vp;
            vector_expression_type v = *vp, v1 = *vp, v2 = *vp;
            size_type i (0);
            value_type t = value_type ();
            // Find (internal?)
            iterator_type it (v.find (i));
            // Beginning of range
            iterator_type it_begin (v.begin ());
            // End of range
            iterator_type it_end (v.end ());
            // Swap
            v1.swap (v2);
            // Beginning of reverse range
            reverse_iterator_type rit_begin (v.rbegin ());
            // End of reverse range
            reverse_iterator_type rit_end (v.rend ());
            // Assignments
            v2 = v1;
            v2.assign (v1);
            v2 += v1;
            v2.plus_assign (v1);
            v2 -= v1;
            v2.minus_assign (v1);
            v *= t;
            ignore_unused_variable_warning (it);
            ignore_unused_variable_warning (it_begin);
            ignore_unused_variable_warning (it_end);
            ignore_unused_variable_warning (rit_begin);
            ignore_unused_variable_warning (rit_end);
        }
    };

    template<class ME>
    struct MatrixExpressionConcept {
        typedef ME matrix_expression_type;
        typedef typename ME::type_category type_category;
        typedef typename ME::size_type size_type;
        typedef typename ME::value_type value_type;
        typedef typename ME::const_iterator1 const_subiterator1_type;
        typedef typename ME::const_iterator2 const_subiterator2_type;
        typedef typename ME::const_reverse_iterator1 const_reverse_subiterator1_type;
        typedef typename ME::const_reverse_iterator2 const_reverse_subiterator2_type;

        void constraints () {
            matrix_expression_type *mp;
            const matrix_expression_type *cmp;
            matrix_expression_type m = *mp;
            const matrix_expression_type cm = *cmp;
            size_type n (0), i (0), j (0);
            value_type t;
            // Find (internal?)
            const_subiterator1_type cit1 (m.find1 (0, i, j));
            const_subiterator2_type cit2 (m.find2 (0, i, j));
            // Beginning of range
            const_subiterator1_type cit1_begin (m.begin1 ());
            const_subiterator2_type cit2_begin (m.begin2 ());
            // End of range
            const_subiterator1_type cit1_end (m.end1 ());
            const_subiterator2_type cit2_end (m.end2 ());
            // Size
            n = m.size1 ();
            n = m.size2 ();
            // Beginning of reverse range
            const_reverse_subiterator1_type crit1_begin (cm.rbegin1 ());
            const_reverse_subiterator2_type crit2_begin (cm.rbegin2 ());
            // End of reverse range
            const_reverse_subiterator1_type crit1_end (cm.rend1 ());
            const_reverse_subiterator2_type crit2_end (cm.rend2 ());
            // Element access
            t = m (i, j);
            ignore_unused_variable_warning (n);
            ignore_unused_variable_warning (cit1);
            ignore_unused_variable_warning (cit2);
            ignore_unused_variable_warning (cit1_begin);
            ignore_unused_variable_warning (cit2_begin);
            ignore_unused_variable_warning (cit1_end);
            ignore_unused_variable_warning (cit2_end);
            ignore_unused_variable_warning (crit1_begin);
            ignore_unused_variable_warning (crit2_begin);
            ignore_unused_variable_warning (crit1_end);
            ignore_unused_variable_warning (crit2_end);
            ignore_unused_variable_warning (t);
        }
    };

    template<class ME>
    struct Mutable_MatrixExpressionConcept {
        typedef ME matrix_expression_type;
        typedef typename ME::size_type size_type;
        typedef typename ME::value_type value_type;
        typedef typename ME::iterator1 subiterator1_type;
        typedef typename ME::iterator2 subiterator2_type;
        typedef typename ME::reverse_iterator1 reverse_subiterator1_type;
        typedef typename ME::reverse_iterator2 reverse_subiterator2_type;

        void constraints () {
            function_requires< AssignableConcept<matrix_expression_type> >();
            function_requires< MatrixExpressionConcept<matrix_expression_type> >();
            matrix_expression_type *mp;
            matrix_expression_type m = *mp, m1 = *mp, m2 = *mp;
            size_type i (0), j (0);
            value_type t = value_type ();
            // Find (internal?)
            subiterator1_type it1 (m.find1 (0, i, j));
            subiterator2_type it2 (m.find2 (0, i, j));
            // Beginning of range
            subiterator1_type it1_begin (m.begin1 ());
            subiterator2_type it2_begin (m.begin2 ());
            // End of range
            subiterator1_type it1_end (m.end1 ());
            subiterator2_type it2_end (m.end2 ());
            // Swap
            m1.swap (m2);
            // Beginning of reverse range
            reverse_subiterator1_type rit1_begin (m.rbegin1 ());
            reverse_subiterator2_type rit2_begin (m.rbegin2 ());
            // End of reverse range
            reverse_subiterator1_type rit1_end (m.rend1 ());
            reverse_subiterator2_type rit2_end (m.rend2 ());
            // Assignments
            m2 = m1;
            m2.assign (m1);
            m2 += m1;
            m2.plus_assign (m1);
            m2 -= m1;
            m2.minus_assign (m1);
            m *= t;
            ignore_unused_variable_warning (it1);
            ignore_unused_variable_warning (it2);
            ignore_unused_variable_warning (it1_begin);
            ignore_unused_variable_warning (it2_begin);
            ignore_unused_variable_warning (it1_end);
            ignore_unused_variable_warning (it2_end);
            ignore_unused_variable_warning (rit1_begin);
            ignore_unused_variable_warning (rit2_begin);
            ignore_unused_variable_warning (rit1_end);
            ignore_unused_variable_warning (rit2_end);
        }
    };

    template<class V>
    struct VectorConcept {
        typedef V vector_type;
        typedef typename V::size_type size_type;
        typedef typename V::value_type value_type;
        typedef const value_type *const_pointer;

        void constraints () {
            function_requires< VectorExpressionConcept<vector_type> >();
            size_type n (0);
            size_type i (0);
            // Sizing constructor
            vector_type v (n);
            // Element support
            const_pointer p = v.find_element (i);

            ignore_unused_variable_warning (p);
        }
    };

    template<class V>
    struct Mutable_VectorConcept {
        typedef V vector_type;
        typedef typename V::size_type size_type;
        typedef typename V::value_type value_type;
        typedef value_type *pointer;

        void constraints () {
            function_requires< VectorConcept<vector_type> >();
            function_requires< DefaultConstructible<vector_type> >();
            function_requires< Mutable_VectorExpressionConcept<vector_type> >();
            size_type n (0);
            value_type t = value_type ();
            size_type i (0);
            vector_type v;
            // Element support
            pointer p = v.find_element (i);
            // Element assignment
            value_type r = v.insert_element (i, t);
            v.insert_element (i, t) = r;
            // Zeroing
            v.clear ();
            // Resize
            v.resize (n);

            ignore_unused_variable_warning (p);
            ignore_unused_variable_warning (r);
        }
    };

    template<class V>
    struct SparseVectorConcept {
        typedef V vector_type;
        typedef typename V::size_type size_type;

        void constraints () {
            function_requires< VectorConcept<vector_type> >();
        }
    };

    template<class V>
    struct Mutable_SparseVectorConcept {
        typedef V vector_type;
        typedef typename V::size_type size_type;
        typedef typename V::value_type value_type;

        void constraints () {
            function_requires< SparseVectorConcept<vector_type> >();
            function_requires< Mutable_VectorConcept<vector_type> >();
            size_type i (0);
            vector_type v;
            // Element erasure
            v.erase_element (i);
        }
    };

    template<class M>
    struct MatrixConcept {
        typedef M matrix_type;
        typedef typename M::size_type size_type;
        typedef typename M::value_type value_type;
        typedef const value_type *const_pointer;

        void constraints () {
            function_requires< MatrixExpressionConcept<matrix_type> >();
            size_type n (0);
            size_type i (0), j (0);
            // Sizing constructor
            matrix_type m (n, n);
            // Element support
#ifndef SKIP_BAD
            const_pointer p = m.find_element (i, j);
#else
            const_pointer p;
            ignore_unused_variable_warning (i);
            ignore_unused_variable_warning (j);
#endif
            ignore_unused_variable_warning (p);
        }
    };

    template<class M>
    struct Mutable_MatrixConcept {
        typedef M matrix_type;
        typedef typename M::size_type size_type;
        typedef typename M::value_type value_type;
        typedef value_type *pointer;

        void constraints () {
            function_requires< MatrixConcept<matrix_type> >();
            function_requires< DefaultConstructible<matrix_type> >();
            function_requires< Mutable_MatrixExpressionConcept<matrix_type> >();
            size_type n (0);
            value_type t = value_type ();
            size_type i (0), j (0);
            matrix_type m;
            // Element support
#ifndef SKIP_BAD
            pointer p = m.find_element (i, j);
            ignore_unused_variable_warning (i);
            ignore_unused_variable_warning (j);
#else
            pointer p;
#endif
            // Element assigment
            value_type r = m.insert_element (i, j, t);
            m.insert_element (i, j, t) = r;
            // Zeroing
            m.clear ();
            // Resize
            m.resize (n, n);
            m.resize (n, n, false);

            ignore_unused_variable_warning (p);
            ignore_unused_variable_warning (r);
        }
    };

    template<class M>
    struct SparseMatrixConcept {
        typedef M matrix_type;
        typedef typename M::size_type size_type;

        void constraints () {
            function_requires< MatrixConcept<matrix_type> >();
        }
    };

    template<class M>
    struct Mutable_SparseMatrixConcept {
        typedef M matrix_type;
        typedef typename M::size_type size_type;
        typedef typename M::value_type value_type;

        void constraints () {
            function_requires< SparseMatrixConcept<matrix_type> >();
            function_requires< Mutable_MatrixConcept<matrix_type> >();
            size_type i (0), j (0);
            matrix_type m;
            // Elemnent erasure
            m.erase_element (i, j);
        }
    };

    /** introduce anonymous namespace to make following functions
     * local to the current compilation unit.
     */
    namespace {

    // Replaced the ZeroElement and OneElement functions with the templated versions
    // because the former where giving warnings with clang
    template<class T>
    T
    ZeroElement (T) {
        return T(0.0);
    }

    template<class T>
    vector<T>
    ZeroElement (vector<T>) {
        return zero_vector<T> ();
    }

    template<class T>
    matrix<T>
    ZeroElement (matrix<T>) {
        return zero_matrix<T> ();
    }

    template<class T>
    T
    OneElement (T) {
        return T(0.0);
    }

    template<class T>
    vector<T>
    OneElement (vector<T>) {
        return zero_vector<T> ();
    }

    template<class T>
    matrix<T>
    OneElement (matrix<T>) {
        return identity_matrix<T> ();
    }

//    template<>
//    float
//    ZeroElement (float) {
//        return 0.f;
//    }
//    template<>
//    double
//    ZeroElement (double) {
//        return 0.;
//    }
//    template<>
//    vector<float>
//    ZeroElement (vector<float>) {
//        return zero_vector<float> ();
//    }
//    template<>
//    vector<double>
//    ZeroElement (vector<double>) {
//        return zero_vector<double> ();
//    }
//    template<>
//    matrix<float>
//    ZeroElement (matrix<float>) {
//        return zero_matrix<float> ();
//    }
//    template<>
//    matrix<double>
//    ZeroElement (matrix<double>) {
//        return zero_matrix<double> ();
//    }
//    template<>
//    std::complex<float>
//    ZeroElement (std::complex<float>) {
//        return std::complex<float> (0.f);
//    }
//    template<>
//    std::complex<double>
//    ZeroElement (std::complex<double>) {
//        return std::complex<double> (0.);
//    }
//    template<>
//    vector<std::complex<float> >
//    ZeroElement (vector<std::complex<float> >) {
//        return zero_vector<std::complex<float> > ();
//    }
//    template<>
//    vector<std::complex<double> >
//    ZeroElement (vector<std::complex<double> >) {
//        return zero_vector<std::complex<double> > ();
//    }
//    template<>
//    matrix<std::complex<float> >
//    ZeroElement (matrix<std::complex<float> >) {
//        return zero_matrix<std::complex<float> > ();
//    }
//    template<>
//    matrix<std::complex<double> >
//    ZeroElement (matrix<std::complex<double> >) {
//        return zero_matrix<std::complex<double> > ();
//    }

//    template<class T>
//    T
//    OneElement (T);
//    template<>
//    float
//    OneElement (float) {
//        return 1.f;
//    }
//    template<>
//    double
//    OneElement (double) {
//        return 1.;
//    }
//    template<>
//    matrix<float>
//    OneElement (matrix<float>) {
//        return identity_matrix<float> ();
//    }
//    template<>
//    matrix<double>
//    OneElement (matrix<double>) {
//        return identity_matrix<double> ();
//    }
//    template<>
//    std::complex<float>
//    OneElement (std::complex<float>) {
//        return std::complex<float> (1.f);
//    }
//    template<>
//    std::complex<double>
//    OneElement (std::complex<double>) {
//        return std::complex<double> (1.);
//    }
//    template<>
//    matrix<std::complex<float> >
//    OneElement (matrix<std::complex<float> >) {
//        return identity_matrix<std::complex<float> > ();
//    }
//    template<>
//    matrix<std::complex<double> >
//    OneElement (matrix<std::complex<double> >) {
//        return identity_matrix<std::complex<double> > ();
//    }

    template<class E1, class E2>
    bool
    operator == (const vector_expression<E1> &e1, const vector_expression<E2> &e2) {
        typedef typename promote_traits<typename E1::value_type,
                                                    typename E2::value_type>::promote_type value_type;
        typedef typename type_traits<value_type>::real_type real_type;
        return norm_inf (e1 - e2) == real_type/*zero*/();
    }
    template<class E1, class E2>
    bool
    operator == (const matrix_expression<E1> &e1, const matrix_expression<E2> &e2) {
        typedef typename promote_traits<typename E1::value_type,
                                                    typename E2::value_type>::promote_type value_type;
        typedef typename type_traits<value_type>::real_type real_type;
        return norm_inf (e1 - e2) == real_type/*zero*/();
    }

    template<class T>
    struct AdditiveAbelianGroupConcept {
        typedef T value_type;

        void constraints () {
            bool r;
            value_type a = value_type (), b = value_type (), c = value_type ();
            r = (a + b) + c == a + (b + c);
            r = ZeroElement (value_type ()) + a == a;
            r = a + ZeroElement (value_type ()) == a;
            r = a + (- a) == ZeroElement (value_type ());
            r = (- a) + a == ZeroElement (value_type ());
            r = a + b == b + a;
            ignore_unused_variable_warning (r);
        }
    };

    template<class T>
    struct MultiplicativeAbelianGroupConcept {
        typedef T value_type;

        void constraints () {
            bool r;
            value_type a = value_type (), b = value_type (), c = value_type ();
            r = (a * b) * c == a * (b * c);
            r = OneElement (value_type ()) * a == a;
            r = a * OneElement (value_type ()) == a;
            r = a * (OneElement (value_type ()) / a) == a;
            r = (OneElement (value_type ()) / a) * a == a;
            r = a * b == b * a;
            ignore_unused_variable_warning (r);
        }
    };

    template<class T>
    struct RingWithIdentityConcept {
        typedef T value_type;

        void constraints () {
            function_requires< AdditiveAbelianGroupConcept<value_type> >();
            bool r;
            value_type a = value_type (), b = value_type (), c = value_type ();
            r = (a * b) * c == a * (b * c);
            r = (a + b) * c == a * c + b * c;
            r = OneElement (value_type ()) * a == a;
            r = a * OneElement (value_type ()) == a;
            ignore_unused_variable_warning (r);
        }
    };

    template<class T>
    struct Prod_RingWithIdentityConcept {
        typedef T value_type;

        void constraints () {
            function_requires< AdditiveAbelianGroupConcept<value_type> >();
            bool r;
            value_type a = value_type (), b = value_type (), c = value_type ();
            r = prod (T (prod (a, b)), c) == prod (a, T (prod (b, c)));
            r = prod (a + b, c) == prod (a, c) + prod (b, c);
            r = prod (OneElement (value_type ()), a) == a;
            r = prod (a, OneElement (value_type ())) == a;
            ignore_unused_variable_warning (r);
        }
    };

    template<class T>
    struct CommutativeRingWithIdentityConcept {
        typedef T value_type;

        void constraints () {
            function_requires< RingWithIdentityConcept<value_type> >();
            bool r;
            value_type a = value_type (), b = value_type ();
            r = a * b == b * a;
            ignore_unused_variable_warning (r);
        }
    };

    template<class T>
    struct FieldConcept {
        typedef T value_type;

        void constraints () {
            function_requires< CommutativeRingWithIdentityConcept<value_type> >();
            bool r;
            value_type a = value_type ();
            r = a == ZeroElement (value_type ()) || a * (OneElement (value_type ()) / a) == a;
            r = a == ZeroElement (value_type ()) || (OneElement (value_type ()) / a) * a == a;
            ignore_unused_variable_warning (r);
        }
    };

    template<class T, class V>
    struct VectorSpaceConcept {
        typedef T value_type;
        typedef V vector_type;

        void constraints () {
            function_requires< FieldConcept<value_type> >();
            function_requires< AdditiveAbelianGroupConcept<vector_type> >();
            bool r;
            value_type alpha = value_type (), beta = value_type ();
            vector_type a = vector_type (), b = vector_type ();
            r = alpha * (a + b) == alpha * a + alpha * b;
            r = (alpha + beta) * a == alpha * a + beta * a;
            r = (alpha * beta) * a == alpha * (beta * a);
            r = OneElement (value_type ()) * a == a;
            ignore_unused_variable_warning (r);
        }
    };

    template<class T, class V, class M>
    struct LinearOperatorConcept {
        typedef T value_type;
        typedef V vector_type;
        typedef M matrix_type;

        void constraints () {
            function_requires< VectorSpaceConcept<value_type, vector_type> >();
            bool r;
            value_type alpha = value_type (), beta = value_type ();
            vector_type a = vector_type (), b = vector_type ();
            matrix_type A = matrix_type ();
            r = prod (A, alpha * a + beta * b) == alpha * prod (A, a) + beta * prod (A, b);
            ignore_unused_variable_warning (r);
        }
    };

inline void concept_checks () {

        // Allow tests to be group to keep down compiler storage requirement
#ifdef INTERAL
#define INTERNAL_STORAGE
#define INTERNAL_VECTOR
#define INTERNAL_MATRIX
#define INTERNAL_SPECIAL
#define INTERNAL_SPARSE
#define INTERNAL_EXPRESSION
#endif

        // TODO enable this for development
        // #define VIEW_CONCEPTS

        // Element value type for tests
        typedef float T;

        // Storage Array
#if defined (INTERNAL_STORAGE) || defined (INTERNAL_STORAGE_DENSE)
        {
            typedef std::vector<T> container_model;
            function_requires< Mutable_StorageArrayConcept<container_model> >();
            function_requires< RandomAccessIteratorConcept<container_model::const_iterator> >();
            function_requires< Mutable_RandomAccessIteratorConcept<container_model::iterator> >();
        }

        {
            typedef bounded_array<T, 1> container_model;
            function_requires< Mutable_StorageArrayConcept<container_model> >();
            function_requires< RandomAccessIteratorConcept<container_model::const_iterator> >();
            function_requires< Mutable_RandomAccessIteratorConcept<container_model::iterator> >();
        }

        {
            typedef unbounded_array<T> container_model;
            function_requires< Mutable_StorageArrayConcept<container_model> >();
            function_requires< RandomAccessIteratorConcept<container_model::const_iterator> >();
            function_requires< Mutable_RandomAccessIteratorConcept<container_model::iterator> >();
        }

/* FIXME array_adaptors are in progress
        {
            typedef array_adaptor<T> container_model;
            function_requires< Mutable_StorageArrayConcept<container_model> >();
            function_requires< RandomAccessIteratorConcept<container_model::const_iterator> >();
            function_requires< Mutable_RandomAccessIteratorConcept<container_model::iterator> >();
        }
*/

        {
            typedef range container_model;
            function_requires< IndexSetConcept<range> >();
            function_requires< RandomAccessIteratorConcept<range::const_iterator> >();
        }

        {
            typedef slice container_model;
            function_requires< IndexSetConcept<range> >();
            function_requires< RandomAccessIteratorConcept<range::const_iterator> >();
        }

        {
            typedef indirect_array<> container_model;
            function_requires< IndexSetConcept<range> >();
            function_requires< RandomAccessIteratorConcept<range::const_iterator> >();
        }
#endif

        // Storage Sparse
#if defined (INTERNAL_STORAGE) || defined (INTERNAL_STORAGE_SPARSE)
        {
           typedef map_array<std::size_t, T> container_model;
           function_requires< Mutable_StorageSparseConcept<container_model> >();
           function_requires< RandomAccessIteratorConcept<container_model::const_iterator> >();
           function_requires< RandomAccessIteratorConcept<container_model::iterator> >();
        }

        {
           typedef std::map<std::size_t, T> container_model;
           function_requires< Mutable_StorageSparseConcept<container_model > >();
           function_requires< BidirectionalIteratorConcept<container_model::const_iterator> >();
           function_requires< BidirectionalIteratorConcept<container_model::iterator> >();
        }
#endif

#ifdef VIEW_CONCEPTS
        // read only vectors
        {
           typedef vector_view<T> container_model;
           function_requires< RandomAccessContainerConcept<container_model> >();
           function_requires< VectorConcept<container_model> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
        }
#endif

        // Vector
#if defined (INTERNAL_VECTOR) || defined (INTERNAL_VECTOR_DENSE)
        {
           typedef vector<T> container_model;
           function_requires< RandomAccessContainerConcept<container_model> >();
           function_requires< Mutable_VectorConcept<container_model> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
           function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
           function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
        }

        {
           typedef zero_vector<T> container_model;
           function_requires< VectorConcept<container_model> >();
           function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_iterator> >();
           function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_reverse_iterator> >();
        }

        {
           typedef unit_vector<T> container_model;
           function_requires< VectorConcept<container_model> >();
           function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_iterator> >();
           function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_reverse_iterator> >();
        }

        {
           typedef scalar_vector<T> container_model;
           function_requires< VectorConcept<container_model> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
        }

        {
           typedef c_vector<T, 1> container_model;
           function_requires< Mutable_VectorConcept<container_model> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
           function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
           function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
        }
#endif

        // Vector Proxies
#if defined (INTERNAL_VECTOR) || defined (INTERNAL_VECTOR_PROXY)
        {
           typedef vector_range<vector<T> > container_model;
           function_requires< Mutable_VectorExpressionConcept<container_model> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
           function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
           function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
        }

        {
           typedef vector_slice<vector<T> > container_model;
           function_requires< Mutable_VectorExpressionConcept<container_model> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
           function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
           function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
        }

        {
           typedef vector_indirect<vector<T> > container_model;
           function_requires< Mutable_VectorExpressionConcept<container_model> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
           function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
           function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
           function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
        }
#endif

        // Sparse Vector
#if defined (INTERNAL_SPARSE) || defined (INTERNAL_VECTOR_SPARSE)
        {
            typedef mapped_vector<T> container_model;
            function_requires< Mutable_SparseVectorConcept<container_model> >();
            function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_iterator> >();
            function_requires< Mutable_IndexedBidirectional1DIteratorConcept<container_model::iterator> >();
            function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_reverse_iterator> >();
            function_requires< Mutable_IndexedBidirectional1DIteratorConcept<container_model::reverse_iterator> >();
        }

        {
            typedef compressed_vector<T> container_model;
            function_requires< Mutable_SparseVectorConcept<container_model> >();
            function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_iterator> >();
            function_requires< Mutable_IndexedBidirectional1DIteratorConcept<container_model::iterator> >();
            function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_reverse_iterator> >();
            function_requires< Mutable_IndexedBidirectional1DIteratorConcept<container_model::reverse_iterator> >();
        }

        {
            typedef coordinate_vector<T> container_model;
            function_requires< Mutable_SparseVectorConcept<container_model> >();
            function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_iterator> >();
            function_requires< Mutable_IndexedBidirectional1DIteratorConcept<container_model::iterator> >();
            function_requires< IndexedBidirectional1DIteratorConcept<container_model::const_reverse_iterator> >();
            function_requires< Mutable_IndexedBidirectional1DIteratorConcept<container_model::reverse_iterator> >();
        }
#endif

        // Matrix
#if defined (INTERNAL_MATRIX) || defined (INTERNAL_MATRIX_DENSE)
        {
            typedef matrix<T> container_model;
            function_requires< Mutable_MatrixConcept<matrix<T> > >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }

        {
            typedef vector_of_vector<T> container_model;
            function_requires< Mutable_MatrixConcept<matrix<T> > >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }

        {
            typedef zero_matrix<T> container_model;
            function_requires< Mutable_MatrixConcept<matrix<T> > >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
        }

        {
            typedef identity_matrix<T> container_model;
            function_requires< Mutable_MatrixConcept<matrix<T> > >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
        }

        {
            typedef scalar_matrix<T> container_model;
            function_requires< Mutable_MatrixConcept<matrix<T> > >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
        }

        {
            typedef c_matrix<T, 1, 1> container_model;
            function_requires< Mutable_MatrixConcept<matrix<T> > >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }
#endif

        // Matrix Proxies
#if defined (INTERNAL_MATRIX) || defined (INTERNAL_MATRIX_PROXY)
        {
            typedef matrix_row<matrix<T> > container_model;
            function_requires< Mutable_VectorExpressionConcept<container_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
            function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
            function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
        }

        {
            typedef matrix_column<matrix<T> > container_model;
            function_requires< Mutable_VectorExpressionConcept<container_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
            function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
            function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
        }

        {
            typedef matrix_vector_range<matrix<T> > container_model;
            function_requires< Mutable_VectorExpressionConcept<container_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
            function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
            function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
        }

        {
            typedef matrix_vector_slice<matrix<T> > container_model;
            function_requires< Mutable_VectorExpressionConcept<container_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
            function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
            function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
        }

        {
            typedef matrix_vector_indirect<matrix<T> > container_model;
            function_requires< Mutable_VectorExpressionConcept<container_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_iterator> >();
            function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<container_model::const_reverse_iterator> >();
            function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<container_model::reverse_iterator> >();
        }

        {
            typedef matrix_range<matrix<T> > container_model;
            function_requires< Mutable_MatrixExpressionConcept<container_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }

        {
            typedef matrix_slice<matrix<T> > container_model;
            function_requires< Mutable_MatrixExpressionConcept<container_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }

        {
            typedef matrix_indirect<matrix<T> > container_model;
            function_requires< Mutable_MatrixExpressionConcept<container_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }
#endif

        // Banded Matrix
#if defined (INTERNAL_SPECIAL) || defined (INTERNAL_BANDED)
        {
            typedef banded_matrix<T> container_model;
            function_requires< Mutable_MatrixConcept<container_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }

        {
            typedef banded_adaptor<matrix<T> > container_model;
            function_requires< Mutable_MatrixExpressionConcept<container_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }
#endif

        // Triangular Matrix
#if defined (INTERNAL_SPECIAL) || defined (INTERNAL_TRIANGULAR)
        {
            typedef triangular_matrix<T> container_model;
            function_requires< Mutable_MatrixConcept<container_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }

        {
            typedef triangular_adaptor<matrix<T> > container_model;
            function_requires< Mutable_MatrixExpressionConcept<container_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }
#endif

        // Symmetric Matrix
#if defined (INTERNA_SPECIAL) || defined (INTERNAL_SYMMETRIC)
        {
            typedef symmetric_matrix<T> container_model;
            function_requires< Mutable_MatrixConcept<container_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }

        {
            typedef banded_adaptor<matrix<T> > container_model;
#ifndef SKIP_BAD
           // const_iterator (iterator) constructor is bad
            function_requires< Mutable_MatrixExpressionConcept<container_model> >();
#endif
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }
#endif

        // Hermitian Matrix
#if defined (INTERNAL_SPECIAL) || defined (INTERNAL_HERMITIAN)
        {
            typedef hermitian_matrix<T> container_model;
            function_requires< Mutable_MatrixConcept<container_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }

        {
            typedef hermitian_adaptor<matrix<T> > container_model;
#ifndef SKIP_BAD
           // const_iterator (iterator) constructor is bad
            function_requires< Mutable_MatrixExpressionConcept<container_model> >();
#endif
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }
#endif

        // Sparse Matrix
#if defined (INTERNAL_SPARSE) || defined (INTERNAL_MATRIX_SPARSE)
        {
            typedef mapped_matrix<T> container_model;
            function_requires< Mutable_SparseMatrixConcept<container_model> >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }
        {
            typedef mapped_vector_of_mapped_vector<T> container_model;
            function_requires< Mutable_SparseMatrixConcept<container_model> >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }
        {
            typedef compressed_matrix<T> container_model;
            function_requires< Mutable_SparseMatrixConcept<container_model> >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }
        {
            typedef coordinate_matrix<T> container_model;
            function_requires< Mutable_SparseMatrixConcept<container_model> >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }
        {
            typedef generalized_vector_of_vector<T, row_major, vector< coordinate_vector<T> > > container_model;
            function_requires< Mutable_SparseMatrixConcept<container_model> >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_iterator1, container_model::const_iterator2> >();
            function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::iterator1, container_model::iterator2> >();
            function_requires< IndexedBidirectional2DIteratorConcept<container_model::const_reverse_iterator1, container_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedBidirectional2DIteratorConcept<container_model::reverse_iterator1, container_model::reverse_iterator2> >();
        }

#endif

        // Scalar Expressions
#if defined (INTERNAL_EXPRESSION) || defined (INTERNAL_VECTOR_EXPRESSION)
        function_requires< ScalarExpressionConcept<scalar_value<T> > >();
        function_requires< ScalarExpressionConcept<scalar_reference<T> > >();

        // Vector Expressions
        {
            typedef vector_reference<vector<T> > expression_model;
            function_requires< VectorExpressionConcept<expression_model> >();
            function_requires< Mutable_VectorExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
            function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<expression_model::iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
            function_requires< Mutable_IndexedRandomAccess1DIteratorConcept<expression_model::reverse_iterator> >();
        }

        {
            typedef vector_unary<vector<T>, scalar_identity<T> > expression_model;
            function_requires< VectorExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
        }

        {
            typedef vector_binary<vector<T>, vector<T>, scalar_plus<T, T> > expression_model;
            function_requires< VectorExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
        }

        {
            typedef vector_binary_scalar1<T, vector<T>, scalar_multiplies<T, T> > expression_model;
            function_requires< VectorExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
        }

        {
            typedef vector_binary_scalar2<vector<T>, scalar_value<T>, scalar_multiplies<T, T> > expression_model;
            function_requires< VectorExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
        }

        {
            typedef vector_binary_scalar1<scalar_value<T>, vector<T>, scalar_multiplies<T, T> > expression_model;
            function_requires< VectorExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
        }

        {
            typedef vector_binary_scalar2<vector<T>, scalar_value<T>, scalar_multiplies<T, T> > expression_model;
            function_requires< VectorExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
        }

        function_requires< ScalarExpressionConcept<vector_scalar_unary<vector<T>, vector_sum<vector<T> > > > >();
        function_requires< ScalarExpressionConcept<vector_scalar_unary<vector<T>, vector_norm_1<vector<T> > > > >();
        function_requires< ScalarExpressionConcept<vector_scalar_unary<vector<T>, vector_norm_2<vector<T> > > > >();
        function_requires< ScalarExpressionConcept<vector_scalar_unary<vector<T>, vector_norm_inf<vector<T> > > > >();

        function_requires< ScalarExpressionConcept<vector_scalar_binary<vector<T>, vector<T>, vector_inner_prod<vector<T>, vector<T>, T> > > >();
#endif

        // Matrix Expressions
#if defined (INTERNAL_EXPRESSION) || defined (INTERNAL_MATRIX_EXPRESSION)
        {
            typedef matrix_reference<matrix<T> > expression_model;
            function_requires< MatrixExpressionConcept<expression_model> >();
            function_requires< Mutable_MatrixExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<expression_model::iterator1, expression_model::iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
            function_requires< Mutable_IndexedRandomAccess2DIteratorConcept<expression_model::reverse_iterator1, expression_model::reverse_iterator2> >();
        }

        {
            typedef vector_matrix_binary<vector<T>, vector<T>, scalar_multiplies<T, T> > expression_model;
            function_requires< MatrixExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
        }

        {
            typedef matrix_unary1<matrix<T>, scalar_identity<T> > expression_model;
            function_requires< MatrixExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
        }

        {
            typedef matrix_unary2<matrix<T>, scalar_identity<T> > expression_model;
            function_requires< MatrixExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
        }

        {
            typedef matrix_binary<matrix<T>, matrix<T>, scalar_plus<T, T> > expression_model;
            function_requires< MatrixExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
        }

        {
            typedef matrix_binary_scalar1<T, matrix<T>, scalar_multiplies<T, T> > expression_model;
            function_requires< MatrixExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
        }

        {
            typedef matrix_binary_scalar2<matrix<T>, T, scalar_multiplies<T, T> > expression_model;
            function_requires< MatrixExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
        }

        {
            typedef matrix_binary_scalar1<scalar_value<T>, matrix<T>, scalar_multiplies<T, T> > expression_model;
            function_requires< MatrixExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
        }

        {
            typedef matrix_binary_scalar2<matrix<T>, scalar_value<T>, scalar_multiplies<T, T> > expression_model;
            function_requires< MatrixExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
        }

        {
            typedef matrix_vector_binary1<matrix<T>, vector<T>, matrix_vector_prod1<matrix<T>, vector<T>, T> > expression_model;
            function_requires< VectorExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
        }

        {
            typedef matrix_vector_binary2<vector<T>, matrix<T>, matrix_vector_prod2<matrix<T>, vector<T>, T > > expression_model;
            function_requires< VectorExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_iterator> >();
            function_requires< IndexedRandomAccess1DIteratorConcept<expression_model::const_reverse_iterator> >();
        }

        {
            typedef matrix_matrix_binary<matrix<T>, matrix<T>, matrix_matrix_prod<matrix<T>, matrix<T>, T > > expression_model;
            function_requires< MatrixExpressionConcept<expression_model> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_iterator1, expression_model::const_iterator2> >();
            function_requires< IndexedRandomAccess2DIteratorConcept<expression_model::const_reverse_iterator1, expression_model::const_reverse_iterator2> >();
        }

        function_requires< ScalarExpressionConcept<matrix_scalar_unary<matrix<T>, matrix_norm_1<vector<T> > > > >();
        function_requires< ScalarExpressionConcept<matrix_scalar_unary<matrix<T>, matrix_norm_frobenius<vector<T> > > > >();
        function_requires< ScalarExpressionConcept<matrix_scalar_unary<matrix<T>, matrix_norm_inf<vector<T> > > > >();
#endif

#ifdef EXTERNAL
        function_requires< AdditiveAbelianGroupConcept<T> >();
        function_requires< CommutativeRingWithIdentityConcept<T> >();
        function_requires< FieldConcept<T> >();
        function_requires< VectorSpaceConcept<T, vector<T> > >();
        function_requires< Prod_RingWithIdentityConcept<matrix<T> > >();
        function_requires< VectorSpaceConcept<T, matrix<T> > >();
        function_requires< LinearOperatorConcept<T, vector<T>, matrix<T> > >();

        function_requires< AdditiveAbelianGroupConcept<std::complex<T> > >();
        function_requires< CommutativeRingWithIdentityConcept<std::complex<T> > >();
        function_requires< FieldConcept<std::complex<T> > >();
        function_requires< VectorSpaceConcept<std::complex<T>, vector<std::complex<T> > > >();
        function_requires< Prod_RingWithIdentityConcept<matrix<std::complex<T> > > >();
        function_requires< VectorSpaceConcept<std::complex<T>, matrix<std::complex<T> > > >();
        function_requires< LinearOperatorConcept<std::complex<T>, vector<std::complex<T> >, matrix<std::complex<T> > > >();
#endif
    }

    } // end of anonymous namespace

}}}

#endif


code:
//
//  Copyright (c) 2000-2002
//  Joerg Walter, Mathias Koch
//
//  Distributed under the Boost Software License, Version 1.0. (See
//  accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt)
//
//  The authors gratefully acknowledge the support of
//  GeNeSys mbH & Co. KG in producing this work.
//

#ifndef _BOOST_UBLAS_CONFIG_
#define _BOOST_UBLAS_CONFIG_

#include <cassert>
#include <cstddef>
#include <algorithm>
#include <limits>

#include <boost/config.hpp>
#include <boost/static_assert.hpp>
#include <boost/noncopyable.hpp>
#include <boost/mpl/if.hpp>
#include <boost/mpl/and.hpp>
#include <boost/type_traits/is_same.hpp>
#include <boost/type_traits/is_convertible.hpp>
#include <boost/type_traits/is_const.hpp>
#include <boost/type_traits/remove_reference.hpp>

// C++11
#if defined(__cplusplus) && __cplusplus >= 201103L

#define BOOST_UBLAS_CPP_GE_2011

#elif BOOST_MSVC >= 1800

#define BOOST_UBLAS_CPP_GE_2011

#else

#undef BOOST_UBLAS_CPP_GE_2011 // Make sure no one defined it

#endif

// Microsoft Visual C++
#if defined (BOOST_MSVC) && ! defined (BOOST_STRICT_CONFIG)

// Version 7.1
#if BOOST_MSVC == 1310
// One of these workarounds is needed for MSVC 7.1 AFAIK
// (thanks to John Maddock and Martin Lauer).
#if !(defined(BOOST_UBLAS_NO_NESTED_CLASS_RELATION) || defined(BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION))
#define BOOST_UBLAS_NO_NESTED_CLASS_RELATION
#endif

#endif

#endif

// GNU Compiler Collection
#if defined (__GNUC__) && ! defined (BOOST_STRICT_CONFIG)

#if __GNUC__ >= 4 || (__GNUC__ >= 3 && __GNUC_MINOR__ >= 4)
// Specified by ABI definition see GCC bug id 9982
#define BOOST_UBLAS_USEFUL_ARRAY_PLACEMENT_NEW
#endif

#if __GNUC__ < 3
#define BOOST_UBLAS_UNSUPPORTED_COMPILER 1
#endif

#endif

// Intel Compiler
#if defined (BOOST_INTEL) && ! defined (BOOST_STRICT_CONFIG)

#if defined (BOOST_INTEL_LINUX) && (BOOST_INTEL_LINUX >= 800)
// By inspection of compiler results
#define BOOST_UBLAS_USEFUL_ARRAY_PLACEMENT_NEW
#endif

#if (BOOST_INTEL < 700)
#define BOOST_UBLAS_UNSUPPORTED_COMPILER 1
#endif

// Define swap for index_pair and triple.
#if (BOOST_INTEL <= 800)
namespace boost { namespace numeric { namespace ublas {
    template<class C, class IC>
    class indexed_iterator;

    template<class V>
    class index_pair;
    template<class M>
    class index_triple;
}}}

namespace std {
    template<class V>
    inline
    void swap (boost::numeric::ublas::index_pair<V> i1, boost::numeric::ublas::index_pair<V> i2) {
        i1.swap (i2);
    }
    template<class M>
    inline
    void swap (boost::numeric::ublas::index_triple<M> i1, boost::numeric::ublas::index_triple<M> i2) {
        i1.swap (i2);
    }
    // iter_swap also needed for ICC on Itanium?
    template<class C, class IC>
    inline
    void iter_swap (boost::numeric::ublas::indexed_iterator<C, IC> it1,
                    boost::numeric::ublas::indexed_iterator<C, IC> it2) {
        swap (*it1, *it2);
    }
}
#endif

#endif

// Comeau compiler - thanks to Kresimir Fresl
#if defined (__COMO__) && ! defined (BOOST_STRICT_CONFIG)

// Missing std::abs overloads for float types in <cmath> are in <cstdlib>
#if defined(__LIBCOMO__) && (__LIBCOMO_VERSION__ <= 31)
#include <cstdlib>
#endif

#endif

// PGI compiler
#ifdef __PGIC__
#define BOOST_UBLAS_UNSUPPORTED_COMPILER 0
#endif

//  HP aCC C++ compiler
#if defined (__HP_aCC) && ! defined (BOOST_STRICT_CONFIG)
#  if (__HP_aCC >= 60000 )
#    define BOOST_UBLAS_USEFUL_ARRAY_PLACEMENT_NEW
#endif
#endif

//  SGI MIPSpro C++ compiler
#if defined (__sgi) && ! defined (BOOST_STRICT_CONFIG)

// Missing std::abs overloads for float types in <cmath> are in <cstdlib>
// This test should be library version specific.
#include <cstdlib>

#if __COMPILER_VERSION >=650
// By inspection of compiler results - thanks to Peter Schmitteckert
#define BOOST_UBLAS_USEFUL_ARRAY_PLACEMENT_NEW
#endif

#endif

// Metrowerks Codewarrior
#if defined (__MWERKS__) && ! defined (BOOST_STRICT_CONFIG)

// 8.x
#if __MWERKS__ <= 0x3003
#define BOOST_UBLAS_UNSUPPORTED_COMPILER 1
#endif

#endif

// Detect other compilers with serious defects - override by defineing BOOST_UBLAS_UNSUPPORTED_COMPILER=0
#ifndef BOOST_UBLAS_UNSUPPORTED_COMPILER
#if defined(BOOST_NO_FUNCTION_TEMPLATE_ORDERING) || defined(BOOST_NO_SFINAE) || defined(BOOST_NO_STDC_NAMESPACE)
#define BOOST_UBLAS_UNSUPPORTED_COMPILER 1
#endif
#endif

// Cannot continue with an unsupported compiler
#if defined(BOOST_UBLAS_UNSUPPORTED_COMPILER) && (BOOST_UBLAS_UNSUPPORTED_COMPILER != 0)
#error Your compiler and/or configuration is unsupported by this verions of uBLAS. Define BOOST_UBLAS_UNSUPPORTED_COMPILER=0 to override this message. Boost 1.32.0 includes uBLAS with support for many older compilers.
#endif

// Enable performance options in RELEASE mode
#if defined (NDEBUG) || defined (BOOST_UBLAS_NDEBUG)

#ifndef BOOST_UBLAS_INLINE
#define BOOST_UBLAS_INLINE inline
#endif

// Do not check sizes!
#define BOOST_UBLAS_USE_FAST_SAME

// NO runtime error checks with BOOST_UBLAS_CHECK macro
#ifndef BOOST_UBLAS_CHECK_ENABLE
#define BOOST_UBLAS_CHECK_ENABLE 0
#endif

// NO type compatibility numeric checks
#ifndef BOOST_UBLAS_TYPE_CHECK
#define BOOST_UBLAS_TYPE_CHECK 0
#endif

// Disable performance options in DEBUG mode
#else

#ifndef BOOST_UBLAS_INLINE
#define BOOST_UBLAS_INLINE
#endif

// Enable runtime error checks with BOOST_UBLAS_CHECK macro. Check bounds etc
#ifndef BOOST_UBLAS_CHECK_ENABLE
#define BOOST_UBLAS_CHECK_ENABLE 1
#endif

// Type compatibiltity numeric checks
#ifndef BOOST_UBLAS_TYPE_CHECK
#define BOOST_UBLAS_TYPE_CHECK 1
#endif

#endif

/*
 * Type compatibility checks
 *  Control type compatibility numeric runtime checks for non dense matrices.
 *  Require additional storage and complexity
 */
#if BOOST_UBLAS_TYPE_CHECK
template <class Dummy>
struct disable_type_check
{
    static bool value;
};
template <class Dummy>
bool disable_type_check<Dummy>::value = false;
#endif
#ifndef BOOST_UBLAS_TYPE_CHECK_EPSILON
#define BOOST_UBLAS_TYPE_CHECK_EPSILON (type_traits<real_type>::type_sqrt (std::numeric_limits<real_type>::epsilon ()))
#endif
#ifndef BOOST_UBLAS_TYPE_CHECK_MIN
#define BOOST_UBLAS_TYPE_CHECK_MIN (type_traits<real_type>::type_sqrt ( (std::numeric_limits<real_type>::min) ()))
#endif

/*
 * General Configuration
 */

// Proxy shortcuts overload the alreadly heavily over used operator ()
//#define BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS

// In order to simplify debugging is is possible to simplify expression template
// so they are restricted to a single operation
// #define BOOST_UBLAS_SIMPLE_ET_DEBUG

// Use invariant hoisting.
// #define BOOST_UBLAS_USE_INVARIANT_HOISTING

// Use Duff's device in element access loops
// #define BOOST_UBLAS_USE_DUFF_DEVICE

// Choose evaluation method for dense vectors and matrices
#if !(defined(BOOST_UBLAS_USE_INDEXING) || defined(BOOST_UBLAS_USE_ITERATING))
#define BOOST_UBLAS_USE_INDEXING
#endif
// #define BOOST_UBLAS_ITERATOR_THRESHOLD 0

// Use indexed iterators - unsupported implementation experiment
// #define BOOST_UBLAS_USE_INDEXED_ITERATOR

// Alignment of bounded_array type
#ifndef BOOST_UBLAS_BOUNDED_ARRAY_ALIGN
#define BOOST_UBLAS_BOUNDED_ARRAY_ALIGN
#endif

// Enable different sparse element proxies
#ifndef BOOST_UBLAS_NO_ELEMENT_PROXIES
// Sparse proxies prevent reference invalidation problems in expressions such as:
// a [1] = a [0] = 1        Thanks to Marc Duflot for spotting this.
// #define BOOST_UBLAS_STRICT_MAP_ARRAY
#define BOOST_UBLAS_STRICT_VECTOR_SPARSE
#define BOOST_UBLAS_STRICT_MATRIX_SPARSE
// Hermitian matrices use element proxies to allow assignment to conjugate triangle
#define BOOST_UBLAS_STRICT_HERMITIAN
#endif

// Define to configure special settings for reference returning members
// #define BOOST_UBLAS_REFERENCE_CONST_MEMBER
// #define BOOST_UBLAS_PROXY_CONST_MEMBER

// Include type declerations and functions
#include <boost/numeric/ublas/fwd.hpp>
#include <boost/numeric/ublas/detail/definitions.hpp>

#endif