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-2007 Stanford University and Jack Middleton.   *
 * Contributors:                                                              *
 *
 * 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 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 <limits.h>
#include "SimTKcommon.h"
#include "SimTKmath.h"
#include "SimTKcommon/internal/BigMatrix.h"
#include "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_