[eigen] Google Summer of Code 2018 - Symmetric Matrices for Eigen

[David A. Tellenbach]

Hello together,

my name is David Tellenbach, I'm currently studying Computer Science at the LMU Munich, Germany and got chosen for the Google Summer of Code 2018 Project "Faster Matrix Algebra for ATLAS", supervised by Dmitry Emeliyanov and Stewart Martin-Haugh.

Google Summer of Code is a global program focused on bringing more student developers into open source software development. Students work with an open source organization on a three month programming project during their break from school.

As you might know, our project's task is to implement support for symmetric matrices for Eigen. A short project description is available via the following link: https://summerofcode.withgoogle.com/projects/#5293950017994752 <https://summerofcode.withgoogle.com/projects/#5293950017994752>

The official coding period hasn't started yet and lasts for three month, from May, 14 until August, 14 2018. The time now is meant to get to know the community and people involved.

As far as I see, Eigen basically provides three types of matrices: Dense, sparse and diagonal matrices. Of these types, the class Eigen::DiagonalMatrix seems to be the one that could be most similar to a possible implementation of a class for symmetric matrices. There is no need for storing all elements (as in the case of dense matrices) neither is a sophisticated mechanism to find the position of scalars in the matrix needed (as in the case of sparse matrices). Therefore I’d like to create the for symmetric matrices by deriving from Eigen::EigenBase (as in the case of diagonal matrices).

Of course one goal is to store only the upper or lower triangular part of the matrix since this already defines it completely. Similar to Eigen::DiagonalMatrix the storage could look something like this:

typedef typename internal::traits<Derived>::SymmetricVectorType SymmetricVectorType;

// Store just one triangular part of the matrix
typedef Matrix<_Scalar, (SizeAtCompileTime * SizeAtCompileTime + SizeAtCompileTime)/2,
1, 0, (MaxSizeAtCompileTime * MaxSizeAtCompileTime + MaxSizeAtCompileTime)/2 ,1> SymmetricVectorType;

We plan to provide constructors which take either matrices of type Eigen::Matrix<
> or Eigen::SelfAdjointView<
>.

What do you think about these broad plans so far? We are happy about any feedback.

Thanks,
David

[Rasmus Munk Larsen]

Hi David,

I support you idea, but would like to suggest that you call the class
HermitianMatrix (or SelfAdjointMatrix) to properly support the extension of
the concept to complex-valued matrices. I will warn you that in my
experience, you may save half the memory, but the speedup may be far less
(or even less than 1 for small to medium sized matrices) due to the less
regular memory access patterns.

Best regards,
Rasmus

Post by David A. Tellenbach
Hello together,
my name is David Tellenbach, I'm currently studying Computer Science at
the LMU Munich, Germany and got chosen for the Google Summer of Code 2018
Project "Faster Matrix Algebra for ATLAS", supervised by Dmitry Emeliyanov
and Stewart Martin-Haugh.
Google Summer of Code is a global program focused on bringing more student
developers into open source software development. Students work with an
open source organization on a three month programming project during their
break from school.
As you might know, our project's task is to implement support for
symmetric matrices for Eigen. A short project description is available via
the following link: https://summerofcode.withgoogle.com/projects/#
5293950017994752
The official coding period hasn't started yet and lasts for three month,
from May, 14 until August, 14 2018. The time now is meant to get to know
the community and people involved.
As far as I see, Eigen basically provides three types of matrices: Dense,
sparse and diagonal matrices. Of these types, the class
Eigen::DiagonalMatrix seems to be the one that could be most similar to a
possible implementation of a class for symmetric matrices. There is no need
for storing all elements (as in the case of dense matrices) neither is a
sophisticated mechanism to find the position of scalars in the matrix
needed (as in the case of sparse matrices). Therefore I’d like to create
the for symmetric matrices by deriving from Eigen::EigenBase (as in the
case of diagonal matrices).
Of course one goal is to store only the upper or lower triangular part of
the matrix since this already defines it completely. Similar to
typedef typename internal::traits<Derived>::SymmetricVectorType SymmetricVectorType;
// Store just one triangular part of the matrix
typedef Matrix<_Scalar, (SizeAtCompileTime * SizeAtCompileTime + SizeAtCompileTime)/2,
1, 0, (MaxSizeAtCompileTime * MaxSizeAtCompileTime +
MaxSizeAtCompileTime)/2 ,1> SymmetricVectorType;
We plan to provide constructors which take either matrices of type
Eigen::Matrix<
> or Eigen::SelfAdjointView<
>.
What do you think about these broad plans so far? We are happy about any feedback.
Thanks,
David

[Márton Danóczy]

Hi David,

