Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Namespace Members | Class Members | File Members | Related Pages

robot/control/kinematics/SVDFullSpaceSolver.cpp

Go to the documentation of this file.
00001 /****************************************************************************
00002   Copyright (C)2002 David Jung <opensim@pobox.com>
00003 
00004   This program/file is free software; you can redistribute it and/or modify
00005   it under the terms of the GNU General Public License as published by
00006   the Free Software Foundation; either version 2 of the License, or
00007   (at your option) any later version.
00008 
00009   This program is distributed in the hope that it will be useful,
00010   but WITHOUT ANY WARRANTY; without even the implied warranty of
00011   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012   GNU General Public License for more details. (http://www.gnu.org)
00013 
00014   You should have received a copy of the GNU General Public License
00015   along with this program; if not, write to the Free Software
00016   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00017 
00018   $Id: SVDFullSpaceSolver.cpp 1080 2004-07-28 19:51:26Z jungd $
00019   $Revision: 1.4 $
00020   $Date: 2004-07-28 15:51:26 -0400 (Wed, 28 Jul 2004) $
00021   $Author: jungd $
00022 
00023 ****************************************************************************/
00024 
00025 #include <robot/control/kinematics/SVDFullSpaceSolver>
00026 
00027 using robot::control::kinematics::SVDFullSpaceSolver;
00028 
00029 #include <base/SVD>
00030 #include <robot/control/kinematics/solution_error>
00031 
00032 //
00033 // Solution coded from Lonnie Love's Matlab code
00034 //
00035 
00036 
00037 using base::IVector;
00038 using base::zeroIVector;
00039 using base::equals;
00040 using base::SVD;
00041 
00042 const Real small = 1.0e-04; // (=def.SMALL)
00043 
00044 
00045 SVDFullSpaceSolver::SVDFullSpaceSolver()
00046  : stopOnIllCondition(false)
00047 {
00048 }
00049 
00050 
00051 
00052 const static Real maxCondition = 10000;
00053 
00054 base::Matrix SVDFullSpaceSolver::solve(const Matrix& A_in, const Vector& b_in,
00055                                        array<Int>& dependentRowsEliminated)
00056 {
00057   const Int N = A_in.rows();
00058   const Int M = A_in.cols();
00059 
00060   const Vector& b(b_in);
00061   Assert( b.size() == N);
00062 
00063   bool nullSpace = b.equals(zeroVector(N));
00064 
00065   // decompose A using SVD
00066   SVD svd(A_in);
00067 
00068   if (!nullSpace && (svd.condition() >= maxCondition)) {
00069     if (stopOnIllCondition)
00070       throw solution_error(Exception("FSP unable to complete solution: A is ill-conditioned (has large SVD condition number="
00071                                             +base::realToString(svd.condition())+")"));
00072     else
00073       Logln("Warning: A is ill-conditioned; singular values are:" << svd.diag() << " condition=" << svd.condition());
00074   }
00075 
00076   const Matrix& U( svd.U() );
00077   const Matrix  S( svd.S() );
00078   const Matrix& V( svd.V() );
00079   Matrix Uinv( transpose(U) );
00080 
00081 
00082   // invert S
00083   Matrix Sinv( zeroMatrix(M,N) );
00084   for(Int i=0; i< Math::minimum(M, N); i++) {
00085     if (Math::abs(S(i,i)) > small)
00086       Sinv(i,i) = 1/S(i,i);
00087     else {
00088       if (stopOnIllCondition)
00089         throw solution_error(Exception(String("FSP unable to complete solution: a singular value of A is near 0 ")
00090                                         +"(S("+base::intToString(i)+","+base::intToString(i)+")="+base::realToString(S(i,i))+")"));
00091       else {
00092 //        Logln("Warning: small singular value, setting 1/S(" << base::intToString(i) << "," << base::intToString(i) << ") to 0");
00093         Sinv(i,i) = 0;
00094       }
00095     }
00096   }
00097 
00098 
00099   Matrix g;
00100 
00101   // special case is b=0 (null-space)
00102   if (nullSpace) {
00103 
00104     // just extract the null-space vectors from V (the last M-N vectors)
00105     g.resize(M, Math::maximum(SInt(M-N),0) );
00106     if (g.cols() > 0)
00107       for(Int i=0; i<M-N; i++)
00108         matrixColumn(g,i) = matrixColumn(V,N+i);
00109 
00110   }
00111   else {
00112 
00113     g.resize(M, Math::maximum(SInt(M-N+1),0) );
00114 
00115     Vector Utb(Uinv * b);
00116 
00117     if (g.cols() > 0) {
00118       matrixColumn(g,0) = V * Sinv * Utb; // particular solution
00119 
00120       // now add homogeneous solution
00121       for(Int i=0; i<M-N; i++)
00122         matrixColumn(g,i+1)=matrixColumn(V,N+i) + matrixColumn(g,0);
00123     }
00124   }
00125 
00126   // return the indices of rows that were eliminated due to dependency (none for this method)
00127   dependentRowsEliminated.clear();
00128 
00129   return g;
00130 }
00131 
00132 
00133 
00134 
00135 
00136 
00137 
00138 
00139 

Generated on Thu Jul 29 15:56:34 2004 for OpenSim by doxygen 1.3.6