Re: ublas & lapack

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