What is the good practice to invert a general blocked BCRSMatrix

[GitPaean]

I have a question about how to invert a BCRSMatrix block matrix in DUNE. The following is the definition of the Matrix type. NumWellEq are usually 2, 3, 4 .. maybe 5 for now.

82         typedef Dune::FieldMatrix<Scalar, numWellEq, numWellEq > DiagMatrixBlockWellType;
83         typedef Dune::BCRSMatrix <DiagMatrixBlockWellType> DiagMatWell;

The the size of the matrix I think will go to maximum hundreds of block size. The matrix is not diagonal, so I can not just invert each diagonal block.

I am wondering what is the good practice to invert this type of matrix with DUNE. Basically, we generate this matrix, and then only use its inversion in a few places.

Thanks in advance.

[dr-robertk]

@GitPaean You'll have to use a direct solver for that. How did you do this before?

[GitPaean]

For Standard Well, the matrix is diagonal, so we can just invert() each diagonal block.

For Multisegment Well, the matrix is not diagonal, and there is no fixed pattern.

[dr-robertk]

@GitPaean : Do you already have a working version based on the old code? Then my question was how did you invert such matrices there. If you don't have that piece of code then of course the question is pointless.

[GitPaean]

For StandardWell, invert() from FieldMatrix is called (same with StandardWellsDense.hpp in master branch now.). For Multisegment Well, I have not reached that point yet. This is one of the few remaining few technique difficulties.

[akva2]

do you need to form the inverted matrix, or do you just need the action of the inverse on a vector? that's the fundamental question here.

[dr-robertk]

I'm not sure if there is a version of a direct solver implemented in dune-istl for that. You might need to rely on UMFPACK or SuperLU with entires being doubles, i.e. FieldMatrix<double,1,1>. In the case of @akva2 is right a standard iterative solver will do the trick.

[akva2]

if you need to form the inverted matrix, you can form it by repeatedly applying the lu decomposition to (1 0 0 0), (0 1 0 0), ... and sticking these as columns in a matrix. if you only need the action, you just use an iterative or a direct solver to solve for the vector.

[dr-robertk]

There is an LU implementation in dune-common that could be reused here. But this might be a bit of work.

[GitPaean]

I think I need both for now, so if we can formulate the inverted matrix, I can reuse the current solution process/framework. I can look into using lu decomposition to formulate the inverted matrix. I will come back with more information when I reach the point. Thank @akva2 and @dr-robertk for the information.

[dr-robertk]

One more thing. If the matrix is a tri-diagonal block matrix then there exist a version of the Thomas algorithm for BCRSMatrix. You might want to check that.

[GitPaean]

One more thing. If the matrix is a tri-diagonal block matrix then there exist a version of the Thomas algorithm for BCRSMatrix.

Thanks. It is good to know. It is possible we will use it when possible.

For a single branched well, the matrix will be a tri-diagonal block matrix. For multi-branched wells, there will be extra entries.

[blattms]

My two cents:

• Both UMFPack and SuperLU bindings in dune-istl can cope with real matrix blocks (FieldMatrix<double,n,n>
• If you do not absolutely have to, then avoid storing the inverse. (Store the LU decomposition or use a direct solver). Otherwise the code will be less numerically stable.

[GitPaean]

A small update. I think I only use the action to calculate the D^-1 x (D is the matrix and x is a vector) without inverting matrix D. Sorry for the misleading information. D^-1 is involved in the middle of the a chain of matrix multiplication. So as @akva2 and @dr-robertk said, the vector can be solved with a iterative linear solver.

[GitPaean]

Hi, again.

The matrix D needs the inversion operation turns out to be a tri-diagonal matrix plus some off-diagonal entries. The pattern is symmetric while the values of the entries are not symmetric. I wrote (modified) the following code to obtain D^-1 * x. So far, it looks like working (have not checked the exact values in details.). Any suggestion to improve it, basically means lighter in term of creation and more efficient in term of solving, since it will happen a few times each iteration. The possible maximum size of the matrix can be hundreds of blocks.

Thanks.

338     // obtain y = D^-1 * x
339     template<typename MatrixType, typename VectorType>
340     VectorType
341     invDX(MatrixType D, VectorType x)
342     {
343         // TODO: checking the problem related to use reference parameter
344         VectorType y(x.size());
345         y = 0.;
346         Dune::MatrixAdapter<MatrixType, VectorType, VectorType> linearOperator(D);
347 
348         // Sequential incomplete LU decomposition as the preconditioner
349         Dune::SeqILU0<MatrixType, VectorType, VectorType> preconditioner(D, 1.0);
350 
351         // Preconditioned BICGSTAB solver
352         Dune::BiCGSTABSolver<VectorType> linsolver(linearOperator,
353                                                    preconditioner,
354                                                    1.e-6, // desired residual reduction factor
355                                                    50, // maximum number of iterations
356                                                    0); // verbosity of the solver
357 
358         // Object storing some statistics about the solving process
359         Dune::InverseOperatorResult statistics ;
360 
361         // Solve
362         linsolver.apply(y, x, statistics );
363 
364         return y;
365     }

[GitPaean]

It looks like the result is generally okay compared with the result from MATLAB/UMFPACK. The result from above function

-1.01072 4.29344e-17 0 0
-0.0393548 4.31181e-14 4.31181e-16 1.08771e+06

The result from MATLAB

         -1.01072221901382
     -4.33680868994202e-19
      8.67361737988404e-19
                         0
       -0.0393546942104009
      2.32240650497066e-27
      2.01948391736579e-28
                   1087710

[GitPaean]

I consider the problem is solved, closing the issue now. Might open a different one when other related issue comes up later.