|
|
Ublas : |
Subject: Re: [ublas] [bindings][lapack] Interface design
From: Karl Meerbergen (Karl.Meerbergen_at_[hidden])
Date: 2008-12-04 03:50:12
Dear all,
I am happy to see a discussion on this topic. I have started (a while
ago) an lapack toolbox in glas, which currently only contains the
functions eig() and qr(), where I am trying to mimic somehow the matlab
interface. Due to lack of time I stopped this a while ago. This
discussion may be a motivation to resume this work.
I quite like
lapack::eig( A, eigenvalues(e), right_eigenvectors(vr),
left_eigenvectors(vl) ) ;
or
lapack::eig( A, B, eigenvalues(e), right_eigenvectors(vr),
left_eigenvectors(vl) ) ;
This is technically not hard to do.
The order of the arguments could then be arbitrary or arguments that
need not be computed could be omitted. The eigenvalues always have to be
computed. Therefore, I think the argument eigenvalues(e) could just be
e. Workspaces and other options can be added in a similar fashion.
However, the bindings traits are not rich enough I believe to handle all
cases that LAPACK supports, e.g. tridiagonal matrices. The bindings
could be enriched for this purpose, i.e. we could add a
tridiagonal_matrix_traits as well.
Best,
Karl
Thomas Klimpel wrote:
> Rutger ter Borg wrote:
>
>> Suppose we take LAPACK's problem type as a base, not their actual driver
>> routines (Thomas' suggestion of 'xxev.hpp' is along this line). See LAPACK
>> User's Guide, especially tables 2.2, 2.3, 2.4, 2.5, etc.. In other words,
>> we might have functions called similar to these problem types (I've
>> included some suggestions after 'or'):
>>
>
> 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.
>
>
>
>> Next, let's take the difficult case mentioned by Thomas, geevx. I think we
>> can achieve an easy and consistent interface by compile-time
>> type-decoration of our variables. Herewith some examples which might not be
>> perfect at their first shot, it's for the idea, suppose we just say
>>
>> lapack::eigen( A, left_eigen_values( vl ) );
>> lapack::eigen( A, right_eigen_values( vl ) );
>> lapack::eigen( A, right_eigen_values( vr ), left_eigen_values( vl ),
>> eigen_vectors( W ) );
>>
>
> The decorations are a nice idea, since they allow to state the role of certain arguments.
>
>
>> lapack::eigen( symmetric_lower( A ), left_eigen_values( vl ) );
>> lapack::eigen( hermitian_upper( A ), right_eigen( v ) );
>>
>
> 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.
>
>
>> Now, why I am happy with this?
>>
>> * it is very readable
>>
>
> 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.
>
>
>> * we might dispatch not only to xgeevx: other algorithms may be selected at
>> compile-time, like the divide-and-conquer algorithm
>>
>
> The selection at compile-time is a good point.
>
>
>> * no expert lapack-knowledge is required by the user (a part of the posts
>> to this list is actually about finding the right lapack driver)
>>
>
> 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.
>
>
>> * there is zero performance penalty in spite of the clear syntax (exploit
>> power of C++)
>>
>
> Except compile times, but C++ programmers are used to it.
>
>
>
>> Now, if a user wants to have more flexibility and/or access to the
>> implementation details, we can offer them the stuff generated by the
>> bindings generator, that may be used for low-level C++ calls to lapack:
>>
>> * templated calls to all driver / computational routines, by calling
>> detail::xxxxx( .... ) (or without detail::, since the actual interface
>> has other names)
>>
>
> without detail please, the user should not try to exploit implementation details.
>
>
>> * detail::optimal_workspace_xxxx( .... );
>> * detail::minimal_workspace_xxxx( .... );
>>
>
> The workspace handling should not be done by different function names, but by a workspace argument which defaults to optimal_workspace.
>
>
>
>> What do you think? Let's figure out a good set of names for decoration...
>>
>
> 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 ...
>
> I'm not really sure whether this is the right way to go, but I tried to enumerate at least some of the potential roles in order to get a better feeling for this.
>
>
> How are you planing to implement this? I could imagine free template functions of the form
>
> struct left_eigen_vectors_tag;
> template <class T> boost::tagged_type<T,left_eigen_vectors_tag> left_eigen_vectors(T& ref) {
> return boost::tagged_type<T,left_eigen_vectors_tag>(ref);
> }
>
> , but I don't know whether there is such a thing as boost::tagged_type.
>
>
> Regards,
> Thomas
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
>