|Re: [eigen] Efficient syntax for Ax=b where A is diagonal/Toeplitz|
[ Thread Index |
| More lists.tuxfamily.org/eigen Archives
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: Re: [eigen] Efficient syntax for Ax=b where A is diagonal/Toeplitz
- From: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
- Date: Fri, 2 Apr 2010 20:21:59 +0200
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :from:date:received:message-id:subject:to:content-type; bh=F/ib0l9+038D2PPkuky+UWwbfE7Tvr2W762a0FjX1UQ=; b=yBZXrtPj3Sg9VzRtNOZtsJJPdsEoK0PkG+0x1VZOt/oywmOIF+uWZXr9SN2E2uuP8N 4JP9BWuavYNGxagAoDCPq8Wr5IeZFWYBZ/LiHecfVmAUDY62wp1Rx7GUv5VsImnfb/zA p9jk4m13A+0mD/p3s8NxzM+RMawr2onuuss6Y=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :content-type; b=fdp/Ts/nPf2uCskyLOg1BtRbeRR9BVixZxQ2fKe3BOkmx0EFolXBWWvoHicniAFbhs eej4lJew3W9Hha1NxbjJCYHjlP1qUTT0fHdtRgjHzFVr5h7yfdmilLidicZVC2fjr2MB FXiV4PjRB2s7eM3B04XEZ3Gu5PV95FCgqr02U=
On Fri, Apr 2, 2010 at 5:03 PM, Manoj Rajagopalan <rmanoj@xxxxxxxxx>
I come across cases where my physics generates different types of
matrices, sometimes diagonal, sometimes Toeplitz depending on the simulation
region. My PDE solution code requires me to solve systems of the form AX=B
where A takes such different patterns. For readability it would be nice to
introduce code that looks likes B = A.solve(X).
we already support such API via a LU, LLT or SVD decomposition or via the triangularView API for triangular matrices.
The Eigen tutorial mentions A.lu().solve(X) for dense matrices but I
couldn't locate similar info for diagonal matrices for which the solution
procedure is very simple. I could explicitly code this in a few lines but I
am just wondering if the above syntactic-sugar is possible with Eigen.
if A is diagonal, "A.inverse() * X" will do the job with full lazy evaluation.
When B,X are dense, column-major, the procedure, in slightly-abused MATLAB
notation, is X[:,j]=B[:,j] ./ A.diagonal()[:]. When B,X are row-major, the
fast procedure is X[i,:] = B[i,:] / A[i,i]. We may be more efficient if A is
triangular or diagonal.
If the convenient syntax doesn't exist, would it be possible to add this? I
don't see any solve() or solveInPlace() methods on classes DiagonalMatrix or
well, for consistency and genericity we could add a solve() function to Diagonal, why not.
Also, efficient storage, arithmetic and linear-system solution (O(n^2))
algorithms are possible when A is (symmetric) Toeplitz (occur frequently in
signal-processing and Fourier-series solution of ODE/PDE). Does the roadmap
consider this case?
I'd say contributions are welcome ;)