Re: [ublas] [bindings][lapack] Interface design

Subject: Re: [ublas] [bindings][lapack] Interface design
From: Thomas Klimpel (Thomas.Klimpel_at_[hidden])
Date: 2008-12-03 17:53:14

Jesse Perla wrote:
> And if there is any way to add in  
> overloaded versions that returns the object as an alternative, that would  
> be amazing as well and make matlab porting even easier. For example:
> matrix<double, column_major> A(3,3);
> vector<double> b(3);
> ....
> vector<double> x = lapack::solve(A,b); //allocates x internally and returns  
> it, consistent with the ublas expression template syntax.... an extra  
> object copy is reasonable for those who want to use this notation since it  
> may not be possible to use the return value optimization

I would prefer to leave this sort of usability features to matrix libraries like GLAS or MTL4. The numeric bindings library would be in a very unfortunate situation when trying to figure out the desired return type, and even if it where able to figure it out, it would still not know how to call its constructors...

> and...
> 
> tuple<vector<double>, vector<double> > my_eigenvalues = eigenvalues(A);  
> //Both Left and right for example, or:
> vector<double> left_eigen = left_eigenvalue(A);

Please don't get me wrong, but the numeric bindings library is currently not focused on people believing that the left eigenvalues of a matrix are different from its right eigenvalues. So it's certainly good when one can learn something from Matlab about making the interface easier to use, but the focus should stay on providing C++ bindings against existing external libraries.

> Another thing is that I will often want to write generic functions that  
> operate on matrices but may not be aware of their structure (other than  
> row/column major), so this could help a great deal if the function looks at  
> the ublas matrix traits/type.
> For example, what I really want to do is (which may be what you meant in  
> your interface suggestion):
> symmetric_matrix<double, lower> A1 (3, 3);
> hermitian_matrix<std::complex<double>, lower> A2(3,3);
> mapped_matrix<double> A3(3, 3, 3 * 3); //For sparse
> 
> lapack::eigen( A1, left_eigen_values( vl ) );
> lapack::eigen( A2, right_eigen( v ) ); //calls the hermitian routine
> lapack::eigen( A3, right_eigen( v ) ); //calls the sparse routine.....

This feature request might be within the scope of the numeric bindings library, but I'm unsure about the potential consequences of implementing this. Maybe we should try and find out.

> Last, is there any way you can throw an exception if the matrix traits  
> doesn't have fortran column-major ordering in the bindings for LAPACK?

I think it would be better to prevent your code from even compiling. At least some bindings have such compile time checks, so it should be easy to add the missing checks consistently to all bindings.

Regards,
Thomas