Hi Lucas,

You should be able to use syntax like

A(i, j) = 5;

to assign values to a FullMatrix - it's what we do in most of the tutorial 
programs. TableIndices are somewhat hard to use (you are right) - they are 
intended to be used as a generic interface to Table<N, T> where N is not known 
until compile time. Once you populate the matrix, you should be able to do what 
you mentioned before with the Gauss-Jordan algorithm to get the inverse.

That being said - do you need the inverse itself, or the action of the inverse? 
If you only ever need the action of the inverse a much better approach would be 
to use LAPACKFullMatrix and factorize it (via 
LAPACKFullMatrix::compute_lu_factorization() and LAPACKFullMatrix::solve()) 
since there is a significant amount of roundoff error accumulation in computing 
inverses.

The Solver classes are more than you need here - they are optimized for large 
sparse problems, not small dense ones.

Best,
David Wells


________________________________
From: dealii@googlegroups.com <dealii@googlegroups.com> on behalf of Lucas 
Myers <lucasmyer...@gmail.com>
Sent: Monday, June 21, 2021 1:50 PM
To: deal.II User Group <dealii@googlegroups.com>
Subject: [deal.II] 5x5 Matrix Inversion Problem

Hi Everyone,

I have a problem wherein one of the terms in my PDE is a 5-component vector 
function \Lambda of the 5-component solution vector Q. However, I only have an 
explicit formula for Q as a function of \Lambda and so I am trying to write a 
numerical inversion class (as a `dealii::Function`) to get \Lambda as a 
function of Q.

I am using Newton's method, and so I need to solve a 5x5 matrix equation. I was 
planning on just inverting the Jacobian which I have implemented as a 
`dealii::Tensor`, but it looks like the `invert` function is only implemented 
up to dimension 3 (?). It looks like the `dealii::FullMatrix` class can do this 
via a Gauss Jordan method, but I was having some trouble using the `FullMatrix` 
class directly. In particular, to read and write individual elements I need to 
create a `dealii::TableIndices` object, which seems tedious and makes me think 
it's meant to be used internally.

Given all this, I am looking for advice at any level: is this a reasonable way 
to numerically get \Lambda as a function of Q? Is it recommended to use a 
`dealii::Tensor` object as opposed to a `dealii::FullMatrix` object? And if I 
am supposed to use the latter, is there a nicer way to index it than with the 
`dealii::TableIndices` object? Lastly is it a good idea to be using an `invert` 
operation to solve the Newton's method linear equation? I know that the `Eigen` 
package recommends against that 
here<https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html>, and 
instead recommends one of its solution methods. I haven't seen any of those 
here besides one of the `Solver` classes, but those seem like a lot of overhead 
for a 5x5 matrix (I'm going to have to do this operation for every degree of 
freedom, so I want too minimize overhead). I also don't really want to use the 
Eigen package, just because (1) I want to avoid unnecessary value-copying and 
(2) I'd like to not have to try to interface multiple packages.

Any help is appreciated, and let me know if you need more information about the 
problem.

Kind regards,
Lucas

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to 
dealii+unsubscr...@googlegroups.com<mailto:dealii+unsubscr...@googlegroups.com>.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/29d3a1c1-6703-4c98-95a6-44ba96add049n%40googlegroups.com<https://groups.google.com/d/msgid/dealii/29d3a1c1-6703-4c98-95a6-44ba96add049n%40googlegroups.com?utm_medium=email&utm_source=footer>.

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/BN7PR03MB4356317B680523F10734A168ED0A9%40BN7PR03MB4356.namprd03.prod.outlook.com.

Reply via email to