Re: [eigen] makeHouseholder() |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] makeHouseholder()*From*: Christoph Hertzberg <chtz@xxxxxxxxxxxxxxxxxxxxxxxx>*Date*: Mon, 19 Sep 2016 17:02:50 +0200

At least in the latter case abs2(c0) should be small as well.

Btw:

Christoph [1] https://github.com/Reference-LAPACK/lapack/blob/master/SRC/clarfgp.f [2] https://github.com/Reference-LAPACK/lapack/blob/master/BLAS/SRC/scnrm2.f [3] https://github.com/Reference-LAPACK/lapack/blob/master/SRC/slapy2.f https://github.com/Reference-LAPACK/lapack/blob/master/SRC/slapy3.f On 2016-09-17 23:19, Gael Guennebaud wrote:

Sorry for late reply, but I don't see how this change can increase stability. Recall that the goal here is to compute a reflection H such that H*v=[x 0 ... 0] with x real. The initial condition simply check whether v was already in the desired form, and if so we set H=I, and we are done. This identity transformation will also be skipped during later computation. In contrast, the new condition won't produce an identity transformation in this case and will thus lost some bit of precision during the useless computations. gael On Fri, Sep 16, 2016 at 2:26 PM, Christoph Hertzberg < chtz@xxxxxxxxxxxxxxxxxxxxxxxx> wrote:Pushed: https://bitbucket.org/eigen/eigen/commits/cc01b1eb81f23a41 On 2016-09-16 14:09, Christoph Hertzberg wrote:On 2016-09-12 20:03, Peter wrote:Dear All, is there are specific reason for if(tailSqNorm <= tol && numext::abs2(numext::imag(c0))<=tol) or could it be changed into if(tailSqNorm <= tol && numext::abs2(c0)<=tol)As discussed at the beginning of the week on IRC, I second that change. I'll be pushing this change (if I don't get any objections within the next minutes). We should probably try to make some unit tests which would fail without this change (preferably not requiring custom types). Cheers, Christoph within makeHouseholder() in Householder.C, [*] see below?I'm asking as I stumbled over the specific case of c0=1.0+0.0i and a 2x2 matrix, where tailSqNorm==0. The current implementation with the use the special case tau = RealScalar(0); beta = numext::real(c0); essential.setZero(); while if(tailSqNorm <= tol && numext::abs2(c0)<=tol) would select the normal path in that case, beta = sqrt(numext::abs2(c0) + tailSqNorm); if (numext::real(c0)>=RealScalar(0)) beta = -beta; essential = tail / (c0 - beta); tau = conj((beta - c0) / beta); which is fine to evaluate. Due to sign choice of beta, (c0 - beta) will even be further away from zero as beta itself, so numerically that shouldn't be a problem. I'm asking because I could pin-point, with the help of Christoph on IRC, this particular if-statement causing problems for my own custom type. Best regards, Peter [*] template<typename Derived> template<typename EssentialPart> void MatrixBase<Derived>::makeHouseholder( EssentialPart& essential, Scalar& tau, RealScalar& beta) const { using std::sqrt; using numext::conj; EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1); RealScalar tailSqNorm = size()==1 ? RealScalar(0) : tail.squaredNorm(); Scalar c0 = coeff(0); const RealScalar tol = (std::numeric_limits<RealScalar>::min)(); if(tailSqNorm <= tol && numext::abs2(numext::imag(c0))<=tol) { tau = RealScalar(0); beta = numext::real(c0); essential.setZero(); } else { beta = sqrt(numext::abs2(c0) + tailSqNorm); if (numext::real(c0)>=RealScalar(0)) beta = -beta; essential = tail / (c0 - beta); tau = conj((beta - c0) / beta); } }-- Dipl. Inf., Dipl. Math. Christoph Hertzberg Universität Bremen FB 3 - Mathematik und Informatik AG Robotik Robert-Hooke-Straße 1 28359 Bremen, Germany Zentrale: +49 421 178 45-6611 Besuchsadresse der Nebengeschäftsstelle: Robert-Hooke-Straße 5 28359 Bremen, Germany Tel.: +49 421 178 45-4021 Empfang: +49 421 178 45-6600 Fax: +49 421 178 45-4150 E-Mail: chtz@xxxxxxxxxxxxxxxxxxxxxxxx Weitere Informationen: http://www.informatik.uni-bremen.de/robotik

-- Dipl. Inf., Dipl. Math. Christoph Hertzberg Universität Bremen FB 3 - Mathematik und Informatik AG Robotik Robert-Hooke-Straße 1 28359 Bremen, Germany Zentrale: +49 421 178 45-6611 Besuchsadresse der Nebengeschäftsstelle: Robert-Hooke-Straße 5 28359 Bremen, Germany Tel.: +49 421 178 45-4021 Empfang: +49 421 178 45-6600 Fax: +49 421 178 45-4150 E-Mail: chtz@xxxxxxxxxxxxxxxxxxxxxxxx Weitere Informationen: http://www.informatik.uni-bremen.de/robotik

**Follow-Ups**:**Re: [eigen] makeHouseholder()***From:*Peter

**References**:**[eigen] makeHouseholder()***From:*Peter

**Re: [eigen] makeHouseholder()***From:*Christoph Hertzberg

**Re: [eigen] makeHouseholder()***From:*Christoph Hertzberg

**Re: [eigen] makeHouseholder()***From:*Gael Guennebaud

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Permutation times Sparse matrix** - Next by Date:
**Re: [eigen] Using Cholmods LDLT factorization with eigen.** - Previous by thread:
**Re: [eigen] makeHouseholder()** - Next by thread:
**Re: [eigen] makeHouseholder()**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |