Boost logo

Boost :

From: jhrwalter (walter_at_[hidden])
Date: 2002-02-01 15:22:02


--- In boost_at_y..., Peter Schmitteckert (boost) <boost_at_s...> wrote:
> Dear all,
>
> I would volunteer to write an interface for ublas to at least some
> Lapack functions.
 
Good idea.
 
> Before implementing them I'd like to have a few
> comments on my appended example, especially:
>
> -> should I use boost::scoped_array for the workspace pointers?
> -> should I avoid dynamic memory allocation within the routines?
i.e.
> one has to provide the workspace as arguments.
 
Yes, IMO you should let the client allocate the workspace (at least
optionally).
 
> -> naming: should I try to match Lapack names ? e.g. I'd prefer
> to implement the various diagonalization routines by overloading
> a 'diagonalize(···)' function instead of implementing
> 'hpevd', 'hevd', ... . Same applies to the various SVD routines.

Yes. This is one step of designing an OO wrapper for LAPACK.
 
> -> Am I correct that this makes only sense for the dense
> numerics::unbounded_array memory shemes.
 
Currently yes.

> -> What's the "canonical" ublas way to specify the start of the
memory block ?
> Is &E[0] fine ?

This could even work. In the long term I would prefer to add a
(somewhat dangerous) accessor to the internal data pointers.
 
> -> How do I get the leading dimension for matrix_ranges ?
 
If I understand your question correctly, you also need access to the
internal matrix reference of a matrix_range.
 
Regards
 
Joerg
 
P.S.: Do you currently try to interface CLAPACK? If yes, you could
get problems with 0 vs. 1 based indexing. Otherwise, you could get
problems with the matrix storage layout ;-)

> Best wishes,
> Peter
>
> -----------------------------------------------------------------
>
> #ifndef INCLUDES_PS_LaPack
> #define INCLUDES_PS_LaPack
>
> #include <complex>
>
> extern "C"
> {
> void zhpevd_( const char& jobz, const char& uplo,
> const int& N, std::complex<double>* pAP, double* W,
> std::complex<double>* pZ, const int& LDZ,
> std::complex<double>* pWork, const int& LWork,
> double* pRWork, const int& LRWork,
> int* pIWork, const int& LIWork, int& Info );
> }
>
> namespace NSp_LaPack // Yes, there should be a better name.
> {
> int Diagonalize(
> numerics::vector< double, numerics::forward,
> numerics::unbounded_array<double> >& E,
> numerics::hermitean_matrix< std::complex<double>,
> numerics::lower, numerics::row_major<>,
> numerics::unbounded_array<
std::complex<double> > >& H,
> numerics::matrix< std::complex<double>, numerics::row_major<>,
>
numerics::unbounded_array<std::complex<double> > >& U )
> {
> const int N = H.size1();
> const int LWork = N*N;
> const int LRWork = 1 + (5 + 2*N)*N;
> const int LIWork = 3 + 5*N;
> int Info;
>
> std::complex<double>* pWork = new std::complex<double>[ LWork];
> double* pRWork = new double[ LRWork ];
> int* pIWork = new int[LIWork];
>
> // 'U' corresponds to 'L' in C++
> zhpevd_( 'V', 'U', N, & H(0,0), & E[0], & U(0,0), N,
> pWork, LWork, pRWork, LRWork, pIWork, LIWork, Info );
>
> delete [] pWork;
> delete [] pRWork;
> delete [] pIWork;
>
> return Info;
> }
> }
>
> #endif


Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk