37 #ifndef VIGRA_SINGULAR_VALUE_DECOMPOSITION_HXX
38 #define VIGRA_SINGULAR_VALUE_DECOMPOSITION_HXX
41 #include "array_vector.hxx"
73 template <
class T,
class C1,
class C2,
class C3,
class C4>
83 vigra_precondition(rows >= cols,
84 "singularValueDecomposition(): Input matrix A must be rectangular with rowCount >= columnCount.");
86 "singularValueDecomposition(): Output S must be column vector with rowCount == columnCount(A).");
88 "singularValueDecomposition(): Output matrix U must have the same dimensions as input matrix A.");
90 "singularValueDecomposition(): Output matrix V must be square with n = columnCount(A).");
111 for (k = 0; k <
std::max(nct,nrt); ++k)
119 for (i = k; i < m; ++i)
121 s(k) =
hypot(s(k), a(i, k));
129 for (i = k; i < m; ++i)
137 for (j = k+1; j < n; ++j)
139 if ((k < nct) && (s(k) != 0.0))
143 for (i = k; i < m; ++i)
145 t += a(i, k)*a(i, j);
148 for (i = k; i < m; ++i)
150 a(i, j) += t*a(i, k);
164 for (i = k; i < m; ++i)
175 for (i = k+1; i < n; ++i)
177 e[k] =
hypot(e[k],e[i]);
185 for (i = k+1; i < n; ++i)
192 if ((k+1 < m) & (e[k] != 0.0))
195 for (i = k+1; i < m; ++i)
199 for (j = k+1; j < n; ++j)
201 for (i = k+1; i < m; ++i)
203 work[i] += e[j]*a(i, j);
206 for (j = k+1; j < n; ++j)
208 Real t(-e[j]/e[k+1]);
209 for (i = k+1; i < m; ++i)
211 a(i, j) += t*work[i];
217 for (i = k+1; i < n; ++i)
229 s(nct) = a(nct, nct);
237 e[nrt] = a(nrt, p-1);
242 for (j = nct; j < nu; ++j)
244 for (i = 0; i < m; ++i)
250 for (k = nct-1; k >= 0; --k)
254 for (j = k+1; j < nu; ++j)
257 for (i = k; i < m; ++i)
259 t += U(i, k)*U(i, j);
262 for (i = k; i < m; ++i)
264 U(i, j) += t*U(i, k);
267 for (i = k; i < m; ++i )
271 U(k, k) = 1.0 + U(k, k);
272 for (i = 0; i < k-1; ++i)
279 for (i = 0; i < m; ++i)
288 for (k = n-1; k >= 0; --k)
290 if ((k < nrt) & (e[k] != 0.0))
292 for (j = k+1; j < nu; ++j)
295 for (i = k+1; i < n; ++i)
297 t += V(i, k)*V(i, j);
300 for (i = k+1; i < n; ++i)
302 V(i, j) += t*V(i, k);
306 for (i = 0; i < n; ++i)
317 Real eps = NumericTraits<Real>::epsilon()*2.0;
335 for (k = p-2; k >= -1; --k)
341 if (
abs(e[k]) <= eps*(
abs(s(k)) +
abs(s(k+1))))
354 for (ks = p-1; ks >= k; --ks)
360 Real t( (ks != p ?
abs(e[ks]) : 0.) +
361 (ks != k+1 ?
abs(e[ks-1]) : 0.));
362 if (
abs(s(ks)) <= eps*t)
392 for (j = p-2; j >= k; --j)
394 Real t(
hypot(s(j),f));
403 for (i = 0; i < n; ++i)
405 t = cs*V(i, j) + sn*V(i, p-1);
406 V(i, p-1) = -sn*V(i, j) + cs*V(i, p-1);
416 for (j = k; j < p; ++j)
418 Real t(
hypot(s(j),f));
424 for (i = 0; i < m; ++i)
426 t = cs*U(i, j) + sn*U(i, k-1);
427 U(i, k-1) = -sn*U(i, j) + cs*U(i, k-1);
436 Real scale =
std::max(std::max(std::max(std::max(
439 Real sp = s(p-1)/scale;
440 Real spm1 = s(p-2)/scale;
441 Real epm1 = e[p-2]/scale;
442 Real sk = s(k)/scale;
443 Real ek = e[k]/scale;
444 Real b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0;
445 Real c = (sp*epm1)*(sp*epm1);
447 if ((b != 0.0) || (c != 0.0))
454 shift = c/(b + shift);
456 Real f = (sk + sp)*(sk - sp) + shift;
460 for (j = k; j < p-1; ++j)
469 f = cs*s(j) + sn*e[j];
470 e[j] = cs*e[j] - sn*s(j);
473 for (i = 0; i < n; ++i)
475 t = cs*V(i, j) + sn*V(i, j+1);
476 V(i, j+1) = -sn*V(i, j) + cs*V(i, j+1);
483 f = cs*e[j] + sn*s(j+1);
484 s(j+1) = -sn*e[j] + cs*s(j+1);
489 for (i = 0; i < m; ++i)
491 t = cs*U(i, j) + sn*U(i, j+1);
492 U(i, j+1) = -sn*U(i, j) + cs*U(i, j+1);
506 s(k) = (s(k) < 0.0 ? -s(k) : 0.0);
507 for (i = 0; i <= pp; ++i)
526 for (i = 0; i < n; ++i)
528 t = V(i, k+1); V(i, k+1) = V(i, k); V(i, k) = t;
533 for (i = 0; i < m; ++i)
535 t = U(i, k+1); U(i, k+1) = U(i, k); U(i, k) = t;
545 vigra_fail(
"vigra::svd(): Internal error.");
549 unsigned int rank = 0;
566 #endif // VIGRA_SINGULAR_VALUE_DECOMPOSITION_HXX