dolfin team mailing list archive

[Matthew Knepley <knepley@xxxxxxxxxxx>]

Anders Logg <logg@xxxxxxxxx> writes:
> I will make a comparison with PETSc (MatSetValuesBlocked(A,...,AK,...))
> so we can see the exact difference. I think it is ok to use this type
> of call within the assembler, but I would like to be able to use a
> clean and simple OO interface elsewhere. (Looks like this is going to
> be added to PETSc.)
    Again, these are just wrappers. They are trivial (except for maintenance).
> 
> We currently add the local element matrix AK to the global matrix A by doing
> 
> for (unsigned int i; ...)
>     for (unsigned int j; ...)
>         A(dof(i), dof(j)) += AK(i, j);
> 
> where dof() is the mapping from local to global degrees of freedom. Each
> time we access the global matrix A on row dof(i), it has to search for
> a nonzero entry in column dof(j) on that row so this is probably slow.
> How do you avoid doing this search in PETSc? Or don't you?
    1) Rob is correct that this interface is inherently limiting. No compiler
       in our lifetime will emit good code for it. That is why some systems
       overload += to allow a matrix operation here (somewhat like completely
       unrolling the loops, but more involved since we have slicing).

    2) For sparse matrices, that kind of general slicing is slow. Therefore,
       you want more assurances about the input from your user.

    3) In PETSc, we have several levels of optimization here. First, the call
       accepts logically dense blocks of values, meaning it is like overloading
       +=. Second, you can indicate that dof(j) is already sorted, which saves
       a lot of work. In the code above, that would not be true. Third, we
       often can use block storage, which means only the first column of a
       block is looked up. Forth, the whole routine is inlined.

       Matt
> 
> /Anders
-- 
"Failure has a thousand explanations. Success doesn't need one" -- Sir Alec Guiness