Revise iteration in PreconditionChebyshev #7703

kronbichler · 2019-02-07T11:35:12Z

The old implementation of PreconditionChebyshev used an iteration that was more complicated than necessary due to a temporary vector representing x^{n-1}-x^{n}. I have rewritten the algorithm directly in terms of the underlying three-term recurrence x^{n+1} = x^{n} + f1 (x^{n}-x^{n-1}) + f2 P^{-1} (b-A x^{n}) where f1 and f2 are some factors. This has the advantage of being slightly faster because we do not need to update the temporary vector and only load from x^{n-1} instead, i.e., we save one vector write, reducing the access in the vector updates from five reads and two writes to five reads (x^n, x^{n-1}, P^{-1}, t=Ax^n, b) and one write to x. Overall, it depends on the cost of the matrix-vector product how much this reduction by around 15% helps in final cost. I measured up to 5% of the solution with a multigrid algorithm, which is not bad.

Furthermore, the new implementation should be more comprehensible to someone familiar with Chebyshev polynomials. I adapted the variable names for the temporary vectors to make things clearer and have also stated the recurrence relation in the class description.

Since the algorithm should be mathematically equivalent to the one we had before, we should not need any test changes (checked on my machine).

kronbichler · 2019-02-07T11:35:24Z

/rebuild

kronbichler · 2019-02-07T11:37:56Z

include/deal.II/lac/precondition.h

-   */
-  void
-  do_transpose_chebyshev_loop(VectorType &dst, const VectorType &src) const;
-


For someone reviewing this part: We had previously put the code for vmult() and step() as well as Tvmult() and Tstep() into these common functions. Given some changes, I would now also have to pass the iteration index. Given that each function only calls two other functions and does an update to a scalar, I have put the code inline because it does not save us any lines in the end.

include/deal.II/lac/precondition.h

jodlbauer · 2019-02-07T15:20:39Z

Probably somehow relevant: in a parallel setting, the initial vector is set to something like 1/norm with 0-th entry explicitly 0.
This destroys the parallel consistency, since entry 0 may not represent the same dof for different number of cores.

include/deal.II/lac/precondition.h

tjhei · 2019-02-07T15:44:27Z

@tcclevenger FYI.

kronbichler · 2019-02-07T16:51:37Z

in a parallel setting, the initial vector is set to something like 1/norm with 0-th entry explicitly 0. This destroys the parallel consistency, since entry 0 may not represent the same dof for different number of cores.

This is something that has bothered me for a while indeed. For deal.II's own vectors, we set something like (-5.5, -4.5, -3.5, .... 3.5, 4.5, 5.5), a vector with zero mean and some offset, that is consistent as long as the numbering of the unknowns is the same. For other vectors we cannot access all entries as easily, which is why we had the norm-type of approach. If the numbering does change different, there is probably little one can do at this abstract level because we have no notion of the numbering or associated DoFHandler.

My advise would be to to pull out the eigenvalue computation to user code. To use this, you should set AdditionalData::eig_cg_n_iterations = 0, AdditionalData::max_eigenvalue to your desired value, and AdditionalData::smoothing_range = max_eigenvalue / min_eigenvalue. I guess we should update the documentation to better describe this procedure, the current version is not really clear.

kronbichler · 2019-02-07T17:01:10Z

One thing we tried to do in the past is to use the right hand side that is given the first time vmult() is called. However, those vectors are often not of sufficiently high frequency, resulting in rather bad eigenvalue estimates. Other libraries use random vectors, but again those are not repeatable at all (in particular not across processor boundaries), so we did not try them.

tamiko · 2019-02-08T23:04:37Z

@kronbichler Please rebase and resolve the merge conflict :-)

kronbichler · 2019-02-09T10:04:55Z

I rebased and also adjusted the Chebyshev degree in a new test that came in in #7696 that was not adjusted in #7708.

tjhei · 2019-02-09T15:28:25Z

@kronbichler , what is this diagonal_matrix.h.2 file thing about?

kronbichler · 2019-02-09T15:38:42Z

what is this diagonal_matrix.h.2 file thing about?

Sorry, left over from some experiments. I removed it.

tjhei · 2019-02-09T16:04:14Z

include/deal.II/lac/precondition.h

+          solution.swap(temp_vector1);
+          solution_old.swap(temp_vector1);
+        }
+      else if (iteration_index > 1)


iteration_index is always positive, which makes this logic somewhat weird (as this check is unnecessary). Is there a reason you don't start with iteration=0? (here and above of course)

wait: that is not true, you also call it with 0. Does it make sense to make an empty

