59 template<
int M,
int N,
class Precision = DefaultPrecision,
bool WANT_U = 1,
bool WANT_V = 1>
85 for(
int i=0; i<N; ++i)
86 inv_diag[i] = (
vDiagonal[i] * condition > dMax) ?
87 static_cast<Precision
>(1)/
vDiagonal[i] : 0;
94 template <
int Rows2,
int Cols2,
typename P2,
typename B2>
107 template <
int Size,
typename P2,
typename B2>
146 template<
int M,
int N,
class Precision,
bool WANT_U,
bool WANT_V>
147 template<
class Precision2,
class Base>
158 template<
int M,
int N,
class Precision,
bool WANT_U,
bool WANT_V>
165 Precision scale = 0.0;
167 for(
int i=0; i<N; ++i)
170 vOffDiagonal[i] = scale * g;
176 for(
int k=i; k<M; ++k)
177 scale +=
abs(mU[k][i]);
180 for(
int k=i; k<M; ++k)
183 s += mU[k][i] * mU[k][i];
185 Precision f = mU[i][i];
186 g = -(f>=0?sqrt(s):-sqrt(s));
187 Precision h = f * g - s;
191 for(
int j=l; j<N; ++j)
194 for(
int k=i; k<M; ++k)
195 s += mU[k][i] * mU[k][j];
197 for(
int k=i; k<M; ++k)
198 mU[k][j] += f * mU[k][i];
201 for(
int k=i; k<M; ++k)
205 vDiagonal[i] = scale * g;
209 if(!(i >= M || i == (N-1)))
211 for(
int k=l; k<N; ++k)
212 scale +=
abs(mU[i][k]);
215 for(
int k=l; k<N; k++)
218 s += mU[i][k] * mU[i][k];
220 Precision f = mU[i][l];
221 g = -(f>=0?sqrt(s):-sqrt(s));
222 Precision h = f * g - s;
224 for(
int k=l; k<N; ++k)
225 vOffDiagonal[k] = mU[i][k] / h;
228 for(
int j=l; j<M; ++j)
231 for(
int k=l; k<N; ++k)
232 s += mU[j][k] * mU[i][k];
233 for(
int k=l; k<N; ++k)
234 mU[j][k] += s * vOffDiagonal[k];
237 for(
int k=l; k<N; ++k)
241 anorm = max(anorm,
abs(vDiagonal[i]) +
abs(vOffDiagonal[i]));
247 template<
int M,
int N,
class Precision,
bool WANT_U,
bool WANT_V>
253 mV[N-1][N-1] =
static_cast<Precision
>(1);
254 Precision g = vOffDiagonal[N-1];
257 for(
int i=N-2; i>=0; --i)
262 for(
int j=l; j<N; ++j)
263 mV[j][i] = (mU[i][j] / mU[i][l]) / g;
264 for(
int j=l; j<N; ++j)
267 for(
int k=l; k<N; ++k)
268 s += mU[i][k] * mV[k][j];
269 for(
int k=l; k<N; ++k)
270 mV[k][j] += s * mV[k][i];
273 for(
int j=l; j<N; ++j)
274 mV[i][j] = mV[j][i] = 0;
275 mV[i][i] =
static_cast<Precision
>(1);
280 template<
int M,
int N,
class Precision,
bool WANT_U,
bool WANT_V>
285 for(
int i=SmallDim-1; i>=0; --i)
288 Precision g = vDiagonal[i];
291 for(
int j=l; j<N; ++j)
294 for(
int j=i; j<M; ++j)
299 Precision inv_g =
static_cast<Precision
>(1) / g;
300 if(i != (SmallDim-1))
302 for(
int j=l; j<N; ++j)
305 for(
int k=l; k<M; ++k)
306 s += mU[k][i] * mU[k][j];
307 Precision f = (s / mU[i][i]) * inv_g;
308 for(
int k=i; k<M; ++k)
309 mU[k][j] += f * mU[k][i];
312 for(
int j=i; j<M; ++j)
315 mU[i][i] +=
static_cast<Precision
>(1);
319 template<
int M,
int N,
class Precision,
bool WANT_U,
bool WANT_V>
323 for(
int k=N-1; k>=0; --k)
327 bool bConverged_Or_Error =
false;
329 bConverged_Or_Error = Diagonalize_SubLoop(k, z);
330 while(!bConverged_Or_Error);
339 for(
int j=0; j<N; ++j)
340 mV[j][k] = -mV[j][k];
346 template<
int M,
int N,
class Precision,
bool WANT_U,
bool WANT_V>
352 for(
int l=k; l>=0; --l)
355 const bool rv1_test = ((
abs(vOffDiagonal[l]) + anorm) == anorm);
356 const bool w_test = ((
abs(vDiagonal[l1]) + anorm) == anorm);
357 if(!rv1_test && w_test)
361 for(
int i=l; i<=k; ++i)
363 Precision f = s * vOffDiagonal[i];
364 vOffDiagonal[i] *= c;
365 if((
abs(f) + anorm) == anorm)
367 Precision g = vDiagonal[i];
368 Precision h = sqrt(f * f + g * g);
373 for(
int j=0; j<M; ++j)
375 Precision y = mU[j][l1];
376 Precision z = mU[j][i];
377 mU[j][l1] = y*c + z*s;
378 mU[j][i] = -y*s + z*c;
382 if(rv1_test || w_test)
388 if(nIterations == 30)
394 Precision x = vDiagonal[l];
395 Precision y = vDiagonal[k1];
396 Precision g = vOffDiagonal[k1];
397 Precision h = vOffDiagonal[k];
398 Precision f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2.0*h*y);
400 Precision signed_g = (f>=0)?g:-g;
401 f = ((x-z)*(x+z) + h*(y/(f + signed_g) - h)) / x;
406 for(
int i1 = l; i1<=k1; ++i1)
414 vOffDiagonal[i1] = z;
422 for(
int j=0; j<N; ++j)
424 Precision xx = mV[j][i1];
425 Precision zz = mV[j][i];
426 mV[j][i1] = xx*c + zz*s;
427 mV[j][i] = -xx*s + zz*c;
439 for(
int j=0; j<M; ++j)
441 Precision yy = mU[j][i1];
442 Precision zz = mU[j][i];
443 mU[j][i1] = yy*c + zz*s;
444 mU[j][i] = -yy*s + zz*c;
460 template<
int M,
int N,
class Precision,
bool WANT_U,
bool WANT_V>
464 Precision d = vDiagonal[0];
465 for(
int i=1; i<N; ++i) d = max(d, vDiagonal[i]);
469 template<
int M,
int N,
class Precision,
bool WANT_U,
bool WANT_V>
473 Precision d = vDiagonal[0];
474 for(
int i=1; i<N; ++i) d = min(d, vDiagonal[i]);
478 template<
int M,
int N,
class Precision,
bool WANT_U,
bool WANT_V>
483 Precision d = vDiagonal[0];
484 for(
int i=1; i<N; ++i)
493 template<
int M,
int N,
class Precision,
bool WANT_U,
bool WANT_V>
496 std::vector<std::pair<Precision, unsigned int> > vSort;
498 for(
unsigned int i=0; i<N; ++i)
499 vSort.push_back(std::make_pair(-vDiagonal[i], i));
500 std::sort(vSort.begin(), vSort.end());
501 for(
unsigned int i=0; i<N; ++i)
502 vDiagonal[i] = -vSort[i].first;
506 for(
unsigned int i=0; i<N; ++i)
507 mU.T()[i] = mU_copy.T()[vSort[i].second];
512 for(
unsigned int i=0; i<N; ++i)
513 mV.T()[i] = mV_copy.T()[vSort[i].second];
Vector< BigDim, Precision > vOffDiagonal
Matrix< N, M, Precision > get_pinv(const Precision condition=1e9)
Get the pseudo-inverse .
static const int SmallDim
Everything lives inside this namespace.
void reorder()
Reorder the components so the singular values are in descending order.
Matrix< N, Cols2, typename Internal::MultiplyType< Precision, P2 >::type > backsub(const Matrix< Rows2, Cols2, P2, B2 > &rhs, const Precision condition=1e9)
Precision get_smallest_singular_value()
void get_inv_diag(Vector< N > &inv_diag, const Precision condition)
const Matrix< N, N, Precision > & get_V()
int get_smallest_singular_value_index()
bool Diagonalize_SubLoop(int k, Precision &z)
Matrix< M, N, Precision > mU
GR_SVD(const Matrix< M, N, Precision2, Base > &A)
const Matrix< M, N, Precision > & get_U()
Matrix< N, N, Precision > mV
Precision get_largest_singular_value()
Vector< N, typename Internal::MultiplyType< Precision, P2 >::type > backsub(const Vector< Size, P2, B2 > &rhs, const Precision condition=1e9)
const Vector< N, Precision > & get_diagonal()
Vector< N, Precision > vDiagonal
Vector< Internal::Sizer< S1, S2 >::size, typename Internal::MultiplyType< P1, P2 >::type > diagmult(const Vector< S1, P1, B1 > &v1, const Vector< S2, P2, B2 > &v2)