5x5 Matrix Inversion Problem

68 views
Skip to first unread message

Lucas Myers

unread,
Jun 21, 2021, 1:50:12 PM6/21/21
to deal.II User Group
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, 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

Wells, David

unread,
Jun 21, 2021, 2:08:50 PM6/21/21
to deal.II User Group
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: dea...@googlegroups.com <dea...@googlegroups.com> on behalf of Lucas Myers <lucasm...@gmail.com>
Sent: Monday, June 21, 2021 1:50 PM
To: deal.II User Group <dea...@googlegroups.com>
Subject: [deal.II] 5x5 Matrix Inversion Problem
 
--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/29d3a1c1-6703-4c98-95a6-44ba96add049n%40googlegroups.com.

Lucas Myers

unread,
Jun 22, 2021, 4:04:58 PM6/22/21
to deal.II User Group
Hi David,

Thanks so much, this is exactly what I was looking for! For some reason I was trying to use the square bracket notation. Is there any reason why the (i, j) syntax is used rather than [i, j]? Is it just a restriction of c++?

A few follow-up questions about the LAPACKFullMatrix class: is there anywhere that I can find more information about the class besides the reference documentation here? In particular, I'm wondering how to deal with the state attribute. It looks like the compute_lu_factorization() method sets the state to lu (reasonable), but I'd like to repopulate the same matrix (meaning it has to be in the matrix state, I think). For now I'm using reinit() because that seems to change the state back to matrix, but I worry that it's expensive to keep reinitializing the matrix every Newton's method iteration. Is there a way to avoid reinitializing it? And is it bad programming practice to not reinitialize it?

Again, thanks so much for the help.

Kind regards,
Lucas

Wells, David

unread,
Jun 23, 2021, 9:15:30 AM6/23/21
to deal.II User Group
Hi Lucas,

C++ doesn't allow us to overload operator[] with multiple arguments so that indexing has to be with single values - so yes, it's a language thing. This matches the C usage of square brackets where they only accept one value at a time.

The state attribute is automatically set by whatever functions you call - you should never need to mess with it (and in debug mode the class should abort if you try to do something that is inconsistent with the current state of the matrix). Internally, LAPACK typically modifies matrices in-place so LAPACKFullMatrix behaves in a similar way. The only other place I would look for documentation is that class's own page:


If you have to do Newton iteration by setting up a Jacobian and solving a matrix then what you did makes sense. No matter what, if you are doing iterations, you have to compute a new Jacobian and solve a new linear system. The reinit() call is not expensive - it just zeros the data, which is O(n^2) for an n x n matrix (as opposed to the LU factorization which is O(n^3)).

If you want to save some time you can reuse Jacobians (and therefore their factorizations) in Newton's method but that will cause things to converge more slowly, but it can be tricky to determine how often to recompute the Jacobian.

Is it bad programming practice to not reinitialize it? It isn't the worst thing you can do, but reinitializing it to guarantee that there aren't any irrelevant values stored is the right choice. In this case I think you have to reinitalize the matrix in order to keep the state of the matrix consistent (it doesn't make sense to modify a matrix containing an LU factorization computed by LAPACK unless you have some knowledge of how LAPACK internally stores the factorization), so you need to do it no matter what.

Best,
David Wells
Sent: Tuesday, June 22, 2021 4:04 PM

To: deal.II User Group <dea...@googlegroups.com>
Subject: Re: [deal.II] 5x5 Matrix Inversion Problem
 
Reply all
Reply to author
Forward
0 new messages