On Fri, Apr 2, 2010 at 5:03 PM, Manoj Rajagopalan
<rmanoj@xxxxxxxxx> wrote:
Hi eigen-users,
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
DiagonalWrapper.
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 ;)
gael
Thanks,
Manoj