cvrSymmetricEigenSystem.h

Go to the documentation of this file.

/*
 * Copyright (C) 1998-2005
 * Lehrstuhl fuer Technische Informatik, RWTH-Aachen, Germany
 *
 *
 * This file is part of the Computer Vision and Robotics Library (CVR-Lib)
 *
 * The CVR-Lib is free software; you can redistribute it and/or
 * modify it under the terms of the BSD License.
 *
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 *
 * 3. Neither the name of the authors nor the names of its contributors may be
 *    used to endorse or promote products derived from this software without
 *    specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */


/**
 * \file   cvrSymmetricEigenSystem.h
 *         Contains the template class symmetricEigenSystem<T> that
 *         uses LAPACK if available.
 * \author Thomas Rusert
 * \author Peter Doerfler
 * \date   16.06.99
 *
 * $Id: cvrSymmetricEigenSystem.h,v 1.4 2006/09/07 13:38:37 doerfler Exp $
 */

#ifndef _CVR_SYMMETRIC_EIGEN_SYSTEM_H_
#define _CVR_SYMMETRIC_EIGEN_SYSTEM_H_

#include "cvrMatrix.h"
#include "cvrVector.h"
#include "cvrLinearAlgebraFunctor.h"

#include "cvrTypeInfo.h"

#ifdef HAVE_LAPACK
#  include "cvrLapackInterface.h"
#endif

