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& 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_