LinearAlgebra.h

Go to the documentation of this file.

#ifndef SimTK_LINEAR_ALGEBRA_H_
#define SimTK_LINEAR_ALGEBRA_H_
 
/* -------------------------------------------------------------------------- *
 *                        Simbody(tm): SimTKmath                              *
 * -------------------------------------------------------------------------- *
 * This is part of the SimTK 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/home/simbody.  *
 *                                                                            *
 * Portions copyright (c) 2006-12 Stanford University and the Authors.        *
 * Authors: Jack Middleton                                                    *
 * Contributors: Michael Sherman                                              *
 *                                                                            *
 * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
 * not use this file except in compliance with the License. You may obtain a  *
 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
 *                                                                            *
 * Unless required by applicable law or agreed to in writing, software        *
 * distributed under the License is distributed on an "AS IS" BASIS,          *
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
 * See the License for the specific language governing permissions and        *
 * limitations under the License.                                             *
 * -------------------------------------------------------------------------- */
 
#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( FactorLU&& ) noexcept;
    FactorLU& operator=(const FactorLU& rhs);
    FactorLU& operator=(FactorLU&& rhs) noexcept;
 
    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:
    ClonePtr<FactorLURepBase> rep;
 
}; // class FactorLU
 
class FactorLLTRepBase;
class SimTK_SIMMATH_EXPORT FactorLLT : public Factor {
   public:
    ~FactorLLT();
 
    FactorLLT();
    FactorLLT(const FactorLLT& c);
    FactorLLT(FactorLLT&& c) noexcept;
    FactorLLT& operator=(const FactorLLT& rhs);
    FactorLLT& operator=(FactorLLT&& rhs) noexcept;
 
    template <class ELT>
    FactorLLT(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 inverse(Matrix_<ELT>& m) const;
 
    template <class ELT>
    void getL(Matrix_<ELT>& l) const;
 
    protected:
    ClonePtr<FactorLLTRepBase> rep;
 
};  // class FactorLLT
 
class FactorQTZRepBase;
class SimTK_SIMMATH_EXPORT FactorQTZ: public Factor {
    public:
 
    ~FactorQTZ();
 
    FactorQTZ();
    FactorQTZ( const FactorQTZ& c );
    FactorQTZ( FactorQTZ&& c ) noexcept;
    FactorQTZ& operator=(const FactorQTZ& rhs);
    FactorQTZ& operator=(FactorQTZ&& rhs) noexcept;
 
    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:
    ClonePtr<FactorQTZRepBase> rep;
}; // class FactorQTZ
 
class EigenRepBase;
class SimTK_SIMMATH_EXPORT Eigen {
    public:
 
    ~Eigen();
 
    Eigen();
    Eigen( const Eigen& c );
    Eigen( Eigen&& c ) noexcept;
    Eigen& operator=(const Eigen& rhs);
    Eigen& operator=(Eigen&& rhs) noexcept;
 
    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:
    ClonePtr<EigenRepBase> rep;
 
}; // class Eigen
 
class FactorSVDRepBase;
class SimTK_SIMMATH_EXPORT FactorSVD: public Factor {
    public:
 
    ~FactorSVD();
    FactorSVD();
    FactorSVD( const FactorSVD& c );
    FactorSVD( FactorSVD&& ) noexcept;
    FactorSVD& operator=(const FactorSVD& rhs);
    FactorSVD& operator=(FactorSVD&&) noexcept;
 
    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:
    ClonePtr<FactorSVDRepBase> rep;
 
}; // class FactorSVD
 
} // namespace SimTK 
 
#endif //SimTK_LINEAR_ALGEBRA_H_