CeedMatrixPseudoinverse utility function

[rezgarshakeri]

See #1250

[jeremylt]

This is looking good. We now need to use this utility in a few places to cut repeated code. I think the two places to update are CeedBasisCreateProjectionMatrices and CeedBasisGetCollocatedGrad.

[rezgarshakeri]

I don't understand the rust failure. How can I fix it?

[jeremylt]

The rust test is saying that the incorrect numerical result was computed, so the pesudoinverse probably isn't right in that case

[rezgarshakeri]

The rust test is saying that the incorrect numerical result was computed, so the pesudoinverse probably isn't right in that case

Hmm, so how the other tests are fine and rust is not? I checked only with make prove in my machine and I got

All tests successful.
Files=232, Tests=2016,  6 wallclock secs ( 1.10 usr  0.34 sys + 44.30 cusr 11.80 csys = 57.54 CPU)
Result: PASS

[jeremylt]

I'd start by seeing if there are any differences between how the Rust test t501 is written vs the C version of the test. It's possible it has a slightly different setup.

The Rust tests allow less memory misbehavior on our part, so I'm inclined to trust them when they say something went wrong.

[rezgarshakeri]

I'd start by seeing if there are any differences between how the Rust test t501 is written vs the C version of the test. It's possible it has a slightly different setup.

The Rust tests allow less memory misbehavior on our part, so I'm inclined to trust them when they say something went wrong.

I see, I checked the pseudoinverse with python. For example t301-basis

  CeedScalar A[12]    = {1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0};
  CeedScalar A_inv[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  CeedCall(CeedMatrixPseudoinverse(ceed, A, 4, 3, A_inv));
  for (CeedInt i =0; i < 12; i++) printf("A_inv[%d]=%f\n", i, A_inv[i]);

gives

A_inv[0]=0.200000
A_inv[1]=0.300000
A_inv[2]=-0.100000
A_inv[3]=0.600000
A_inv[4]=-0.050000
A_inv[5]=0.050000
A_inv[6]=0.150000
A_inv[7]=-0.150000
A_inv[8]=0.125000
A_inv[9]=-0.125000
A_inv[10]=0.125000
A_inv[11]=-0.125000

which is exactly the same as python

array([[ 0.2  ,  0.3  , -0.1  ,  0.6  ],
       [-0.05 ,  0.05 ,  0.15 , -0.15 ],
       [ 0.125, -0.125,  0.125, -0.125]])

so maybe we have memory/size issue in the middle of the code but result is correct.

[rezgarshakeri]

Now the error is clear

Incorrect interval length computed. Expected: 2.0, Found: 1.999999999999999, Error: 1.110223024625e-15

[jeremylt]

Incorrect interval length computed. Expected: 2.0, Found: 1.999999999999999, Error: 1.110223024625e-15

So this gave us a mild degradation in accuracy on this test case, but not enough to be worried. I'd say we bump the tolerance to 5e-15 for this test.

[rezgarshakeri]

Incorrect interval length computed. Expected: 2.0, Found: 1.999999999999999, Error: 1.110223024625e-15

So this gave us a mild degradation in accuracy on this test case, but not enough to be worried. I'd say we bump the tolerance to 5e-15 for this test.

Yeah. I'll update that.

[rezgarshakeri]

@jeremylt Do I need to update

    /// // Check
    /// let sum: Scalar = v.view()?.iter().sum();
    /// assert!(
    ///     (sum - 2.0).abs() < 10.0 * libceed::EPSILON,
    ///     "Incorrect interval length computed"
    /// );

in operator.rs to

    /// // Check
    /// let sum: Scalar = v.view()?.iter().sum();
    /// assert!(
    ///     (sum - 2.0).abs() < 50.0 * libceed::EPSILON,
    ///     "Incorrect interval length computed. Expected: 2.0, Found: {}, Error: {:.12e}"
    /// );

?

[jeremylt]

Probably. And the t501 in main should probably use CeedEpsilon too

[rezgarshakeri]

Probably. And the t501 in main should probably use CeedEpsilon too

I'll make small PR for that.

[jeremylt]

This is needed but may not fix the Rust issue. There shouldn't be any problem with any slice from raw parts and I can't reproduce that locally

$ git diff
diff --git a/examples/rust/ex1-volume/src/main.rs b/examples/rust/ex1-volume/src/main.rs
index ca755438..ed22b9fb 100644
--- a/examples/rust/ex1-volume/src/main.rs
+++ b/examples/rust/ex1-volume/src/main.rs
@@ -247,7 +247,7 @@ fn example_1(options: opt::Opt) -> libceed::Result<()> {
         );
     }
     let tolerance = match dim {
-        1 => 100.0 * libceed::EPSILON,
+        1 => 500.0 * libceed::EPSILON,
         _ => 1E-5,
     };
     let error = (volume - exact_volume).abs();

[jeremylt]

Well, Rust is doing some crappy misbehavior I haven't been able to identify yet

[jeremylt]

Its happening on a different branch too, so the Rust issue isn't related to this branch and is something that happened to the nightly toolchain

[jrwrigh]

Possibly another solution (though I'm not sure how feasible it would be) is to detect whether the projection should result in an identity matrix and, if it should, just return an identity matrix. Determining that ahead of time while not limiting ourselves to nodal basis functions might be a bit challenging. The other option would be verifying identity after matrix-matrix multiply, but that might be just as ad-hoc as forcing small values to zero.

[jrwrigh]

Basically, I discovered that the entries of the projection basis interp matrix were different from main by ~e-16 or so,

I'm a bit surprised by that. Aren't those matrices created by CPU code, not the GPU? Or is the PseudoInverse code also used for generating those projection matrices?

[rezgarshakeri]

This is needed but may not fix the Rust issue. There shouldn't be any problem with any slice from raw parts and I can't reproduce that locally

$ git diff
diff --git a/examples/rust/ex1-volume/src/main.rs b/examples/rust/ex1-volume/src/main.rs
index ca755438..ed22b9fb 100644
--- a/examples/rust/ex1-volume/src/main.rs
+++ b/examples/rust/ex1-volume/src/main.rs
@@ -247,7 +247,7 @@ fn example_1(options: opt::Opt) -> libceed::Result<()> {
         );
     }
     let tolerance = match dim {
-        1 => 100.0 * libceed::EPSILON,
+        1 => 500.0 * libceed::EPSILON,
         _ => 1E-5,
     };
     let error = (volume - exact_volume).abs();

Do I still need this change?

[jeremylt]

The pesudoinverse is use in the core of the library to create the projection matrices. The new computations are in this weird place where they are not different enough to cause appreciable differences on the CPU but they do create differences on the GPU, likely due to slightly different algebra on the machine to apply the basis matrices

[jeremylt]

Do I still need this change?

I think so, yea

[rezgarshakeri]

I don't know why codecov breaks.

[jeremylt]

Its a minor change and codecov is flaky, so I think we can ignore it

[rezgarshakeri]

Cool, so I merge it.

[jrwrigh]

The pesudoinverse is use in the core of the library to create the projection matrices.

Oh, I misread the earlier message to be that the interpolation matrices themselves had slightly different values. This makes more sense.