ScaLAPACKMatrix: add functions to perform addition, multiplication and row/column scaling #5869

BenBrands · 2018-02-06T16:01:13Z

No description provided.

davydden

LGTM, thanks for improvements 👍 few minor comments. Just wanted to raise the discussion on naming, so don't rush in re-naming anything.

p.s. also add a changelog for the new functionality

davydden · 2018-02-06T16:18:08Z

include/deal.II/lac/scalapack.h

+   *
+   * The matrices <tt>A</tt> and <tt>B</tt> must have the same process grid
+   *
+   * Following alignment conditions have to be fulfilled: MB_A = NB_B and NB_A = MB_B


in class description we have

...into $MB$ by $NB$ blocks...

not sure how it renders, but probably you should do $MB_A$ and alike here and below in this PR to keep it consistent.

I will change it at the one point in the class description, if that is fine with you.

i would actually prefer the opposite, i think it will render better as latex formula. Otherwise you have it as a text MB_A

p.s. that's also how we do it all other places inside ScaLAPACK, including the most recent PR with some docu changes #5875

davydden · 2018-02-06T16:18:51Z

include/deal.II/lac/scalapack.h

+  /**
+   * Transposing assignment: <i>A = B<sup>T</sup></i>
+   *
+   * The matrices <tt>A</tt> and <tt>B</tt> must have the same process grid


davydden · 2018-02-06T16:19:07Z

include/deal.II/lac/scalapack.h

+   *
+   * else  <i>op(B) = B</i>
+   *
+   * The matrices <tt>A</tt> and <tt>B</tt> must have the same process grid


davydden · 2018-02-06T16:23:05Z

include/deal.II/lac/scalapack.h

+   * |   false     | MB_A = MB_B <br> NB_A = NB_B |
+   * |   true      | MB_A = NB_B <br> NB_A = MB_B |
+   */
+  void add(const ScaLAPACKMatrix<NumberType> &B,


now i see what you wanted to implement, elsewhere (for vector classes) we call this:

sadd (const Number s, const Number a, const VectorSpaceVector< Number > & V)

since other matrix classes don't support this operation, i would rather stick with

sadd (const .. s, const ... a,const ... B, const bool transpose_B = false)

@dealii/dealii any opinions on naming of this matrix operation?

davydden · 2018-02-06T16:29:20Z

include/deal.II/lac/scalapack.h

+   *
+   * The matrices <tt>A</tt>, <tt>B</tt> and <tt>C</tt> must have the same process grid.
+   */
+  void mult(ScaLAPACKMatrix<NumberType> &C,


@dealii/dealii this is kind of smult (scale matrix and add multiplication). For consistency with vectors, i would probably move two numbers in front i.e. (const NumberType s, const NumberType a, const Matrix....).
But then it would be inconsistent with FullMatrix

void mmult (FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const;

hm...

We used to have similar functions in some of the vector classes (e.g here), so I would say that the factors and matrices should be ordered according to the operation along with the input matrices, even if that means that we lose the default values. Does a default value of 0 for b make sense anyway? This then just becomes an assignment operation.

davydden · 2018-02-06T16:42:30Z

source/lac/scalapack.cc

+      Assert(this->row_block_size==B.column_block_size,ExcDimensionMismatch(this->column_block_size,B.row_block_size));
+      Assert(B.row_block_size==C.column_block_size,ExcDimensionMismatch(B.column_block_size,C.column_block_size));
+    }
+  Threads::Mutex::ScopedLock lock (mutex);


we need mutex only when mutable member variables like work array are being used. I think you don't need it here.

davydden · 2018-02-06T16:42:57Z

source/lac/scalapack.cc

+  if (this->grid->mpi_process_is_active)
+    {
+      for (int i=0; i<n_local_rows; ++i)
+        {


don't need these brackets

davydden · 2018-02-06T16:44:38Z

source/lac/scalapack.cc

+    {
+      for (int i=0; i<n_local_rows; ++i)
+        {
+          for (int j=0; j<n_local_columns; ++j)


you should probably loop over columns first, take value and multiply all row elements.

davydden · 2018-02-06T16:45:42Z

source/lac/scalapack.cc

+          const int glob_i = global_row(i);
+          for (int j=0; j<n_local_columns; ++j)
+            {
+              local_el(i,j) *= factors[glob_i];


just do const NumberType & s = factors[global_row(i)];

davydden · 2018-02-06T16:46:03Z

source/lac/scalapack.cc

+        {
+          const int glob_i = global_row(i);
+          for (int j=0; j<n_local_columns; ++j)
+            {


and remove these brackets

davydden · 2018-02-06T18:14:35Z

/run-tests

davydden

LGTM

davydden · 2018-02-12T06:12:32Z

@dealii/dealii can someone please have a look at the documentation here and the chosen names/arguments for member functions?

bangerth · 2018-02-13T21:14:24Z

include/deal.II/lac/scalapack.h

+   *
+   * The matrices <tt>A</tt> and <tt>B</tt> must have the same process grid.
+   *
+   * Following alignment conditions have to be fulfilled: $MB_A$ = $NB_B$ and $NB_A$ = $MB_B$


Say "The following ...". I had to think about whether "following" is a verb or adjective. Same in a few places below.

Also end the sentence with a period.

You beat me to it! I was also going to make the suggestion that the style of all of the documentation be made more consistent with that of the rest of the class. The conditions can also be put in a @note.

bangerth · 2018-02-13T21:15:27Z

include/deal.II/lac/scalapack.h

+  /**
+   * Matrix-addition:
+   *
+   * <i>A = a A + b op(B)</i>


remove the empty line so that this reads back to back

think about whether you could do the same in the following lines

bangerth · 2018-02-13T21:31:19Z

include/deal.II/lac/scalapack.h

+   *
+   * Copies of @p factors have to available on all processes of the underlying MPI communicator
+   */
+  void scale_columns(std::vector<NumberType> &factors);


make the argument const

bangerth · 2018-02-13T21:31:27Z

include/deal.II/lac/scalapack.h

+   *
+   * Copies of @p factors have to available on all processes of the underlying MPI communicator
+   */
+  void scale_rows(std::vector<NumberType> &factors);


btw, we should probably use Vector<NumberType> instead here. It won't make much of a difference at the end of the day, but is probably more consistent representation of a diagonal matrix used here as pre/post multiplication.

Why is Vector<NumberType> more consistent than std::vector<NumberType>?
We do not use Vector<NumberType> in ScaLAPACKMatrix unlike std::vector<NumberType>, which is used at three places.
We could provide two public interfaces and use one private function (working with pointers) internally.

we may seek just a few eigenvalues or all eigenvalues of a matrix. In other places (like Arpack for sparse large matrices) we also use std::vector<> which is problably ok.

Here, however, this vector has to have as many elements as matrix rows/columns. So in case of a square matrix, this is simply a diagonal matrix represented by a vector. With other matrix classes AFAIK we never use std::vector<> to represent a diagonal, i.e. https://www.dealii.org/developer/doxygen/deal.II/classSparseMatrix.html#a07ba942c5b7466b3895fed3104fbcc26

bangerth · 2018-02-13T21:32:12Z

source/lac/scalapack.cc

+  // https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lgemm.htm
+  if (!transpose_A && !transpose_B)
+    {
+      Assert(this->n_columns==B.n_rows,ExcDimensionMismatch(this->n_columns,B.n_rows));


thanks for adding all of these asserts!

jppelteret · 2018-02-13T21:33:17Z

@BenBrands I'll take a look through the documentation tomorrow!

jppelteret · 2018-02-13T21:36:28Z

Ok, so @bangerth happened post an initial few comment just commented, but I'll stand by my initial commitment and will give it a proper look tomorrow.

jppelteret · 2018-02-14T07:22:19Z

include/deal.II/lac/scalapack.h

+   *
+   * The matrices <tt>A</tt> and <tt>B</tt> must have the same process grid.
+   *
+   * Following alignment conditions have to be fulfilled: $MB_A$ = $NB_B$ and $NB_A$ = $MB_B$


I think that for all of the functions that you're documented, the style could be made a bit more similar to that of the existing document. So here, for example, you could consider writing it as follows:

Copy the transpose of the contents of the distributed matrix @p B into this one. This transpose assignment operation modifies the calling matrix @p A such that $A = B^{T}$. @note The matrices @p A and @p B must have the same process grid. The following alignment conditions also have to be fulfilled: $MB_A = NB_B$ and $NB_A = MB_B$.

Note also that $MB_A$ = $NB_B$ should be $MB_A = NB_B$ so that the equals is also rendered using LaTeX. Same applies below.

I discussed with @davydden the documentation of ScaLAPACKMatrix and we agree that it has to be reworked to make it uniform. I will do that in a different PR based on the style used for the doc in this PR.

jppelteret · 2018-02-14T07:22:49Z

include/deal.II/lac/scalapack.h

+  /**
+   * Transposing assignment: <i>A = B<sup>T</sup></i>
+   *
+   * The matrices <tt>A</tt> and <tt>B</tt> must have the same process grid.


Why have you used HTML as opposed to LaTeX formatting (both here and below)?

You should use @p to format the parameters of a function, i.e. @p A and @p B, when inline with other text.

Same applies below.

I am also in favour of LaTeX, but we have a mixture of the two, apparently #5897

jppelteret · 2018-02-14T07:30:01Z

include/deal.II/lac/scalapack.h

+   *
+   * if(transpose_B)  <i>op(B) = B<sup>T</sup></i>
+   *
+   * else  <i>op(B) = B</i>


I'd prefer for this to be expressed in words rather than as pseudo-code. i.e.

If the @p transpose_B is set then the transpose of @p B is used. That is to say that the operation $op(B)$ corresponds to $B^{T}$. If this flag is set to false, then $op(B)$ applies @p B with its default ordering.

This also applies below.

I would favor the pseudo-code as it is short and precise.
But I can adapt it if wanted.

I suppose that the pseudo-code is consistent with LapackFullMatrix and elements of the Vector class. Although I'm not a huge fan of it, I'm not going to push the issue since there's precedent for it.

jppelteret · 2018-02-14T07:44:44Z

include/deal.II/lac/scalapack.h

+   *
+   * The matrices <tt>A</tt>, <tt>B</tt> and <tt>C</tt> must have the same process grid.
+   */
+  void mult(ScaLAPACKMatrix<NumberType> &C,


We used to have similar functions in some of the vector classes (e.g here), so I would say that the factors and matrices should be ordered according to the operation along with the input matrices, even if that means that we lose the default values. Does a default value of 0 for b make sense anyway? This then just becomes an assignment operation.

jppelteret · 2018-02-14T07:45:08Z

include/deal.II/lac/scalapack.h

+  /**
+   * Matrix-matrix-multiplication.
+   *
+   * The optional parameter <tt>adding</tt> determines, whether the result is


jppelteret · 2018-02-14T07:46:52Z

include/deal.II/lac/scalapack.h

+   *
+   * else <i>C = A<sup>T</sup>*B</i>
+   *
+   * Assumes that <tt>A</tt> and <tt>B</tt> have compatible sizes and that


It is assumed that...

jppelteret · 2018-02-14T07:49:07Z

include/deal.II/lac/scalapack.h

@@ -389,6 +546,24 @@ class ScaLAPACKMatrix : protected TransposeTable<NumberType>
   */
  NumberType &local_el(const unsigned int loc_row, const unsigned int loc_column);

+  /**
+   * Scale the columns of the distributed matrix by scalars in array @p factors.


by the scalars provided in the array

jppelteret · 2018-02-14T07:49:27Z

include/deal.II/lac/scalapack.h

+   *
+   * The array @p factors must have as many entries as the matrix columns.
+   *
+   * Copies of @p factors have to available on all processes of the underlying MPI communicator


have to be; fullstop at end of sentence.

jppelteret · 2018-02-14T07:49:53Z

include/deal.II/lac/scalapack.h

+   *
+   * The array @p factors must have as many entries as the matrix rows.
+   *
+   * Copies of @p factors have to available on all processes of the underlying MPI communicator


have to be; fullstop at end of sentence.

jppelteret · 2018-02-14T07:50:04Z

include/deal.II/lac/scalapack.h

+  void scale_columns(std::vector<NumberType> &factors);
+
+  /**
+   * Scale the rows of the distributed matrix by scalars in array @p factors.


by the scalars provided in the array

jppelteret · 2018-02-14T07:52:09Z

Thanks for taking the time to add this functionality @BenBrands! Since @davydden checked the implementation already, lets clean up the documentation a little bit and then I think that this is good to go!

BenBrands · 2018-02-14T14:09:12Z

Thank you all for your comments.
The documentation for this class needs an update, I would like to do so in a different PR (based on the doc in this PR and #5898).

BenBrands · 2018-02-14T20:17:21Z

The functions scale_columns() and scale_rows() take const ArrayView<const NumberType> as input!

jppelteret · 2018-02-15T05:57:11Z

@BenBrands Thanks again for the contribution. I've created an issue for the proposed documentation enhancements.

jppelteret · 2018-02-15T07:36:29Z

@davydden @BenBrands Lets discuss the changes made to scale_columns and scale_rows before we merge this. This did not work out quite as I thought it would. I guess that we would need to more closely follow the idiom used here and here to achieve it.

davydden · 2018-02-15T11:50:43Z

include/deal.II/lac/scalapack.h

+   *
+   * | transpose_B |          Block Sizes         |                    Operation                  |
+   * | :---------: | :--------------------------: | :-------------------------------------------: |
+   * |   false     | $MB_A=MB_B$ <br> $NB_A=NB_B$ |  $\mathbf{A} = a \mathbf{A} + b \mathbf{B}$   |


in another PR #5899 (comment) @kronbichler raised the point whether we shall use bold symbols (\mathbf) for matrices and vectors.

Personally, I don't mind sticking with plain LaTeX A, B, c for linear algebra, probably with \cdots to denote multiplication. Downside is that we won't distinguish between scalars and vectors in this notation.

@dealii/dealii thoughts?

jppelteret · 2018-02-15T12:38:44Z

include/deal.II/base/array_view.h

+template <typename Number>
+inline
+ArrayView<const Number>
+make_array_view (const ArrayView<Number> &array_view)


@masterleinad I think that you added some of the other make_array_view functions previously. Are you happy with these two additions (they help instantiation of some functions)?

Sure. In the end, I prefer to have external make_array_functions instead of modifying the class ArrayView itself. If this helps with compatibility, I am totally fine with it.

Great, thanks for the confirmation!

jppelteret · 2018-02-15T18:31:06Z

Ok, so I think that this is ready to merge but since I've contributed directly to this PR it would be more appropriate for someone else to make the final call.

davydden · 2018-02-17T21:20:42Z

tests/scalapack/scalapack_13.cc

+
+  const std::vector<unsigned int> sizes = {{400,500}};
+
+  // test A = alpha A + beta B


could you please split this test into 2 (each does different operation), previous into 4 (one operation in each test) and the next one into two (rows/cols). If/when something goes south, it's much easier to debug a small tests as opposed to a huge one which covers all the cases .

davydden · 2018-02-18T21:46:07Z

source/lac/scalapack.cc

+void ScaLAPACKMatrix<NumberType>::scale_columns(const InputVector &factors)
+{
+  internal::scale_columns(*this, make_array_view(factors),
+                          this->grid->mpi_process_is_active);


minor nitpick, i think it's easier to read as

if (this->grid->mpi_process_is_active) internal::scale_columns(...);

that is, when you move that if outside of internal function which actually does the job.

davydden added the Linear Algebra label Feb 6, 2018

davydden reviewed Feb 6, 2018

View reviewed changes

BenBrands force-pushed the scalapack branch from 25dfca1 to c31ecad Compare February 7, 2018 18:38

davydden approved these changes Feb 7, 2018

View reviewed changes

davydden added the ready for review label Feb 7, 2018

bangerth reviewed Feb 13, 2018

View reviewed changes

jppelteret requested changes Feb 14, 2018

View reviewed changes

BenBrands force-pushed the scalapack branch from c31ecad to 0135313 Compare February 14, 2018 14:05

BenBrands added 5 commits February 14, 2018 21:10

routines to scale rows/columns of ScaLAPACKMatrix + tests

8be5a8d

add routines to add ScaLAPACKMatrices + tests

be1b887

add ScaLAPACKMatrix.copy_transposed()

5a41381

add routines for multiplication of ScaLAPACKMatrices + tests

3572a1e

avoid compiler warning -unused-variable

ba17171

BenBrands force-pushed the scalapack branch from 0135313 to ba17171 Compare February 14, 2018 20:10

jppelteret mentioned this pull request Feb 15, 2018

Improve documentation of LAPACKFullMatrix and ScaLAPACKFullMatrix #5903

Open

jppelteret approved these changes Feb 15, 2018

View reviewed changes

Generalise Scalapack scale_columns and scale_rows functions

f5660de

BenBrands force-pushed the scalapack branch from e67b91a to f5660de Compare February 15, 2018 11:45

davydden reviewed Feb 15, 2018

View reviewed changes

davydden approved these changes Feb 15, 2018

View reviewed changes

jppelteret reviewed Feb 15, 2018

View reviewed changes

davydden mentioned this pull request Feb 15, 2018

LAPACKFullMatrix - add another Tmmult() #5898

Merged

davydden reviewed Feb 17, 2018

View reviewed changes

davydden reviewed Feb 18, 2018

View reviewed changes

Change scale rows/columns and split tests

51601c9

davydden added the Reviewed and ready to merge label Feb 19, 2018

davydden merged commit abd4921 into dealii:master Feb 19, 2018


		const std::vector<unsigned int> sizes = {{400,500}};

		// test A = alpha A + beta B

ScaLAPACKMatrix: add functions to perform addition, multiplication and row/column scaling #5869

ScaLAPACKMatrix: add functions to perform addition, multiplication and row/column scaling #5869

Conversation

BenBrands commented Feb 6, 2018

davydden left a comment • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

davydden commented Feb 6, 2018

davydden left a comment

Choose a reason for hiding this comment

davydden commented Feb 12, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

davydden Feb 14, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jppelteret commented Feb 13, 2018

jppelteret commented Feb 13, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jppelteret commented Feb 14, 2018

BenBrands commented Feb 14, 2018

BenBrands commented Feb 14, 2018

jppelteret commented Feb 15, 2018

jppelteret commented Feb 15, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jppelteret commented Feb 15, 2018

davydden Feb 17, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

davydden left a comment •

edited

davydden Feb 14, 2018 •

edited

davydden Feb 17, 2018 •

edited