get the partition of the system matrix A associated with the unconstrained dofs

[Simon]

Dear all,

some degrees of freedom of my solution vector are constrained.

When solsing the linear system A*x=b, I use the first approach as described in "Constraints on degrees of freedom".

That is, I call

constraints.distribute_local_to_global(...,false)

and, after solving the linear system,

constraints.distribute(...).

My question pertains to the following:

Say, I have in total 20 dofs and the dof with global dof index Zero is constrained.

As a consequence, the first component of the rhs b is set to Zero as well as

the first row and column of the system matrix A (except the diagonal value).

In a postprecessing step, I have to solve another linear system, however, with the system matrix being only the 19x19 matrix associated with the 19 unconstrained dofs; I do not need the rhs anymore.

Is there a way to get this portion of the system matrix based on the existing system matrix (sparsity pattern)?

Thank you,

Simon

[Wolfgang Bangerth]

On 8/18/22 10:38, Simon wrote:
> Say, I have in total 20 dofs and the dof with global dof index Zero is
> constrained.
> As a consequence, the first component of the rhs b is set to Zero as well as
> the first row and column of the system matrix A (except the diagonal value).
> In a postprecessing step, I have to solve another linear system, however, with
> the system matrix being only the 19x19 matrix associated with the 19
> unconstrained dofs; I do not need the rhs anymore.
>
> Is there a way to get this portion of the system matrix based on the existing
> system matrix (sparsity pattern)?

Both the sparsity pattern and the sparse matrix have iterators that allow you
to iterate over all entries of a matrix. You can do that and just filter out
those rows and columns you're not interested in, and then copy the rest into
output objects with translated indices.

But the easier approach may be to use the same 20x20 matrix and just copy the
new rhs you want to solve with into a vector with size 20, leaving the entries
of the rhs vector that correspond to constrained DoFs zero (or, in fact,
whatever you want -- the value there doesn't matter). By zeroing out rows and
columns, you are in essence solving a linear system with only the 19x19
matrix. You then copy the 19 DoFs you care about out of the solution vector
into an output vector.

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Simon Wiesheier]