there is very interesting prior art on this: The LAPACK Rectangular Full
Packed format (see http://www.netlib.org/lapack/lawnspdf/lawn199.pdf) is
able to take advantage of GEMM when multiplying Hermitian matrices. This is
a great improvement over the conventional packed format (which you plan to
implement as far as I can tell), which does not support Level-3 BLAS. You
should definitely look at the paper for inspiration.

Best
Márton

Post by David A. Tellenbach
Hello together,
my name is David Tellenbach, I'm currently studying Computer Science at
the LMU Munich, Germany and got chosen for the Google Summer of Code 2018
Project "Faster Matrix Algebra for ATLAS", supervised by Dmitry Emeliyanov
and Stewart Martin-Haugh.
Google Summer of Code is a global program focused on bringing more student
developers into open source software development. Students work with an
open source organization on a three month programming project during their
break from school.
As you might know, our project's task is to implement support for
symmetric matrices for Eigen. A short project description is available via
the following link: https://summerofcode.withgoogle.com/projects/#
5293950017994752
The official coding period hasn't started yet and lasts for three month,
from May, 14 until August, 14 2018. The time now is meant to get to know
the community and people involved.
As far as I see, Eigen basically provides three types of matrices: Dense,
sparse and diagonal matrices. Of these types, the class
Eigen::DiagonalMatrix seems to be the one that could be most similar to a
possible implementation of a class for symmetric matrices. There is no need
for storing all elements (as in the case of dense matrices) neither is a
sophisticated mechanism to find the position of scalars in the matrix
needed (as in the case of sparse matrices). Therefore I’d like to create
the for symmetric matrices by deriving from Eigen::EigenBase (as in the
case of diagonal matrices).
Of course one goal is to store only the upper or lower triangular part of
the matrix since this already defines it completely. Similar to
typedef typename internal::traits<Derived>::SymmetricVectorType SymmetricVectorType;
// Store just one triangular part of the matrix
typedef Matrix<_Scalar, (SizeAtCompileTime * SizeAtCompileTime + SizeAtCompileTime)/2,
1, 0, (MaxSizeAtCompileTime * MaxSizeAtCompileTime +
MaxSizeAtCompileTime)/2 ,1> SymmetricVectorType;
We plan to provide constructors which take either matrices of type
Eigen::Matrix<
> or Eigen::SelfAdjointView<
>.
What do you think about these broad plans so far? We are happy about any feedback.
Thanks,
David

[David A. Tellenbach]

Hi Márton,

thanks for the suggestion. I just had a short look at it but the idea of supporting Level-3 BLAS by design sounds promising.

I'll keep you updated.

Thanks,
David

Post by Rasmus Munk Larsen
Hi David,
there is very interesting prior art on this: The LAPACK Rectangular Full Packed format (see http://www.netlib.org/lapack/lawnspdf/lawn199.pdf <http://www.netlib.org/lapack/lawnspdf/lawn199.pdf>) is able to take advantage of GEMM when multiplying Hermitian matrices. This is a great improvement over the conventional packed format (which you plan to implement as far as I can tell), which does not support Level-3 BLAS. You should definitely look at the paper for inspiration.
Best
Márton
Hello together,
my name is David Tellenbach, I'm currently studying Computer Science at the LMU Munich, Germany and got chosen for the Google Summer of Code 2018 Project "Faster Matrix Algebra for ATLAS", supervised by Dmitry Emeliyanov and Stewart Martin-Haugh.
Google Summer of Code is a global program focused on bringing more student developers into open source software development. Students work with an open source organization on a three month programming project during their break from school.
As you might know, our project's task is to implement support for symmetric matrices for Eigen. A short project description is available via the following link: https://summerofcode.withgoogle.com/projects/#5293950017994752 <https://summerofcode.withgoogle.com/projects/#5293950017994752>
The official coding period hasn't started yet and lasts for three month, from May, 14 until August, 14 2018. The time now is meant to get to know the community and people involved.
As far as I see, Eigen basically provides three types of matrices: Dense, sparse and diagonal matrices. Of these types, the class Eigen::DiagonalMatrix seems to be the one that could be most similar to a possible implementation of a class for symmetric matrices. There is no need for storing all elements (as in the case of dense matrices) neither is a sophisticated mechanism to find the position of scalars in the matrix needed (as in the case of sparse matrices). Therefore I’d like to create the for symmetric matrices by deriving from Eigen::EigenBase (as in the case of diagonal matrices).
typedef typename internal::traits<Derived>::SymmetricVectorType SymmetricVectorType;
// Store just one triangular part of the matrix
typedef Matrix<_Scalar, (SizeAtCompileTime * SizeAtCompileTime + SizeAtCompileTime)/2,
1, 0, (MaxSizeAtCompileTime * MaxSizeAtCompileTime + MaxSizeAtCompileTime)/2 ,1> SymmetricVectorType;
We plan to provide constructors which take either matrices of type Eigen::Matrix<
> or Eigen::SelfAdjointView<
>.
What do you think about these broad plans so far? We are happy about any feedback.
Thanks,
David

[Christoph Hertzberg]

Hi,

mostly repeating what I told David privately already:

I suggest implementing a compact _triangular_ matrix which can be used
together with SelfAdjointView to represent selfadjoint matrices as well.

And then I'd like to point anyone interested in this discussion to this
(very old) feature request:
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=42
This also mentions the RFPF suggested by Márton -- but we probably
should support both and start with the one which is easier to implement.


Cheers,
Christoph

Post by David A. Tellenbach
Hello together,
my name is David Tellenbach, I'm currently studying Computer Science at the LMU Munich, Germany and got chosen for the Google Summer of Code 2018 Project "Faster Matrix Algebra for ATLAS", supervised by Dmitry Emeliyanov and Stewart Martin-Haugh.
Google Summer of Code is a global program focused on bringing more student developers into open source software development. Students work with an open source organization on a three month programming project during their break from school.
As you might know, our project's task is to implement support for symmetric matrices for Eigen. A short project description is available via the following link: https://summerofcode.withgoogle.com/projects/#5293950017994752 <https://summerofcode.withgoogle.com/projects/#5293950017994752>
The official coding period hasn't started yet and lasts for three month, from May, 14 until August, 14 2018. The time now is meant to get to know the community and people involved.
As far as I see, Eigen basically provides three types of matrices: Dense, sparse and diagonal matrices. Of these types, the class Eigen::DiagonalMatrix seems to be the one that could be most similar to a possible implementation of a class for symmetric matrices. There is no need for storing all elements (as in the case of dense matrices) neither is a sophisticated mechanism to find the position of scalars in the matrix needed (as in the case of sparse matrices). Therefore I’d like to create the for symmetric matrices by deriving from Eigen::EigenBase (as in the case of diagonal matrices).
typedef typename internal::traits<Derived>::SymmetricVectorType SymmetricVectorType;
// Store just one triangular part of the matrix
typedef Matrix<_Scalar, (SizeAtCompileTime * SizeAtCompileTime + SizeAtCompileTime)/2,
1, 0, (MaxSizeAtCompileTime * MaxSizeAtCompileTime + MaxSizeAtCompileTime)/2 ,1> SymmetricVectorType;
We plan to provide constructors which take either matrices of type Eigen::Matrix<…> or Eigen::SelfAdjointView<…>.
What do you think about these broad plans so far? We are happy about any feedback.
Thanks,
David

--
Dr.-Ing. Christoph Hertzberg

Besuchsadresse der Nebengeschäftsstelle:
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 5
28359 Bremen, Germany

Postadresse der Hauptgeschäftsstelle Standort Bremen:
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 1
28359 Bremen, Germany

Tel.: [phone]
Zentrale: [phone]
E-Mail: ***@dfki.de

Weitere Informationen: http://www.dfki.de/robotik
-----------------------------------------------------------------------
Deutsches Forschungszentrum fuer Kuenstliche Intelligenz GmbH
Firmensitz: Trippstadter Straße 122, D-67663 Kaiserslautern
Geschaeftsfuehrung: Prof. Dr. Dr. h.c. mult. Wolfgang Wahlster
(Vorsitzender) Dr. Walter Olthoff
Vorsitzender des Aufsichtsrats: Prof. Dr. h.c. Hans A. Aukes
Amtsgericht Kaiserslautern, HRB 2313
Sitz der Gesellschaft: Kaiserslautern (HRB 2313)
USt-Id.Nr.: DE 148646973
Steuernummer: 19/672/50006
-----------------------------------------------------------------------

[Márton Danóczy]

Hi Christoph,

supporting the "Packed Format" is not worth it, even LAPACK have dropped
support because of the FP routines' subpar performance (see the second
comment in this thread:
https://github.com/Reference-LAPACK/lapack/issues/245). The RFP format is
and will keep being supported by LAPACK (same Github thread).

Another reference:
https://software.intel.com/en-us/mkl-developer-reference-c-matrix-storage-schemes-for-lapack-routines

Best
Márton

Post by Christoph Hertzberg
Hi,
I suggest implementing a compact _triangular_ matrix which can be used
together with SelfAdjointView to represent selfadjoint matrices as well.
And then I'd like to point anyone interested in this discussion to this
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=42
This also mentions the RFPF suggested by Márton -- but we probably should
support both and start with the one which is easier to implement.
Cheers,
Christoph

Post by David A. Tellenbach
Hello together,
my name is David Tellenbach, I'm currently studying Computer Science at
the LMU Munich, Germany and got chosen for the Google Summer of Code 2018
Project "Faster Matrix Algebra for ATLAS", supervised by Dmitry Emeliyanov
and Stewart Martin-Haugh.
Google Summer of Code is a global program focused on bringing more
student developers into open source software development. Students work
with an open source organization on a three month programming project
during their break from school.
As you might know, our project's task is to implement support for
symmetric matrices for Eigen. A short project description is available via
the following link: https://summerofcode.withgoogl
e.com/projects/#5293950017994752 <https://summerofcode.withgoog
le.com/projects/#5293950017994752>
The official coding period hasn't started yet and lasts for three month,
from May, 14 until August, 14 2018. The time now is meant to get to know
the community and people involved.
As far as I see, Eigen basically provides three types of matrices: Dense,
sparse and diagonal matrices. Of these types, the class
Eigen::DiagonalMatrix seems to be the one that could be most similar to a
possible implementation of a class for symmetric matrices. There is no need
for storing all elements (as in the case of dense matrices) neither is a
sophisticated mechanism to find the position of scalars in the matrix
needed (as in the case of sparse matrices). Therefore I’d like to create
the for symmetric matrices by deriving from Eigen::EigenBase (as in the
case of diagonal matrices).
Of course one goal is to store only the upper or lower triangular part of
the matrix since this already defines it completely. Similar to
typedef typename internal::traits<Derived>::SymmetricVectorType SymmetricVectorType;
// Store just one triangular part of the matrix
typedef Matrix<_Scalar, (SizeAtCompileTime * SizeAtCompileTime + SizeAtCompileTime)/2,
1, 0, (MaxSizeAtCompileTime * MaxSizeAtCompileTime +
MaxSizeAtCompileTime)/2 ,1> SymmetricVectorType;
We plan to provide constructors which take either matrices of type
Eigen::Matrix<
> or Eigen::SelfAdjointView<
>.
What do you think about these broad plans so far? We are happy about any feedback.
Thanks,
David

--
Dr.-Ing. Christoph Hertzberg
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 5
28359 Bremen, Germany
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 1
28359 Bremen, Germany
Tel.: [phone]
Zentrale: [phone]
Weitere Informationen: http://www.dfki.de/robotik
-----------------------------------------------------------------------
Deutsches Forschungszentrum fuer Kuenstliche Intellig
<https://maps.google.com/?q=ches+Forschungszentrum+fuer+Kuenstliche+Intellig&entry=gmail&source=g>enz
GmbH
Firmensitz: Trippstadter Straße 122, D-67663 Kaiserslautern
Geschaeftsfuehrung: Prof. Dr. Dr. h.c. mult. Wolfgang Wahlster
(Vorsitzender) Dr. Walter Olthoff
Vorsitzender des Aufsichtsrats: Prof. Dr. h.c. Hans A. Aukes
Amtsgericht Kaiserslautern, HRB 2313
Sitz der Gesellschaft: Kaiserslautern (HRB 2313)
USt-Id.Nr.: DE 148646973
Steuernummer: 19/672/50006
-----------------------------------------------------------------------

[Christoph Hertzberg]

I'd leave the choice what to implement to David. Of course, before this
will be merged, there should be some benchmarking to make sure it does
not totally kill performance (and it would be wise to start benchmarking
early on). Slightly subpar performance may be ok in certain situations
(e.g., if someone is using this format for storage or data transfer
already and can't change that). It may also depend on the size of the
matrices which format is better (and on the used CPU architecture). I
doubt that RFPF is superior to simple packed format for 4x4 matrices (I
did not benchmark this).

Cheers,
Christoph

Post by Márton Danóczy
Hi Christoph,
supporting the "Packed Format" is not worth it, even LAPACK have dropped
support because of the FP routines' subpar performance (see the second
https://github.com/Reference-LAPACK/lapack/issues/245). The RFP format is
and will keep being supported by LAPACK (same Github thread).
https://software.intel.com/en-us/mkl-developer-reference-c-matrix-storage-schemes-for-lapack-routines
Best
Márton

Post by Christoph Hertzberg
Hi,
I suggest implementing a compact _triangular_ matrix which can be used
together with SelfAdjointView to represent selfadjoint matrices as well.
And then I'd like to point anyone interested in this discussion to this
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=42
This also mentions the RFPF suggested by Márton -- but we probably should
support both and start with the one which is easier to implement.
Cheers,
Christoph

Post by David A. Tellenbach
Hello together,
my name is David Tellenbach, I'm currently studying Computer Science at
the LMU Munich, Germany and got chosen for the Google Summer of Code 2018
Project "Faster Matrix Algebra for ATLAS", supervised by Dmitry Emeliyanov
and Stewart Martin-Haugh.
Google Summer of Code is a global program focused on bringing more
student developers into open source software development. Students work
with an open source organization on a three month programming project
during their break from school.
As you might know, our project's task is to implement support for
symmetric matrices for Eigen. A short project description is available via
the following link: https://summerofcode.withgoogl
e.com/projects/#5293950017994752 <https://summerofcode.withgoog
le.com/projects/#5293950017994752>
The official coding period hasn't started yet and lasts for three month,
from May, 14 until August, 14 2018. The time now is meant to get to know
the community and people involved.
As far as I see, Eigen basically provides three types of matrices: Dense,
sparse and diagonal matrices. Of these types, the class
Eigen::DiagonalMatrix seems to be the one that could be most similar to a
possible implementation of a class for symmetric matrices. There is no need
for storing all elements (as in the case of dense matrices) neither is a
sophisticated mechanism to find the position of scalars in the matrix
needed (as in the case of sparse matrices). Therefore I’d like to create
the for symmetric matrices by deriving from Eigen::EigenBase (as in the
case of diagonal matrices).
Of course one goal is to store only the upper or lower triangular part of
the matrix since this already defines it completely. Similar to
typedef typename internal::traits<Derived>::SymmetricVectorType
SymmetricVectorType;
// Store just one triangular part of the matrix
typedef Matrix<_Scalar, (SizeAtCompileTime * SizeAtCompileTime +
SizeAtCompileTime)/2,
1, 0, (MaxSizeAtCompileTime * MaxSizeAtCompileTime +
MaxSizeAtCompileTime)/2 ,1> SymmetricVectorType;
We plan to provide constructors which take either matrices of type
Eigen::Matrix<…> or Eigen::SelfAdjointView<…>.
What do you think about these broad plans so far? We are happy about any feedback.
Thanks,
David

--
Dr.-Ing. Christoph Hertzberg
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 5
28359 Bremen, Germany
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 1
28359 Bremen, Germany
Tel.: [phone]
Zentrale: [phone]
Weitere Informationen: http://www.dfki.de/robotik
-----------------------------------------------------------------------
Deutsches Forschungszentrum fuer Kuenstliche Intellig
<https://maps.google.com/?q=ches+Forschungszentrum+fuer+Kuenstliche+Intellig&entry=gmail&source=g>enz
GmbH
Firmensitz: Trippstadter Straße 122, D-67663 Kaiserslautern
Geschaeftsfuehrung: Prof. Dr. Dr. h.c. mult. Wolfgang Wahlster
(Vorsitzender) Dr. Walter Olthoff
Vorsitzender des Aufsichtsrats: Prof. Dr. h.c. Hans A. Aukes
Amtsgericht Kaiserslautern, HRB 2313
Sitz der Gesellschaft: Kaiserslautern (HRB 2313)
USt-Id.Nr.: DE 148646973
Steuernummer: 19/672/50006
-----------------------------------------------------------------------

--
Dr.-Ing. Christoph Hertzberg

Besuchsadresse der Nebengeschäftsstelle:
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 5
28359 Bremen, Germany

Postadresse der Hauptgeschäftsstelle Standort Bremen:
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 1
28359 Bremen, Germany

Tel.: [phone]
Zentrale: [phone]
E-Mail: ***@dfki.de

Weitere Informationen: http://www.dfki.de/robotik
-----------------------------------------------------------------------
Deutsches Forschungszentrum fuer Kuenstliche Intelligenz GmbH
Firmensitz: Trippstadter Straße 122, D-67663 Kaiserslautern
Geschaeftsfuehrung: Prof. Dr. Dr. h.c. mult. Wolfgang Wahlster
(Vorsitzender) Dr. Walter Olthoff
Vorsitzender des Aufsichtsrats: Prof. Dr. h.c. Hans A. Aukes
Amtsgericht Kaiserslautern, HRB 2313
Sitz der Gesellschaft: Kaiserslautern (HRB 2313)
USt-Id.Nr.: DE 148646973
Steuernummer: 19/672/50006
-----------------------------------------------------------------------

[Peter]

Dear All,

a small remark on the packed symmetric matrices.

The packing A11, A21, A22, A31, A32, A33 ...

has one advantage that I make regularly use of, namely it can easily be resized,
especially if one reserves the memory at the beginning.

In addition, a nice feature would be the rotation of a symmetric/hermitian matrix which will be symmetric again,
provided the rotation matrix U is unitary, A = U * H * U^+ .

Best regards,
Peter

[Christoph Hertzberg]

Post by Peter
Dear All,
a small remark on the packed symmetric matrices.
The packing A11, A21, A22, A31, A32, A33 ...
has one advantage that I make regularly use of, namely it can easily be resized,
especially if one reserves the memory at the beginning.
In addition, a nice feature would be the rotation of a symmetric/hermitian matrix which will be symmetric again,
provided the rotation matrix U is unitary, A = U * H * U^+ .

`A` being hermitian is true for any `U` of course.
Does there exist an efficient (inplace?, without temporaries?) algorithm
for unitarian `U`? And does it depend on how `A`/`H` are stored?

Cheers,
Christoph

Post by Peter
Best regards,
Peter

--
Dr.-Ing. Christoph Hertzberg

Besuchsadresse der Nebengeschäftsstelle:
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 5
28359 Bremen, Germany

Postadresse der Hauptgeschäftsstelle Standort Bremen:
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 1
28359 Bremen, Germany

Tel.: [phone]
Zentrale: [phone]
E-Mail: ***@dfki.de

Weitere Informationen: http://www.dfki.de/robotik
-----------------------------------------------------------------------
Deutsches Forschungszentrum fuer Kuenstliche Intelligenz GmbH
Firmensitz: Trippstadter Straße 122, D-67663 Kaiserslautern
Geschaeftsfuehrung: Prof. Dr. Dr. h.c. mult. Wolfgang Wahlster
(Vorsitzender) Dr. Walter Olthoff
Vorsitzender des Aufsichtsrats: Prof. Dr. h.c. Hans A. Aukes
Amtsgericht Kaiserslautern, HRB 2313
Sitz der Gesellschaft: Kaiserslautern (HRB 2313)
USt-Id.Nr.: DE 148646973
Steuernummer: 19/672/50006
-----------------------------------------------------------------------

[Rasmus Munk Larsen]

I agree with Marton. It doesn' seem worthwhile to complicate the code and
increase the API surface for something that is know not to perform well on
modern architectures. (It was a different story on Crays with (tiny) solid
state memories :-) )

Post by Márton Danóczy
Hi Christoph,
supporting the "Packed Format" is not worth it, even LAPACK have dropped
support because of the FP routines' subpar performance (see the second
comment in this thread: https://github.com/Reference-
LAPACK/lapack/issues/245). The RFP format is and will keep being
supported by LAPACK (same Github thread).
Another reference: https://software.intel.com/en-us/mkl-developer-
reference-c-matrix-storage-schemes-for-lapack-routines
Best
Márton

Post by Christoph Hertzberg
Hi,
I suggest implementing a compact _triangular_ matrix which can be used
together with SelfAdjointView to represent selfadjoint matrices as well.
And then I'd like to point anyone interested in this discussion to this
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=42
This also mentions the RFPF suggested by Márton -- but we probably should
support both and start with the one which is easier to implement.
Cheers,
Christoph

Post by David A. Tellenbach
Hello together,
my name is David Tellenbach, I'm currently studying Computer Science at
the LMU Munich, Germany and got chosen for the Google Summer of Code 2018
Project "Faster Matrix Algebra for ATLAS", supervised by Dmitry Emeliyanov
and Stewart Martin-Haugh.
Google Summer of Code is a global program focused on bringing more
student developers into open source software development. Students work
with an open source organization on a three month programming project
during their break from school.
As you might know, our project's task is to implement support for
symmetric matrices for Eigen. A short project description is available via
the following link: https://summerofcode.withgoogl
e.com/projects/#5293950017994752 <https://summerofcode.withgoog
le.com/projects/#5293950017994752>
The official coding period hasn't started yet and lasts for three month,
from May, 14 until August, 14 2018. The time now is meant to get to know
the community and people involved.
Dense, sparse and diagonal matrices. Of these types, the class
Eigen::DiagonalMatrix seems to be the one that could be most similar to a
possible implementation of a class for symmetric matrices. There is no need
for storing all elements (as in the case of dense matrices) neither is a
sophisticated mechanism to find the position of scalars in the matrix
needed (as in the case of sparse matrices). Therefore I’d like to create
the for symmetric matrices by deriving from Eigen::EigenBase (as in the
case of diagonal matrices).
Of course one goal is to store only the upper or lower triangular part
of the matrix since this already defines it completely. Similar to
typedef typename internal::traits<Derived>::SymmetricVectorType
SymmetricVectorType;
// Store just one triangular part of the matrix
typedef Matrix<_Scalar, (SizeAtCompileTime * SizeAtCompileTime +
SizeAtCompileTime)/2,
1, 0, (MaxSizeAtCompileTime * MaxSizeAtCompileTime +
MaxSizeAtCompileTime)/2 ,1> SymmetricVectorType;
We plan to provide constructors which take either matrices of type
Eigen::Matrix<
> or Eigen::SelfAdjointView<
>.
What do you think about these broad plans so far? We are happy about any feedback.
Thanks,
David

--
Dr.-Ing. Christoph Hertzberg
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 5
<https://maps.google.com/?q=Robert-Hooke-Stra%C3%9Fe+5+%0D%0A%C2%A028359+Bremen,+Germany&entry=gmail&source=g>
28359 Bremen, Germany
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 1
<https://maps.google.com/?q=Robert-Hooke-Stra%C3%9Fe+1+%0D%0A%C2%A028359+Bremen,+Germany&entry=gmail&source=g>
28359 Bremen, Germany
Tel.: [phone]
Zentrale: [phone]
Weitere Informationen: http://www.dfki.de/robotik
-----------------------------------------------------------------------
Deutsches Forschungszentrum fuer Kuenstliche Intellig
<https://maps.google.com/?q=ches+Forschungszentrum+fuer+Kuenstliche+Intellig&entry=gmail&source=g>enz
GmbH
Firmensitz: Trippstadter Straße 122, D-67663 Kaiserslautern
<https://maps.google.com/?q=Trippstadter+Stra%C3%9Fe+122,+D-67663+Kaiserslautern&entry=gmail&source=g>
Geschaeftsfuehrung: Prof. Dr. Dr. h.c. mult. Wolfgang Wahlster
(Vorsitzender) Dr. Walter Olthoff
Vorsitzender des Aufsichtsrats: Prof. Dr. h.c. Hans A. Aukes
Amtsgericht Kaiserslautern, HRB 2313
Sitz der Gesellschaft: Kaiserslautern (HRB 2313)
USt-Id.Nr.: DE 148646973
Steuernummer: 19/672/50006
-----------------------------------------------------------------------

[David A. Tellenbach]

Hello together,

first of all: Thank you again for supporting the idea of implementing a class
for hermitian matrices in Eigen.

I would be happy if you could give me some feedback on what I have implemented
so far:

https://bitbucket.org/tellenbach/eigen/src/HermitianMatrix/

I've decided to name the class HermitianMatrix and to store either the upper or
lower triangular part in rectangular full packed format as suggested by some of
you. For further information about this storage format see

http://www.netlib.org/lapack/lawnspdf/lawn199.pdf

The current implementation consists of the following functionality:

1. Constructors to build a HermitianMatrix by passing a dense Matrix. The
elements are stored in a dense Matrix in either row or column major.
2. Core evaluators for providing the methods coeff(), coeffRef() and the
operator ()
3. Sum and difference assignment operators += and -=
4. CwiseBinaryOp for supporting addition and subtraction of two hermitian
matrices.
5. Assignment of CwiseBinaryOp to HermitianMatrix.

All methods try to use transform the problems of hermitian matrices to problems
of the underlying dense matrix that stores the matrix elements.

A problem I'm currently facing is the following: In the case of complex scalars,
I want to support real hermitian and not just symmetric matrices. That means when
the upper triangular part of a matrix is stored in RFPF and the method coeff() or
coeffRef() is called for an element in the lower triangular part, the complex
conjugated should be returned. Any idea how this could be handled?

I have benchmarked (using Google benchmark) some functionality I've implemented
and it looks very good. As an attachment you find some of the benchmark results
for addition (including evaluation) for fixed and dynamically sized matrices.

I would be really happy about some feedback.

All the best,
David

Hi all,

I doubt that RFPF is superior to simple packed format for 4x4 matrices

at least in the case of simple operations I agree (even for bigger matrices). If we additionally consider Peter's argument that a packed format is easier to resize and Christoph's argument that someone could use packed storage and cannot change that, I think we shouldn't just drop the idea of implementing it completely.
Nevertheless the RFPF performance speaks for itself. Without deciding to never implement packed storage, the implementation of RFPF should have priority, I guess.
Thanks,
David

I agree with Marton. It doesn' seem worthwhile to complicate the code and increase the API surface for something that is know not to perform well on modern architectures. (It was a different story on Crays with (tiny) solid state memories :-) )
Hi Christoph,
supporting the "Packed Format" is not worth it, even LAPACK have dropped support because of the FP routines' subpar performance (see the second comment in this thread: https://github.com/Reference-LAPACK/lapack/issues/245 <https://github.com/Reference-LAPACK/lapack/issues/245>). The RFP format is and will keep being supported by LAPACK (same Github thread).
Another reference: https://software.intel.com/en-us/mkl-developer-reference-c-matrix-storage-schemes-for-lapack-routines <https://software.intel.com/en-us/mkl-developer-reference-c-matrix-storage-schemes-for-lapack-routines>
Best
Márton
Hi,
I suggest implementing a compact _triangular_ matrix which can be used together with SelfAdjointView to represent selfadjoint matrices as well.
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=42 <http://eigen.tuxfamily.org/bz/show_bug.cgi?id=42>
This also mentions the RFPF suggested by Márton -- but we probably should support both and start with the one which is easier to implement.
Cheers,
Christoph
Hello together,
my name is David Tellenbach, I'm currently studying Computer Science at the LMU Munich, Germany and got chosen for the Google Summer of Code 2018 Project "Faster Matrix Algebra for ATLAS", supervised by Dmitry Emeliyanov and Stewart Martin-Haugh.
Google Summer of Code is a global program focused on bringing more student developers into open source software development. Students work with an open source organization on a three month programming project during their break from school.
As you might know, our project's task is to implement support for symmetric matrices for Eigen. A short project description is available via the following link: https://summerofcode.withgoogle.com/projects/#5293950017994752 <https://summerofcode.withgoogle.com/projects/#5293950017994752> <https://summerofcode.withgoogle.com/projects/#5293950017994752 <https://summerofcode.withgoogle.com/projects/#5293950017994752>>
The official coding period hasn't started yet and lasts for three month, from May, 14 until August, 14 2018. The time now is meant to get to know the community and people involved.
As far as I see, Eigen basically provides three types of matrices: Dense, sparse and diagonal matrices. Of these types, the class Eigen::DiagonalMatrix seems to be the one that could be most similar to a possible implementation of a class for symmetric matrices. There is no need for storing all elements (as in the case of dense matrices) neither is a sophisticated mechanism to find the position of scalars in the matrix needed (as in the case of sparse matrices). Therefore I’d like to create the for symmetric matrices by deriving from Eigen::EigenBase (as in the case of diagonal matrices).
typedef typename internal::traits<Derived>::SymmetricVectorType SymmetricVectorType;
// Store just one triangular part of the matrix
typedef Matrix<_Scalar, (SizeAtCompileTime * SizeAtCompileTime + SizeAtCompileTime)/2,
1, 0, (MaxSizeAtCompileTime * MaxSizeAtCompileTime + MaxSizeAtCompileTime)/2 ,1> SymmetricVectorType;
We plan to provide constructors which take either matrices of type Eigen::Matrix<
> or Eigen::SelfAdjointView<
>.
What do you think about these broad plans so far? We are happy about any feedback.
Thanks,
David
--
Dr.-Ing. Christoph Hertzberg
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 5 <https://maps.google.com/?q=Robert-Hooke-Stra%C3%9Fe+5+%0D%0A%C2%A028359+Bremen,+Germany&entry=gmail&source=g>
28359 Bremen, Germany
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 1 <https://maps.google.com/?q=Robert-Hooke-Stra%C3%9Fe+1+%0D%0A%C2%A028359+Bremen,+Germany&entry=gmail&source=g>
28359 Bremen, Germany
Tel.: [phone]
Zentrale: [phone]
Weitere Informationen: http://www.dfki.de/robotik <http://www.dfki.de/robotik>
-----------------------------------------------------------------------
Deutsches Forschungszentrum fuer Kuenstliche Intellig <https://maps.google.com/?q=ches+Forschungszentrum+fuer+Kuenstliche+Intellig&entry=gmail&source=g>enz GmbH
Firmensitz: Trippstadter Straße 122, D-67663 Kaiserslautern <https://maps.google.com/?q=Trippstadter+Stra%C3%9Fe+122,+D-67663+Kaiserslautern&entry=gmail&source=g>
Geschaeftsfuehrung: Prof. Dr. Dr. h.c. mult. Wolfgang Wahlster
(Vorsitzender) Dr. Walter Olthoff
Vorsitzender des Aufsichtsrats: Prof. Dr. h.c. Hans A. Aukes
Amtsgericht Kaiserslautern, HRB 2313
Sitz der Gesellschaft: Kaiserslautern (HRB 2313)
USt-Id.Nr.: DE 148646973
Steuernummer: 19/672/50006
-----------------------------------------------------------------------

[Rasmus Munk Larsen]

Hi David,

The performance numbers look impressive. Nice work. Looking forward to
seeing matrix multiply numbers :-)