if (iteration_index==0) { // nothing to do here, because .... }

?

Yes, it would make sense to add this part to make it similar to the other if statements cases.

tjhei · 2019-02-09T18:06:52Z

tests/matrix_free/parallel_multigrid_adaptive_08.cc

@@ -342,7 +342,7 @@ do_test(const DoFHandler<dim> &dof)
       ++level)
    {
      smoother_data[level].smoothing_range     = 15.;
-      smoother_data[level].degree              = 5;
+      smoother_data[level].degree              = 6;


this test change does not change the output?

This is necessary after #7708 with the change in what the Chebyshev degree means. That PR was open while the test, added in #7696, was also open. So this simply fixes an only remotely related issue, see also here:
https://cdash.kyomu.43-1.org/viewTest.php?onlydelta&buildid=7762

kronbichler · 2019-02-09T18:32:08Z

/rebuild

kronbichler · 2019-02-12T07:53:10Z

I rebased to re-start the testers that were failing spuriously on a base/timer and base/thread_validity test. Should be ready now.

masterleinad · 2019-02-17T18:04:43Z

include/deal.II/lac/precondition.h

@@ -976,7 +1018,7 @@ class PreconditionChebyshev : public Subscriptor
    /**
     * Constructor.
     */
-    AdditionalData(const unsigned int degree              = 0,
+    AdditionalData(const unsigned int degree              = 1,


This is a follow-up change to #7708, isn't it?

Yes, I had forgotten that change there and one of the exceptions I introduced here revealed this issue in a test.

kronbichler added Linear Algebra ready to test labels Feb 7, 2019

kronbichler commented Feb 7, 2019

View reviewed changes

davydden reviewed Feb 7, 2019

View reviewed changes

include/deal.II/lac/precondition.h Show resolved Hide resolved

tjhei reviewed Feb 7, 2019

View reviewed changes

include/deal.II/lac/precondition.h Show resolved Hide resolved

kronbichler force-pushed the chebyshev_cleanup branch 2 times, most recently from c16245f to 9b51ace Compare February 8, 2019 07:19

kronbichler mentioned this pull request Feb 8, 2019

Adjust degree of PreconditionChebyshev to convention in literature #7708

Merged

kronbichler force-pushed the chebyshev_cleanup branch from 9b51ace to e82b695 Compare February 9, 2019 10:02

kronbichler force-pushed the chebyshev_cleanup branch from e82b695 to b3dc62d Compare February 9, 2019 15:36

tjhei reviewed Feb 9, 2019

View reviewed changes

kronbichler force-pushed the chebyshev_cleanup branch 2 times, most recently from 90901ac to 8bb4172 Compare February 9, 2019 18:17

kronbichler force-pushed the chebyshev_cleanup branch from 8bb4172 to 4845efd Compare February 10, 2019 09:51

kronbichler added 5 commits February 12, 2019 08:51

Clean up iteration of PreconditionChebyshev.

8a6f81c

Add changelog.

63f7956

Augment documentation about eigenvalue computation.

1b3760a

Adjust test to definition of Chebyshev degree.

d3b96e1

Add a comment about iteration 0.

75eec85

kronbichler force-pushed the chebyshev_cleanup branch from 4845efd to 75eec85 Compare February 12, 2019 07:51

masterleinad reviewed Feb 17, 2019

View reviewed changes

masterleinad approved these changes Feb 17, 2019

View reviewed changes

masterleinad merged commit 53adee8 into dealii:master Feb 18, 2019

kronbichler deleted the chebyshev_cleanup branch June 3, 2019 13:32

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Revise iteration in PreconditionChebyshev #7703

Revise iteration in PreconditionChebyshev #7703

kronbichler commented Feb 7, 2019

kronbichler commented Feb 7, 2019

kronbichler Feb 7, 2019

jodlbauer commented Feb 7, 2019

tjhei commented Feb 7, 2019

kronbichler commented Feb 7, 2019

kronbichler commented Feb 7, 2019

tamiko commented Feb 8, 2019

kronbichler commented Feb 9, 2019

tjhei commented Feb 9, 2019

kronbichler commented Feb 9, 2019

tjhei Feb 9, 2019

tjhei Feb 9, 2019 •

edited

kronbichler Feb 9, 2019

tjhei Feb 9, 2019

kronbichler Feb 9, 2019

kronbichler commented Feb 9, 2019

kronbichler commented Feb 12, 2019

masterleinad Feb 17, 2019

kronbichler Feb 17, 2019

Revise iteration in PreconditionChebyshev #7703

Revise iteration in PreconditionChebyshev #7703

Conversation

kronbichler commented Feb 7, 2019

kronbichler commented Feb 7, 2019

kronbichler Feb 7, 2019

Choose a reason for hiding this comment

jodlbauer commented Feb 7, 2019

tjhei commented Feb 7, 2019

kronbichler commented Feb 7, 2019

kronbichler commented Feb 7, 2019

tamiko commented Feb 8, 2019

kronbichler commented Feb 9, 2019

tjhei commented Feb 9, 2019

kronbichler commented Feb 9, 2019

tjhei Feb 9, 2019

Choose a reason for hiding this comment

tjhei Feb 9, 2019 • edited

Choose a reason for hiding this comment

kronbichler Feb 9, 2019

Choose a reason for hiding this comment

tjhei Feb 9, 2019

Choose a reason for hiding this comment

kronbichler Feb 9, 2019

Choose a reason for hiding this comment

kronbichler commented Feb 9, 2019

kronbichler commented Feb 12, 2019

masterleinad Feb 17, 2019

Choose a reason for hiding this comment

kronbichler Feb 17, 2019

Choose a reason for hiding this comment

tjhei Feb 9, 2019 •

edited