5 template <
typename Scalar>
16 typedef DenseIndex
Index;
20 Scalar sum, temp, alpha, bnorm;
25 const Scalar epsmch = NumTraits<Scalar>::epsilon();
26 const Index n = qrfac.cols();
27 eigen_assert(n==qtb.size());
28 eigen_assert(n==x.size());
29 eigen_assert(n==diag.size());
33 for (j = n-1; j >=0; --j) {
36 temp = epsmch * qrfac.col(j).head(j+1).maxCoeff();
43 x[j] = (qtb[j] - qrfac.row(j).tail(n-j-1).dot(x.tail(n-j-1))) / temp;
47 qnorm = diag.cwiseProduct(x).stableNorm();
57 for (j = 0; j < n; ++j) {
58 wa1.tail(n-j) += qrfac.row(j).tail(n-j) * qtb[j];
64 gnorm = wa1.stableNorm();
66 alpha = delta / qnorm;
72 wa1.array() /= (diag*gnorm).array();
75 for (j = 0; j < n; ++j) {
77 for (i = j; i < n; ++i) {
78 sum += qrfac(j,i) * wa1[i];
82 temp = wa2.stableNorm();
83 sgnorm = gnorm / temp / temp;
93 bnorm = qtb.stableNorm();
94 temp = bnorm / gnorm * (bnorm / qnorm) * (sgnorm / delta);
95 temp = temp - delta / qnorm * numext::abs2(sgnorm / delta) +
sqrt(numext::abs2(temp - delta / qnorm) + (1.-numext::abs2(delta / qnorm)) * (1.-numext::abs2(sgnorm / delta)));
96 alpha = delta / qnorm * (1. - numext::abs2(sgnorm / delta)) / temp;
101 temp = (1.-alpha) * (std::min)(sgnorm,delta);
102 x = temp * wa1 + alpha * x;
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Namespace containing all symbols from the Eigen library.
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index