Cheers,
Rasmus

Post by David A. Tellenbach
Hello together,
first of all: Thank you again for supporting the idea of implementing a class
for hermitian matrices in Eigen.
I would be happy if you could give me some feedback on what I have implemented
https://bitbucket.org/tellenbach/eigen/src/HermitianMatrix/
I've decided to name the class HermitianMatrix and to store either the upper or
lower triangular part in rectangular full packed format as suggested by some of
you. For further information about this storage format see
http://www.netlib.org/lapack/lawnspdf/lawn199.pdf
1. Constructors to build a HermitianMatrix by passing a dense Matrix. The
elements are stored in a dense Matrix in either row or column major.
2. Core evaluators for providing the methods coeff(), coeffRef() and the
operator ()
3. Sum and difference assignment operators += and -=
4. CwiseBinaryOp for supporting addition and subtraction of two hermitian
matrices.
5. Assignment of CwiseBinaryOp to HermitianMatrix.
All methods try to use transform the problems of hermitian matrices to problems
of the underlying dense matrix that stores the matrix elements.
A problem I'm currently facing is the following: In the case of complex scalars,
I want to support real hermitian and not just symmetric matrices. That means when
the upper triangular part of a matrix is stored in RFPF and the method coeff() or
coeffRef() is called for an element in the lower triangular part, the complex
conjugated should be returned. Any idea how this could be handled?
I have benchmarked (using Google benchmark) some functionality I've implemented
and it looks very good. As an attachment you find some of the benchmark results
for addition (including evaluation) for fixed and dynamically sized matrices.
I would be really happy about some feedback.
All the best,
David
Am 04.05.2018 um 20:28 schrieb David A. Tellenbach <
Hi all,
I doubt that RFPF is superior to simple packed format for 4x4 matrices
at least in the case of simple operations I agree (even for bigger
matrices). If we additionally consider Peter's argument that a packed
format is easier to resize and Christoph's argument that someone could use
packed storage and cannot change that, I think we shouldn't just drop the
idea of implementing it completely.
Nevertheless the RFPF performance speaks for itself. Without deciding to
never implement packed storage, the implementation of RFPF should have
priority, I guess.
Thanks,
David
I agree with Marton. It doesn' seem worthwhile to complicate the code and
increase the API surface for something that is know not to perform well on
modern architectures. (It was a different story on Crays with (tiny) solid
state memories :-) )

Post by Márton Danóczy
Hi Christoph,
supporting the "Packed Format" is not worth it, even LAPACK have dropped
support because of the FP routines' subpar performance (see the second
https://github.com/Reference-LAPACK/lapack/issues/245). The RFP format
is and will keep being supported by LAPACK (same Github thread).
https://software.intel.com/en-us/mkl-developer-reference-c-matrix-storage-schemes-for-lapack-routines
Best
Márton
On 3 May 2018 at 19:28, Christoph Hertzberg <

Post by Christoph Hertzberg
Hi,
I suggest implementing a compact _triangular_ matrix which can be used
together with SelfAdjointView to represent selfadjoint matrices as well.
And then I'd like to point anyone interested in this discussion to this
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=42
This also mentions the RFPF suggested by Márton -- but we probably
should support both and start with the one which is easier to implement.
Cheers,
Christoph

Post by David A. Tellenbach
Hello together,
my name is David Tellenbach, I'm currently studying Computer Science at
the LMU Munich, Germany and got chosen for the Google Summer of Code 2018
Project "Faster Matrix Algebra for ATLAS", supervised by Dmitry Emeliyanov
and Stewart Martin-Haugh.
Google Summer of Code is a global program focused on bringing more
student developers into open source software development. Students work
with an open source organization on a three month programming project
during their break from school.
As you might know, our project's task is to implement support for
symmetric matrices for Eigen. A short project description is available via
https://summerofcode.withgoogle.com/projects/#5293950017994752 <
https://summerofcode.withgoogle.com/projects/#5293950017994752>
The official coding period hasn't started yet and lasts for three
month, from May, 14 until August, 14 2018. The time now is meant to get to
know the community and people involved.
Dense, sparse and diagonal matrices. Of these types, the class
Eigen::DiagonalMatrix seems to be the one that could be most similar to a
possible implementation of a class for symmetric matrices. There is no need
for storing all elements (as in the case of dense matrices) neither is a
sophisticated mechanism to find the position of scalars in the matrix
needed (as in the case of sparse matrices). Therefore I’d like to create
the for symmetric matrices by deriving from Eigen::EigenBase (as in the
case of diagonal matrices).
Of course one goal is to store only the upper or lower triangular part
of the matrix since this already defines it completely. Similar to
typedef typename internal::traits<Derived>::SymmetricVectorType
SymmetricVectorType;
// Store just one triangular part of the matrix
typedef Matrix<_Scalar, (SizeAtCompileTime * SizeAtCompileTime +
SizeAtCompileTime)/2,
1, 0, (MaxSizeAtCompileTime * MaxSizeAtCompileTime +
MaxSizeAtCompileTime)/2 ,1> SymmetricVectorType;
We plan to provide constructors which take either matrices of type
Eigen::Matrix<
> or Eigen::SelfAdjointView<
>.
What do you think about these broad plans so far? We are happy about any feedback.
Thanks,
David

