1 #ifndef SimTK_SIMMATH_FACTORSVD_REP_H_ 2 #define SimTK_SIMMATH_FACTORSVD_REP_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) 2007-12 Stanford University and the Authors. * 13 * Authors: Jack Middleton * 14 * Contributors: * 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 #include "SimTKmath.h" 28 #include "WorkSpace.h" 29 30 namespace SimTK { 31 32 class FactorSVDRepBase { 33 public: 34 ~FactorSVDRepBase()35 virtual ~FactorSVDRepBase(){}; 36 clone()37 virtual FactorSVDRepBase* clone() const { return 0; }; 38 39 getSingularValues(Vector_<float> & values)40 virtual void getSingularValues( Vector_<float>& values ){ 41 checkIfFactored( "getSingularValues" ); 42 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","getSingularValues", 43 "getSingularValues( float) called with types that are inconsistent with the original matrix \n"); 44 } getSingularValues(Vector_<double> & values)45 virtual void getSingularValues( Vector_<double>& values ){ 46 checkIfFactored( "getSingularValues" ); 47 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","getSingularValues", 48 "getSingularValues( double) called with types that are inconsistent with the original matrix \n"); 49 } getSingularValuesAndVectors(Vector_<float> & values,Matrix_<float> & leftVectors,Matrix_<float> & rightVectors)50 virtual void getSingularValuesAndVectors( Vector_<float>& values, Matrix_<float>& leftVectors, Matrix_<float>& rightVectors ){ 51 checkIfFactored( "getSingularValuesAndVectors" ); 52 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","getSingularValuesAndVectors", 53 "getSingularValuesAndVectors( float, float, float) called with types that are inconsistent with the original matrix \n"); 54 } getSingularValuesAndVectors(Vector_<double> & values,Matrix_<double> & leftVectors,Matrix_<double> & rightVectors)55 virtual void getSingularValuesAndVectors( Vector_<double>& values, Matrix_<double>& leftVectors, Matrix_<double>& rightVectors ){ 56 checkIfFactored( "getSingularValuesAndVectors" ); 57 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","getSingularValuesAndVectors", 58 "getSingularValuesAndVectors( double, double, double) called with types that are inconsistent with the original matrix \n"); 59 } 60 getSingularValuesAndVectors(Vector_<float> & values,Matrix_<std::complex<float>> & leftVectors,Matrix_<std::complex<float>> & rightVectors)61 virtual void getSingularValuesAndVectors( Vector_<float>& values, Matrix_<std::complex<float> >& leftVectors, Matrix_<std::complex<float> >& rightVectors ){ 62 checkIfFactored( "getSingularValuesAndVectors" ); 63 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","getSingularValuesAndVectors", 64 "getSingularValuesAndVectors( float, std::complex<float> , std::complex<float> ) called with types that are inconsistent with the original matrix \n"); 65 } 66 getSingularValuesAndVectors(Vector_<double> & values,Matrix_<std::complex<double>> & leftVectors,Matrix_<std::complex<double>> & rightVectors)67 virtual void getSingularValuesAndVectors( Vector_<double>& values, Matrix_<std::complex<double> >& leftVectors, Matrix_<std::complex<double> >& rightVectors ){ 68 checkIfFactored( "getSingularValuesAndVectors" ); 69 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","getSingularValuesAndVectors", 70 "getSingularValuesAndVectors( double, std::complex<double> , std::complex<double> ) called with types that are inconsistent with the original matrix \n"); 71 } solve(const Vector_<float> & b,Vector_<float> & x)72 virtual void solve( const Vector_<float>& b, Vector_<float>& x ) { 73 checkIfFactored("solve"); 74 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","solve", 75 "solve called with rhs of type <float> which does not match type of original linear system \n"); 76 } solve(const Vector_<double> & b,Vector_<double> & x)77 virtual void solve( const Vector_<double>& b, Vector_<double>& x ){ 78 checkIfFactored("solve"); 79 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","solve", 80 "solve called with rhs of type <double> which does not match type of original linear system \n"); 81 } solve(const Vector_<std::complex<float>> & b,Vector_<std::complex<float>> & x)82 virtual void solve( const Vector_<std::complex<float> >& b, Vector_<std::complex<float> >& x ) { 83 checkIfFactored("solve"); 84 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","solve", 85 "solve called with rhs of type complex<float> which does not match type of original linear system \n"); 86 } solve(const Vector_<std::complex<double>> & b,Vector_<std::complex<double>> & x)87 virtual void solve( const Vector_<std::complex<double> >& b, Vector_<std::complex<double> >& x ) { 88 checkIfFactored("solve"); 89 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","solve", 90 "solve called with rhs of type complex<double> which does not match type of original linear system \n"); 91 } solve(const Matrix_<float> & b,Matrix_<float> & x)92 virtual void solve( const Matrix_<float>& b, Matrix_<float>& x ) { 93 checkIfFactored("solve"); 94 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","solve", 95 "solve called with rhs of type <float> which does not match type of original linear system \n"); 96 } solve(const Matrix_<double> & b,Matrix_<double> & x)97 virtual void solve( const Matrix_<double>& b, Matrix_<double>& x ) { 98 checkIfFactored("solve"); 99 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","solve", 100 "solve called with rhs of type <double> which does not match type of original linear system \n"); 101 } solve(const Matrix_<std::complex<float>> & b,Matrix_<std::complex<float>> & x)102 virtual void solve( const Matrix_<std::complex<float> >& b, Matrix_<std::complex<float> >& x ) { 103 checkIfFactored("solve"); 104 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","solve", 105 "solve called with rhs of type complex<float> which does not match type of original linear system \n"); 106 } solve(const Matrix_<std::complex<double>> & b,Matrix_<std::complex<double>> & x)107 virtual void solve ( const Matrix_<std::complex<double> >& b, Matrix_<std::complex<double> >& x ) { 108 checkIfFactored("solve"); 109 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","solve", 110 "solve called with rhs of type complex<double> which does not match type of original linear system \n"); 111 } inverse(Matrix_<double> & inverse)112 virtual void inverse( Matrix_<double>& inverse ){ 113 checkIfFactored( "inverse" ); 114 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","inverse", 115 "inverse( <double> ) called with type that is inconsistent with the original matrix \n"); 116 } inverse(Matrix_<float> & inverse)117 virtual void inverse( Matrix_<float>& inverse ){ 118 checkIfFactored( "inverse" ); 119 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","inverse", 120 "inverse( <float> ) called with type that is inconsistent with the original matrix \n"); 121 } inverse(Matrix_<std::complex<float>> & inverse)122 virtual void inverse( Matrix_<std::complex<float> >& inverse ){ 123 checkIfFactored( "inverse" ); 124 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","inverse", 125 "inverse( std::complex<float> ) called with type that is inconsistent with the original matrix \n"); 126 } inverse(Matrix_<std::complex<double>> & inverse)127 virtual void inverse( Matrix_<std::complex<double> >& inverse ){ 128 checkIfFactored( "inverse" ); 129 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD","inverse", 130 "inverse( std::complex<double> ) called with type that is inconsistent with the original matrix \n"); 131 } getRank()132 virtual int getRank() const { 133 checkIfFactored( "getRank" ); 134 return(0); 135 } 136 checkIfFactored(const char * function_name)137 void checkIfFactored(const char* function_name) const { 138 if( !isFactored ) { 139 SimTK_APIARGCHECK_ALWAYS(false,"FactorSVD",function_name, 140 "matrix has not been specified\n"); 141 } 142 return; 143 } 144 bool isFactored; 145 146 147 148 149 }; // class FactorSVDRepBase 150 151 class FactorSVDDefault : public FactorSVDRepBase { 152 public: 153 FactorSVDDefault(); 154 FactorSVDRepBase* clone() const override; 155 }; 156 157 158 template <typename T> 159 class FactorSVDRep : public FactorSVDRepBase { 160 public: 161 template <class ELT> FactorSVDRep( const Matrix_<ELT>&, typename CNT<T>::TReal ); 162 163 ~FactorSVDRep(); 164 FactorSVDRepBase* clone() const override; 165 166 typedef typename CNT<T>::TReal RType; 167 168 void getSingularValuesAndVectors( Vector_<RType>& values, Matrix_<T>& leftVectors, Matrix_<T>& rightVectors ) override; 169 void getSingularValues( Vector_<RType>& values ) override; 170 int getRank(); 171 void solve( const Vector_<T>& b, Vector_<T>& x ) override; 172 void solve( const Matrix_<T>& b, Matrix_<T>& x ) override; 173 174 175 private: 176 177 void computeSVD( bool, RType*, T*, T* ); 178 void doSolve( Matrix_<T>& b, Matrix_<T>& x ); 179 void inverse( Matrix_<T>& b ) override; 180 181 int nCol; // number of columns in original matrix 182 int nRow; // number of rows in original matrix 183 int mn; // min(m,n) 184 int maxmn; // max(m,n) 185 int rank; 186 RType rcond; // reciprocol condition number 187 RType abstol; 188 MatrixStructure structure; 189 TypedWorkSpace<T> inputMatrix; 190 TypedWorkSpace<RType> singularValues; 191 192 }; // end class FactorSVDRep 193 } // namespace SimTK 194 #endif // SimTK_SIMMATH_FACTORSVD_REP_H_ 195