Assembling sparse matrix from Matrix-free vmult with constrains

[Michał Wichrowski]

Dear  all,

I've got matrix-free multigrid solver for Stokes problem. The main bottleneck is solution of coarse problem, so I tried to assemble the regular sparse matrix and use direct solver. Since the coarse problem is (relatively) small, I used vmults by unit vector to obtain columns of the matrix.  This is my code:

    setup_matrix(); // generates sparsity patters and reinits matrices

LinearAlgebra::distributed::BlockVector<LevelNumber>  dst;

LinearAlgebra::distributed::BlockVector<LevelNumber>  src;

src.reinit(2);

    for (unsigned int b = 0; b < 2; ++b)

    stokes_matrix.get_matrix_free()->initialize_dof_vector(

    src.block(b), b);

    src.collect_sizes();

    src =0;

    dst.reinit(2);

    dst.block(0).reinit(owned_dofs_u, relevant_dofs_u, MPI_COMM_WORLD);

    dst.block(1).reinit(owned_dofs_p, relevant_dofs_p, MPI_COMM_WORLD);

    dst.collect_sizes();

for(types::global_dof_index i =0; i< owned_dofs_u.size(); ++i){

src=0;

dst=0;

if(owned_dofs_u.is_element(i) )

src.block(0)(i)=1;

src.compress(VectorOperation::insert);

stokes_matrix.vmult(dst, src);

dst.update_ghost_values();

for(IndexSet::ElementIterator  index_iter =owned_dofs_u.begin();

index_iter != owned_dofs_u.end();

++index_iter){

if(dst.block(0)(*index_iter)!=0 )

block_A->set(*index_iter,i, dst.block(0)(*index_iter) );

}

}

    block_A->compress(VectorOperation::unknown);

Without constrains the matrix matches the matrix-free operator, but with constrains present it does not. What is the proper way to assemble the matrix with vmult?

Best,

Michał Wichrowski

[peterrum]

Hi Michal,

any chance that you post or send me a small runnable test program.

By the way, there is an open PR (https://github.com/dealii/dealii/pull/9343) for the computation of the diagonal in a matrix-free manner. Once this is merged, I will work the matrix-free assembly of sparse matrices.

Thanks,
Peter

[Michał Wichrowski]

It cames out it is not so easy to provide simplest working example. I've tried to remake step-37, but surprisingly I was not able to reproduce the problem. I'll take a closer look what is causing a problem. 

[Wolfgang Bangerth]

On 3/9/20 8:26 AM, Michał Wichrowski wrote:
> I've got matrix-free multigrid solver for Stokes problem. The main
> bottleneck is solution of coarse problem, so I tried to assemble the
> regular sparse matrix and use direct solver. Since the coarse problem is
> (relatively) small, I used vmults by unit vector to obtain columns of
> the matrix.

This is unrelated to your original question, but worth mentioning
anyway: Conceptually, you are multiplying a unit vector by a sparse
matrix, and what you should get is a vector with nonzero entries in only
those rows for which the matrix has a nonzero entry in the column
corresponding to the unit vector.

This is inefficient, because you will have to do N such matrix-vector
products. But if you know which entries of the matrix are non-zero
(which you do, because you know which DoFs couple with each other), then
you can multiply with a sum of unit vectors, get an output vector, and
you will know which entries of the output vector correspond to the
product with each of these unit vectors. You should be able to get the
number of matrix-vector products down from N to N/100 or N/1000 for
sufficiently large problems -- i.e., *much much faster*.

(In fact, I think that we should at one point add a function that takes
an operator and that computes a matrix representation. If we would also
give that function a SparsityPattern, it could do the algorithm outlined
above automatically. As always, any help would be welcome!)

Best
W.

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

[Michał Wichrowski]

W dniu wtorek, 10 marca 2020 17:33:43 UTC+1 użytkownik Wolfgang Bangerth napisał:

On 3/9/20 8:26 AM, Michał Wichrowski wrote:
> I've got matrix-free multigrid solver for Stokes problem. The main
> bottleneck is solution of coarse problem, so I tried to assemble the
> regular sparse matrix and use direct solver. Since the coarse problem is
> (relatively) small, I used vmults by unit vector to obtain columns of
> the matrix.

This is unrelated to your original question, but worth mentioning
anyway: Conceptually, you are multiplying a unit vector by a sparse
matrix, and what you should get is a vector with nonzero entries in only
those rows for which the matrix has a nonzero entry in the column
corresponding to the unit vector.

This is inefficient, because you will have to do N such matrix-vector
products. But if you know which entries of the matrix are non-zero
(which you do, because you know which DoFs couple with each other), then
you can multiply with a sum of unit vectors, get an output vector, and
you will know which entries of the output vector correspond to the
product with each of these unit vectors. You should be able to get the
number of matrix-vector products down from N to N/100 or N/1000 for
sufficiently large problems -- i.e., *much much faster*.

(In fact, I think that we should at one point add a function that takes
an operator and that computes a matrix representation. If we would also
give that function a SparsityPattern, it could do the algorithm outlined
above automatically. As always, any help would be welcome!)

Yes, I thought about that, but I my case is not worth an effort. I'm doing it to assemble coarse matrix, so the problem is relatively small (N_dof around 500). The assemble time is  unnoticeable comparing to solution time, even in case of small problems. Moreover, I need to apply inverse of some-kind Schur complement matrix (on coarse level only), and for that it would be extremely hard to guess sparsity pattern (it is still sparse). Fortunately, this matrix is even smaller. 

Your strategy (if implemented correctly) would require constant number of matrix-vector multiplications so the complexity will be linear. The number of vmults will depend on the mesh connectivity and type of FE. So the overall complexity will  be linear with the problem size (instead (O(N^2), as I understand from Your email). But indeed, it will be much faster. 

Best,
Michał

[Michał Wichrowski]

It looks like it was hidden bug inside my code. I've wrote the coarse grid direct solver for matrix-free operators. It ignores sparsity pattern and assembles the matrix by doing N vmults. The interface may be extended for optimalization. I think it might be an useful tool. Usage, for step-37:

constructor needs vector partitioning:

  CoarseDirectSolver<double> coarse_solver(owned_dofs_u,relevant_dofs_u );

This  

  coarse_solver.initialize(mg_matrices[0]/*, dof_handler*/);

or

  coarse_solver.initialize(mg_matrices[0], dof_handler);

assembles the marix. In second case the matrix is initilized with sparsity pattern, it could be optimized to reduce the number of vmults.

Michał