Effective UBLAS

Effective uBLAS

For uBLAS Users and Developers

This document is intended both to facilitate discussion and to document the current status of uBLAS. Its positioning in the Boost Wiki should make it easier for new comers to avoid problems and understand the many complex issues uBLAS has to deal with.

If you are interested, you can join the ublas mailing list:

      (new) http://lists.boost.org/mailman/listinfo.cgi/ublas

and have a look at the main page of the boost project.

There is a Sourceforge site where development took place

      (old) http://sourceforge.net/projects/ublas/ 

The old mailing list is still available at

      (old) http://groups.yahoo.com/group/ublas-dev/

but no new members will be accepted. All discussion should run via the new list.

Boost 1.33.x

First uBLAS users and developers meeting

Further sources of information

Frequently asked Questions using uBLAS mostly performance issues.

Extensions and new usage patterns in uBLAS.

Examples - How to extend uBLAS.

Gunter Winkler's tricks and sparse matrix information (and a lot more)

     http://www.guwi17.de/ublas/index.html

Linear Algebra with uBLAS

LU Matrix Inversion - matrix inversion using uBLAS' LU-decomposition functions.

Effective UBLAS/Matrix Inversion.

Making Efficient use of uBLAS

Make sure you define NDEBUG

You can either rely on your build system to do this for you. Boost.build is highly recommended. Or you must add this yourself to your compiler options.

Controlling the complexity of nested products

 R = prod(A, prod(B,C)); 

uBLAS evaluates expressions without matrix of vector temporaries and honours the bracketing structure. However avoiding temporaies for tested product unnecessarly increases the complexity. Conversly by explictly using temporary matrices the complexity of a nested product can be reduced.

uBLAS will give a compile time error using the syntax above, and instead, it provides 3 alternative syntaxes for this purpose:

 temp_type T = prod(B,C); R = prod(A,T);   // Preferable if T is preallocated

 prod(A, temp_type(prod(B,C));

 prod(A, prod<temp_type>(B,C));

The 'temp_type' is important. Given A,B,C are all of the same type. Say matrix<float>, the choice is easy. However if the value_type is mixed (int with float or double) or the matrix type is mixed (sparse with symmetric) the best solution is not so obvious. It is up to you! It depends on numerical properties of A and the result of the prod(B,C).

Properties of symmetric_matrix

When assigning from a non another matrix type is only valid when its value is numerically (upper identical to lower) symmetric. This is checked as a precondition in the debug build. This check can be avoided (ignoring the other half of the matrix by using a symmetric_adaptor.

Fixed size vector or matrices

 namespace ublas = boost::numeric::ublas;
 class vector3 : public ublas::vector<double, ublas::bounded_array<double,3> >
 {
	typedef ublas::vector<double, ublas::bounded_array<double,3> > Base_vector;
 public:
	// Default construction
	vector3 () : Base_vector(3)
	{}
	// Construction and assignment from a uBLAS vector expression or copy assignment
	template <class R> vector3 (const ublas::vector_expression<R>& r) : Base_vector(r)
	{}
	template <class R> void operator=(const ublas::vector_expression<R>& r)
	{
		Base_vector::operator=(r);
	}
	template <class R> void operator=(const Base_vector& r)
	{
		Base_vector::operator=(r);
	}
 };

It should be noted that this only changes the storage uBLAS uses for the vector3. uBLAS will still use all the same algorithm (which assume a variable size) to manipulate the vector3. In practice this seems to have no negative impact on speed. The above runs just as quickly as a hand crafted vector3 which does not use uBLAS. The only negative impact is that the vector3 always store a "size" member which in this case is redundant.

uBLAS now has a bounded_vector type that effectively wraps the above.

 ublas::bounded_vector<double, 3>

'rank' of an Iterator

It is not the rank (number of independent row or columns) of the matrix!

Joerg's answer: It's something like the distinction over which rank the iterators have to traverse: 0 means it is an outer loop iterator, which is usually not able to determine (easily), if a row/column is empty and 1 means it is in inner loop iterator, which has the task to determine the next/previous element in a row or column.

edit: Thoughts about matrix iterators, see UBLAS matrix iterators.

Debugging in Visual Studio

It is possible to get vectors to display cleanly in VS debugger. Add the following in the autoexp.dat file in Common7/Packages?/Debugger?/ under [Visualizer] where it says DO NOT MODIFY.

 boost::numeric::ublas::bounded_vector<*,*>{
	children
	(
	    #array
	    (
			expr :	($c.data_.data_)[$i],
			size :	$c.data_.size_
		)
	)

    preview
    (
        #(
			"[",
            $c.data_.size_ ,
            "](",

            #array
            (
				expr : 	($c.data_.data_)[$i],
				size : 	$c.data_.size_
			),
			")"
		)
	)
 }

 boost::numeric::ublas::vector<*,*>{
	children
	(
	    #array
	    (
			expr :	($c.data_.data_)[$i],
			size :	$c.data_.size_
		)
	)

    preview
    (
        #(
			"[",
            $c.data_.size_ ,
            "](",

            #array
            (
				expr : 	($c.data_.data_)[$i],
				size : 	$c.data_.size_
			),
			")"
		)
	)
 }

