Lapack.hh

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/


/** @file Lapack.hh
 *  @ingroup g_mathalgo 
    BIAS c++ wrapper class for lapack fortran77 functions
    @author Daniel Grest, Jan Woetzel, Felix Woelk
    @status tested
*/
#ifndef BIASLAPACK_H
#define BIASLAPACK_H

#include "bias_config.h" 

#include <Base/Math/Matrix.hh>
#include <Base/Math/Vector.hh>

#include <vector>


/** @addtogroup g_mathalgo
    \brief C++ Wrappers for important Lapack fortran77 functions

    Functions from Lapack, e.g. solving systems of linear equations
    or computing eigenvalues
    For further information have a look here 
    <a href="http://www.netlib.org/lapack/lug/">
    http://www.netlib.org/lapack/lug/</a>.

    @{
*/

/** solve linear equations using LU factorization */
BIASMathAlgo_EXPORT
BIAS::Vector<double> 
Lapack_LU_linear_solve(const BIAS::Matrix<double> &A, 
                       const BIAS::Vector<double> &b);

/** linear least squares solves |Ax-b|=min 
    res=0 on success, anything else 
    rudimentary tested (woelk)
    @author woetzel  */
BIASMathAlgo_EXPORT
BIAS::Vector<double> 
Lapack_LLS_QR_linear_solve(const BIAS::Matrix<double> &A, 
                           const BIAS::Vector<double> &b, int &res);

/** as above but ignores success info (res) */
BIASMathAlgo_EXPORT
inline BIAS::Vector<double> 
Lapack_LLS_QR_linear_solve(const BIAS::Matrix<double> &A, 
                           const BIAS::Vector<double> &b)
{ int res; return Lapack_LLS_QR_linear_solve(A, b, res); }


/** weighted linear least squares
    solves |B^-1(b-Ax)|=min
    res=0 on success, anything else 
    rudimentary tested
    @author woelk 04 2003 */
BIASMathAlgo_EXPORT
BIAS::Vector<double> 
Lapack_WLLS_solve(const BIAS::Matrix<double> &A, 
                  const BIAS::Vector<double> &b, 
                  const BIAS::Matrix<double> &B, int &res);


/** @brief Coputes the Cholesky decomposition of a symmetric positive definit
    matrix A. 
    
    When ul='U', the factorization has the form 
       A = UL^T * UL 
    where UL is a upper triangular matrix. 
    When ul='L',the factorization has the form 
       A = UL * UL^T 
    where UL is a lower triangular matrix.  
    Wrapper for lapacks dpotf2.
    

    @author woelk 09/2007 */
BIASMathAlgo_EXPORT
int Lapack_Cholesky_SymmetricPositiveDefinit(const BIAS::Matrix<double> &A,
                                             BIAS::Matrix<double> &UL,
                         const char ul='U');


/** @brief computes Mahalanobis distance V^T Sigma^-1 V efficiently using 
 *         cholesky decomposition and exploitation of symmetry 
 *
 *         Sigma must be square and symmetric positive definite, same dim as V
 *
 *         About 10 times faster (for 6x6) than V^T svd.invert(Sigma) V,
 *         but it still has potential for optimization, see doc in code.
 *
 *  @author koeser 11/2007 */    
BIASMathAlgo_EXPORT
double SquaredMahalanobisDistance(const BIAS::Matrix<double> &Sigma,
                                  const BIAS::Vector<double> &d);

/** @brief computes squared Mahalanobis distance V^T Sigma^-1 V efficiently 
 *         using cholesky decomposition and exploitation of symmetry 
 *
 *         Sigma must be square and symmetric positive definite, same dim as V
 *
 *         About 10 times faster (for 6x6) than V^T svd.invert(Sigma) V,
 *         but it still has potential for optimization, see doc in code.
 *
 *  @author koeser 11/2007 */    
BIASMathAlgo_EXPORT
double MahalanobisDistance(const BIAS::Matrix<double> &Sigma,
                           const BIAS::Vector<double> &d);



// *********************** Eigenvalue problems *******************

/** @brief Solve symmetric eigenvalue problem (eigenvalues only)
    @note Wrapper for LAPACK routine dsyev.
    @status untested
    @author grest
*/
BIASMathAlgo_EXPORT
BIAS::Vector<double> 
Upper_symmetric_eigenvalue_solve(const BIAS::Matrix<double> &A);

/** @brief Compute eigenvalues and eigenvectors of a symmetric 
           real square matrix.
    @note Wrapper for LAPACK routine dsyevd.
    @author woelk 11/2006
*/
BIASMathAlgo_EXPORT
int Upper_symmetric_eigenvalue_solve(const BIAS::Matrix<double> &A,
                                     BIAS::Vector<double>& eigenvalues,
                                     std::vector<BIAS::Vector<double> >& eigenvectors);