--
Dr.-Ing. Christoph Hertzberg
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 5
<https://maps.google.com/?q=Robert-Hooke-Stra%C3%9Fe+5+%0D%0A%C2%A028359+Bremen,+Germany&entry=gmail&source=g>
28359 Bremen, Germany
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 1
<https://maps.google.com/?q=Robert-Hooke-Stra%C3%9Fe+1+%0D%0A%C2%A028359+Bremen,+Germany&entry=gmail&source=g>
28359 Bremen, Germany
Tel.: [phone]
Zentrale: [phone]
Weitere Informationen: http://www.dfki.de/robotik
-----------------------------------------------------------------------
Deutsches Forschungszentrum fuer Kuenstliche Intellig
<https://maps.google.com/?q=ches+Forschungszentrum+fuer+Kuenstliche+Intellig&entry=gmail&source=g>enz
GmbH
Firmensitz: Trippstadter Straße 122, D-67663 Kaiserslautern
<https://maps.google.com/?q=Trippstadter+Stra%C3%9Fe+122,+D-67663+Kaiserslautern&entry=gmail&source=g>
Geschaeftsfuehrung: Prof. Dr. Dr. h.c. mult. Wolfgang Wahlster
(Vorsitzender) Dr. Walter Olthoff
Vorsitzender des Aufsichtsrats: Prof. Dr. h.c. Hans A. Aukes
Amtsgericht Kaiserslautern, HRB 2313
Sitz der Gesellschaft: Kaiserslautern (HRB 2313)
USt-Id.Nr.: DE 148646973
Steuernummer: 19/672/50006
-----------------------------------------------------------------------