Re: [ublas] [bindings] Workspaces

Subject: Re: [ublas] [bindings] Workspaces
From: Thomas Klimpel (Thomas.Klimpel_at_[hidden])
Date: 2008-12-10 15:45:02

Rutger ter Borg wrote:
> I also see is that dgelsd.f defines two workspace arrays,
> which seems to be contradicting the assumption that the
> number of workspaces is a 1:1 connection to the value_type of the
underlying. 
>
> Are these correct observations?

Yes and No. I have a bit the impression that you think I should have
rewritten Jesse Manning's gelsd bindings to more closely match the style
of existing bindings like heevd before including it into the subversion
repository. My reasons for not doing this is that I want to go on with
the 64-bit porting, and I considered it a good idea to add Jesse
Manning's gels bindings before I change the integer types.

> Could it be that the option for 3 workspace arguments is missing?
> E.g., AFAICT, cgelsd could be called with three workspace arguments
> (WORK, RWORK, IWORK). At the moment,
> this is not supported through workspace.hpp.

This is fine, because n_workspace_args<value_type>::value and
detail::workspace1, detail::workspace2 refers to the distinction between
the cases with a single real workspace and a workspace for both real and
complex. The integer workspace is independent of this and the same for
the real and complex case.

Take a look at heevd for an example how to handle this case:

template <typename T, typename R>
void operator() (..., optimal_workspace, ...)

template <typename T, typename R, typename W, typename RW, typename WI>
void operator() (..., std::pair<detail::workspace2<W,RW>,
detail::workspace1<WI> > work, ...)

> Perhaps it's safer to select the number of workspaces by the number of
> workspaces defined in the Fortran code :-).

Not really.

> Another question: are there cases where the formulae for the
> minimum_workspace requirements differ for each value_type? (I.e., 4
cases)?

No, the double precision versions of the lapack routines are generated
automatically from the single precision versions.

> Or is two formulae the maximum (real vs. complex)? I remember I saw
> something along the lines of template specialisations based on
> real/complex/mixed scenarios, which would make it three, or is this
based
> on something else?

I guess you are referring to geev. The mixed case is when the matrix A
is real, but the user requests complex eigenvectors (for convenience
reasons...).

Regards,
Thomas