pinv.hpp
Go to the documentation of this file.00001 /*************************************************************************** 00002 pinv.h - description 00003 ------------------------- 00004 begin : June 2006 00005 copyright : (C) 2006 Erwin Aertbelien 00006 email : firstname.lastname@mech.kuleuven.ac.be 00007 00008 History (only major changes)( AUTHOR-Description ) : 00009 00010 *************************************************************************** 00011 * This library is free software; you can redistribute it and/or * 00012 * modify it under the terms of the GNU Lesser General Public * 00013 * License as published by the Free Software Foundation; either * 00014 * version 2.1 of the License, or (at your option) any later version. * 00015 * * 00016 * This library is distributed in the hope that it will be useful, * 00017 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 00018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * 00019 * Lesser General Public License for more details. * 00020 * * 00021 * You should have received a copy of the GNU Lesser General Public * 00022 * License along with this library; if not, write to the Free Software * 00023 * Foundation, Inc., 59 Temple Place, * 00024 * Suite 330, Boston, MA 02111-1307 USA * 00025 * * 00026 ***************************************************************************/ 00027 00028 00029 #ifndef PINV_H 00030 #define PINV_H 00031 00032 #include <math.h> 00033 /* 00034 inline void MatrixOutput(std::ostream& os,double* a,int stridea,int rows,int cols) { 00035 int i,j; 00036 for (i=0;i<rows;++i) { 00037 for (j=0;j<cols;++j) { 00038 os << a[stridea*i+j]; 00039 if (j!=cols-1) os << "\t"; 00040 } 00041 if (i!=rows-1) os << "\n"; 00042 } 00043 } 00044 */ 00045 00049 class PseudoInverse { 00050 int maxsize; 00051 double * tmp; 00052 double* y2; 00053 public: 00054 00058 const double* A; 00059 00063 int strideA; 00064 00068 const double* Lx; 00069 00073 const double* Ly; 00074 00078 int m; 00079 00083 int n; 00084 00088 bool prepared; 00094 double * U; 00098 int strideU; 00104 double * V; 00108 int strideV; 00112 double * S; 00118 int *ndx; 00119 00120 public: 00127 PseudoInverse(int _maxsize); 00128 00135 int smallestSV(double eps); 00136 00137 /*int pinv(double* A,int strideA,int m,int n,double *Ai,int strideAi, 00138 int nrOfIts,double eps); 00139 ** 00140 * Weighted PseudoInverse 00141 * - t = A*q 00142 * - The norm of t is given by t'*Wt*t 00143 * - The norm of q is given by q'*Wq*q 00144 * - Wt is POSITIVE SEMI-DEFINITE diagonal matrix 00145 * - Wq is POSITIVE DEFINITE diagonal matrix 00146 * - Lt is a diagonal matrix of size m x m, such that Wt = Lt'*Lt 00147 * - Lq is a diagonal matrix of size n x n, such that Wq = Lq'*Lq 00148 * 00149 * \todo 00150 * check Lq for POS DEF. COND 00151 * 00152 void wpinv( 00153 double* A,int strideA,int m,int n,double *Awi,int strideAwi, 00154 double* Lt,double* Lq,int nrOfIts,double eps); 00155 */ 00156 ~PseudoInverse(); 00157 00178 int prepare( const double* _A,int _strideA,int _m,int _n, 00179 const double *_Ly,const double* _Lx, 00180 int nrOfIts); 00181 00182 00197 void inverseWithNullSpace( 00198 const double* y, 00199 double* x, 00200 const double* xd, 00201 double eps); 00214 void inverse( 00215 const double* y, 00216 double* x, 00217 double eps); 00218 00229 double derivsv( 00230 const double* Adot, 00231 int strideAdot, 00232 int svnr); 00233 00234 }; 00235 00236 00237 00238 #endif 00239