inverse matrix of matrix with dealii::BlockSparseMatrix<double>

31 views
Skip to first unread message

Lance Zhang

unread,
Jun 22, 2023, 8:35:52 AM6/22/23
to deal.II User Group
Hello,

Currently, I met one problem when I calculate the inverse matrix with datatype of class dealii::BlockSparseMatrix<double>.

I searched the website of dealii for the method,but it seems there is no method to get the inverse matrix in BlockSparseMatrix<double>.

May I know whether there is any other method to change the matrix with data type BlockSparseMatrix<double> into a matrix with BlockSparseMatrix<double>?

Best regards
Lance

Lance Zhang

unread,
Jun 22, 2023, 8:54:13 AM6/22/23
to deal.II User Group
Hello,

The code is related to the file cook_membrane.cc.


I added the content below in main function
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
        double **get_matrix=new double*[get_rows];
        for(int i=0;i<get_rows;i++){
           get_matrix[i]=new double[get_cols];


        }

         for(int i=0;i<get_rows;i++){
            for(int j=0;j<get_cols;j++){
                 get_matrix[i][j]=double(solid_3d.tangent_matrix(i,j));

             }
          }
<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
But I got error by using gdb 
:[Thread 0x7fffd67fc640 (LWP 114719) exited]

Thread 1 "cook_membrane" received signal SIGSEGV, Segmentation fault.
0x00005555555be580 in dealii::BlockMatrixBase<dealii::SparseMatrix<double> >::operator() (j=24, i=0, this=0x7fffffffdde8) at /usr/include/deal.II/lac/block_matrix_base.h:2112
2112                                                  col_index.second);
(gdb) exit
A debugging session is active.

Inferior 1 [process 113432] will be killed.


But from the process information it seems the compiling is successful  and some result can be displayed.

In code,I wrote the content 2d array to obtain the value from solid_3d.tangent_matrix for inverse matrix:
<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

      std::cout<< "L1 norm:"<<solid_3d.system_rhs.l1_norm()<<std::endl;
         std::cout<< "L2 norm:"<<solid_3d.system_rhs.l2_norm()<<std::endl;

         int get_rows=solid_3d.tangent_matrix.m();
         int get_cols=solid_3d.tangent_matrix.n();


        double **get_matrix=new double*[get_rows];
        for(int i=0;i<get_rows;i++){
           get_matrix[i]=new double[get_cols];


        }

         for(int i=0;i<get_rows;i++){
            for(int j=0;j<get_cols;j++){
                 get_matrix[i][j]=double(solid_3d.tangent_matrix(i,j));

             }
          }
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

Could anyone provide any suggestion or hint?

Thanks in advance!
Lance

Wolfgang Bangerth

unread,
Jun 22, 2023, 1:29:18 PM6/22/23
to dea...@googlegroups.com
On 6/22/23 06:35, Lance Zhang wrote:
>
> May I know whether there is any other method to change the matrix with data
> type BlockSparseMatrix<double> into a matrix with BlockSparseMatrix<double>?

No. The inverse of a sparse matrix is, in general, not sparse. For all
problems of a realistic size, even if you can store a matrix A, you will not
be able to store A^{-1}. As a consequence, we never use the inverse of a
sparse matrix in deal.II, and there are no functions to compute inverses of
sparse matrices.

Best
W.

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


Lance Zhang

unread,
Jun 22, 2023, 4:00:44 PM6/22/23
to deal.II User Group
Hello  Wolfgang,

thanks for your message.

May I know if there is any method which can be used for  adjoint sensitivity method?

Best regards
Lance

Wolfgang Bangerth

unread,
Jun 22, 2023, 5:26:17 PM6/22/23
to dea...@googlegroups.com
On 6/22/23 14:00, Lance Zhang wrote:
>
> May I know if there is any method which can be used for  adjoint
> sensitivity method?

What is the connection between adjoint sensitivities and inverse matrices?

In general, if you ever find that you want to compute something like
y = A^{-1} x
you never do that by explicitly computing A^{-1}. Rather, you solve a
linear system
A y = x
with any of the iterative solvers such as CG, GMRES, ...
Reply all
Reply to author
Forward
0 new messages