linear sum input order #660

maggul · 2025-02-04T18:10:04Z

An N_VLinearSum input/output issue arises when N_VLinearCombination is not provided by the user or when N_VLinearSum follows a different path than the one followed by SUNDIALS. The output is allowed as an input only in the first component for external library compatibility.

To that end, we have changed the order of inputs if the output vector is in the second component.

N_VLinearSum(a, x, b, z, z) →N_VLinearSum(b, z, a, x, z)

drreynolds · 2025-02-05T04:46:34Z

benchmarks/nvector/test_nvector_performance.c

@@ -2344,7 +2344,7 @@ int Test_N_VScaleAddMultiVectorArray(N_Vector V, sunindextype local_length,
    {
      for (j = 0; j < nsums; j++)
      {
-        N_VLinearSum(c[j], X[k], ONE, Y[j][k], Z[j][k]);
+        N_VLinearSum(ONE, Y[j][k], c[j], X[k], Z[j][k]);


I don't see why this one was changed, since the output Z[j][k] is not one of the inputs.

I made changes whenever I saw (x, z, z). However, it might actually appear as (x, y, z), but y and z could be pointing to the same memory for some other reason. Since I didn’t have the chance to verify this everywhere* while adjusting the code, I wanted to at least reduce the likelihood of the second argument and the output being the same by replacing Y with X, which holds only one vector. If necessary, I can undo the changes and track each vector in the array to see if they ever become identical.

*After making the adjustment, I tested whether (x, z, z) was still hidden somewhere or remained unchanged by forcing "N_VLinearSum" to return without performing its operation and checking the results on DEV_TEST. Since none of the tests failed, this strongly indicates that we haven't returned from "N_VLinearSum" without executing its task.

drreynolds · 2025-02-05T04:56:44Z

src/idas/idas.c

@@ -6934,9 +6934,9 @@ static int IDASensTestError(IDAMem IDA_mem, sunrealtype ck, sunrealtype* err_k,
      if (IDA_mem->ida_kk > 2)
      {
        /* Update error at order k-2 */
-        retval = N_VLinearSumVectorArray(IDA_mem->ida_Ns, ONE,
+        retval = N_VLinearSumVectorArray(IDA_mem->ida_Ns, ONE, tempv, ONE,


So the issue was with both LinearSum and LinearSumVectorArray (since you adjusted the order for N_VLinearSumVectorArray both here and on line 7038)?

Functions like N_VLinearSumVectorArray and several others use N_VLinearSum under the hood. For example, when nvec = 1, the execution directly enters N_VLinearSum. Therefore, this change was necessary to cover such cases.

drreynolds · 2025-02-05T04:58:40Z

src/sundials/sundials_nvector.c

@@ -584,7 +584,7 @@ SUNErrCode N_VScaleAddMulti(int nvec, sunrealtype* a, N_Vector x, N_Vector* Y,
  {
    for (i = 0; i < nvec; i++)
    {
-      x->ops->nvlinearsum(a[i], x, SUN_RCONST(1.0), Y[i], Z[i]);
+      x->ops->nvlinearsum(SUN_RCONST(1.0), Y[i], a[i], x, Z[i]);


Why did this one need to be modified (the output Z[i] is not one of the inputs)?

Similarly, to reduce the probability. Again, I can revert it if necessary and run the DEV_TEST.

drreynolds · 2025-02-05T04:59:01Z

src/sundials/sundials_nvector.c

@@ -773,7 +773,7 @@ SUNErrCode N_VScaleAddMultiVectorArray(int nvec, int nsum, sunrealtype* a,
    {
      for (j = 0; j < nsum; j++)
      {
-        X[0]->ops->nvlinearsum(a[j], X[i], SUN_RCONST(1.0), Y[j][i], Z[j][i]);
+        X[0]->ops->nvlinearsum(SUN_RCONST(1.0), Y[j][i], a[j], X[i], Z[j][i]);


Why did this one need to be modified (the output Z[j][i] is not one of the inputs)?

I made changes whenever I saw (x, z, z). However, it might actually appear as (x, y, z), but y and z could be pointing to the same memory for some other reason. Since I didn’t have the chance to verify this everywhere* while adjusting the code, I wanted to at least reduce the likelihood of the second argument and the output being the same by replacing Y with X, which holds only one vector. If necessary, I can undo the changes and track each vector in the array to see if they ever become identical.

*After making the adjustment, I tested whether (x, z, z) was still hidden somewhere or remained unchanged by forcing "N_VLinearSum" to return without performing its operation and checking the results on DEV_TEST. Since none of the tests failed, this strongly indicates that we haven't returned from "N_VLinearSum" without executing its task.

drreynolds · 2025-02-05T05:01:19Z

test/unit_tests/nvector/test_nvector.c

@@ -2756,7 +2756,7 @@ int Test_N_VLinearSumVectorArray(N_Vector V, sunindextype local_length, int myid
  N_VConst(NEG_ONE, Y[2]);

  start_time = get_time();
-  ierr       = N_VLinearSumVectorArray(3, ONE, X, ONE, Y, Y);
+  ierr       = N_VLinearSumVectorArray(3, ONE, Y, ONE, X, Y);


Again, was there also an issue with N_VLinearSumVectorArray (since you changed it in this file on this line, as well as lines 2802 and 2845)?

I'll note that if there was also an issue with N_VLinearSumVectorArray then you should edit the documentation for that function in doc/shared/nvectors/NVector_Operations.rst to also mention the requirement that the output can equal only the first vector input. If there was not an issue with N_VLinearSumVectorArray then perhaps the edits to those calls in this PR should be reverted?

I think it is really an issue since when nvec = 1, the execution directly enters N_VLinearSum. If we document this convention, then we should also be consistent with the order herein as well, even if no chance is necessary here since nvec = 3.

Steven-Roberts

You mentioned that some external libraries do not support certain orderings of inputs and output. Could you add a comment in the nvector implementations that have this constraint?

Also the changelog will need to be update.

Steven-Roberts · 2025-02-05T20:54:32Z

doc/shared/nvectors/NVector_Operations.rst

@@ -204,7 +204,7 @@ operations below.
   .. math::
      z_i = a x_i + b y_i, \quad i=0,\ldots,n-1.

-   The output vector *z* can be the same as either of the input vectors (*x* or *y*).
+   The output vector *z* can be the same as the input vector (*x*).


Suggested change

The output vector *z* can be the same as the input vector (*x*).

The output vector *z* can be the same as the input vector *x* but not *y*.

@Steven-Roberts: Regarding your comment above, do you think we should add anything else beyond including but not *y* here? I'll also update the changelog to mention this modification. Is there anything else we should consider beyond these two changes?

Can you also add comments in the code indicating why this change was made. It's not entirely clear to me yet. You said external library compatibility, but what external library is incompatible with y==z? Is it one of the NVector implementations, e.g., LAPACK, RAJA, etc? If so, can you add a comment in the relevant operation.

In SUNDIALS, the linear sum operation is defined as follows:
N_VLinearSum(a, x, b, z, z) := (z = a*x + b*z).

However, an external library might implement it differently:
EX_LinearSum(a, x, b, z, z) := (z = a*x, then z += b*z).

After the first step, the second component z no longer refers to the original input value, leading to incorrect results. To prevent this issue, we restrict the output variable z to be an input only in the first component. This way external libraries that follow the alternative approach can keep their existing linear sum.

Thanks for explaining. That's clearer now. I suppose the external library could do some checks to avoid this, but it does seem simpler to just avoid it altogether. It offers some potential for optimization with restrict keyword too. An analogous restriction on SUNMatScaleAdd would be useful for a separate task I'm working on.

balos1 · 2025-02-06T00:53:42Z

examples/arkode/CXX_parallel/ark_heat2D_lsrk_p.cpp

@@ -1418,7 +1418,7 @@ static int SolutionError(sunrealtype t, N_Vector u, N_Vector e, UserData* udata)
  if (flag != 0) { return -1; }

  // Compute absolute error
-  N_VLinearSum(ONE, u, -ONE, e, e);
+  N_VLinearSum(-ONE, e, ONE, u, e);


To be clear, with the new assumption, the old code would still work fine with any existing N_Vector implementation (our own or someone else's), right?

That is, we are imposing a new assumption on how we call N_VLinearSum from our own code, but users could violate the assumption and things would still work.

In SUNDIALS implementations, this modification only changes the order of summation. Instead of computing:
Z[i] = a * X[i] + b * Z[i]
it will now compute:
Z[i] = b * Z[i] + a * X[i]
This change should not affect existing working implementations. Additionally, in its current state, using the output variable as the second input argument does not cause any issues, as long as users rely on N_VLinearSum or follow our approach of performing the linear sum in a single step (adding all terms at once).

To answer @balos1: yes. All existing NVector implementations (either provided by SUNDIALS or user-supplied) will still work. All this PR does is relax one assumption made by how SUNDIALS uses N_VLinearSum, allowing new user-supplied NVectors to be constructed more easily than at present.

balos1 · 2025-02-06T01:00:12Z

CHANGELOG.md

@@ -6,6 +6,10 @@

 ### New Features and Enhancements

+The `N_VLinearSum` and `N_VLinearSumVectorArray` operations now assume that the output 
+can be used as input only in the first component, ensuring compatibility with external 


This seems sort of like we are breaking the API, which would require a new major version. The reason I say "sort of" is because the existing vectors would still work since our change is in the assumption not the implementations. However, the change in the assumption now instructs us to create new vectors that would possibly not work with old code.

I disagree, since nothing about this PR breaks any existing/working code.

The only thing that would "break" is if a user creates an N_Vector that works with SUNDIALS v7.2.1, and then tries to use that same vector with an older version of SUNDIALS. If anything, the documentation update could add a .. versionchanged comment to note that previous SUNDIALS versions required that all vectors support the case where either input could also be the output.

I disagree, since nothing about this PR breaks any existing/working code.

That why I said it was "sort of" a breaking change. It changes the contract of the API, therefore pedantically speaking it is a breaking change to the API. However, in practice, its impact is more like a "new feature".

drreynolds · 2025-02-06T02:39:10Z

You mentioned that some external libraries do not support certain orderings of inputs and output. Could you add a comment in the nvector implementations that have this constraint?

Also the changelog will need to be update.

Mustafa was not referring to external libraries with NVector implementations that are included in SUNDIALS; this comment is relevant to user-supplied NVector interfaces to external libraries (in this case, Gkeyll). I do not see any need to add comments to any NVector implementations in the SUNDIALS repository, only to make the assumptions on the SUNDIALS end clear.

drreynolds

Here are a few suggestions to resolve @balos1's concern about this PR being a kind-of breaking change.

drreynolds · 2025-02-06T02:49:37Z

doc/shared/nvectors/NVector_Operations.rst

@@ -204,7 +204,7 @@ operations below.
   .. math::
      z_i = a x_i + b y_i, \quad i=0,\ldots,n-1.

-   The output vector *z* can be the same as either of the input vectors (*x* or *y*).
+   The output vector *z* can be the same as the input vector *x* but not *y*.


Suggested change

The output vector *z* can be the same as the input vector *x* but not *y*.

The output vector *z* can be the same as the input vector *x* but not *y*.

.. versionchanged:: x.y.z

Prior versions of SUNDIALS assumed that the output vector *z* could be used as

**either** of the input vectors (*x* or *y*).

I think a user could easily interpret this statement as a breaking change. I.e., they may think "oh, to update to sundials X.Y.Z I need to go update my code now to ensure that z is not the same as y". I suggest:

Prior versions of SUNDIALS assumed that the output vector *z* could be used as **either** of the input vectors (*x* or *y*). While all of the `N_Vector` implementations included with SUNDIALS still allow the output vector *z* to be the same as the input vector *x* or *y*, the API no longer requires it.

drreynolds · 2025-02-06T02:52:26Z

doc/shared/nvectors/NVector_Operations.rst

@@ -606,6 +606,8 @@ of its mathematical operations below.
   where *a* and *b* are scalars, :math:`x_j` and :math:`y_j` are
   vectors in the vector arrays *X* and *Y* respectively, and
   :math:`z_j` is a vector in the output vector array *Z*. The operation returns a :c:type:`SUNErrCode`.
+
+   The output vector array *Z* can be the same as the input vector array *X* but not *Y*.


Suggested change

The output vector array *Z* can be the same as the input vector array *X* but not *Y*.

The output vector array *Z* can be the same as the input vector array *X* but not *Y*.

.. versionchanged:: x.y.z

Prior versions of SUNDIALS assumed that the input vector :math:`x` could be the same

as one of the output vectors :math:`z_j`.

doc/shared/nvectors/NVector_Operations.rst

linear sum order update

e477552

maggul marked this pull request as ready for review February 4, 2025 20:20

maggul requested review from balos1, gardner48, cswoodward and drreynolds as code owners February 4, 2025 20:20

drreynolds reviewed Feb 5, 2025

View reviewed changes

Steven-Roberts reviewed Feb 5, 2025

View reviewed changes

changelog and doc update

61b725e

balos1 reviewed Feb 6, 2025

View reviewed changes

drreynolds reviewed Feb 6, 2025

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

linear sum input order #660

linear sum input order #660

maggul commented Feb 4, 2025

drreynolds Feb 5, 2025

maggul Feb 5, 2025

drreynolds Feb 5, 2025

maggul Feb 5, 2025

drreynolds Feb 5, 2025

maggul Feb 5, 2025

drreynolds Feb 5, 2025

maggul Feb 5, 2025

drreynolds Feb 5, 2025

drreynolds Feb 5, 2025

maggul Feb 5, 2025

Steven-Roberts left a comment

Steven-Roberts Feb 5, 2025

maggul Feb 5, 2025

Steven-Roberts Feb 5, 2025 •

edited

Loading

maggul Feb 5, 2025

Steven-Roberts Feb 6, 2025

balos1 Feb 6, 2025 •

edited

Loading

maggul Feb 6, 2025

drreynolds Feb 6, 2025

balos1 Feb 6, 2025

drreynolds Feb 6, 2025

balos1 Feb 6, 2025

drreynolds commented Feb 6, 2025

drreynolds left a comment

drreynolds Feb 6, 2025 •

edited

Loading

balos1 Feb 6, 2025

drreynolds Feb 6, 2025 •

edited

Loading

	The output vector z can be the same as the input vector (x).
	The output vector z can be the same as the input vector x but not y.

-   The output vector *z* can be the same as the input vector *x* but not *y*.
+   The output vector *z* can be the same as the input vector *x* but not *y*.
+   .. versionchanged:: x.y.z
+      Prior versions of SUNDIALS assumed that the output vector *z* could be used as
+      **either** of the input vectors (*x* or *y*).

linear sum input order #660

Are you sure you want to change the base?

linear sum input order #660

Conversation

maggul commented Feb 4, 2025

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

Steven-Roberts left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Steven-Roberts Feb 5, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

balos1 Feb 6, 2025 • edited Loading

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

drreynolds commented Feb 6, 2025

drreynolds left a comment

Choose a reason for hiding this comment

drreynolds Feb 6, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

drreynolds Feb 6, 2025 • edited Loading

Choose a reason for hiding this comment

Steven-Roberts Feb 5, 2025 •

edited

Loading

balos1 Feb 6, 2025 •

edited

Loading

drreynolds Feb 6, 2025 •

edited

Loading

drreynolds Feb 6, 2025 •

edited

Loading