Solve block Trilinos matrix directly

[yy.wayne]

I've post a similar question before(Solve distributed BlockVector with direct solver), where I had a Trilinos SparseMatrix(the system has 2 components) and LinearAlgebra::distributed::BlockVector as dst and src for that SparseMatirx. I choose  LinearAlgebra::distributed::BlockVector because the bigger problem is intend to solve by MatrixFree operators, and use direct solver on coarsest grid.

However, by copy to and from between LinearAlgebra::distributed::BlockVector and LinearAlgebra::distributed::Vector, it produces wrong result. Though I cannot confirm the wrong result comes from transfer between Vector and BlockVector(I tested simple x=inv(A)*b and result is right, but treating it as the coarse solver for bigger problem the result is wrong), I reassembled the matrix in Trilinos BlockSparseMatrix and solve it with a iterative solver and the result is correct. While solving coarse problem iteratively works for very small matrix, my system is indefinite so I have to use a direct solver.

So given a BlockSparseMatrix and BlockVector src and dst, is there a recommended strategy to solve it directly(or block diagonal precondition)? The matrix has form [A, B; -B^T A]

[Wolfgang Bangerth]

On 11/21/22 07:45, 'yy.wayne' via deal.II User Group wrote:
>
> However, by copy to and from between LinearAlgebra::distributed::BlockVector
> and LinearAlgebra::distributed::Vector, it produces wrong result. Though I
> cannot confirm the wrong result comes from transfer between Vector and
> BlockVector(I tested simple x=inv(A)*b and result is right, but treating it as
> the coarse solver for bigger problem the result is wrong), I reassembled the
> matrix in Trilinos BlockSparseMatrix and solve it with a iterative solver and
> the result is correct. While solving coarse problem iteratively works for very
> small matrix, my system is indefinite so I have to use a direct solver.
>
> So given a BlockSparseMatrix and BlockVector src and dst, is there a
> recommended strategy to solve it directly(or block diagonal precondition)? The
> matrix has form [A, B; -B^T A]

No, not without copying data back and forth.

But you can test whether the copying introduces the error. Let's say you have
copied the data into a nonblock system, solved it, and copied it back to the
block system. Let's say your block system consists of a matrix A, right hand
side b, and the solution vector x you obtained by copying the solution of the
nonblock version into x. All of these objects consist of blocks.

Then if x solved the nonblock version, you would expect that the residual,
b-Ax
is small. You can test that via the BlockSparseMatrix::residual() function.
Specifically, you would expect that
||b-Ax|| / ||b||
is, say, 1e-12. You can compute this via
BlockVector r(...);
double relative_residual = A.residual (r,x,b) / b.l2_norm();

If the relative residual is not small, then something must have gone wrong
during the copying step. If it is small, then the solution of the nonblock
system, copied into the block vector x, is also the solution of the block system.

Best
W.

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

[yy.wayne]

Hi Prof. Bangerth,

By comparing residuals between normal matrix, block matrix, and direct and iterative solvers as you said, I believe the mistake is somewhere else. Thank you.

Best,

Wayne