In This Package:

EigenSystems.cpp

Go to the documentation of this file.

// $Id: EigenSystems.cpp,v 1.4 2008/07/28 08:11:22 truf Exp $
// ============================================================================
// CVS tag $Name: v3r3 $ 
// ============================================================================
// $Log: EigenSystems.cpp,v $
// Revision 1.4  2008/07/28 08:11:22  truf
// just to please cmt or cvs or whatever
//
// Revision 1.3  2006/03/09 16:48:15  odescham
// v2r1 - migrated to LHCb v20r0 - to be completed
//
// Revision 1.2  2004/09/10 14:06:58  ibelyaev
//  fix a stupid bug!
//
// ============================================================================
// Include files
// ============================================================================
// local
// ============================================================================
#include "SoUtils/EigenSystems.h"
// ============================================================================
// use GSL
// ============================================================================
// GSL
// ============================================================================
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_eigen.h"
// ============================================================================

// ============================================================================
// ============================================================================

// ============================================================================
// ============================================================================
StatusCode SoUtils::eigensystem 
( const Gaudi::SymMatrix3x3&      matrix       ,
  Gaudi::Vector3&                  eigenvalues  ,
  std::vector<Gaudi::Vector3> &  eigenvectors ) 
{
  // dimension of matrix 
  //  const size_t num = matrix.num_row() ;
  const size_t num = 3 ;
  // reset the output values 
  eigenvalues  = Gaudi::Vector3 (0,0,0);
  eigenvectors = std::vector<Gaudi::Vector3>( num , Gaudi::Vector3( 0, 0 , 0 ) ) ;
  if( 0 == num ) { return StatusCode::SUCCESS ; }                    // RETURN
  
  // allocate GSL matrices 
  gsl_matrix* mtrx = gsl_matrix_alloc( num , num ) ;
  gsl_matrix* evec = gsl_matrix_alloc( num , num ) ;
  if ( 0 == mtrx  ) 
  { GSL_ERROR_VAL("SoCaloUtils::eigenvalues: the matrix(1) is not allocated", 
                  1000 , StatusCode( 1000 ) ) ; }                  // RETURN
  if ( 0 == evec  ) 
  { GSL_ERROR_VAL("SoCaloUtils::eigenvalues: the matrix(2) is not allocated", 
                  1000 , StatusCode( 1000 ) ) ; }                  // RETURN
  
  // allocate working space 
  gsl_eigen_symmv_workspace* wspace = gsl_eigen_symmv_alloc( num ) ;
  if ( 0 == wspace  ) 
  { GSL_ERROR_VAL("SoCaloUtil::eigenvalues: workspace is not allocated",
                  1001 , StatusCode ( 1001 ) ) ; }               // RETURN
  
  // allocate the vector of eigenvalues 
  gsl_vector* eval  = gsl_vector_alloc ( num ) ;
  gsl_vector* vcol  = gsl_vector_alloc ( num ) ;
  if ( 0 == eval    ) 
  { GSL_ERROR_VAL("SoCaloUtil::eigenvalues: vector(1) is not allocated",
                  1002 , StatusCode ( 1002 ) ) ; }             // RETURN  
  if ( 0 == vcol  ) 
  { GSL_ERROR_VAL("SoCaloUtil::eigenvalues: vector(2) is not allocated",
                  1002 , StatusCode ( 1002 ) ) ; }             // RETURN  
  
  // copy  symmetric matrix into GSL matrix 
  { for( size_t i = 0 ; i < num ; ++i ) 
  { for( size_t j = 0 ; j < num ; ++j ) 
  { gsl_matrix_set ( mtrx , i , j , matrix( i + 1 , j + 1 ) ) ; } } }
  
  int ierror = gsl_eigen_symmv ( mtrx , eval , evec , wspace ) ;
  if ( ierror      )
  { GSL_ERROR_VAL("SoCaloUtil::eigenvalues: values are not computed",
                  1010 + ierror , StatusCode ( 1010 + ierror ) ) ; } 
  
  // decode  GSL values into Vector
  { for( size_t i = 0 ; i < num ; ++i ) 
  { eigenvalues( i + 1 ) = gsl_vector_get( eval , i ) ; } }
  
  
  { // get eigenvectors 
    for( size_t i = 0 ; i < num ; ++i ) 
    {
      ierror = gsl_matrix_get_col( vcol , evec , i ) ;
      if ( ierror ) 
      { GSL_ERROR_VAL("SoCaloUtil::eigenvalues: invalid matrix operation",
                      1100 + ierror , StatusCode ( 1100 + ierror ) ) ; } 
      // copy eigenvector 
      for( size_t j = 0 ; j < num ; ++j ) 
      { eigenvectors[i]( j + 1 ) = gsl_vector_get( vcol , j ) ; }
    } 
  }
  
  // free everything at the end 
  if ( 0 != mtrx   ) { gsl_matrix_free      ( mtrx   ) ; }
  if ( 0 != evec   ) { gsl_matrix_free      ( evec   ) ; }
  if ( 0 != wspace ) { gsl_eigen_symmv_free ( wspace ) ; }
  if ( 0 != eval   ) { gsl_vector_free      ( eval   ) ; }
  if ( 0 != vcol   ) { gsl_vector_free      ( vcol   ) ; }
  
  return StatusCode::SUCCESS ;
  
};