namespace cvr {
  /**
   * Computes the Eigenvectors and Eigenvalues of a symmetric,
   * positive definite, real matrix.
   *
   * The most common source of such matrices is a covariance matrix
   * (see cvr::secondOrderStatistics).
   *
   * This class uses the LAPACK functions by default if they are
   * installed on your system. Otherwise a simple Jacobi algorithm is
   * used. This can also be chosen deliberately by setting
   * parameters::useLapack to false.
   *
   * The eigenvectors can be sorted according to their eigenvalues by
   * setting parameters::sort. A fixed number of eigenvectors can be
   * set via parameters::dimension. For LAPACK this also reduces
   * computation time.
   *
   * Please note that the eigenvector matrices will contain the
   * eigenvectors in the COLUMNS and not in the rows, as could be
   * expected. This avoids the requirement of transposing matrices in
   * eigenvector-based transformations like PCA!
   *
   * @ingroup gLinearAlgebra, gStatistics
   */
  template <typename T>
  class symmetricEigenSystem : public linearAlgebraFunctor
#ifdef HAVE_LAPACK
    , public lapackInterface
#endif
  {
  public:

    /**
     * eigenSystem parameter class
     */
    class parameters : public linearAlgebraFunctor::parameters {
    public:
      /**
       * default constructor
       */
      parameters();

      /**
       * copy constructor
       */
      parameters(const parameters& other);

      /**
       * copy member.
       */
      parameters& copy(const parameters& other);

      /**
       * Returns the name of this class.
       */
      const std::string& name() const;

      /**
       * returns a pointer to a clone of the parameters.
       */
      virtual parameters* clone() const;

      /**
       * returns a pointer to a clone of the parameters.
       */
      virtual parameters* newInstance() const;

      /**
       * write the parameters in the given ioHandler
       * \param handler the ioHandler to be used
       * \param complete if true (the default) the enclosing begin/end will
       *        be also written, otherwise only the data block will be written.
       * \return true if write was successful
       */
      virtual bool write(ioHandler& handler, const bool complete=true) const;

      /**
       * read the parameters from the given ioHandler
       * \param handler the ioHandler to be used
       * \param complete if true (the default) the enclosing begin/end will
       *        be also written, otherwise only the data block will be written.
       * \return true if write was successful
       */
      virtual bool read(ioHandler& handler, const bool complete=true);


      //****************************
      // the parameters
      //****************************

      /**
       * If true, the eigenvalues and corresponding eigenvectors will be
       * sorted in decreasing order of the eigenvalues.
       *
       * For complex-valued matrices a complex number will be considered
       * "greater than" another one if its real part is greater, or, if both
       * real components are identical, if its imaginary part is greater.
       * 
       * Default value: false
       */
      bool sort;

      /**
       * Mumber of dimensions calculated.
       *
       * If set to zero, then all eigenvectors and eigenvalues are calculated.
       *
       * This option is only provided for compatibility with the
       * fastEigenSystem functor based on LAPACK, but it does not make (yet)
       * the computation of the eigensolution any faster.  It just will cut
       * the already computed complete solution to the desired size!
       *
       * Default value: 0 (implying all eigenvalues will be computed)
       */
      int dimensions;
    };


    /**
     * Default constructor
     */
    symmetricEigenSystem();

    /**
     * Default constructor with parameters
     */
    symmetricEigenSystem(const parameters& params);


    /**
     * Calculate the Eigensystem with \a eigenvalues and \a eigenvalues
     * from the given data in \a theMatrix (Samples in rows). Due to the
     * symmetricity of the data it suffices if \a theMatrix is upper
     * triangular.
     *
     * @param theMatrix  Real symmetric matrix to be transformed.
     * Only the diagonal and above diagonal elements have
     * to be set.
     * @param eigenvalues  Elements will contain the eigenvalues.
     * @param eigenvectors Columns will contain the eigenvectors corresponding
     * to the eigenvalues.
     * returns true if successful, false otherwise.
     */
    virtual bool apply(const matrix<T>& theMatrix,
                       vector<T>& eigenvalues,
                       matrix<T>& eigenvectors) const;


    /**
     * Convenience function that calculates the Eigensystem like
     * apply, but uses the given argument \a dimensions instead of
     * those given in parameters::dimensions.
     *
     * @param theMatrix  Real symmetric matrix to be transformed.
     * Only the diagonal and above diagonal elements have
     * to be set.
     * @param eigenvalues  Elements will contain the eigenvalues.
     * @param eigenvectors Columns will contain the eigenvectors corresponding
     * to the eigenvalues.
     * @param dimensions number of dimensions in eigenspace
     * returns true if successful, false otherwise.
     */
    virtual bool reducedEigenSystem(const matrix<T>& theMatrix,
                                    vector<T>& eigenvalues,
                                    matrix<T>& eigenvectors,
                                    const int dimensions) const;

    /**
     * Returns the name of this class.
     */
    const std::string& name() const;

    /**
     * Returns a pointer to a clone of this functor.
     */
    virtual symmetricEigenSystem* clone() const;

    /**
     * Returns a pointer to a new instance of this functor.
     */
    virtual symmetricEigenSystem* newInstance() const;

    /**
     * returns the current parameters.
     */
    const parameters& getParameters() const;

  protected:
#ifdef HAVE_LAPACK
  private:
    /**
     * Contained type in T
     * 
     * This is a trick of generic programming to be an alias of type T if it is
     * a normal type (double, float) or the contained type A if T is
     * complex<A>.
     */
    typedef typename typeInfo<T>::value_type value_type;

    inline void rotateT(T& g,T& h,matrix<T>& a,
                        const int i,const int j,const int k,const int l,
                        const T s,const T tau) const;

    inline void rotate(T& g,T& h,matrix<T>& a,
                       const int i,const int j,const int k,const int l,
                       const T s,const T tau) const;


    /**
     * Apply the LAPACK function with the correct parameters
     * 
     * The matrix \c theMatrix will be modified, since LAPACK uses the unused
      * triangular submatrix as buffer.
     */
    bool applyLapack(matrix<T>& theMatrix,
                     vector<T>& eigenvalues,
                     matrix<T>& eigenvectors,
                     const int dimensions) const;
    
    /**
     * Calls the lapack function of eigensystem of symmetric matrix
     * specialized for float and double
     */
    int evr(char* jobz, char* range, char* uplo,
            integer* n, T* a, integer* lda,
            value_type* vl, value_type* vu,
            integer* il, integer* iu,
            value_type* abstol,
            integer* m, value_type* w,
            T* z, integer* ldz, integer* isuppz,
            T* work, integer* lwork,
            value_type* rwork,integer* lrwork,
            integer* iwork, integer* liwork,
            integer* info) const;
#endif

  };
}


#endif