本文整理汇总了C++中orthotope::copy方法的典型用法代码示例。如果您正苦于以下问题：C++ orthotope::copy方法的具体用法？C++ orthotope::copy怎么用？C++ orthotope::copy使用的例子？那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类orthotope的用法示例。
在下文中一共展示了orthotope::copy方法的3个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: check_QR
void check_QR(
orthotope<T> const& A
, orthotope<T> const& Q
, orthotope<T> const& R
)
{ // {{{
BOOST_ASSERT(2 == A.order());
BOOST_ASSERT(2 == Q.order());
BOOST_ASSERT(2 == R.order());
BOOST_ASSERT(A.hypercube());
BOOST_ASSERT(Q.hypercube());
BOOST_ASSERT(R.hypercube());
std::size_t const n = A.extent(0);
BOOST_ASSERT(n == Q.extent(0));
BOOST_ASSERT(n == R.extent(0));
///////////////////////////////////////////////////////////////////////////
/// Make sure Q * R equals A.
orthotope<T> QR = matrix_multiply(Q, R);
for (std::size_t l = 0; l < n; ++l)
{
for (std::size_t i = 0; i < n; ++i)
{
if (!compare_floating(A(l, i), QR(l, i), 1e-6))
std::cout << "WARNING: QR[" << l << "][" << i << "] (value "
<< QR(l, i) << ") is not equal to A[" << l << "]["
<< i << "] (value " << A(l, i) << ")\n";
}
}
///////////////////////////////////////////////////////////////////////////
/// Make sure R is an upper triangular matrix.
for (std::size_t l = 0; l < (n - 1); ++l)
{
for (std::size_t i = l + 1; i < n; ++i)
{
if (!compare_floating(0.0, R(i, l), 1e-6))
std::cout << "WARNING: R[" << i << "][" << l << "] is not 0 "
"(value is " << R(i, l) << "), R is not an upper "
"triangular matrix\n";
}
}
///////////////////////////////////////////////////////////////////////////
/// Make sure Q is orthogonal. A matrix is orthogonal if its transpose is
/// equal to its inverse:
///
/// Q^T = Q^-1
///
/// This implies that:
///
/// Q^T * Q = Q * Q^T = I
///
/// We use the above formula to verify Q's orthogonality.
orthotope<T> QT = Q.copy();
// Transpose QT.
for (std::size_t l = 0; l < (n - 1); ++l)
for (std::size_t i = l + 1; i < n; ++i)
std::swap(QT(l, i), QT(i, l));
// Compute Q^T * Q and store the result in QT.
QT = matrix_multiply(Q, QT);
for (std::size_t l = 0; l < n; ++l)
{
for (std::size_t i = 0; i < n; ++i)
{
// Diagonals should be 1.
if (l == i)
{
if (!compare_floating(1.0, QT(l, i), 1e-6))
std::cout << "WARNING: (Q^T * Q)[" << l << "][" << i << "] "
"is not 1 (value is " << QT(l, i) << "), Q is "
"not an orthogonal matrix\n";
}
// All other entries should be 0.
else
{
if (!compare_floating(0.0, QT(l, i), 1e-6))
std::cout << "WARNING: (Q^T * Q)[" << l << "][" << i << "] "
"is not 0 (value is " << QT(l, i) << "), Q is "
"not an orthogonal matrix\n";
}
}
}
} // }}}示例2: householders
void householders(
orthotope<T> const& A
)
{
BOOST_ASSERT(2 == A.order());
BOOST_ASSERT(A.hypercube());
std::size_t const n = A.extent(0);
orthotope<T> R = A.copy();
orthotope<T> Q({n, n});
for (std::size_t l = 0; l < n; ++l)
Q(l, l) = 1.0;
for (std::size_t l = 0; l < (n - 1); ++l)
{
T const sigma = compute_sigma(R, n, l);
boost::int16_t const sign = compute_sign(R(l, l));
#if defined(HPXLA_DEBUG_HOUSEHOLDERS)
std::cout << (std::string(80, '#') + "\n")
<< "ROUND " << l << "\n\n";
print(sigma, "sigma");
print(sign, "sign");
#endif
orthotope<T> w({n});
w(l) = R(l, l) + sign * sigma;
for (std::size_t i = (l + 1); i < w.extent(0); ++i)
w(i) = R(i, l);
#if defined(HPXLA_DEBUG_HOUSEHOLDERS)
print(w, "u");
#endif
T const w_norm = euclidean_norm(w);
for (std::size_t i = l; i < n; ++i)
w(i) /= w_norm;
#if defined(HPXLA_DEBUG_HOUSEHOLDERS)
print(w, "v");
#endif
orthotope<T> H = compute_H(w);
#if defined(HPXLA_DEBUG_HOUSEHOLDERS)
print(H, "H");
#endif
R = matrix_multiply(H, R);
Q = matrix_multiply(Q, H);
for (std::size_t i = l + 1; i < n; ++i)
R(i, l) = 0;
}
#if defined(HPXLA_DEBUG_HOUSEHOLDERS)
std::cout << std::string(80, '#') << "\n";
#endif
print(A, "A");
print(Q, "Q");
print(R, "R");
#if defined(HPXLA_DEBUG_HOUSEHOLDERS)
check_QR(A, Q, R);
#endif
}示例3: evs
std::vector<std::complex<T> > qr_eigenvalue(
orthotope<T> const& A
, std::size_t max_iterations
, std::size_t block_size
, T const& tolerance = 1.0
)
{
BOOST_ASSERT(2 == A.order());
BOOST_ASSERT(A.hypercube());
std::size_t const n = A.extent(0);
/*
std::vector<std::complex<T> > evs;
evs.reserve(n);
*/
std::complex<T> const nan_(std::numeric_limits<T>::quiet_NaN()
, std::numeric_limits<T>::quiet_NaN());
std::vector<std::complex<T> > evs(n, nan_), old(n, nan_);
orthotope<T> Ak = A.copy(), R, Q;
householders_tri_factor(Ak, block_size);
write_matrix_to_octave_file(Ak, "hess_A0");
std::size_t iterations = 0;
while (true)
{
/*
T const mu = Ak(n, n);
if (0 != iterations)
{
for (std::size_t i = 0; i < (n - 1); ++i)
Ak(i, i) -= mu;
Ak(n, n) = 0.0;
}
*/
householders_qr_factor(Ak, Q, R, block_size);
Ak = blocked_matrix_multiply(R, Q, block_size);
/*
if (0 != iterations)
for (std::size_t i = 0; i < n; ++i)
Ak(i, i) += mu;
*/
/*
bool pseudo_upper_triangular = true;
for (std::size_t j = 0; j < (n - 1); ++j)
{
// Make sure we're in Hessenberg form.
for (std::size_t i = j + 2; i < n; ++i)
{
if (!compare_floating(0.0, Ak(i, j), 1e-6))
pseudo_upper_triangular = false;
}
/// Check for convergence. Either we converge to 2x2 complex
/// conjugates eigenvalues, which take the form:
///
/// [ a b ]
/// [ c a ]
///
/// Where b * c < 0. Or, we converge to real eigenvalues which take
/// the form:
///
/// [ a ]
/// [ 0 ]
///
// Determine if we've failed to converge to a real eigenvalue.
if (!compare_floating(0.0, Ak(j, j + 1), 1e-6))
{
// Determine if we've failed to converge to a pair of complex
// eigenvalues.
if (!compare_floating(Ak(j, j), Ak(j + 1, j + 1), 1e-6))
pseudo_upper_triangular = false;
}
}
*/
bool converged = true;
// std::cout << "ITERATION " << iterations << "\n";
for (std::size_t j = 0; j < n; ++j)
{
if (j != n)
{
// Check for complex eigenvalues.
if (!compare_floating(0.0, Ak(j + 1, j), 1e-6))
//.........这里部分代码省略.........