// ============================================================================
// ============================================================================
StatusCode SoUtils::eigenvalues
( const Gaudi::SymMatrix3x3&     matrix       ,
  Gaudi::Vector3&              eigenvalues  )  
{
  // dimension of matrix 
  //  const size_t num = matrix.num_row() ;
  const size_t num = 3 ;
  // reset the output values 
  eigenvalues = Gaudi::Vector3( 0, 0 , 0 );
  if ( 0 == num ) { return StatusCode::SUCCESS ; }             // RETURN
  
  // alocate GSL matrix 
  gsl_matrix* mtrx = gsl_matrix_alloc( num , num ) ;
  if ( 0 == mtrx    ) 
  { GSL_ERROR_VAL("SoCaloUtils::eigenvalues: the matrix is not allocated", 
                  1000 , StatusCode( 1000 ) ) ; }               // RETURN
  
  // allocate working space 
  gsl_eigen_symm_workspace* wspace = gsl_eigen_symm_alloc( num ) ;
  if ( 0 == wspace  ) 
  { GSL_ERROR_VAL("SoCaloUtil::eigenvalues: workspace is not allocated",
                  1001 , StatusCode ( 1001 ) ) ; }             // RETURN  
  
  // allocate the vector of eigenvalues 
  gsl_vector* eval  = gsl_vector_alloc ( num ) ;
  if ( 0 == eval    ) 
  { GSL_ERROR_VAL("SoCaloUtil::eigenvalues: vector is not allocated",
                  1002 , StatusCode ( 1002 ) ) ; }             // RETURN  
  
  // copy CLHEP symmetric matrix into GSL matrix 
  { for( size_t i = 0 ; i < num ; ++i ) 
  { for( size_t j = 0 ; j < num ; ++j ) 
  { gsl_matrix_set ( mtrx , i , j , matrix( i + 1 , j + 1 ) ) ; } } }
  
  // use GSL 
  int ierror = gsl_eigen_symm( mtrx , eval , wspace ) ;
  if ( ierror      ) 
  { GSL_ERROR_VAL("SoCaloUtil::eigenvalues: values are not computed",
                  1010 + ierror , StatusCode ( 1010 + ierror ) ) ; }
  
  // decode  GSL values into Vector
  { for( size_t i = 0 ; i < num ; ++i ) 
  { eigenvalues( i + 1 ) = gsl_vector_get( eval , i ) ; } }
  
  // free everything at the end 
  if ( 0 != mtrx   ) { gsl_matrix_free     ( mtrx   ) ; }
  if ( 0 != wspace ) { gsl_eigen_symm_free ( wspace ) ; }
  if ( 0 != eval   ) { gsl_vector_free     ( eval   ) ; }
  
  return StatusCode::SUCCESS ;
};

// ============================================================================
// The End 
// ============================================================================

---