57 m_(inMatrix.
shape().rows),
58 n_(inMatrix.
shape().cols),
62 eps_(std::numeric_limits<double>::epsilon())
115 if (inInput.
size() != m_)
124 tsh_ = (inThresh >= 0. ? inThresh : 0.5 *
sqrt(m_ + n_ + 1.) * s_.
front() * eps_);
131 for (
uint32 i = 0; i < m_; i++)
133 ss += u_(i,
j) * inInput[i];
143 for (
uint32 jj = 0; jj < n_; jj++)
145 ss += v_(
j, jj) * tmp[jj];
164 static double SIGN(
double inA,
double inB) noexcept
166 return inB >= 0 ? (inA >= 0 ? inA : -inA) : (inA >= 0 ? -inA : inA);
197 for (i = 0; i < n_; ++i)
201 g = ss = scale = 0.0;
205 for (k = i; k < m_; ++k)
212 for (k = i; k < m_; ++k)
215 ss += u_(k, i) * u_(k, i);
223 for (
j = l - 1;
j < n_; ++
j)
225 for (ss = 0.0, k = i; k < m_; ++k)
227 ss += u_(k, i) * u_(k,
j);
232 for (k = i; k < m_; ++k)
234 u_(k,
j) +=
f * u_(k, i);
238 for (k = i; k < m_; ++k)
246 g = ss = scale = 0.0;
248 if (i + 1 <= m_ && i + 1 != n_)
250 for (k = l - 1; k < n_; ++k)
257 for (k = l - 1; k < n_; ++k)
260 ss += u_(i, k) * u_(i, k);
266 u_(i, l - 1) =
f - g;
268 for (k = l - 1; k < n_; ++k)
270 rv1[k] = u_(i, k) / h;
273 for (
j = l - 1;
j < m_; ++
j)
275 for (ss = 0.0, k = l - 1; k < n_; ++k)
277 ss += u_(
j, k) * u_(i, k);
280 for (k = l - 1; k < n_; ++k)
282 u_(
j, k) += ss * rv1[k];
286 for (k = l - 1; k < n_; ++k)
296 for (i = n_ - 1; i !=
static_cast<uint32>(-1); --i)
302 for (
j = l;
j < n_; ++
j)
304 v_(
j, i) = (u_(i,
j) / u_(i, l)) / g;
307 for (
j = l;
j < n_; ++
j)
309 for (ss = 0.0, k = l; k < n_; ++k)
311 ss += u_(i, k) * v_(k,
j);
314 for (k = l; k < n_; ++k)
316 v_(k,
j) += ss * v_(k, i);
321 for (
j = l;
j < n_; ++
j)
323 v_(i,
j) = v_(
j, i) = 0.0;
337 for (
j = l;
j < n_; ++
j)
346 for (
j = l;
j < n_; ++
j)
348 for (ss = 0.0, k = l; k < m_; ++k)
350 ss += u_(k, i) * u_(k,
j);
353 f = (ss / u_(i, i)) * g;
355 for (k = i; k < m_; ++k)
357 u_(k,
j) +=
f * u_(k, i);
361 for (
j = i;
j < m_; ++
j)
368 for (
j = i;
j < m_; ++
j)
377 for (k = n_ - 1; k !=
static_cast<uint32>(-1); --k)
379 for (its = 0; its < 30; ++its)
382 for (l = k; l !=
static_cast<uint32>(-1); --l)
385 if (l == 0 ||
std::abs(rv1[l]) <= eps_ * anorm)
391 if (
std::abs(s_[nm]) <= eps_ * anorm)
401 for (i = l; i < k + 1; ++i)
418 for (
j = 0;
j < m_; ++
j)
422 u_(
j, nm) = y *
c + z * ss;
423 u_(
j, i) = z *
c - y * ss;
434 for (
j = 0;
j < n_; ++
j)
436 v_(
j, k) = -v_(
j, k);
452 f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
454 f = ((x - z) * (x + z) + h * ((y / (
f + SIGN(g,
f))) - h)) / x;
457 for (
j = l;
j <= nm;
j++)
473 for (jj = 0; jj < n_; ++jj)
477 v_(jj,
j) = x *
c + z * ss;
478 v_(jj, i) = z *
c - x * ss;
494 for (jj = 0; jj < m_; ++jj)
498 u_(jj,
j) = y *
c + z * ss;
499 u_(jj, i) = z *
c - y * ss;
522 NdArray<double> su(m_, 1);
523 NdArray<double> sv(n_, 1);
534 for (i = inc; i < n_; ++i)
537 for (k = 0; k < m_; ++k)
542 for (k = 0; k < n_; ++k)
548 while (s_[
j - inc] < sw)
552 for (k = 0; k < m_; ++k)
554 u_(k,
j) = u_(k,
j - inc);
557 for (k = 0; k < n_; ++k)
559 v_(k,
j) = v_(k,
j - inc);
572 for (k = 0; k < m_; ++k)
577 for (k = 0; k < n_; ++k)
584 for (k = 0; k < n_; ++k)
588 for (i = 0; i < m_; i++)
596 for (
j = 0;
j < n_; ++
j)
604 if (ss > (m_ + n_) / 2)
606 for (i = 0; i < m_; ++i)
608 u_(i, k) = -u_(i, k);
611 for (
j = 0;
j < n_; ++
j)
613 v_(
j, k) = -v_(
j, k);
628 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:4289
const_reference front() const noexcept
Definition: NdArrayCore.hpp:2764
Definition: SVDClass.hpp:48
const NdArray< double > & u() noexcept
Definition: SVDClass.hpp:75
const NdArray< double > & s() noexcept
Definition: SVDClass.hpp:97
const NdArray< double > & v() noexcept
Definition: SVDClass.hpp:86
NdArray< double > solve(const NdArray< double > &inInput, double inThresh=-1.0)
Definition: SVDClass.hpp:111
SVD(const NdArray< double > &inMatrix)
Definition: SVDClass.hpp:56
constexpr auto j
Definition: Constants.hpp:45
constexpr double c
speed of light
Definition: Constants.hpp:40
dtype f(GeneratorType &generator, dtype inDofN, dtype inDofD)
Definition: f.hpp:58
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:44
auto abs(dtype inValue) noexcept
Definition: abs.hpp:49
NdArray< dtype > min(const NdArray< dtype > &inArray, Axis inAxis=Axis::NONE)
Definition: min.hpp:44
auto sqrt(dtype inValue) noexcept
Definition: sqrt.hpp:48
Shape shape(const NdArray< dtype > &inArray) noexcept
Definition: Functions/Shape.hpp:42
std::uint32_t uint32
Definition: Types.hpp:40