/** @brief Solve symmetric eigenvalue problem for matrix in packed storage.
    @param[in] N Number of rows/cols of N x N matrix A
    @param[in] A Data vector of A in packed storage of length N*(N+1)/2.
               If upper triangle matrix is given, vector contains matrix rows
               above the diagonal (i.e. A11, A12, A13, ..., A1N, A22, A23, ...),
               otherweise vector contains matrix columns below the diagonal
               (i.e. A11, A21, A31, ..., AN1, A22, A32, ...).
    @param[out] eigVals Returns eigenvalues in ascending order
    @param[out] eigVecs Return eigenvectors for eigenvalues
    @param[in] upperTriangle Specifies if packed vector for matrix A is given
                             as upper (row-wise) or lower (column-wise) triangle.
    @note Wrapper for LAPACK routine dspev.
    @author esquivel 08/2013
 */
BIASMathAlgo_EXPORT
int Packed_symmetric_eigenvalue_solve(long int N,
                                      const BIAS::Vector<double> &A,
                                      BIAS::Vector<double>& eigVals,
                                      std::vector<BIAS::Vector<double> >& eigVecs,
                                      const bool upperTriangle = true);

/** solve unsymmetric eigenvalue problems (for a rectangular matrix)
    @status untested, grest
*/
BIASMathAlgo_EXPORT
int
Eigenvalue_solve(const BIAS::Matrix<double> &A,
                 BIAS::Vector<double> &wr, BIAS::Vector<double> &wi);

/** solve general eigenproblem for a general quadratix (n x n)  matrix.
    computes eigenvalues and eigenvectors of A.
    (Jan Woetzel), 03/2002
**/
BIASMathAlgo_EXPORT
int
Eigenproblem_quadratic_matrix(const BIAS::Matrix<double> &A, 
                  BIAS::Vector<double> &ret_EigenValuesReal, 
                  BIAS::Vector<double> &ret_EigenValuesImag,
                  BIAS::Matrix<double> &eigenVectors);

/** solve general eigenproblem.
    computes singular value decomposition of a rectangular m x n matrix A 
    calls extern liblapack routine. 
    better use SVD.hh
    @params ret_S vector with singular values of A with s(i) >= s(i+1) 
    @params ret_U m x m orthogonal/unitary matrix U 
    @params ret_VT n x n orthogonal/unitary matrx transpose(V)
    Jan Woetzel 03/2002, changed for Matrix by Daniel Grest, July 2002
    woelk: added jobu and jobvt
    jobu:  'A' - calculate U
           'N' - do not calculate U (much faster)
    not tested
    see man dgesvd*/
BIASMathAlgo_EXPORT
int
General_singular_value_decomposition(char jobu, char jobvt,
                                     const BIAS::Matrix<double> &A, 
                                     BIAS::Vector<double> &ret_S, 
                                     BIAS::Matrix<double> &ret_U, 
                                     BIAS::Matrix<double> &ret_VT);

BIASMathAlgo_EXPORT
inline int
General_singular_value_decomposition(const BIAS::Matrix<double> &A, 
                                     BIAS::Vector<double> &ret_S, 
                                     BIAS::Matrix<double> &ret_U, 
                                     BIAS::Matrix<double> &ret_VT)
{ return General_singular_value_decomposition('A', 'A', A, ret_S, ret_U, 
                                              ret_VT); }

BIASMathAlgo_EXPORT
inline int
General_singular_value_decomposition(const BIAS::Matrix<double> &A, 
                                     BIAS::Vector<double> &ret_S, 
                                     BIAS::Matrix<double> &ret_VT)
{ BIAS::Matrix<double> ret_U;
  return General_singular_value_decomposition('N', 'A', A, ret_S, ret_U, 
                                              ret_VT); }


/** converts the argument Fortran_Matrix to a Matrix and returns it.
    it mainly tranposes the data array because Fortran-Matrices 
    different row/col-order
    @author Jan Woetzel
    @status bug fixed by Daniel Grest (04/2003)  */
BIASMathAlgo_EXPORT
BIAS::Matrix<double> 
Fortran_Matrix_to_Matrix(const TNT::Fortran_Matrix<double>&FMat);

/** converts the argument Fortran_Matrix to a Matrix and returns it.
    it mainly tranposes the data array because Fortran-Matrices 
    different row/col-order, same as above but now copying necessary
    @author Daniel Grest
    @status beta April 2003*/
BIASMathAlgo_EXPORT
void
Fortran_Matrix_to_Matrix(const TNT::Fortran_Matrix<double> &FMat, 
                         BIAS::Matrix<double>  &res);


BIASMathAlgo_EXPORT
BIAS::Matrix<double> generalised_eigenvalue_matrix_solve(
        BIAS::Matrix<double> &A_, /* m x m matrix */
        BIAS::Matrix<double> &B_ /* m x m matrx */ );

//@}
 

#endif
// BIASLAPACK_H