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.
> 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?
> one has to provide the workspace as arguments.
Yes, IMO you should let the client allocate the workspace (at least
> -> 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.
> -> What's the "canonical" ublas way to specify the start of the
memory block ?
> Is &E 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.
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,
> #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<>,
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, & U(0,0), N,
> pWork, LWork, pRWork, LRWork, pIWork, LIWork, Info );
> delete  pWork;
> delete  pRWork;
> delete  pIWork;
> return Info;
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk