Kratos: Porting to AMatrix

[pooyan-dadvand]

Dear @KratosMultiphysics/team-maintainers This issue is a preparation point for the porting to AMatrix library.

The idea is to remove the dependency to ublas library which is not as active as it was (lack of c+11 improvement supports for instance) also aligned with our strategic plan to be independent of the boost. ublas has a wide variety of features which are not currently used in Kratos but causes very large compilation time and not giving the best in class performance due to its complexity and design restrictions. Nevertheless the AMatrix is tailored to have only row based matrix with few used expressions to be lightweight and fast.

In this issue I would discuss some design decisions and eventually would suggest some preparation steps before the full migrations

[pooyan-dadvand]

@msandre I have noticed that you have used the boost indirect array in your implementation (in math_utils.h). Is it crucial or can be avoided? Consider that my intention was to avoid implementing the array view (with elementwise operations) in AMatrix.

[loumalouomega]

My main concern how backward compatible it will be. In https://github.com/KratosMultiphysics/Kratos/blob/master/applications/ContactStructuralMechanicsApplication/custom_linear_solvers/mixedulm_linear_solver.h and specially in https://github.com/KratosMultiphysics/Kratos/blob/master/applications/ContactStructuralMechanicsApplication/custom_utilities/sparse_matrix_multiplication_utility.h I manipulate the compressed_matrix manually for performance reasons.

Do I need to do any special modification before the merge?

[pooyan-dadvand]

In a very first step we will only substitute the dense matrix and vectors. The sparse matrices will be changes in a second step.

[jcotela]

@pooyan-dadvand Maybe I am misunderstanding you, but I am now a bit worried about your comment about the "array view". Shape functions are implemented as Matrix (Num Nodes x Num Gauss Points) in all geometries. I use row views in virtually all fluid elements every time we access the shape functions "at gauss point i".

Also, what is the final decision regarding aliasing? Do we still need the noalias access? (I fully support assuming no-aliasing by default, by the way, the interface is much cleaner).

[RiccardoRossi]

@jcotela do you think that we need the "indirect" access?

i mean, one thing is to have row(A,1) to give back a vector containing the values of row 1 of A (this is definitely needed) a different one is to be able to modify A by modifying the row object. I am not saying that such feature isn't useful...just that i don't see it as vital.

in any case my point of view is that we should do this change in steps:
step1 - try to be as compatible as possible (mantain the syntax, and throw errors when the behaviour is different)
step2 - eventually change the syntax.

with this i mean, that i think we should maintain the "noalias" access, at least as first step. To make an example

   Matrix A;
   A = ZeroMatrix(3,3) //here matrix A has size 3*3 (currently)

   Matrix B;
   noalias(B) = ZeroMatrix(3,3) //here matrix B has size 0*0 (currently, i mean) [this may well be a memory error though...]

we need to be careful in changing this, since it may change the behaviour of valid code.
in my opinion what we should do as step1 is to put a lot of KRATOS_DEBUG checks to catch the behaviour that may be different. For example we could put that aliasing check is always done (in debug mode) so that a lot of common errors could be detected.

Just for the records, my biggest concern right now is with the "resize" function.

the point is that CURRENTLY

  Matrix A(3,3)
  A.resize(2,2) //the top left 2*2 part of the matrix is preserved in the new A
  A.resize(2,2,false) //would discard everything and allocate a new empty matrix of size 2*2

While i also would be in favor of not mantaining such feature (that is, to have the second behaviour without the false) argument, by changing this we would silently kill a lot of valid code.
My proposal would be thus to implement resize as

   resize(size1,size2, preserve = true)
   {
          if(preserve == true) KRATOS_ERROR ...
   }

once everything works, we could then do a second sweep, and get rid of all of the third arguments. This would be definitely boring, but not risky.

[jcotela]

@RiccardoRossi I don't need to be able to modify the values, but I'd like to still be able to pass a row of the matrix as a (read-only) vector. The alternative is either making a copy or leting all underlying functions know about the gauss index (so that Ni = N[g,i]). Neither is impossible, but it would be a non-trivial amount of work and either a performance penalty or a degradation of the interfaces.

I don't have a strong opinion regarding the resize, so I'll just follow what is decided on that point.

[msandre]

@pooyan-dadvand I don't see it as a requirement for the functions in math_utils.h. But we will need to work with copies in that case. For performance this is a bigger issue for the elements that use such features as @jcotela mentioned. There are tests so we should know when these stop working and we can add the missing operations then. I wasn't aware of the AMatrix lib. It would be great if it could support std array and vector and initializer lists.

[pooyan-dadvand]

Now I see the situation, let me check how I would add this to AMatrix

Edit: typo

[philbucher]

@pooyan-dadvand do you still plan to change to AMatrix for the next release?
I am asking because I think this will be quite some effort for the developers to adapt to the changes

Are you planning to do it in the same way as with pybind, i.e. in a separate branch where we commit the ported code to?

[pooyan-dadvand]

I don't think that we can manage to have this change for release 6 even with a major effort and more I dig into it I see more problems to be handled. As you said the ublas is much more intrusive in our code than any other library rather than STL.

The plan would be more or less the same as the pybind but at least I should achieve a full core compilation (in my local branch) before starting any serious migration.

Still there are several features are missing in AMatrix (LU solver, range, slice, row, column, etc.) to be used in Kratos. Saying that, any help in improving the AMatrix is more than welcomed! :sos:

[RiccardoRossi]

@pooyan-dadvand just to tell that i don't think you should support slicing. You do need to support "rows" and "cols" operator, however that's a easier one.

[pooyan-dadvand]

Now the PR #1931 has been merged and I would suggest all @KratosMultiphysics/team-maintainers to proceed with changes suggested on that PR for their applications. I would try to compile the Kratos in daily bases to see the evolution of the migration and keep you informed here in this issue.

[pooyan-dadvand]

I have added the PR #1939 which adds the amatrix interface and migration macro.

[pooyan-dadvand]

I have merged the PR #1939 but still far from having the core compiled....

[roigcarlo]

Won't be ready for 6.0.0 so I am removing it from the milestone

[philbucher]

@pooyan-dadvand I asked a while ago if you will introduce an alias option (the opposite of noalias) in case there is aliasing.
It seems we need quite some auxiliar objects if we don't have it, see in #2923 => in the svd-utils

If you plan to add it then I wait with porting StructuralMechanics

FYI @jcotela @loumalouomega I think you had some problems with this

[pooyan-dadvand]

alias option would be to copy it in temp matrix.

We can define it in Kratos side:

template <typename TMatrixType>
class AliasExpression
    : public AMatrix::MatrixExpression<AliasExpression<TMatrixType>, AMatrix::row_major_access> {
    TMatrixType& _original_matrix;
   public:
    AliasExpression() = delete;

    AliasExpression(TMatrixType& OriginalMatrix) : _original_matrix(OriginalMatrix){}

    template<TExpressionType>
    inline TMatrixType& operator=(TExpressionType OtherExpression) { 
        TMatrixType temp(OtherExpression)
        _original_expression = temp; 
        return _original_expression; 
    }

    bool check_aliasing(const data_type* From, const data_type* To) const {
        return false;
    }
};

However, IMO using an explicit temp matrix is more clear.

[philbucher]

thanks @pooyan-dadvand , then I will work with temp objects

[pooyan-dadvand]

temp objects or the assign operators like += and -= if they are available...