58 m_(inMatrix.
shape().rows),
59 n_(inMatrix.
shape().cols),
63 eps_(std::numeric_limits<double>::epsilon())
120 if (inInput.
size() != m_)
129 tsh_ = (inThresh >= 0. ? inThresh : 0.5 *
sqrt(m_ + n_ + 1.) * s_.
front() * eps_);
136 for (
uint32 i = 0; i < m_; i++)
138 ss += u_(i,
j) * inInput[i];
148 for (
uint32 jj = 0; jj < n_; jj++)
150 ss += v_(
j, jj) * tmp[jj];
170 static double SIGN(
double inA,
double inB) noexcept
172 return inB >= 0 ? (inA >= 0 ? inA : -inA) : (inA >= 0 ? -inA : inA);
203 for (i = 0; i < n_; ++i)
207 g = ss = scale = 0.0;
211 for (k = i; k < m_; ++k)
218 for (k = i; k < m_; ++k)
221 ss += u_(k, i) * u_(k, i);
229 for (
j = l - 1;
j < n_; ++
j)
231 for (ss = 0.0, k = i; k < m_; ++k)
233 ss += u_(k, i) * u_(k,
j);
238 for (k = i; k < m_; ++k)
240 u_(k,
j) +=
f * u_(k, i);
244 for (k = i; k < m_; ++k)
252 g = ss = scale = 0.0;
254 if (i + 1 <= m_ && i + 1 != n_)
256 for (k = l - 1; k < n_; ++k)
263 for (k = l - 1; k < n_; ++k)
266 ss += u_(i, k) * u_(i, k);
272 u_(i, l - 1) =
f - g;
274 for (k = l - 1; k < n_; ++k)
276 rv1[k] = u_(i, k) / h;
279 for (
j = l - 1;
j < m_; ++
j)
281 for (ss = 0.0, k = l - 1; k < n_; ++k)
283 ss += u_(
j, k) * u_(i, k);
286 for (k = l - 1; k < n_; ++k)
288 u_(
j, k) += ss * rv1[k];
292 for (k = l - 1; k < n_; ++k)
302 for (i = n_ - 1; i !=
static_cast<uint32>(-1); --i)
308 for (
j = l;
j < n_; ++
j)
310 v_(
j, i) = (u_(i,
j) / u_(i, l)) / g;
313 for (
j = l;
j < n_; ++
j)
315 for (ss = 0.0, k = l; k < n_; ++k)
317 ss += u_(i, k) * v_(k,
j);
320 for (k = l; k < n_; ++k)
322 v_(k,
j) += ss * v_(k, i);
327 for (
j = l;
j < n_; ++
j)
329 v_(i,
j) = v_(
j, i) = 0.0;
343 for (
j = l;
j < n_; ++
j)
352 for (
j = l;
j < n_; ++
j)
354 for (ss = 0.0, k = l; k < m_; ++k)
356 ss += u_(k, i) * u_(k,
j);
359 f = (ss / u_(i, i)) * g;
361 for (k = i; k < m_; ++k)
363 u_(k,
j) +=
f * u_(k, i);
367 for (
j = i;
j < m_; ++
j)
375 for (
j = i;
j < m_; ++
j)
384 for (k = n_ - 1; k !=
static_cast<uint32>(-1); --k)
386 for (its = 0; its < 30; ++its)
389 for (l = k; l !=
static_cast<uint32>(-1); --l)
392 if (l == 0 ||
std::abs(rv1[l]) <= eps_ * anorm)
398 if (
std::abs(s_[nm]) <= eps_ * anorm)
408 for (i = l; i < k + 1; ++i)
425 for (
j = 0;
j < m_; ++
j)
429 u_(
j, nm) = y *
c + z * ss;
430 u_(
j, i) = z *
c - y * ss;
441 for (
j = 0;
j < n_; ++
j)
443 v_(
j, k) = -v_(
j, k);
459 f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
461 f = ((x - z) * (x + z) + h * ((y / (
f + SIGN(g,
f))) - h)) / x;
464 for (
j = l;
j <= nm;
j++)
480 for (jj = 0; jj < n_; ++jj)
484 v_(jj,
j) = x *
c + z * ss;
485 v_(jj, i) = z *
c - x * ss;
501 for (jj = 0; jj < m_; ++jj)
505 u_(jj,
j) = y *
c + z * ss;
506 u_(jj, i) = z *
c - y * ss;
529 NdArray<double> su(m_, 1);
530 NdArray<double> sv(n_, 1);
541 for (i = inc; i < n_; ++i)
544 for (k = 0; k < m_; ++k)
549 for (k = 0; k < n_; ++k)
555 while (s_[
j - inc] < sw)
559 for (k = 0; k < m_; ++k)
561 u_(k,
j) = u_(k,
j - inc);
564 for (k = 0; k < n_; ++k)
566 v_(k,
j) = v_(k,
j - inc);
579 for (k = 0; k < m_; ++k)
584 for (k = 0; k < n_; ++k)
592 for (k = 0; k < n_; ++k)
596 for (i = 0; i < m_; i++)
604 for (
j = 0;
j < n_; ++
j)
612 if (ss > (m_ + n_) / 2)
614 for (i = 0; i < m_; ++i)
616 u_(i, k) = -u_(i, k);
619 for (
j = 0;
j < n_; ++
j)
621 v_(
j, k) = -v_(
j, k);
637 static double pythag(
double inA,
double inB) noexcept
#define THROW_INVALID_ARGUMENT_ERROR(msg)
Definition: Error.hpp:36
size_type size() const noexcept
Definition: NdArrayCore.hpp:4497
const_reference front() const noexcept
Definition: NdArrayCore.hpp:2933
Definition: SVDClass.hpp:48
const NdArray< double > & u() noexcept
Definition: SVDClass.hpp:77
const NdArray< double > & s() noexcept
Definition: SVDClass.hpp:101
const NdArray< double > & v() noexcept
Definition: SVDClass.hpp:89
NdArray< double > solve(const NdArray< double > &inInput, double inThresh=-1.0)
Definition: SVDClass.hpp:116
SVD(const NdArray< double > &inMatrix)
Definition: SVDClass.hpp:57
constexpr auto j
Definition: Constants.hpp:45
constexpr double c
speed of light
Definition: Constants.hpp:40
dtype f(dtype inDofN, dtype inDofD)
Definition: f.hpp:56
constexpr dtype sqr(dtype inValue) noexcept
Definition: sqr.hpp:44
Definition: Coordinate.hpp:45
NdArray< dtype > max(const NdArray< dtype > &inArray, Axis inAxis=Axis::NONE)
Definition: max.hpp:45
auto abs(dtype inValue) noexcept
Definition: abs.hpp:51
NdArray< dtype > min(const NdArray< dtype > &inArray, Axis inAxis=Axis::NONE)
Definition: min.hpp:45
auto sqrt(dtype inValue) noexcept
Definition: sqrt.hpp:50
Shape shape(const NdArray< dtype > &inArray) noexcept
Definition: Functions/Shape.hpp:44
std::uint32_t uint32
Definition: Types.hpp:40