LinearAlgebra.h

Go to the documentation of this file.

#ifndef SimTK_LINEAR_ALGEBRA_H_
#define SimTK_LINEAR_ALGEBRA_H_
/* -------------------------------------------------------------------------- *
 *                      SimTK Core: SimTK Simmath(tm)                         *
 * -------------------------------------------------------------------------- *
 * This is part of the SimTK Core biosimulation toolkit originating from      *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
 * Portions copyright (c) 2006-10 Stanford University and the Authors.        *
 * Authors: Jack Middleton                                                    *
 * Contributors: Michael Sherman                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */


#include "SimTKcommon.h"
#include "simmath/internal/common.h"


namespace SimTK {

//  default for reciprocal of the condition number
// TODO: sherm 080128 I changed this from 0.01 to a more reasonable
// value but it is still wrong because the default should depend
// on the matrix size, something like max(m,n)*eps^(7/8) where
// eps is machine precision for float or double as appropriate.
static const double DefaultRecpCondition = 1e-12;

class SimTK_SIMMATH_EXPORT Factor {
public:

  Factor() {}
  template <class ELT> Factor( Matrix_<ELT> m );
  template <class ELT> void solve( const Vector_<ELT>& b, Vector_<ELT>& x ) const;
  template <class ELT> void solve( const Matrix_<ELT>& b, Matrix_<ELT>& x ) const;
  
}; // class Factor

class FactorLURepBase;

class SimTK_SIMMATH_EXPORT FactorLU: public Factor {
    public:

    ~FactorLU();

    FactorLU();
    FactorLU( const FactorLU& c );
    FactorLU& operator=(const FactorLU& rhs);

    template <class ELT> FactorLU( const Matrix_<ELT>& m );
    template <class ELT> void factor( const Matrix_<ELT>& m );
    template <class ELT> void solve( const Vector_<ELT>& b, Vector_<ELT>& x ) const;
    template <class ELT> void solve( const Matrix_<ELT>& b, Matrix_<ELT>& x ) const;

    template <class ELT> void getL( Matrix_<ELT>& l ) const;
    template <class ELT> void getU( Matrix_<ELT>& u ) const;
    template < class ELT > void inverse(  Matrix_<ELT>& m ) const;

    bool isSingular() const;
    int getSingularIndex() const;


    protected:
    class FactorLURepBase *rep;

}; // class FactorLU


class FactorQTZRepBase;
class SimTK_SIMMATH_EXPORT FactorQTZ: public Factor {
    public:

    ~FactorQTZ();

    FactorQTZ();
    FactorQTZ( const FactorQTZ& c );
    FactorQTZ& operator=(const FactorQTZ& rhs);
    template <typename ELT> FactorQTZ( const Matrix_<ELT>& m);
    template <typename ELT> FactorQTZ( const Matrix_<ELT>& m, double rcond );
    template <typename ELT> FactorQTZ( const Matrix_<ELT>& m, float rcond );
    template <typename ELT> void factor( const Matrix_<ELT>& m);
    template <typename ELT> void factor( const Matrix_<ELT>& m, float rcond );
    template <typename ELT> void factor( const Matrix_<ELT>& m, double rcond );
    template <typename ELT> void solve( const Vector_<ELT>& b, Vector_<ELT>& x ) const;
    template <typename ELT> void solve( const Matrix_<ELT>& b, Matrix_<ELT>& x ) const;

    template < class ELT > void inverse(  Matrix_<ELT>& m ) const;

    int getRank() const;
    double getRCondEstimate() const;
//    void setRank(int rank); TBD

    protected:
    class FactorQTZRepBase *rep;
}; // class FactorQTZ
class SimTK_SIMMATH_EXPORT Eigen {
    public:

    ~Eigen();

    Eigen();
    Eigen( const Eigen& c );
    Eigen& operator=(const Eigen& rhs);

    template <class ELT> Eigen( const Matrix_<ELT>& m );
    template <class ELT> void factor( const Matrix_<ELT>& m );
    template <class VAL, class VEC> void getAllEigenValuesAndVectors( Vector_<VAL>& values, Matrix_<VEC>& vectors);
    template <class T> void getAllEigenValues( Vector_<T>& values);

    template <class VAL, class VEC> void getFewEigenValuesAndVectors( Vector_<VAL>& values, Matrix_<VEC>& vectors, int ilow, int ihi);
    template <class T> void getFewEigenVectors( Matrix_<T>& vectors, int ilow, int ihi );
    template <class T> void getFewEigenValues( Vector_<T>& values, int ilow, int ihi );

    template <class VAL, class VEC> void getFewEigenValuesAndVectors( Vector_<VAL>& values, Matrix_<VEC>& vectors, typename CNT<VAL>::TReal rlow, typename CNT<VAL>::TReal rhi);
    template <class T> void getFewEigenVectors( Matrix_<T>& vectors, typename CNT<T>::TReal rlow, typename CNT<T>::TReal rhi );
    template <class T> void getFewEigenValues( Vector_<T>& values, typename CNT<T>::TReal rlow, typename CNT<T>::TReal rhi );

     
    protected:
    class EigenRepBase *rep;

}; // class Eigen
class SimTK_SIMMATH_EXPORT FactorSVD: public Factor {
    public:

    ~FactorSVD();
    FactorSVD();
    FactorSVD( const FactorSVD& c );
    FactorSVD& operator=(const FactorSVD& rhs);

    template < class ELT > FactorSVD( const Matrix_<ELT>& m );
    template < class ELT > FactorSVD( const Matrix_<ELT>& m, float rcond );
    template < class ELT > FactorSVD( const Matrix_<ELT>& m, double rcond );
    template < class ELT > void factor( const Matrix_<ELT>& m );
    template < class ELT > void factor( const Matrix_<ELT>& m, float rcond );
    template < class ELT > void factor( const Matrix_<ELT>& m, double rcond );

    template < class T > void getSingularValuesAndVectors( Vector_<typename CNT<T>::TReal>& values, 
                              Matrix_<T>& leftVectors,  Matrix_<T>& rightVectors );
    template < class T > void getSingularValues( Vector_<T>& values);

    int getRank();
    template < class ELT > void inverse(  Matrix_<ELT>& m );
    template <class ELT> void solve( const Vector_<ELT>& b, Vector_<ELT>& x );
    template <class ELT> void solve( const Matrix_<ELT>& b, Matrix_<ELT>& x );

    protected:
    class FactorSVDRepBase *rep;

}; // class FactorSVD

} // namespace SimTK 

#endif //SimTK_LINEAR_ALGEBRA_H_