37#ifndef VIGRA_SINGULAR_VALUE_DECOMPOSITION_HXX
38#define VIGRA_SINGULAR_VALUE_DECOMPOSITION_HXX
41#include "array_vector.hxx"
73template <
class T,
class C1,
class C2,
class C3,
class C4>
83 "singularValueDecomposition(): Input matrix A must be rectangular with rowCount >= columnCount.");
85 "singularValueDecomposition(): Output S must be column vector with rowCount == columnCount(A).");
87 "singularValueDecomposition(): Output matrix U must have the same dimensions as input matrix A.");
89 "singularValueDecomposition(): Output matrix V must be square with n = columnCount(A).");
118 for (
i =
k;
i <
m; ++
i)
128 for (
i =
k;
i <
m; ++
i)
136 for (
j =
k+1;
j < n; ++
j)
138 if ((
k <
nct) && (s(
k) != 0.0))
142 for (
i =
k;
i <
m; ++
i)
144 t += a(
i,
k)*a(
i,
j);
147 for (
i =
k;
i <
m; ++
i)
149 a(
i,
j) += t*a(
i,
k);
163 for (
i =
k;
i <
m; ++
i)
174 for (
i =
k+1;
i < n; ++
i)
184 for (
i =
k+1;
i < n; ++
i)
191 if ((
k+1 <
m) & (e[
k] != 0.0))
194 for (
i =
k+1;
i <
m; ++
i)
198 for (
j =
k+1;
j < n; ++
j)
200 for (
i =
k+1;
i <
m; ++
i)
205 for (
j =
k+1;
j < n; ++
j)
207 Real t(-e[
j]/e[
k+1]);
208 for (
i =
k+1;
i <
m; ++
i)
216 for (
i =
k+1;
i < n; ++
i)
243 for (
i = 0;
i <
m; ++
i)
249 for (
k =
nct-1;
k >= 0; --
k)
253 for (
j =
k+1;
j <
nu; ++
j)
256 for (
i =
k;
i <
m; ++
i)
261 for (
i =
k;
i <
m; ++
i)
266 for (
i =
k;
i <
m; ++
i )
271 for (
i = 0;
i <
k-1; ++
i)
278 for (
i = 0;
i <
m; ++
i)
287 for (
k = n-1;
k >= 0; --
k)
289 if ((
k <
nrt) & (e[
k] != 0.0))
291 for (
j =
k+1;
j <
nu; ++
j)
294 for (
i =
k+1;
i < n; ++
i)
296 t += V(
i,
k)*V(
i,
j);
299 for (
i =
k+1;
i < n; ++
i)
301 V(
i,
j) += t*V(
i,
k);
305 for (
i = 0;
i < n; ++
i)
316 Real
eps = NumericTraits<Real>::epsilon()*2.0;
334 for (
k = p-2;
k >= -1; --
k)
340 if (abs(e[
k]) <=
eps*(abs(s(
k)) + abs(s(
k+1))))
359 Real t( (
ks != p ? abs(e[
ks]) : 0.) +
360 (
ks !=
k+1 ? abs(e[
ks-1]) : 0.));
361 if (abs(s(
ks)) <=
eps*t)
391 for (
j = p-2;
j >=
k; --
j)
402 for (
i = 0;
i < n; ++
i)
405 V(
i, p-1) = -
sn*V(
i,
j) +
cs*V(
i, p-1);
415 for (
j =
k;
j < p; ++
j)
423 for (
i = 0;
i <
m; ++
i)
435 Real scale = std::max(std::max(std::max(std::max(
436 abs(s(p-1)),abs(s(p-2))),abs(e[p-2])),
437 abs(s(
k))),abs(e[
k]));
438 Real
sp = s(p-1)/scale;
439 Real
spm1 = s(p-2)/scale;
440 Real
epm1 = e[p-2]/scale;
441 Real
sk = s(
k)/scale;
442 Real
ek = e[
k]/scale;
446 if ((b != 0.0) || (c != 0.0))
448 shift = VIGRA_CSTD::sqrt(b*b + c);
453 shift = c/(b + shift);
455 Real f = (
sk +
sp)*(
sk -
sp) + shift;
459 for (
j =
k;
j < p-1; ++
j)
472 for (
i = 0;
i < n; ++
i)
488 for (
i = 0;
i <
m; ++
i)
505 s(
k) = (s(
k) < 0.0 ? -s(
k) : 0.0);
506 for (
i = 0;
i <=
pp; ++
i)
525 for (
i = 0;
i < n; ++
i)
527 t = V(
i,
k+1); V(
i,
k+1) = V(
i,
k); V(
i,
k) = t;
532 for (
i = 0;
i <
m; ++
i)
544 vigra_fail(
"vigra::svd(): Internal error.");
547 Real
tol = std::max(
m,n)*s(0)*
eps;
548 unsigned int rank = 0;
549 for (
i = 0;
i < n; ++
i)
Class for a single RGB value.
Definition rgbvalue.hxx:128
void init(Iterator i, Iterator end)
Definition tinyvector.hxx:708
unsigned int singularValueDecomposition(MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > &U, MultiArrayView< 2, T, C3 > &S, MultiArrayView< 2, T, C4 > &V)
Definition singular_value_decomposition.hxx:75
MultiArrayIndex columnCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:684
MultiArrayIndex rowCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:671
std::ptrdiff_t MultiArrayIndex
Definition multi_fwd.hxx:60
FixedPoint16< IntBits, OverflowHandling > hypot(FixedPoint16< IntBits, OverflowHandling > v1, FixedPoint16< IntBits, OverflowHandling > v2)
Length of hypotenuse.
Definition fixedpoint.hxx:1640