[eigen] solve API |

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

*To*: eigen <eigen@xxxxxxxxxxxxxxxxxxx>*Subject*: [eigen] solve API*From*: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>*Date*: Mon, 28 Jun 2010 11:42:49 +0200*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:received:mime-version:received:from:date :message-id:subject:to:content-type; bh=MI9CSgSJJSVmwiH0fJWl8QQBk7CktLvt0+ybL6KuhT8=; b=g+aOwteMkujP1Z4MEV14nxr2iNbgoqBdaDPzZy6nfFqBB+wQ6q8SCOI7S3Fr8tkjEA jpXbw1I2uxUaHQFoNDHOECo+hMPGvtXsjxmfuPmK+34yJS5TebIqcDU4GFfVOkcHsLyW UVDPKH0BUv85dVzMJw+7oRMULsHLLQUE4N/ls=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:from:date:message-id:subject:to:content-type; b=d/Vp3H2sTs4VI11GQQnkkANeHre6afSW0ScN/Ekr3WOaZPakHvM/Z/i7+Gmwr6gU4s fZuLC0sQyKa9RGIdHLJTOMRMK1zAOd038ST7ma49hGXb4zbHfH+rQ3E9SLeqaXQoejRP YrIgdX+SU4TwT2+UTMVjbE8LUQwucijcEyjWY=

Hi, I was thinking that we did not discussed the solved API for a while, at so it would be perfect time to discuss it again ;) Seriously, currently our solvers only allows you to compute x = A^1 b via some decompositions of A. However, it is sometime needed to also solve for x = A^-T b or x = b A^1. A naive way to add support for that would be to add template options to the solve functions, as we currently do for triangular solving: M.triangularView<Lower>().solveInPlace<OnTheRight>(b); <=> b = b * lower(M)^-1; For "b = lower(M)^-* * b" you do: M.triangularView<Lower>().adjoint().solveInPlace(b); (the "OnTheLeft" is implicit) To extend this API to the solver we could naively add a "Transpose" or "Adjoint" option: x = m.lu().solve<OnTheRight | Adjoint>(b); which is ... well, ... pure Lapack style ;) This API is really confusing and not clear at all. Note that "OnTheRight/Left" is also used in function names such as: a.applyOnTheRight(b); Here it is pretty clear that's a which is on the right, and so we are doing: b = b * a; However, x = m.lu().solve<OnTheRight>(b); makes no clear sense. It is not clear whether it is b which is on the right, or the inverse of m. Recall that one goal is to be able to write code as you write math. So here is the proposal: x = m.lu().inverse() * b; x = b * m.lu().inverse(); x = m.lu().adjoint().inverse() * b; x = b * m.lu().adjoint().inverse(); ..inverse() and .adjoint() would return simple proxy with operator * overloaded to return a "solve" proxy that we already have (so implementation-wise that's not much). I know that computer scientists used to other linear algebra packages might be shocked by such an API and will look forever for a function having "solve" in its name, but this is pure Eigen style IMO and it is consistent with the rest of the API (diagonal, permutation, transposition already use this API). It remains the problem of convenient "inplace" solving. Currently, when doing: v = m.lu().solve(v); we can detect that we have the same object for the "x" and the "b" via one runtime pointer comparison, and then skip the matrix copy of b to x before the inplace solving path. So the same mechanism would still work: b = b * m.lu().adjoint().inverse(); => inplace solving, no temporary, no matrix copy (except for the decomposition itself but that's another story). The pb, is when "b" is a complex expression: b.block(i+s, j+block_size-1, rows-r, std::min(blocksize,cols)) =m.lu().inverse() * b.block(i+s, j+block_size-1, rows-r, std::min(blocksize,cols)); You don't want to write it twice!! Here a solveInPlace() function could make sense: m.lu().solveInPlace(b.block(i+s, j+block_size-1, rows-r, std::min(blocksize,cols))); better ? not so sure: 1) what about "ontheright" and transposed variants ? 2) and even worse, this does not work because block(...) returns a temporary object which has to be passed by const reference, so you have to name it anyway: Block<...> b_block(n,i+s, j+block_size-1, rows-r, std::min(blocksize,cols)); and then repeating twice b_block is fine. As some of you probably guessed, I did not slept well this night.... cheers, gael

**Follow-Ups**:**Re: [eigen] solve API***From:*Benoit Jacob

**Re: [eigen] solve API***From:*Manoj Rajagopalan

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] block of dense-block should again be block** - Next by Date:
**[eigen] Tutorial example does not complie** - Previous by thread:
**Re: [eigen] (sparse) matrix multiplication** - Next by thread:
**Re: [eigen] solve API**

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