|
|
Blitz Support : |
From: carlos.rega_at_[hidden]
Date: 2004-06-01 10:15:16
Hi Thomas
[...]
** For my problem, I need SVD, FFTs, QMRES, PLU, QR, ....*
Blitz++ does not currently provide any of these. However, there are
numerous C++ and C packages out there which do, and it is easy to move
data back and forth between Blitz++ and other libraries. See these
terms in the index: creating an array from pre-existing data, `data()',
`stride()', `extent()', `fortranArray'. For a list of other numerical
C++ libraries, see the Object Oriented Numerics Page at
`http://oonumerics.org/oon/'.
The FAQ answer, plus a trawl through the blitz-support archives, leads
me to believe that interfacing Blitz++ to LAPACK is probably straight
foward, however, I am a novice when it comes to both Blitz++ and LAPACK.
[My Question:] Has anyone built a Blitz++-LAPACK interface? i.e. a C++
wrapper around LAPACK that can be called from Blitz++?
Failing that, could anyone please post/email a small example of a
Blitz++ access to LAPACK, i.e. .cpp example file(s) with a Makefile?
[...]
I am working on a class interface to LAPACK (specifically the Intel
implementation, MKL) using blitz++ as the matrix library. As you say the
implementation is not terribly complicated. Here is an outline of how I
propose tackling the problem.
For the factorisations (SVD, QR ...) you will have something like:
template<class T> Factorisation
as a superclass that provides the interface, mainly a pure virtual method
(DoFactor) that will perform the factorisation.
For instance, the QR factorisation would be something like:
template <class T> class QRFactorisation : virtual public Factorisation<T>{
public:
/// default constructor
CQRFactorisation(void):mTau(Array<T,1>(FortranArray<1>())){;}
/// performs the factorisation
virtual void DoFactor(const Array<T,2>&);
/// returns the Q matrix
Array<T,2> Q(void);
/// returns the R matrix
Array<T,2> R(void);
private:
/// contains the tau vector
Array<T,1> mTau;
};
The FortranArray thing makes sure that the data is stored in the order
that LAPACK expects. Now, in principle, DoFactor only has to take a
blitz++ object, use the data() method to get a pointer to pass to the
appropriate LAPACK function and you are done. There is a little snag,
though. If you want to use the C++ complex<T> types you will have to cast
the pointer you get from data() to the correct pointer type for LAPACK (in
the MKL case it is MKL_Complex16 for T = double and MKL_Complex8 for T =
float) and a bit more magic because of having to mix pointers to the
complex type and to the real type of the correct precision for the
workspace arrays. For that I create an auxiliary structure that returns
the correct type for what I want:
template<class T> struct MKL_Traits {
typedef T mType;
typedef T mPrec;
};
and specialise it for complex<float> and complex<double>
template<> struct MKL_Traits<complex<float> > {
typedef MKL_Complex8 mType;
typedef float mPrec;
};
template<> struct MKL_Traits<complex<double> > {
typedef MKL_Complex16 mType;
typedef double mPrec;
};
then you can use
MKL_Traits<T>::mType or MKL_Traits<T>::mPrec to obtain the correct type
for the pointers.
I haven't checked the standard LAPACK release yet, but I presume that they
will use some equivalent complex type, so this approach should work with
that as well.
Though we have a working example of this approach, the code is still
evolving very fast, and is not ready for release. However, when we have a
working version we intend to make it available for general use.
Hope this helps
Regards
Carlos A. Rega
Development Scientist
Malvern Instruments Ltd
Grovewood Rd, Malvern
WR14 1XZ
Tel: +44 (0)1684 581 304
Fax: +44 (0)1684 892 789
e-mail: carlos.rega_at_[hidden]
--------------------------------------------------------------------------------------------------------------------
This email and any files transmitted with it are confidential and
intended solely for the use of the individual or entity to whom
they are addressed.
If you have received this email in error please notify the
originator of the message.
Any views expressed in this message are those of the individual
sender, except where the sender specifies and with authority,
states them to be the views of Malvern Instruments.