Re: [ublas] CLAPACK and GEEV

From: Peter Schmitteckert (peter_at_[hidden])
Date: 2004-12-17 06:20:58

Salut,

On Fri, 17 Dec 2004, Daniel Z wrote:

> body{font:12px Arial;margin:3px;overflow-y:auto;overflow-x:auto}p{margin:0px;}blockquote, ol, ul{margin-top:0px;margin-bottom:0px;} Hello,
>
> i intend to use a eigenvalue decompostion, but i dont manage to get it work, and since it is not documented it is very hard to me. Could
> someone explain me how to get geev work?
>
> Thanks.

as an implementation hint you can look at the following function:

/*
       SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
           CHARACTER     JOBVL, JOBVR
           INTEGER       INFO, LDA, LDVL, LDVR, LWORK, N
           DOUBLE        PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ), WI( * ), WORK( * ), WR( * )
*/

int Diagonalize
  (  boost::numeric::ublas::vector< std::complex<double>, boost::numeric::ublas::unbounded_array< std::complex<double> > >& E
  ,  boost::numeric::ublas::matrix< double, boost::numeric::ublas::row_major, boost::numeric::ublas::unbounded_array< double > >& A
  ,  boost::numeric::ublas::matrix< double, boost::numeric::ublas::row_major, boost::numeric::ublas::unbounded_array< double > >& UL
  ,  boost::numeric::ublas::matrix< double, boost::numeric::ublas::row_major, boost::numeric::ublas::unbounded_array< double > >& UR
  )
{
    const int Nx = A.size1();
    const int Ny = A.size2();
    int Info;

    double *pER = new double[Nx];
    double *pEI = new double[Ny];
    double dummy;
    int LWork;

#ifndef D_USE_OLDLAPACK
   dgeev( 'V', 'V', Nx, & A(0,0), Nx, pER, pEI, &UL(0,0), Nx, & UR(0,0), Nx, &dummy, -1, Info );
   LWork = int( dummy);
#else
   LWork = 8*Nx;
#endif

  double* pWork = new double[ LWork];

   dgeev( 'V', 'V', Nx, & A(0,0), Nx, pER, pEI, &UL(0,0), Nx, &UR(0,0), Nx, pWork, LWork, Info );

   std::cerr << "# *** GEEV-Info: " << Info << " LWork: " << LWork << std::endl;

   for( int i=0; i<Nx; ++i)
   {
      E(i) = std::complex<double>( pER[i],  pEI[i] );
   }

   delete [] pWork;
   delete [] pER;
   delete [] pEI;

  return Info;
}

Best wishes,
Peter