#1669 Release of dune 3.0

Added discussion task and removed bug labels

mentioned in issue core/dune-common#11 (closed)

@christoph.gersbacher: Does 21) refer to all three versions of istl_assign_to_fmatrix? As I've laid out in core/dune-common#11 (closed), the C-style matrix version is broken. There are also versions for the ScaledIdentityMatrix and DiagonalMatrix classes, however, without which neither

  Dune::FieldMatrix<double, 3, 3> const m1 =
      Dune::DiagonalMatrix<double, 3>(Dune::FieldVector<double, 3>{1, 2, 3});

nor

  Dune::FieldMatrix<double, 3, 3> const m2 =
    Dune::ScaledIdentityMatrix<double, 3>(7);

would run. Is this supposed to be accounted for somehow?

Hi Elias. Yes, we intended to get rid of the method istl_assign_to_fmatrix everywhere. The idea was to replace each definition of this function by a template specialization of

template<class DenseMatrix, class RHS>
struct DenseMatrixAssigner;

which is declared in dune/common/densematrix.hh. I guess the problems you described are related to the default implementation of this class: in case no template specialization is found the class tries to call the legacy function istl_assign_to_fmatrix (or as a last resort a cast operator). In an ideal world, the solution would be

  template<class DenseMatrix, class Field, int N>
  struct DenseMatrixAssigner<DenseMatrix, DiagonalMatrix<Field, N>>
  {
    static void apply (DenseMatrix &denseMatrix, const DiagonalMatrix<Field, N>& rhs)
    {
      // body of istl_assign_to_fmatrix goes here
    }
  };

mentioned in merge request core/dune-istl!23 (closed)

mentioned in merge request core/dune-common!31 (closed)

Alright. That doesn't seem all too difficult. So I've carried out a few changes while making sure that

  Dune::FieldMatrix<double, 3, 3> const m1 =
      Dune::DiagonalMatrix<double, 3>(Dune::FieldVector<double, 3>{1, 2, 3});
  Dune::FieldMatrix<double, 3, 3> const m2 =
      Dune::ScaledIdentityMatrix<double, 3>(7);
  Dune::FieldMatrix<double, 3, 3> const m3 =
      Dune::IdentityMatrix<double, 3>();

continues to work as expected. These currently work because

There's istl_assign_to_fmatrix for the DiagonalMatrix class
There's istl_assign_to_fmatrix for the ScaledIdentityMatrix class
IdentityMatrix provides a conversion operator to a FieldMatrix (but not other DenseMatrix implentations).

After my changes, they work because DenseMatrixAssigner has been specialised for each (whether this is the right thing to do to the IdentityMatrix class might deserve a debate).

This obsoletes istl_assign_to_fmatrix completely, so that I've removed the corresponding logic from DenseMatrixAssignerImplementation as well.

Please see core/dune-common!31 (closed) and core/dune-istl!23 (closed).

I am tempted to start lobbying for the removal of IdentityMatrix again -- I still don't see its point.

@oliver.sander I see your point. Nothing is specialised to handle instances of the IdentityMatrix class so that in almost all cases, it will be converted to a FieldMatrix, and just like that the saving in arithmetic operations is gone...

