mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-01-18 14:34:17 +08:00
Adaptions from .lazy() towards .noalias().
Added missing casts.
This commit is contained in:
parent
bc7aec0ef5
commit
ab6eb6a1a4
@ -218,8 +218,8 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
|
||||
|
||||
int endSize = size - j - 1;
|
||||
if (endSize > 0) {
|
||||
_temporary.end(endSize) = ( m_matrix.block(j+1,0, endSize, j)
|
||||
* m_matrix.col(j).start(j).conjugate() ).lazy();
|
||||
_temporary.end(endSize).noalias() = m_matrix.block(j+1,0, endSize, j)
|
||||
* m_matrix.col(j).start(j).conjugate();
|
||||
|
||||
m_matrix.row(j).end(endSize) = m_matrix.row(j).end(endSize).conjugate()
|
||||
- _temporary.end(endSize).transpose();
|
||||
|
@ -141,7 +141,7 @@ template<> struct ei_llt_inplace<LowerTriangular>
|
||||
if (x<=RealScalar(0))
|
||||
return false;
|
||||
mat.coeffRef(k,k) = x = ei_sqrt(x);
|
||||
if (k>0 && rs>0) A21 -= (A20 * A10.adjoint()).lazy();
|
||||
if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
|
||||
if (rs>0) A21 *= RealScalar(1)/x;
|
||||
}
|
||||
return true;
|
||||
|
@ -244,11 +244,11 @@ void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort)
|
||||
// Apply Householder similarity transformation
|
||||
// H = (I-u*u'/h)*H*(I-u*u')/h)
|
||||
int bSize = high-m+1;
|
||||
matH.block(m, m, bSize, n-m) -= ((ort.segment(m, bSize)/h)
|
||||
* (ort.segment(m, bSize).transpose() * matH.block(m, m, bSize, n-m))).lazy();
|
||||
matH.block(m, m, bSize, n-m).noalias() -= ((ort.segment(m, bSize)/h)
|
||||
* (ort.segment(m, bSize).transpose() * matH.block(m, m, bSize, n-m)));
|
||||
|
||||
matH.block(0, m, high+1, bSize) -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize))
|
||||
* (ort.segment(m, bSize)/h).transpose()).lazy();
|
||||
matH.block(0, m, high+1, bSize).noalias() -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize))
|
||||
* (ort.segment(m, bSize)/h).transpose());
|
||||
|
||||
ort.coeffRef(m) = scale*ort.coeff(m);
|
||||
matH.coeffRef(m,m-1) = scale*g;
|
||||
@ -265,8 +265,8 @@ void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort)
|
||||
ort.segment(m+1, high-m) = matH.col(m-1).segment(m+1, high-m);
|
||||
|
||||
int bSize = high-m+1;
|
||||
m_eivec.block(m, m, bSize, bSize) += ( (ort.segment(m, bSize) / (matH.coeff(m,m-1) * ort.coeff(m)))
|
||||
* (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)) ).lazy();
|
||||
m_eivec.block(m, m, bSize, bSize).noalias() += ( (ort.segment(m, bSize) / (matH.coeff(m,m-1) * ort.coeff(m)))
|
||||
* (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)) );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
@ -44,21 +44,21 @@ template<typename MatrixType> void bandmatrix(const MatrixType& _m)
|
||||
dm1.diagonal().setConstant(123);
|
||||
for (int i=1; i<=m.supers();++i)
|
||||
{
|
||||
m.diagonal(i).setConstant(i);
|
||||
dm1.diagonal(i).setConstant(i);
|
||||
m.diagonal(i).setConstant(static_cast<RealScalar>(i));
|
||||
dm1.diagonal(i).setConstant(static_cast<RealScalar>(i));
|
||||
}
|
||||
for (int i=1; i<=m.subs();++i)
|
||||
{
|
||||
m.diagonal(-i).setConstant(-i);
|
||||
dm1.diagonal(-i).setConstant(-i);
|
||||
m.diagonal(-i).setConstant(-static_cast<RealScalar>(i));
|
||||
dm1.diagonal(-i).setConstant(-static_cast<RealScalar>(i));
|
||||
}
|
||||
//std::cerr << m.m_data << "\n\n" << m.toDense() << "\n\n" << dm1 << "\n\n\n\n";
|
||||
VERIFY_IS_APPROX(dm1,m.toDense());
|
||||
|
||||
for (int i=0; i<cols; ++i)
|
||||
{
|
||||
m.col(i).setConstant(i+1);
|
||||
dm1.col(i).setConstant(i+1);
|
||||
m.col(i).setConstant(static_cast<RealScalar>(i+1));
|
||||
dm1.col(i).setConstant(static_cast<RealScalar>(i+1));
|
||||
}
|
||||
int d = std::min(rows,cols);
|
||||
int a = std::max(0,cols-d-supers);
|
||||
|
@ -81,7 +81,7 @@ template<typename MatrixType> void product(const MatrixType& m)
|
||||
m3 = m1;
|
||||
m3 *= m1.transpose() * m2;
|
||||
VERIFY_IS_APPROX(m3, m1 * (m1.transpose()*m2));
|
||||
VERIFY_IS_APPROX(m3, m1.lazy() * (m1.transpose()*m2));
|
||||
VERIFY_IS_APPROX(m3, m1 * (m1.transpose()*m2));
|
||||
|
||||
// continue testing Product.h: distributivity
|
||||
VERIFY_IS_APPROX(square*(m1 + m2), square*m1+square*m2);
|
||||
@ -109,26 +109,26 @@ template<typename MatrixType> void product(const MatrixType& m)
|
||||
|
||||
// test optimized operator+= path
|
||||
res = square;
|
||||
res += (m1 * m2.transpose()).lazy();
|
||||
res.noalias() += m1 * m2.transpose();
|
||||
VERIFY_IS_APPROX(res, square + m1 * m2.transpose());
|
||||
if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1)
|
||||
{
|
||||
VERIFY(areNotApprox(res,square + m2 * m1.transpose()));
|
||||
}
|
||||
vcres = vc2;
|
||||
vcres += (m1.transpose() * v1).lazy();
|
||||
vcres.noalias() += m1.transpose() * v1;
|
||||
VERIFY_IS_APPROX(vcres, vc2 + m1.transpose() * v1);
|
||||
|
||||
// test optimized operator-= path
|
||||
res = square;
|
||||
res -= (m1 * m2.transpose()).lazy();
|
||||
res.noalias() -= m1 * m2.transpose();
|
||||
VERIFY_IS_APPROX(res, square - (m1 * m2.transpose()));
|
||||
if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1)
|
||||
{
|
||||
VERIFY(areNotApprox(res,square - m2 * m1.transpose()));
|
||||
}
|
||||
vcres = vc2;
|
||||
vcres -= (m1.transpose() * v1).lazy();
|
||||
vcres.noalias() -= m1.transpose() * v1;
|
||||
VERIFY_IS_APPROX(vcres, vc2 - m1.transpose() * v1);
|
||||
|
||||
tm1 = m1;
|
||||
@ -145,7 +145,7 @@ template<typename MatrixType> void product(const MatrixType& m)
|
||||
VERIFY_IS_APPROX(res, m1 * m2.transpose());
|
||||
|
||||
res2 = square2;
|
||||
res2 += (m1.transpose() * m2).lazy();
|
||||
res2.noalias() += m1.transpose() * m2;
|
||||
VERIFY_IS_APPROX(res2, square2 + m1.transpose() * m2);
|
||||
if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1)
|
||||
{
|
||||
|
@ -96,7 +96,7 @@ template<typename Scalar, int Size, int OtherSize> void symm(int size = Size, in
|
||||
|
||||
|
||||
m2 = m1.template triangularView<UpperTriangular>(); rhs13 = rhs12;
|
||||
VERIFY_IS_APPROX(rhs12 += (s1 * ((m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate())).lazy(),
|
||||
VERIFY_IS_APPROX(rhs12.noalias() += s1 * ((m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate()),
|
||||
rhs13 += (s1*m1.adjoint()) * (s2*rhs3).conjugate());
|
||||
|
||||
// test matrix * selfadjoint
|
||||
|
@ -90,9 +90,9 @@ namespace MatrixExponentialInternal {
|
||||
{
|
||||
typedef typename ei_traits<MatrixType>::Scalar Scalar;
|
||||
const Scalar b[] = {120., 60., 12., 1.};
|
||||
M2 = (M * M).lazy();
|
||||
M2.noalias() = M * M;
|
||||
tmp = b[3]*M2 + b[1]*Id;
|
||||
U = (M * tmp).lazy();
|
||||
U.noalias() = M * tmp;
|
||||
V = b[2]*M2 + b[0]*Id;
|
||||
}
|
||||
|
||||
@ -115,10 +115,10 @@ namespace MatrixExponentialInternal {
|
||||
{
|
||||
typedef typename ei_traits<MatrixType>::Scalar Scalar;
|
||||
const Scalar b[] = {30240., 15120., 3360., 420., 30., 1.};
|
||||
M2 = (M * M).lazy();
|
||||
MatrixType M4 = (M2 * M2).lazy();
|
||||
M2.noalias() = M * M;
|
||||
MatrixType M4 = M2 * M2;
|
||||
tmp = b[5]*M4 + b[3]*M2 + b[1]*Id;
|
||||
U = (M * tmp).lazy();
|
||||
U.noalias() = M * tmp;
|
||||
V = b[4]*M4 + b[2]*M2 + b[0]*Id;
|
||||
}
|
||||
|
||||
@ -141,11 +141,11 @@ namespace MatrixExponentialInternal {
|
||||
{
|
||||
typedef typename ei_traits<MatrixType>::Scalar Scalar;
|
||||
const Scalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.};
|
||||
M2 = (M * M).lazy();
|
||||
MatrixType M4 = (M2 * M2).lazy();
|
||||
MatrixType M6 = (M4 * M2).lazy();
|
||||
M2.noalias() = M * M;
|
||||
MatrixType M4 = M2 * M2;
|
||||
MatrixType M6 = M4 * M2;
|
||||
tmp = b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id;
|
||||
U = (M * tmp).lazy();
|
||||
U.noalias() = M * tmp;
|
||||
V = b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id;
|
||||
}
|
||||
|
||||
@ -169,12 +169,12 @@ namespace MatrixExponentialInternal {
|
||||
typedef typename ei_traits<MatrixType>::Scalar Scalar;
|
||||
const Scalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240.,
|
||||
2162160., 110880., 3960., 90., 1.};
|
||||
M2 = (M * M).lazy();
|
||||
MatrixType M4 = (M2 * M2).lazy();
|
||||
MatrixType M6 = (M4 * M2).lazy();
|
||||
MatrixType M8 = (M6 * M2).lazy();
|
||||
M2.noalias() = M * M;
|
||||
MatrixType M4 = M2 * M2;
|
||||
MatrixType M6 = M4 * M2;
|
||||
MatrixType M8 = M6 * M2;
|
||||
tmp = b[9]*M8 + b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id;
|
||||
U = (M * tmp).lazy();
|
||||
U.noalias() = M * tmp;
|
||||
V = b[8]*M8 + b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id;
|
||||
}
|
||||
|
||||
@ -199,15 +199,15 @@ namespace MatrixExponentialInternal {
|
||||
const Scalar b[] = {64764752532480000., 32382376266240000., 7771770303897600.,
|
||||
1187353796428800., 129060195264000., 10559470521600., 670442572800.,
|
||||
33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.};
|
||||
M2 = (M * M).lazy();
|
||||
MatrixType M4 = (M2 * M2).lazy();
|
||||
MatrixType M6 = (M4 * M2).lazy();
|
||||
M2.noalias() = M * M;
|
||||
MatrixType M4 = M2 * M2;
|
||||
MatrixType M6 = M4 * M2;
|
||||
V = b[13]*M6 + b[11]*M4 + b[9]*M2;
|
||||
tmp = (M6 * V).lazy();
|
||||
tmp.noalias() = M6 * V;
|
||||
tmp += b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id;
|
||||
U = (M * tmp).lazy();
|
||||
U.noalias() = M * tmp;
|
||||
tmp = b[12]*M6 + b[10]*M4 + b[8]*M2;
|
||||
V = (M6 * tmp).lazy();
|
||||
V.noalias() = M6 * tmp;
|
||||
V += b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id;
|
||||
}
|
||||
|
||||
@ -252,7 +252,7 @@ namespace MatrixExponentialInternal {
|
||||
} else if (l1norm < 1.880152677804762e+000) {
|
||||
pade5(M, Id, tmp1, tmp2, U, V);
|
||||
} else {
|
||||
const float maxnorm = 3.925724783138660;
|
||||
const float maxnorm = 3.925724783138660f;
|
||||
*squarings = std::max(0, (int)ceil(log2(l1norm / maxnorm)));
|
||||
MatrixType A = M / std::pow(typename ei_traits<MatrixType>::Scalar(2), *squarings);
|
||||
pade7(A, Id, tmp1, tmp2, U, V);
|
||||
@ -294,7 +294,7 @@ namespace MatrixExponentialInternal {
|
||||
{
|
||||
MatrixType num, den, U, V;
|
||||
MatrixType Id = MatrixType::Identity(M.rows(), M.cols());
|
||||
float l1norm = M.cwise().abs().colwise().sum().maxCoeff();
|
||||
float l1norm = static_cast<float>(M.cwise().abs().colwise().sum().maxCoeff());
|
||||
int squarings;
|
||||
computeUV_selector<MatrixType>::run(M, Id, num, den, U, V, l1norm, &squarings);
|
||||
num = U + V; // numerator of Pade approximant
|
||||
|
@ -43,10 +43,10 @@ void test2dRotation(double tol)
|
||||
A << 0, 1, -1, 0;
|
||||
for (int i=0; i<=20; i++)
|
||||
{
|
||||
angle = pow(10, i / 5. - 2);
|
||||
angle = static_cast<T>(pow(10, i / 5. - 2));
|
||||
B << cos(angle), sin(angle), -sin(angle), cos(angle);
|
||||
ei_matrix_exponential(angle*A, &C);
|
||||
VERIFY(C.isApprox(B, tol));
|
||||
VERIFY(C.isApprox(B, static_cast<T>(tol)));
|
||||
}
|
||||
}
|
||||
|
||||
@ -59,13 +59,13 @@ void test2dHyperbolicRotation(double tol)
|
||||
|
||||
for (int i=0; i<=20; i++)
|
||||
{
|
||||
angle = (i-10) / 2.0;
|
||||
angle = static_cast<T>((i-10) / 2.0);
|
||||
ch = std::cosh(angle);
|
||||
sh = std::sinh(angle);
|
||||
A << 0, angle*imagUnit, -angle*imagUnit, 0;
|
||||
B << ch, sh*imagUnit, -sh*imagUnit, ch;
|
||||
ei_matrix_exponential(A, &C);
|
||||
VERIFY(C.isApprox(B, tol));
|
||||
VERIFY(C.isApprox(B, static_cast<T>(tol)));
|
||||
}
|
||||
}
|
||||
|
||||
@ -77,13 +77,13 @@ void testPascal(double tol)
|
||||
Matrix<T,Dynamic,Dynamic> A(size,size), B(size,size), C(size,size);
|
||||
A.setZero();
|
||||
for (int i=0; i<size-1; i++)
|
||||
A(i+1,i) = i+1;
|
||||
A(i+1,i) = static_cast<T>(i+1);
|
||||
B.setZero();
|
||||
for (int i=0; i<size; i++)
|
||||
for (int j=0; j<=i; j++)
|
||||
B(i,j) = binom(i,j);
|
||||
B(i,j) = static_cast<T>(binom(i,j));
|
||||
ei_matrix_exponential(A, &C);
|
||||
VERIFY(C.isApprox(B, tol));
|
||||
VERIFY(C.isApprox(B, static_cast<T>(tol)));
|
||||
}
|
||||
}
|
||||
|
||||
@ -98,11 +98,13 @@ void randomTest(const MatrixType& m, double tol)
|
||||
MatrixType m1(rows, cols), m2(rows, cols), m3(rows, cols),
|
||||
identity = MatrixType::Identity(rows, rows);
|
||||
|
||||
typedef typename NumTraits<typename ei_traits<MatrixType>::Scalar>::Real RealScalar;
|
||||
|
||||
for(int i = 0; i < g_repeat; i++) {
|
||||
m1 = MatrixType::Random(rows, cols);
|
||||
ei_matrix_exponential(m1, &m2);
|
||||
ei_matrix_exponential(-m1, &m3);
|
||||
VERIFY(identity.isApprox(m2 * m3, tol));
|
||||
VERIFY(identity.isApprox(m2 * m3, static_cast<RealScalar>(tol)));
|
||||
}
|
||||
}
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user