65 tsh_ = 0.5 * std::sqrt(m_ + n_ + 1.) * s_.
front() * eps_;
125 for (
uint32 j = 0; j < n_; j++)
130 for (
uint32 i = 0; i < m_; i++)
139 for (
uint32 j = 0; j < n_; j++)
163 static double SIGN(
double inA,
double inB)
noexcept
194 NdArray<double>
rv1(n_, 1);
196 for (i = 0; i < n_; ++i)
204 for (k = i; k < m_; ++k)
206 scale += std::abs(u_(k, i));
211 for (k = i; k < m_; ++k)
214 ss += u_(k, i) * u_(k, i);
218 g = -SIGN(std::sqrt(
ss), f);
222 for (j =
l - 1;
j < n_; ++
j)
224 for (
ss = 0., k = i; k < m_; ++k)
226 ss += u_(k, i) * u_(k, j);
231 for (k = i; k < m_; ++k)
233 u_(k, j) +=
f * u_(k, i);
237 for (k = i; k < m_; ++k)
247 if (i + 1 <= m_ && i + 1 != n_)
249 for (k =
l - 1; k < n_; ++k)
251 scale += std::abs(u_(i, k));
256 for (k =
l - 1; k < n_; ++k)
259 ss += u_(i, k) * u_(i, k);
263 g = -SIGN(std::sqrt(
ss), f);
265 u_(i,
l - 1) =
f -
g;
267 for (k =
l - 1; k < n_; ++k)
269 rv1[k] = u_(i, k) /
h;
272 for (j =
l - 1;
j < m_; ++
j)
274 for (
ss = 0., k =
l - 1; k < n_; ++k)
276 ss += u_(j, k) * u_(i, k);
279 for (k =
l - 1; k < n_; ++k)
281 u_(j, k) +=
ss *
rv1[k];
285 for (k =
l - 1; k < n_; ++k)
292 anorm = std::max(
anorm, (std::abs(s_[i]) + std::abs(
rv1[i])));
295 for (i = n_ - 1; i !=
static_cast<uint32>(-1); --i)
301 for (j =
l;
j < n_; ++
j)
303 v_(j, i) = (u_(i, j) / u_(i,
l)) /
g;
306 for (j =
l;
j < n_; ++
j)
308 for (
ss = 0., k =
l; k < n_; ++k)
310 ss += u_(i, k) * v_(k, j);
313 for (k =
l; k < n_; ++k)
315 v_(k, j) +=
ss * v_(k, i);
320 for (j =
l;
j < n_; ++
j)
322 v_(i, j) = v_(j, i) = 0.;
331 for (i = std::min(m_, n_) - 1; i !=
static_cast<uint32>(-1); --i)
336 for (j =
l;
j < n_; ++
j)
345 for (j =
l;
j < n_; ++
j)
347 for (
ss = 0., k =
l; k < m_; ++k)
349 ss += u_(k, i) * u_(k, j);
352 f = (
ss / u_(i, i)) *
g;
354 for (k = i; k < m_; ++k)
356 u_(k, j) +=
f * u_(k, i);
360 for (j = i;
j < m_; ++
j)
367 for (j = i;
j < m_; ++
j)
376 for (k = n_ - 1; k !=
static_cast<uint32>(-1); --k)
381 for (
l = k;
l !=
static_cast<uint32>(-1); --
l)
384 if (
l == 0 || std::abs(
rv1[
l]) <= eps_ *
anorm)
390 if (std::abs(s_[
nm]) <= eps_ *
anorm)
400 for (i =
l; i < k + 1; ++i)
405 if (std::abs(f) <= eps_ *
anorm)
417 for (j = 0;
j < m_; ++
j)
421 u_(j,
nm) = y *
c + z *
ss;
422 u_(j, i) = z *
c - y *
ss;
433 for (j = 0;
j < n_; ++
j)
435 v_(j, k) = -v_(j, k);
451 f = ((y - z) * (y + z) + (
g -
h) * (
g +
h)) / (2. *
h * y);
453 f = ((x - z) * (x + z) +
h * ((y / (
f + SIGN(
g, f))) -
h)) / x;
456 for (j =
l;
j <=
nm;
j++)
472 for (
jj = 0;
jj < n_; ++
jj)
476 v_(
jj, j) = x *
c + z *
ss;
477 v_(
jj, i) = z *
c - x *
ss;
493 for (
jj = 0;
jj < m_; ++
jj)
497 u_(
jj, j) = y *
c + z *
ss;
498 u_(
jj, i) = z *
c - y *
ss;
521 NdArray<double>
su(m_, 1);
522 NdArray<double>
sv(n_, 1);
533 for (i =
inc; i < n_; ++i)
536 for (k = 0; k < m_; ++k)
541 for (k = 0; k < n_; ++k)
547 while (s_[j -
inc] <
sw)
551 for (k = 0; k < m_; ++k)
553 u_(k, j) = u_(k, j -
inc);
556 for (k = 0; k < n_; ++k)
558 v_(k, j) = v_(k, j -
inc);
571 for (k = 0; k < m_; ++k)
576 for (k = 0; k < n_; ++k)
583 for (k = 0; k < n_; ++k)
587 for (i = 0; i < m_; i++)
595 for (j = 0;
j < n_; ++
j)
603 if (
ss > (m_ + n_) / 2)
605 for (i = 0; i < m_; ++i)
607 u_(i, k) = -u_(i, k);
610 for (j = 0;
j < n_; ++
j)
612 v_(j, k) = -v_(j, k);
627 static double pythag(
double inA,
double inB)
noexcept
629 const double absa = std::abs(
inA);
630 const double absb = std::abs(
inB);
640 NdArray<double> u_{};
641 NdArray<double> v_{};
642 NdArray<double> s_{};
#define THROW_INVALID_ARGUMENT_ERROR(msg)
Definition Error.hpp:37
Holds 1D and 2D arrays, the main work horse of the NumCpp library.
Definition NdArrayCore.hpp:139
const_reference front() const noexcept
Definition NdArrayCore.hpp:2936
Definition SVDClass.hpp:47
NdArray< double > solve(const NdArray< double > &inInput, double inThresh=-1.)
Definition SVDClass.hpp:110
const NdArray< double > & s() noexcept
Definition SVDClass.hpp:96
const NdArray< double > & v() noexcept
Definition SVDClass.hpp:85
SVD(const NdArray< double > &inMatrix)
Definition SVDClass.hpp:55
const NdArray< double > & u() noexcept
Definition SVDClass.hpp:74
constexpr auto j
Definition Core/Constants.hpp:42
constexpr double c
speed of light
Definition Core/Constants.hpp:36
Definition cholesky.hpp:41
dtype f(GeneratorType &generator, dtype inDofN, dtype inDofD)
Definition f.hpp:56
bool essentiallyEqual(dtype inValue1, dtype inValue2) noexcept
Definition essentiallyEqual.hpp:49
constexpr dtype sqr(dtype inValue) noexcept
Definition sqr.hpp:42
NdArray< dtype > arange(dtype inStart, dtype inStop, dtype inStep=1)
Definition arange.hpp:59
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