@pipping No, my point is: if you really ever do need a super-efficient IdentityMatrix (I can't think of such a case, but hey...), then ScaledIdentityMatrix is at your service. No need for a separate class.

@oliver.sander Well, yes, I expected you to be hinting at ScaledIdentityMatrix. Which is nearly as space-efficient (I guess the single scalar typically isn't worth discussing). I was thinking more about computational efficiency beyond .mv(), though -- things like DenseMatrix::rightmultiply. For an IdentityMatrix that could be turned into no-op and for a ScaledIdentityMatrix into x *= y.scalar(), neither of which is currently happening, though, no? (indeed, I'm not even sure if a ScaledIdentityMatrix is accepted as an argument.

Speaking of which: Should such a function really be a method (of DenseMatrix)?

I expect ScaledIdMatrix to be as efficient as IdentityMatrix in all regards, if you add a sufficient numbers of constepxrs in the right places.

To begin with, I don't care much about the fate of the identity matrix. I would like to remark that there are linear mappings I'd rather not store as dense matrix objects. As Elias already said the the application of the identity operator to a vector has a complexity of O(1) against O(n) in case of multiplication with the scaled identity matrix. A while ago we agreed on a condensed interface for 'constant dense matrices'. These matrices provide the usual mappings but do not grant access to rows or individual matrix entries. If needed these objects must be cast into a dense matrix first. The identity matrix is an example of such an object; I use a number of similar objects in my applications, e.g., for outer products or rotation matrices, and I find them extremely helpful. I think that DUNE's diagonal matrix or the scaled identity matrix suffer from being implemented as dense matrices. They go at great length to provide the [i][j]-operator for accessing matrix entries only to throw an exception in case of i != j (I may however call the non-interface method exists(i,j)).

Please remember that we never agreed on introducing a 'const dense matrix' interface. This is also reflected in the fact that we don't require individual entry access. For some reason someone added the term 'dense' afterwards.

Exactly. Still it's the name of the interface class that I forgot to mention can be found in the feature/ConstDenseMatrix branch of dune-common.

Christoph, nobody says that having special implementations for special types of matrices is not useful. But I do not see the use for a hard-coded identity.

First of all, I think this discussion on the identity matrix should long have gone into a separate issue.

But to the discussion: Funnily enough, it's the ScaledIdentityMatrix that is counterintuitive to me. At least to me, it makes perfect sense to scale an arbitrary matrix (e.g., a rotation) by a scalar. Consequently, I'd prefer a ScaledMatrix< IdentityMatrix > over a ScaledIdentityMatrix. We could even specialize the operator* to return a ScaledMatrix to have an expression template.

Moreover, I think this discussion is less about the extra 'double' but rather about possible template specialization. A ScaledIdentityMatrix already hides information at run time that might be known at compile time. But in the end, it might just be a matter of taste.

In the end, @oliver is right. We should not have both, a ScaledIdentityMatrix and an IdentityMatrix. Maybe we should simply vote on:

remove IdentityMatrix in favor of the ScaledIdentityMatrix.
remove ScaledIdentityMatrix in favor of ScaledMatrix< IdentityMatrix >.

Since rotation matrices have been brought up repeatedly now: Is there already such a class implementation?

Here's how these specialised matrices appear to waste computational efficiency:

#include <config.h>

#include <dune/common/diagonalmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/identitymatrix.hh>
#include <dune/istl/scaledidmatrix.hh>

int main() {
  size_t constexpr m = 2;
  size_t constexpr n = 3;
  Dune::FieldMatrix<double, m, n> fieldMatrix = {{1, 2, 3}, {4, 5, 6}};
  Dune::IdentityMatrix<double, n> identityMatrix;
  Dune::ScaledIdentityMatrix<double, n> scaledIdentityMatrix(2);
  Dune::DiagonalMatrix<double, n> diagonalMatrix{1, 2, 3};

  // Q: Can we expect any of the following to be optimised so that                     
  // zero-entries do not lead to a multiplication so that the order is                 
  // effectively lower?                                                                

  /// IdentityMatrix                                                                   
  // Fails to compile                                                                  
  //fieldMatrix.rightmultiply(identityMatrix);                                         
  // Compiles; O(m*n^2) instead of O(0);                                               
  fieldMatrix.rightmultiply(Dune::FieldMatrix<double, n, n>(identityMatrix));

  /// ScaledIdentityMatrix                                                             
  // Fails to compile                                                                  
  //fieldMatrix.rightmultiply(scaledIdentityMatrix);                                   
  // Compiles. O(m*n^2) instead of O(m*n)                                              
  fieldMatrix.rightmultiply(Dune::FieldMatrix<double, n, n>(scaledIdentityMatrix));

  /// DiagonalMatrix                                                                   
  // Fails to compile                                                                  
  //fieldMatrix.rightmultiply(diagonalMatrix);                                         
  // Compiles. O(m*n^2) instead of O(m*n)                                              
  fieldMatrix.rightmultiply(Dune::FieldMatrix<double, n, n>(diagonalMatrix));
}

Of course efficiency is wasted if you convert to FieldMatrix first. Are you implying that these conversions are done implicitly?

While these compilation failures are certainly bugs, I am not really shocked by them: I don't think I have ever used rightmultiply anywhere. mv and umv is much more important (to me).

Also note that you may have to mark certain methods (in particular, the constructors) of the Matrix classes as constexpr to make the whole constexpr mechanism really work.

No, they don't happen implicitly (even though I would've expected this for the IdentityMatrix class because it has a conversion operator to FieldMatrix). The sense in which efficiency is wasted is: If I want to be able to write

  SomeMatrix m1; // 3 x 3
  Dune::FieldMatrix<double, 2, 3> m2;

  m2.rightmultiply(m1);

then this will generally not work. I need to write

  SomeMatrix m1; // 3 x 3
  Dune::FieldMatrix<double, 2, 3> m2;

  Dune::FieldMatrix<double, 3, 3> fatM1(m1);
  m2.rightmultiply(fatM1);

and thereby waste efficiency.

I would say this is a bug.

Implementing those methods efficiently for ScaledIdentityMatrix and DiagonalMatrix should be easy: One just has to use iterators instead of random access operators.

@martin.nolte: This vote is not an option unless we have a fully functional implementation of IdentityMatrix. This would require to have iterators. Otherwise users can only do those operations that we have foreseen.

BTW: I propose to not do any fundamental quick changes to our matrices unless we have carefully rethought the design which currently has so many flaws. And @martin.nolte is completely right: This should go into a separate issue.

mentioned in issue core/dune-istl#7

I've now opened core/dune-istl#7 for this purpose. I'll keep the summary at the top updated. Please yell and make yourself heard!

mentioned in commit pipping/dune-common@f2e3a98c

mentioned in commit pipping/dune-common@d255f7ef

I believe that in contrast to the plans in core/dune-istl#7 which have far-reaching consequences, core/dune-common!31 (closed) is rather uncontroversial. I'm trying to highlight it here again so that it doesn't get tossed aside. Don't let this opportunity to kill istl_assign_to_fmatrix once and for all go to waste ;)

mentioned in commit core/dune-istl@49be9867

mentioned in commit core/dune-common@1d7528ea

mentioned in commit core/dune-common@9e176a0d

mentioned in commit pipping/dune-common@1d7528ea

mentioned in commit pipping/dune-common@9e176a0d

The deprecation in 11) has been implemented, 21) has been fixed by @pipping 'a patches. 24) is considered in a separate issue (#1435 (closed)) and 17) is now covered by core/dune-common#26 (closed). Since there won't be a 3.0 soon and the mentioned issues either have been done or are covered elsewhere I'll close this.

Status changed to closed

Mentioned in commit core/dune-common@1d7528ea

Mentioned in commit core/dune-common@9e176a0d

Property	Value
Reported by	Carsten Gräser (graeser@math.fu-berlin.de)
Reported at	Jun 16, 2015 17:50
Type	Bug Report
Version	2.3
Operating System	Unspecified / All
Last edited by	Carsten Gräser (graeser@math.fu-berlin.de)
Last edited at	Jun 16, 2015 17:52

#1669 Release of dune 3.0

Metadata

Description

Designs

Child items ...

Activity