Replace matrix with Matrix in various cyclic inversions by d7919 · Pull Request #900 · boutproject/BOUT-dev

d7919 · 2018-02-23T22:12:36Z

Helps enable bounds checking for these various inversion classes built on cyclic_reduction.hxx.

Typically we use new to create a new cyclic_reduction instance. Now all of the allocation is associated with Matrix or Array so should be fast. We may be able to avoid using new and instead just create the solver at the scope in which it is needed. Similar could be done with the storage arrays used in these solvers which are typically currently created in the constructor.

Instead of matrix (pointer to pointer). This is the first step in allowing the cyclic solvers to replace all static/class level memory.

One remains unchanged (see #895 for why)

in cyclic_reduction.hxx

This ensures bounds checking when compiled with appropriate CHECK

Could now move all storage from class level to method scope

Also replaces new with Array in laplacexz cyclic

bendudson · 2018-02-24T06:02:58Z

include/cyclic_reduction.hxx

  /// Solve a single triadiagonal system
  /// 
-  void solve(T rhs[], T x[]) {
+  void solve(Matrix<T> &rhs, Matrix<T> &x) {


I think this should be Array, rather than Matrix. This function solves a single system (Ax=b), so x and rhs here should be Arrays. On line 128 the address of rhs and x were used, to convert T* to T** . Is there a way to convert an Array to a single-column (or single-row) Matrix?

There's not currently a way to convert an Array to a Matrix, but it should be easy enough to add a constructor to Matrix that does this. In general there isn't a unique way to do this, but we could always make the choice that either the matrix is 1xN or Nx1.

I've been thinking about the possibility of making Matrix::operator[] return an Array slice of the data (either by creating a new Array and copying in values, or making Matrix::data actually be a vector of Array). WIth this I think it would make sense to get a 1xN shape Matrix when constructing from an Array of size N.

This interface to solve isn't actually used anywhere at the moment (at least not in any of the framework or the tests) so we could remove it and revisit if/when this is needed/the Matrix<-->Array conversion is sorted?

I should note there's a similar argument for setcoefs.

I've added interfaces to setCoefs and solve that take Array instead of Matrix. These then copy into a Matrix and pass to the main routine before copying back into the Array (where needed).

bendudson · 2018-02-24T06:04:36Z

include/cyclic_reduction.hxx

  /// Solve a set of tridiagonal systems
  /// 
-  void solve(int nrhs, T **rhs, T **x) {
+  void solve(int nrhs, Matrix<T> &rhs, Matrix<T> &x) {


Can the nrhs argument be removed? This should be the number of columns in rhs? I guess it is bounds checked anyway, but shouldn't be needed.

Now removed.

Can get this information from the passed objects.

ZedThree

Some petty things, some minor things, a couple of larger questions:

This changes the signatures of a few functions. Is this going to break backwards-compatibility? Or more precisely, will people need to update physics models? If so, can we just deprecate the old signatures instead of removing them?
This is somewhat orthogonal to this PR, but we appear to always new a CyclicReduce<>. I'm not sure there's any need for this, can we replace with a plain value?

ZedThree · 2018-02-26T09:32:23Z

include/cyclic_reduction.hxx

-
  /// Set the entries in the matrix to be inverted
  ///
  /// @param[in] nsys   The number of independent matrices to be solved


Doxygen comment is now out of date

ZedThree · 2018-02-26T09:33:19Z

include/cyclic_reduction.hxx

  void setPeriodic(bool p=true) {periodic=p;}

-  /// Set up a single equation to invert
-  void setCoefs(T a[], T b[], T c[]) {


Should this be deprecated instead of removed? Or is the signature compatible with Matrix<T>?

Yes I guess this should be deprecated -- I've replaced all uses in the repository but I suppose someone could be using it in their private model. My main issue was how to then construct a Matrix from this to pass through to the actual routine.

The same for the solve interfaces I guess.

ZedThree · 2018-02-26T09:34:10Z

include/cyclic_reduction.hxx

+  // void solve(Matrix<T> &rhs, Matrix<T> &x) {
+  //   // Solving single system
+  //   solve(1, &rhs, &x);
+  // }


Commented-out code should be removed instead

I thought I'd not staged this, oops.

ZedThree · 2018-02-26T09:35:14Z

include/cyclic_reduction.hxx

+  // }

  /// Solve a set of tridiagonal systems
  /// 


Could the doxygen comment be expanded here?

ZedThree · 2018-02-26T09:37:40Z

include/cyclic_reduction.hxx

-	ic[2] *= -beta;
+        ifc(j, 1) = co(j, 4*i + 1) - beta * ifc(j, 0);
+        ifc(j, 0) = co(j, 4*i);
+	ifc(j, 2) *= -beta;


Whitespace issue

ZedThree · 2018-02-26T09:37:54Z

include/cyclic_reduction.hxx

-	ic[0] *= -alpha;
-        ic[1] = c[4*i + 1] - alpha*ic[2];
-        ic[2] = c[4*i + 2];
+	ifc(j, 4+0) *= -alpha;


ZedThree · 2018-02-26T09:46:50Z

src/invert/parderiv/impls/cyclic/cyclic.cxx

-	c[k][y+y0] =            bcoef + 0.5*Im*kwave*ccoef                  +0.5*ecoef;
-	
-	rhsk[k][y+y0] = rhs[y+y0][k]; // Transpose
+        a(k, y + y0) = bcoef - 0.5 * Im * kwave * ccoef - 0.5 * ecoef;


The whitespace changes have made the preceding comment not very useful. Should we either preserve the whitespace, or delete/modify the comment?

I'll revert the whitespace.

ZedThree · 2018-02-26T09:48:33Z

include/cyclic_reduction.hxx

  }

  /// Free all memory arrays allocated by allocMemory()
  void freeMemory() {


Is this needed at all now? Can it be deprecated/removed?

It currently also resets the problem size. I don't think we really need this as it should be cheap enough to fully tear down and rebuild a new instance should we need to. May be worth deprecating?

I've marked as deprecated, but this is private so I think we should probably be free to just remove it if we decide that we don't need to reset the problem size or if we want to do this a different way.

Restore whitespace in parderiv/cyclic to keep alignment with terms

May be helpful in the situation where the is only one system to solve at a time. Adds a few ASSERTs in places to further guard.

…/BOUT-dev into remove_matrix_cyclic

Routine is private so removing uses in framework sufficient to remove entirely without breaking compatibility

All comments addressed except the 2nd larger question which may be outside the scope of this PR

ZedThree · 2018-02-28T17:48:08Z

include/cyclic_reduction.hxx

          output << "Expecting receive from " << p << " of size " << len << endl;
 #endif
-	  MPI_Irecv(recvbuffer[p],
+	  MPI_Irecv(&recvbuffer(p, 0),


Can we use std::begin here, similar to how ifp gets passed below? Or is that not possible due to the subscripting?

Yes, I think the indexing prevents this from working with std::begin, as recvbuffer(p,0) refers to a single value I'm not sure what begin would do here.

Most of the places I take the address or use std::begin with all the changes related to arrays, matrices and tensors are places I would like to change the interface of the functions to take the bounds-checked class rather than reverting to passing pointers. I've not made such changes yet as it would touch a lot of signatures and there are a couple of things that we might want to change in the design of Matrix and similar which might change the best route forward.

d7919 added 13 commits February 22, 2018 11:17

Change cyclic_reduction public routines to take Matrix

e146c9c

Instead of matrix (pointer to pointer). This is the first step in allowing the cyclic solvers to replace all static/class level memory.

Replace matrix with Matrix in LaplaceXY

e43eed0

Replace a number of matrices with Matrices in cylcic_reduction

cb21629

One remains unchanged (see #895 for why)

Make sure output is available

49f9037

Replace some raw blocks with Array

81288f3

in cyclic_reduction.hxx

Replace some direct access to Matrix data with () access

46f728e

This ensures bounds checking when compiled with appropriate CHECK

Remove deprecated uses of Matrix indexing in cyclic Laplace

ed27e64

Could now move all storage from class level to method scope

Remove deprecated Matrix indexing in cyclic parderiv

74ff62f

Merge branch 'origin-next' into remove_matrix_cyclic

f36544d

Fix merge

b5b18f8

Replace last matrix with Matrix in cyclic reduction

a6b3595

Remove more deprecated Matrix accesses in cyclic-reduction

c1ebed3

Remove deprecated Matrix access in laplacexz cyclic

6da400a

Also replaces new with Array in laplacexz cyclic

d7919 mentioned this pull request Feb 24, 2018

Remove matrix remainder #901

Merged

bendudson reviewed Feb 24, 2018

View reviewed changes

No longer need to pass in problem size for cyclic reduction

7ef7435

Can get this information from the passed objects.

d7919 mentioned this pull request Feb 24, 2018

Replace new with array #902

Merged

ZedThree previously requested changes Feb 26, 2018

View reviewed changes

d7919 added 7 commits February 26, 2018 15:09

Restore old interfaces in cyclic_reduction but mark as deprecated

87dac6a

Restore whitespace in parderiv/cyclic to keep alignment with terms

Add interface to setCoefs/solve that take Arrays instead of Matrices

50c3314

May be helpful in the situation where the is only one system to solve at a time. Adds a few ASSERTs in places to further guard.

Add interface to setCoefs/solve that take Arrays instead of Matrices

bf9ee50

May be helpful in the situation where the is only one system to solve at a time. Adds a few ASSERTs in places to further guard.

Merge branch 'remove_matrix_cyclic' of https://github.com/boutproject…

0c457f6

…/BOUT-dev into remove_matrix_cyclic

Fix bad merge/commit ammend

70f3760

Fix ASSERTs in cyclic_reduction

41bf552

Remove uses of deprecated freeMemory in cyclic_reduction.

7ee5c81

Routine is private so removing uses in framework sufficient to remove entirely without breaking compatibility

ZedThree reviewed Feb 28, 2018

View reviewed changes

bendudson approved these changes Mar 2, 2018

View reviewed changes

d7919 merged commit 30c83fb into next Mar 2, 2018

d7919 deleted the remove_matrix_cyclic branch March 2, 2018 11:37

Conversation

d7919 commented Feb 23, 2018

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ZedThree left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants