cvrLuDecomposition.h

Go to the documentation of this file.

/*
 * Copyright (C) 1998
 * 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   cvrLuDecomposition.h
 *         Compute the LU decomposition of a given matrix.
 * \author Thomas Rusert
 * \author Gustavo Quiros
 * \date   23.06.1999
 *
 * $Id: cvrLuDecomposition.h,v 1.8 2006/09/07 13:38:37 doerfler Exp $
 */

#ifndef _CVR_LU_DECOMPOSITON_H_
#define _CVR_LU_DECOMPOSITON_H_

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

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

namespace cvr {

#ifdef HAVE_LAPACK
  namespace internal {

    /**
     * Offers a direct interface to the LAPACK function getrf for
     * computing an LU decomposition with partial pivoting.
     *
     * This function is needed by other classes that use it as part of
     * the LAPACK processing chain.
     */
    template <typename T>
    class luDecompositionLapackInterface : public lapackInterface {
    public:
      /**
       *  GETRF computes an LU factorization of a general M-by-N matrix A
       *  using partial pivoting with row interchanges.
       *
       *  The factorization has the form
       *  \code
       *     A = P * L * U
       *  \endcode
       *  where P is a permutation matrix, L is lower triangular with unit
       *  diagonal elements (lower trapezoidal if m > n), and U is upper
       *  triangular (upper trapezoidal if m < n).
       */
      int getrf(integer* rows, integer* cols, T* a,
                integer* lda, integer *ipiv, integer* info) const;
    };
  }
#endif


  /**
   * LU decomposition functor.
   * Computes the LU decomposition of a square matrix.
   *
   * This class uses LAPACK if it is available.  Note that if LAPACK
   * is not used or not available, A <b>must</b> be of full rank!
   *
   * @ingroup gLinearAlgebra
   */
  template<typename T>
  class luDecomposition
    : public linearAlgebraFunctor
#ifdef HAVE_LAPACK
      , protected internal::luDecompositionLapackInterface<T>
#endif
  {

  public:

    /**
     * The parameters for the class luDecomposition
     */
    class parameters : public linearAlgebraFunctor::parameters {
    public:

      /**
       * Default constructor
       */
      parameters();

      /**
       * Copy constructor
       * @param other the parameters object to be copied
       */
      parameters(const parameters& other);

      /**
       * Destructor
       */
      ~parameters();

      /**
       * Copy the contents of a parameters object
       * @param other the parameters object to be copied
       * @return a reference to this parameters object
       */
      parameters& copy(const parameters& other);

      /**
       * Copy the contents of a parameters object
       * @param other the parameters object to be copied
       * @return a reference to this parameters object
       */
      parameters& operator=(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 new instance 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
      // ------------------------------------------------

    };

    /**
     * Default constructor
     */
    luDecomposition();

    /**
     * Construct a functor using the given parameters
     */
    luDecomposition(const parameters& par);

    /**
     * Copy constructor
     * @param other the object to be copied
     */
    luDecomposition(const luDecomposition& other);

    /**
     * Destructor
     */
    virtual ~luDecomposition();

    /**
     * onPlace version of apply.
     *
     * Given a matrix a[0..n-1][0..n-1], this routine replaces it by
     * the LU decomposition of a rowwise permutation of itself.
     * permutation[0..n-1] is an output vector that records the row
     * permutation effected by the partial pivoting. 'pivot' is
     * +/-1 depending on whether the number of row interchanges was even or
     * odd, respectively.
     * @return true, if the decomposition could be computed, false
     *          otherwise (typically, the matrix was singular).
     */
    bool apply(matrix<T>& a, vector<integer>& permutation, int& pivot) const;

    /**
     * onCopy version of apply.
     *
     * Given a matrix a[0..n-1][0..n-1], this routine returns a matrix
     * decomposition[0..n-1][0..n-1] which contains the LU decomposition of a
     * rowwise permutation of theMatrix. permutation[0..n-1] is an output
     * vector that records the row permutation effected by the partial
     * pivoting. 'pivot' is +/-1 depending on whether the number of
     * row interchanges was even or odd, respectively.
     * @return true, if the decomposition could be computed, false
     *         otherwise (typically, the matrix was singular).
     */
    bool apply(const matrix<T>& a, matrix<T>& lu,
         vector<integer>& permutation, int& pivot) const;

    /**
     * Copy data of "other" functor.
     */
    luDecomposition<T>& copy(const luDecomposition<T>& other);

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

    /**
     * Returns a pointer to a clone of the functor.
     */
    virtual luDecomposition<T>* clone() const;

    /**
     * Returns a pointer to a new instance of the functor.
     */
    virtual luDecomposition<T>* newInstance() const;

    /**
     * Returns used parameters
     */
    const parameters& getParameters() const;

    /**
     * Returns a new Matrix which contains the "L" part of the given
     * LU decomposition
     */
    matrix<T> extractL(const matrix<T>& lu) const;

    /**
     * Returns a new Matrix which contains the "U" part of the given
     * LU decomposition
     */
    matrix<T> extractU(const matrix<T>& lu) const;

    /**
     * Returns a new Matrix which contains the "L" part of the given
     * LU decomposition.  This is much faster than the on copy
     * version.
     */
    void extractL(const matrix<T>& lu, matrix<T>& l) const;

    /**
     * Returns a new Matrix which contains the "U" part of the given
     * LU decomposition.  This is much faster than the on copy
     * version.
     */
    void extractU(const matrix<T>& lu, matrix<T>& u) const;

    /**
     * Builds the permutation matrix for the given permutation vector
     *
     * @param pv The row permutation vector on input
     * @param pm The row permutation matrix on output
     */
    void buildPermutationMatrix(const vector<integer>& pv,
                                matrix<T>& pm) const;

  protected:

    static const T epsilon;

#ifdef HAVE_LAPACK
    bool applyLapack(matrix<T>& theMatrix,
                     vector<integer>& permutation,
                     int& pivot) const;
#endif

  };
}

#endif