ublas & lapack

From: boost (boost_at_[hidden])
Date: 2002-01-30 17:58:27

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? i.e.
   one has to provide the workspace as arguments.
-> 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.
-> 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[0] fine ? 
-> How do I get the leading dimension for  matrix_ranges ?

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