Configuring uBLAS with Macros

• Automatically set by 'boost/numeric/ublas/config.hpp' based on NDEBUG. Makes the distinction between debug (safe) and release (fast) mode. Similar to STLport
  • Release mode (NDEBUG defined)
    • BOOST_UBLAS_INLINE Compiler dependant definition to control function inlining.
    • BOOST_UBLAS_USE_FAST_SAME
  • Debug mode
    • BOOST_UBLAS_CHECK_ENABLE Enable checking of indexs, iterators and parameters. Prevents out of bound access etc.
    • BOOST_UBLAS_TYPE_CHECK Enable additional checks for the results of expressions using non dense types. Picks up runtime error such as the assignment of a numerically non-symmetric matrix to symmertic_matrix. Use #define BOOST_UBLAS_TYPE_CHECK 0 to disable expensive numeric type checks. (Note: "structure check" would be a much better name.)
    • BOOST_UBLAS_TYPE_CHECK_EPSILON default: sqrt(epsilon), controls how large the difference between the expected result and the computed result may become. Increase this value if you are going to use near singular or badly scaled matrices. Please, refer to detail/matrix_assign.hpp for implementation of these type checks.

• Automatically set by 'boost/numeric/ublas/config.hpp' based on compiler and boost/config.hpp macro's. Augments the compiler deficiency workarounds already supplied by boost/config.hpp
  • BOOST_UBLAS_NO_NESTED_CLASS_RELATION A particularly nasty problem with VC7.1 Requires that uBLAS and the user use begin(it) rather then it.begin()
  • BOOST_UBLAS_NO_SMART_PROXIES Disable the automatic propagation of 'constantness' to proxies. Smart proxies automatically determine if the underling container they reference is constant or not. They adjust there definition of iterators and container access to reflect this constantness.

• For use by uBLAS authors to test implementation methods. Preset in config.hpp
  • BOOST_UBLAS_USE_INVARIANT_HOISTING
  • BOOST_UBLAS_USE_INDEXING
  • BOOST_UBLAS_USE_INDEXED_ITERATOR
  • BOOST_UBLAS_NON_CONFORMANT_PROXIES Gappy containers may be non-conformant, that is contain elements at different indices. Assigning between proxies (vector ranges for example) of these containers is difficult as the LHS may need insert new elements. This is slow.
  • BOOST_UBLAS_USE_DUFF_DEVICE Near useless on all platforms (see GCC's -funroll-loops)

• User options. Can be predefined by user before including any uBLAS headers. They may also be automatically defined for some compilers to work around compile bugs.
  • BOOST_UBLAS_USE_LONG_DOUBLE Enable uBLAS expressions that involve containers of 'long double'
  • BOOST_UBLAS_USE_INTERVAL Enable uBLAS expressions that involve containers of 'boost::numeric::interval' types
  • BOOST_UBLAS_SIMPLE_ET_DEBUG In order to simplify debugging is is possible to simplify expression templateso they are restricted to a single operation
  • BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
  • BOOST_UBLAS_NO_ELEMENT_PROXIES Disables the use of element proxies for gappy types.
  • The Gappy types (sparse, coordinate, compressed) store non-zero elements in their own containers. When new non-zero elements are assigned they must rearrange these containers. This invalidates references, iterators or pointers to these elements. This can happen at some surprising times such as the expression "a [1] = a [0] = 1;". Element proxies guarantee all such expressions will work as expected. However they bring their own restrictions and efficiency problems. For example as of Boost 1.30.0 they prevent the assignment of elements between different types.

• BOOST_UBLAS_REFERENCE_CONST_MEMBER Enable to allow refernces to be returned to fixed (zero or one) elements of triangular or banded matrices
• BOOST_UBLAS_NO_EXCEPTIONS Disable the use exceptions of uBLAS internal checks and error conditions. BOOST_NO_EXCEPTIONS has same effect.
• BOOST_UBLAS_SINGULAR_CHECK Check the for singularity in triangular solve() functions

uBLAS Current Issues

Bug reports

      http://sourceforge.net/tracker/?group_id=94014&atid=606391

warning C4224: nonstandard extension used : formal parameter 'lower' was previously defined as a type

Index differences may overflow

 int diff = a.index() - b.index();

size_t not consistently used

sparse_matrix::iterator not sparse for rank==0 (slow index)

Conformant sparse proxies