27 #ifndef CDPL_MATH_SVDECOMPOSITION_HPP
28 #define CDPL_MATH_SVDECOMPOSITION_HPP
69 template <
typename A,
typename W,
typename V>
73 typename V::ValueType>::Type ValueType;
76 typename V::SizeType>::Type SizeType;
78 SizeType
m = a().getSize1();
79 SizeType n = a().getSize2();
81 CDPL_MATH_CHECK(SizeType(w().getSize()) >= n && SizeType(v().getSize1()) >= n && SizeType(v().getSize2()) >= n,
86 ValueType g, scale, anorm,
s;
88 g = scale = anorm = ValueType(0);
90 for (SizeType i = 0; i < n; i++) {
93 g =
s = scale = ValueType(0);
96 for (SizeType k = i; k <
m; k++)
99 if (scale != ValueType(0)) {
100 for (SizeType k = i; k <
m; k++) {
102 s += a()(k, i) * a()(k, i);
105 ValueType f = a()(i, i);
109 ValueType h = f * g -
s;
113 for (SizeType j = l - 1; j < n; j++) {
116 for (SizeType k = i; k <
m; k++)
117 s += a()(k, i) * a()(k, j);
121 for (SizeType k = i; k <
m; k++)
122 a()(k, j) += f * a()(k, i);
125 for (SizeType k = i; k <
m; k++)
131 g =
s = scale = ValueType(0);
133 if (i + 1 <=
m && i + 1 != n) {
134 for (SizeType k = l - 1; k < n; k++)
137 if (scale != ValueType(0)) {
138 for (SizeType k = l - 1; k < n; k++) {
140 s += a()(i, k) * a()(i, k);
143 ValueType f = a()(i, l - 1);
147 ValueType h = f * g -
s;
149 a()(i, l - 1) = f - g;
151 for (SizeType k = l - 1; k < n; k++)
152 rv1(k) = a()(i, k) / h;
154 for (SizeType j = l - 1; j <
m; j++) {
157 for (SizeType k = l - 1; k < n; k++)
158 s += a()(j, k) * a()(i, k);
160 for (SizeType k = l - 1; k < n; k++)
161 a()(j, k) +=
s * rv1(k);
164 for (SizeType k = l - 1; k < n; k++)
172 for (SizeType ic = n; ic > 0; ic--) {
176 if (g != ValueType(0)) {
177 for (SizeType j = l; j < n; j++)
178 v()(j, i) = (a()(i, j) / a()(i, l)) / g;
180 for (SizeType j = l; j < n; j++) {
183 for (SizeType k = l; k < n; k++)
184 s += a()(i, k) * v()(k, j);
186 for (SizeType k = l; k < n; k++)
187 v()(k, j) +=
s * v()(k, i);
191 for (SizeType j = l; j < n; j++)
192 v()(i, j) = v()(j, i) = ValueType(0);
195 v()(i, i) = ValueType(1);
200 for (SizeType ic = std::min(
m, n); ic > 0; ic--) {
206 for (SizeType j = l; j < n; j++)
207 a()(i, j) = ValueType(0);
209 if (g != ValueType(0)) {
210 g = ValueType(1) / g;
212 for (SizeType j = l; j < n; j++) {
215 for (SizeType k = l; k <
m; k++)
216 s += a()(k, i) * a()(k, j);
218 ValueType f = (
s / a()(i, i)) * g;
220 for (SizeType k = i; k <
m; k++)
221 a()(k, j) += f * a()(k, i);
224 for (SizeType j = i; j <
m; j++)
228 for (SizeType j = i; j <
m; j++)
229 a()(j, i) = ValueType(0);
234 for (SizeType kc = n; kc > 0; kc--) {
238 for (std::size_t its = 0; max_iter == 0 || its < max_iter; its++) {
241 for (SizeType lc = k + 1; lc > 0; lc--) {
255 ValueType c = ValueType(0);
259 for (SizeType i = l; i < k + 1; i++) {
260 ValueType f =
s * rv1(i);
269 ValueType h =
pythag(f, g);
272 h = ValueType(1) / h;
276 for (SizeType j = 0; j <
m; j++) {
277 ValueType y = a()(j, nm);
278 ValueType z = a()(j, i);
280 a()(j, nm) = y * c + z *
s;
281 a()(j, i) = z * c - y *
s;
286 ValueType z = w()(k);
289 if (z < ValueType(0)) {
292 for (SizeType j = 0; j < n; j++)
293 v()(j, k) = -v()(j, k);
298 if (max_iter > 0 && its >= max_iter - 1)
301 ValueType x = w()(l);
302 ValueType c = ValueType(1);
306 ValueType y = w()(nm);
310 ValueType h = rv1(k);
311 ValueType f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
313 g =
pythag(f, ValueType(1));
314 f = ((x - z) * (x + z) + h * ((y / (f +
sign(g, f))) - h)) / x;
317 for (SizeType j = l; j <= nm; j++) {
333 for (SizeType jj = 0; jj < n; jj++) {
336 v()(jj, j) = x * c + z *
s;
337 v()(jj, i) = z * c - x *
s;
344 z = ValueType(1) / z;
352 for (SizeType jj = 0; jj <
m; jj++) {
355 a()(jj, j) = y * c + z *
s;
356 a()(jj, i) = z * c - y *
s;
360 rv1(l) = ValueType(0);
381 for (SizeType i = inc; i < n; i++) {
382 ValueType sw = w()(i);
384 for (SizeType k = 0; k <
m; k++)
387 for (SizeType k = 0; k < n; k++)
392 while (w()(j - inc) < sw) {
393 w()(j) = w()(j - inc);
395 for (SizeType k = 0; k <
m; k++)
396 a()(k, j) = a()(k, j - inc);
398 for (SizeType k = 0; k < n; k++)
399 v()(k, j) = v()(k, j - inc);
409 for (SizeType k = 0; k <
m; k++)
412 for (SizeType k = 0; k < n; k++)
418 for (SizeType k = 0; k < n; k++) {
421 for (SizeType i = 0; i <
m; i++)
422 if (a()(i, k) <
typename A::ValueType(0))
425 for (SizeType j = 0; j < n; j++)
426 if (v()(j, k) <
typename V::ValueType(0))
429 if (
s > (
m + n) / 2) {
430 for (SizeType i = 0; i <
m; i++)
431 a()(i, k) = -a()(i, k);
433 for (SizeType j = 0; j < n; j++)
434 v()(j, k) = -v()(j, k);
465 template <
typename U,
typename W,
typename V,
typename B,
typename X>
470 typename W::ValueType>::Type,
471 typename V::ValueType>::Type,
472 typename B::ValueType>::Type,
473 typename X::ValueType>::Type ValueType;
477 SizeType
m = u().getSize1();
478 SizeType n = u().getSize2();
480 CDPL_MATH_CHECK(w().getSize() == n && v().getSize1() == n && v().getSize2() == n &&
481 b().getSize() ==
m && x().getSize() == n,
485 ValueType thresh = 0.5 * std::sqrt(
m + n + 1.0) * w()(0) * std::numeric_limits<ValueType>::epsilon();
487 for (SizeType i = 0; i < n; i++) {
488 ValueType
s = ValueType(0);
496 x().assign(
prod(v, tmp));
524 template <
typename U,
typename W,
typename V,
typename B,
typename X>
530 typename W::ValueType>::Type,
531 typename V::ValueType>::Type,
532 typename B::ValueType>::Type,
533 typename X::ValueType>::Type ValueType;
536 typename W::SizeType>::Type,
537 typename B::SizeType>::Type,
538 typename X::SizeType>::Type SizeType;
540 SizeType
m = u().getSize1();
541 SizeType n = u().getSize2();
542 SizeType
p = b().getSize2();
544 CDPL_MATH_CHECK(SizeType(w().getSize()) == n && SizeType(v().getSize1()) == n && SizeType(v().getSize2()) == n &&
545 SizeType(b().getSize1()) ==
m && SizeType(x().getSize1()) == n && SizeType(x().getSize2()) ==
p,
549 ValueType thresh = 0.5 * std::sqrt(
m + n + 1.0) * w()(0) * std::numeric_limits<ValueType>::epsilon();
551 for (SizeType i = 0; i <
p; i++) {
554 for (SizeType j = 0; j < n; j++) {
555 ValueType
s = ValueType(0);
569 #endif // CDPL_MATH_SVDECOMPOSITION_HPP