" But the easier approach may be to use the same 20x20 matrix and just copy the
new rhs you want to solve with into a vector with size 20, leaving the entries
of the rhs vector that correspond to constrained DoFs zero (or, in fact,
whatever you want -- the value there doesn't matter). By zeroing out rows and
columns, you are in essence solving a linear system with only the 19x19
matrix. You then copy the 19 DoFs you care about out of the solution vector

into an output vector."

I will assemble the new rhs as vector with size 20.

If the values of the new rhs correspondong to constrained DoFs are irrlevant,

I do not have to modify the new rhs as well as the 20x20 system matrix at all. 

But let me make a short example with the following linear 3x3 linear system A*x = b:

A =

[K00   K01     0  ;

 K10   K11   K12;

   0     K21   K22;]

x=[x0; x1; x2] 

b=[b1, b2; b3]

Say x0=c=const.

The linear system looks like this after the assembly routine:

A =

[K00     0        0  ;

   0     K11   K12;

   0     K21   K22;]

x=[0 ; x1; x2] 

b=[0 ; b2-c*K10; b3]

This system boilds down to a 2x2 system for x1 and x2 with x0=0.

This is exactly what I want to compute, but without having -c*K10 substracted.

(Because the new rhs comes from a different problem and has nothing to do with the constrainted dofs - I just need the 2x2 portion from the original problem)

Currently I only deal with homogeneous constraints (x0=0).

In this case constraints.distribute(...) does not change the unconstrained entries of the rhs, right?

But I will also have to deal with inhomogeneous constraints, in which case it makes a difference.

That said, is there a way to leave the components of the rhs corresponding to unconstrained DoFs unchanged?

Best

Simon

--
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/17c15884-e251-1c38-f563-16bec6ad37a7%40colostate.edu.

[Wolfgang Bangerth]

On 8/19/22 05:46, Simon Wiesheier wrote:
>
> This system boilds down to a 2x2 system for x1 and x2 with x0=0.
> This is exactly what I want to compute, but without having -c*K10 substracted.
> (Because the new rhs comes from a different problem and has nothing to do with
> the constrainted dofs - I just need the 2x2 portion from the original problem)

Right. But c*K10 is subtracted from the right hand side when you call
constraints.distribute_local_to_global (...)
It is not magically stored somewhere. When you solve a new linear system with
the matrix, that linear system knows nothing about what happened when you
first built the matrix and the original right hand side.

[Simon Wiesheier]

" When you solve a new linear system with
the matrix, that linear system knows nothing about what happened when you

first built the matrix and the original right hand side."

Yes, but I have to call constraints.distribute_local_to_global (...)

also when building the new linear system. But I observed that there is also

a constraints.distribute_local_to_global () function which writes only into

the global matrix.

I also need the system matrix A for a second purpose, namely

to compute a matrix multiplication:

res = A^{-1} * B ,

where B is another matrix.

-To be more precise, I need the inverse of the 19x19 submatrix

corresponding to the unconstrained DoFs only -- not the inverse of the full system matrix..

I could not find a function which computes the inverse of a sparse matrix directly (without solving a linear system).

What I tried is,

LAPACKFullMatrix<double> new_matrix = my_system_matrix ,

thence calling the invert function.

But I am not sure if this is the right way to go.

-Also, after calling constraints.distribute_local_to_global(),

does it make sense at all to compute an inverse matrix, given that some rows and columns were set to zero?

Thank you again!

Best

Simon

--
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/66de0c9f-8dce-48a3-a7fb-6a047bc4f1c7%40colostate.edu.

[Wolfgang Bangerth]

On 8/19/22 13:14, Simon Wiesheier wrote:
>
> I also need the system matrix A for a second purpose, namely
> to compute a matrix multiplication:
> res = A^{-1} * B ,
> where B is another matrix.
> -To be more precise, I need the inverse of the 19x19 submatrix
> corresponding to the unconstrained DoFs only -- not the inverse of the
> full system matrix..

Right. But the inverse of the 19x19 matrix is the 19x19 subblock of the
20x20 big matrix. That's because after zeroing out the row and column,
you have a block diagonal matrix where the inverse of the matrix
consists of the inverses of the individual blocks.

> I could not find a function which computes the inverse of a sparse
> matrix directly (without solving a linear system).
> What I tried is,
> LAPACKFullMatrix<double> new_matrix = my_system_matrix ,
> thence calling the invert function.
> But I am not sure if this is the right way to go.

That's one way to go. FullMatrix::gauss_jordan() also computes the
inverse of a matrix.

> -Also, after calling constraints.distribute_local_to_global(),
> does it make sense at all to compute an inverse matrix, given that some
> rows and columns were set to zero?

Yes -- as mentioned above, you should consider the resulting matrix as
block diagonal.

[Simon Wiesheier]

Thanks for your input.

In the meantime, I replaced the matrix multiplication

res = A^_{-1}*B

by solving 'p' linear systems

A*res[p] = B[p],

where p is the number of columns of the matrix B.

" That's one way to go. FullMatrix::gauss_jordan() also computes the
inverse of a matrix."

As stated, what I tried is to use the operator= according to

LAPACKFullMatrix<double> new_matrix = my_system_matrix .

However, there is an error message

"error: conversion from ‘dealii::SparseMatrix<double>’ to non-scalar type ‘dealii::LAPACKFullMatrix<double>’ requested
   LAPACKFullMatrix<double> new_matrix = tangent_matrix"

How can I fix this?

Best

Simon

--
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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/4Swu5xeNU3U/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/593cd135-fa19-1572-23a9-53bcfb34aad9%40colostate.edu.

[Wolfgang Bangerth]

On 8/22/22 10:08, Simon Wiesheier wrote:
>
> As stated, what I tried is to use the operator= according to
> LAPACKFullMatrix<double> new_matrix = my_system_matrix .
> However, there is an error message
> "error: conversion from ‘dealii::SparseMatrix<double>’ to non-scalar
> type ‘dealii::LAPACKFullMatrix<double>’ requested
>    LAPACKFullMatrix<double> new_matrix = tangent_matrix"

There doesn't appear a copy-constructor from FullMatrix (which is what
the compiler is looking for here), but there is a copy operator. Just write
LAPACKFullMatrix<double> new_matrix;
new_matrix = tangent_matrix

Even better, of course, if you wrote a patch to add the copy constructor!