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

Subject: Re: [ublas] [bindings][lapack] Interface design
From: Rutger ter Borg (rutger_at_[hidden])
Date: 2008-12-04 04:23:34

Thomas Klimpel wrote:
> We should try to ensure that it's both easy to find the name of the
> binding for a given lapack routine, and easy to find the name of the
> lapack routine for a given binding. Using the LAPACK User's Guide for that
> purpose is sufficiently easy, as long as it's stated explicitly enough in
> the numeric bindings documentation where to look in the LAPACK User's
> Guide in order to deduce the names.

Yes, I agree, usability is most important. Also good to note that the number
of problem types is countable :-).

> symmetric_lower and hermitian_upper fall into a different category. These
> are some sort of type-refinement. I thought about these before in the
> context of interpreting a vector as a matrix, because a vector might be
> interpreted as a (1,n) or as a (n,1) matrix, and this sort of
> type-refinement allows to distinguish the two cases.

Ok. Sorry about the left/right eigenvalues, I meant eigenvectors of course.

> Indeed, but it also raises some questions. The boost parameter library
> also allows to specify the role (/name) of the arguments at the calling
> site. For me, this raises the question whether the role should always be
> specified, or only in the cases where it is necessary.

I'll look into the details of the Boost.Parameter Library. Perhaps it's a
bit overkill for what I meant.

> I don't like this one. From my point of view, the numeric bindings library
> is concerned with providing good C++ interfaces to existing external
> numerical libraries. So it should enable its users the application of
> typical C++ programming techniques like generic programming, but it should
> limit itself to remove the restrictions of the original interfaces caused
> by fortran or C. LAPACK is differnt from Matlab, and the numeric bindings
> library should not try to turn LAPACK into Matlab. The user of the numeric
> bindings library should still be willing to read the original
> documentation of the LAPACK routine he uses, in case he runs into
> problems.

I am not very familiar with Matlab, so I am not really sure how close the
syntax is. I agree that easing the LAPACK interface is the primary
objective, and this is quite an ambitious goal already :-). For me, this
means that it should be possible for a user to use (almost all of) LAPACK
in a convenient and error-preventing way. To elaborate a bit further, what
I am thinking of, is

 * a normal interface, based on LAPACK's problem types, which works fine for
at least 80% of the cases of everyday use. Here we have lapack::eigen etc..
This is handcrafted code.
 * a low-level interface, based on LAPACK's routines, the automatically
generated ones, with similar but probably less possibilities of the current
bindings, being able to call each and every LAPACK routine. Stuff like
lapack::geev etc.. Includes stuff to calculate the workspace.

The normal interface makes compile-time decisions to call which function
from the low-level interface. If someone really wants to call that slower
routine for the same job, or have an explicit answer on the optimal
workspace size, our answer would be to use the low-level interface.

> symmetric_lower, symmetric_upper, hermitian_lower, hermitian_upper,
> matrix_1_n, matrix_n_1 + other refined types for cases with insufficient
> information in the available types.
> 
> left_eigen_vectors, right_eigen_vectors, left_singular_vectors,
> right_singular_vectors, eigen_values, singular_values, right_hand_sides,
> permutation_matrix, index_range, value_range, abs_tolerance, transposed,
> reflectors ...

Perhaps we could also consider the possibility that decoration to be only
necessary and/or non-deductable from the argument types, and/or deviates
from the default? I.e., if one writes

// A, VR = matrices, lambda = vector
lapack::eigen( A, VR );
lapack::eigen( A, lambda );

he/she will just end up with the right eigenvectors of A in VR, and
eigenvalues in lamdba, whereas

lapack::eigen( A, left_eigenvector( VR ) );

will give the non-default?

I think eigenvector, eigenmatrix, eigenspace, etc. are written as one word,
and I also think there is a tendency to write words in a singular form.

right_eigenvector
left_eigenvector
eigenvalue

left_singular_vector
right_singular_vector
singular_value

Here we also have some overlap in functionality (could be reduced to
right_vector for eigen/svd cases), but I guess the current wording makes
everything clearer.

> How are you planing to implement this? I could imagine free template
> functions of the form

I am thinking of wrapping it in another type, and have its traits dispatched
to its containing type. See an example below. In this way, the decorated
arguments can be passed directly to the low-level C++ functions, so we can
focus on compile-time decisions in the dispatching part. 

template< typename T >
struct left_eigenvector_impl {
  typedef T containing_type;
  left_eigenvector_impl( T &ref ): m_ref( boost::ref( ref ) );
  T& operator() {
    return m_ref;
  }
  boost::reference_wrapper<T> m_ref;
};

template< typename T >
struct vector_traits< left_eigenvector_impl< T > >: unwrap_traits<
left_eigenvector_impl< T > >{};

template< typename T >
struct unwrap_traits: vector_traits< typename T::containing_type > {
  typedef vector_traits< typename T::containing_type >::pointer pointer;
  inline pointer storage( T &v ) {
    // pass reference
    return vector_traits< T::containing_type >::storage( v() ); 
  }
  // same for size and stride
};

template< typename T >
left_eigenvector_impl<T> left_eigenvector( T &ref ) {
  return left_eigenvector_impl<T>(ref);
}

I'm curious to hear what you guys think. 

Cheers,

Rutger