1 #ifndef SimTK_LINEAR_ALGEBRA_H_ 2 #define SimTK_LINEAR_ALGEBRA_H_ 3 4 /* -------------------------------------------------------------------------- * 5 * Simbody(tm): SimTKmath * 6 * -------------------------------------------------------------------------- * 7 * This is part of the SimTK biosimulation toolkit originating from * 8 * Simbios, the NIH National Center for Physics-Based Simulation of * 9 * Biological Structures at Stanford, funded under the NIH Roadmap for * 10 * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. * 11 * * 12 * Portions copyright (c) 2006-12 Stanford University and the Authors. * 13 * Authors: Jack Middleton * 14 * Contributors: Michael Sherman * 15 * * 16 * Licensed under the Apache License, Version 2.0 (the "License"); you may * 17 * not use this file except in compliance with the License. You may obtain a * 18 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. * 19 * * 20 * Unless required by applicable law or agreed to in writing, software * 21 * distributed under the License is distributed on an "AS IS" BASIS, * 22 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * 23 * See the License for the specific language governing permissions and * 24 * limitations under the License. * 25 * -------------------------------------------------------------------------- */ 26 27 /** @file 28 * This is the header file that user code should include to pick up the 29 * SimTK Simmath linear algebra tools. 30 */ 31 32 33 #include "SimTKcommon.h" 34 #include "simmath/internal/common.h" 35 36 37 namespace SimTK { 38 39 // default for reciprocal of the condition number 40 // TODO: sherm 080128 I changed this from 0.01 to a more reasonable 41 // value but it is still wrong because the default should depend 42 // on the matrix size, something like max(m,n)*eps^(7/8) where 43 // eps is machine precision for float or double as appropriate. 44 static const double DefaultRecpCondition = 1e-12; 45 46 /** 47 * Base class for the various matrix factorizations. 48 */ 49 class SimTK_SIMMATH_EXPORT Factor { 50 public: 51 Factor()52 Factor() {} 53 /// creates an factorization of a matrix 54 template <class ELT> Factor( Matrix_<ELT> m ); 55 /// solves a single right hand side using a factorization 56 template <class ELT> void solve( const Vector_<ELT>& b, Vector_<ELT>& x ) const; 57 /// solves multiple right hand sides using a factorization 58 template <class ELT> void solve( const Matrix_<ELT>& b, Matrix_<ELT>& x ) const; 59 60 }; // class Factor 61 62 class FactorLURepBase; 63 64 /** 65 * Class for performing LU matrix factorizations 66 */ 67 class SimTK_SIMMATH_EXPORT FactorLU: public Factor { 68 public: 69 70 ~FactorLU(); 71 72 FactorLU(); 73 FactorLU( const FactorLU& c ); 74 FactorLU& operator=(const FactorLU& rhs); 75 76 template <class ELT> FactorLU( const Matrix_<ELT>& m ); 77 /// factors a matrix 78 template <class ELT> void factor( const Matrix_<ELT>& m ); 79 /// solves a single right hand side 80 template <class ELT> void solve( const Vector_<ELT>& b, Vector_<ELT>& x ) const; 81 /// solves multiple right hand sides 82 template <class ELT> void solve( const Matrix_<ELT>& b, Matrix_<ELT>& x ) const; 83 84 /// returns the lower triangle of an LU factorization 85 template <class ELT> void getL( Matrix_<ELT>& l ) const; 86 /// returns the upper triangle of an LU factorization 87 template <class ELT> void getU( Matrix_<ELT>& u ) const; 88 /// returns the inverse of a matrix using an LU factorization 89 template < class ELT > void inverse( Matrix_<ELT>& m ) const; 90 91 /// returns true if matrix was singular 92 bool isSingular() const; 93 /// returns the first diagonal which was found to be singular 94 int getSingularIndex() const; 95 96 97 protected: 98 class FactorLURepBase *rep; 99 100 }; // class FactorLU 101 102 103 class FactorQTZRepBase; 104 /** 105 * Class to perform a QTZ (linear least squares) factorization 106 */ 107 class SimTK_SIMMATH_EXPORT FactorQTZ: public Factor { 108 public: 109 110 ~FactorQTZ(); 111 112 FactorQTZ(); 113 FactorQTZ( const FactorQTZ& c ); 114 FactorQTZ& operator=(const FactorQTZ& rhs); 115 /// do QTZ factorization of a matrix 116 template <typename ELT> FactorQTZ( const Matrix_<ELT>& m); 117 /// do QTZ factorization of a matrix for a given reciprocal condition number 118 template <typename ELT> FactorQTZ( const Matrix_<ELT>& m, double rcond ); 119 /// do QTZ factorization of a matrix for a given reciprocal condition number 120 template <typename ELT> FactorQTZ( const Matrix_<ELT>& m, float rcond ); 121 /// do QTZ factorization of a matrix 122 template <typename ELT> void factor( const Matrix_<ELT>& m); 123 /// do QTZ factorization of a matrix for a given reciprocal condition number 124 template <typename ELT> void factor( const Matrix_<ELT>& m, float rcond ); 125 /// do QTZ factorization of a matrix for a given reciprocal condition number 126 template <typename ELT> void factor( const Matrix_<ELT>& m, double rcond ); 127 /// solve for a vector x given a right hand side vector b 128 template <typename ELT> void solve( const Vector_<ELT>& b, Vector_<ELT>& x ) const; 129 /// solve for an array of vectors given multiple right hand sides 130 template <typename ELT> void solve( const Matrix_<ELT>& b, Matrix_<ELT>& x ) const; 131 132 template < class ELT > void inverse( Matrix_<ELT>& m ) const; 133 134 /// returns the rank of the matrix 135 int getRank() const; 136 /// returns the actual reciprocal condition number at this rank 137 double getRCondEstimate() const; 138 // void setRank(int rank); TBD 139 140 protected: 141 class FactorQTZRepBase *rep; 142 }; // class FactorQTZ 143 /** 144 * Class to compute Eigen values and Eigen vectors of a matrix 145 */ 146 class SimTK_SIMMATH_EXPORT Eigen { 147 public: 148 149 ~Eigen(); 150 151 Eigen(); 152 Eigen( const Eigen& c ); 153 Eigen& operator=(const Eigen& rhs); 154 155 /// create a default eigen class 156 template <class ELT> Eigen( const Matrix_<ELT>& m ); 157 /// supply matrix which eigen values will be computed for 158 template <class ELT> void factor( const Matrix_<ELT>& m ); 159 /// get all the eigen values and eigen vectors of a matrix 160 template <class VAL, class VEC> void getAllEigenValuesAndVectors( Vector_<VAL>& values, Matrix_<VEC>& vectors); 161 /// get all the eigen values of a matrix 162 template <class T> void getAllEigenValues( Vector_<T>& values); 163 164 /// get a few eigen values and eigen vectors of a symmetric matrix which are within a range of indices 165 template <class VAL, class VEC> void getFewEigenValuesAndVectors( Vector_<VAL>& values, Matrix_<VEC>& vectors, int ilow, int ihi); 166 /// get a few eigen vectors of a symmetric matrix which are within a range of indices 167 template <class T> void getFewEigenVectors( Matrix_<T>& vectors, int ilow, int ihi ); 168 /// get a few eigen values of a symmetric matrix which are within a range of indices 169 template <class T> void getFewEigenValues( Vector_<T>& values, int ilow, int ihi ); 170 171 /// get a few eigen values and eigen vectors of a symmetric matrix which are within a range of eigen values 172 template <class VAL, class VEC> void getFewEigenValuesAndVectors( Vector_<VAL>& values, Matrix_<VEC>& vectors, typename CNT<VAL>::TReal rlow, typename CNT<VAL>::TReal rhi); 173 /// get a few eigen vectors of a symmetric matrix which are within a range of eigen values 174 template <class T> void getFewEigenVectors( Matrix_<T>& vectors, typename CNT<T>::TReal rlow, typename CNT<T>::TReal rhi ); 175 /// get a few eigen values of a symmetric matrix which are within a range of eigen values 176 template <class T> void getFewEigenValues( Vector_<T>& values, typename CNT<T>::TReal rlow, typename CNT<T>::TReal rhi ); 177 178 179 protected: 180 class EigenRepBase *rep; 181 182 }; // class Eigen 183 /** 184 * Class to compute a singular value decomposition of a matrix 185 */ 186 class SimTK_SIMMATH_EXPORT FactorSVD: public Factor { 187 public: 188 189 ~FactorSVD(); 190 /// default constructor 191 FactorSVD(); 192 /// copy constructor 193 FactorSVD( const FactorSVD& c ); 194 /// copy assign 195 FactorSVD& operator=(const FactorSVD& rhs); 196 197 /// constructor 198 template < class ELT > FactorSVD( const Matrix_<ELT>& m ); 199 /// singular value decomposition of a matrix using the specified reciprocal of the condition 200 /// number rcond 201 template < class ELT > FactorSVD( const Matrix_<ELT>& m, float rcond ); 202 /// singular value decomposition of a matrix using the specified reciprocal of the condition 203 /// number rcond 204 template < class ELT > FactorSVD( const Matrix_<ELT>& m, double rcond ); 205 /// supply the matrix to do a singular value decomposition 206 template < class ELT > void factor( const Matrix_<ELT>& m ); 207 /// supply the matrix to do a singular value decomposition using the specified 208 /// reciprocal of the condition number rcond 209 template < class ELT > void factor( const Matrix_<ELT>& m, float rcond ); 210 /// supply the matrix to do a singular value decomposition using the specified reciprocal of the condition 211 /// reciprocal of the condition number rcond 212 template < class ELT > void factor( const Matrix_<ELT>& m, double rcond ); 213 214 /// get the singular values and singular vectors of the matrix 215 template < class T > void getSingularValuesAndVectors( Vector_<typename CNT<T>::TReal>& values, 216 Matrix_<T>& leftVectors, Matrix_<T>& rightVectors ); 217 /// get just the singular values of the matrix 218 template < class T > void getSingularValues( Vector_<T>& values); 219 220 /// get rank of the matrix 221 int getRank(); 222 /// get inverse of the matrix using singular value decomposition (sometimes called the pseudo inverse) 223 template < class ELT > void inverse( Matrix_<ELT>& m ); 224 /// solve for x given a right hand side vector using the singular value decomposition 225 template <class ELT> void solve( const Vector_<ELT>& b, Vector_<ELT>& x ); 226 /// solve for a set of x vectors given multiple right hand side vectors 227 /// using the singular value decomposition 228 template <class ELT> void solve( const Matrix_<ELT>& b, Matrix_<ELT>& x ); 229 230 protected: 231 class FactorSVDRepBase *rep; 232 233 }; // class FactorSVD 234 235 } // namespace SimTK 236 237 #endif //SimTK_LINEAR_ALGEBRA_H_ 238