mls_interpolation with SVD solver #167
Conversation
jacobmerson
left a comment
jacobmerson
left a comment
There was a problem hiding this comment.
long term, I'd like for the internal functions to take Rank1View and Rank2View instead of Kokkos or Omega_h specific types (Read, View, etc). This will allow user to pass in data from application layer easier.
This change does not need to be part of this pull request.
| // the value of 'a' is higher if the data is localized | ||
| // the value of 'a' is smaller if the data is farther | ||
| int a = 20; | ||
| int a = 5; |
There was a problem hiding this comment.
This should be a member of the RBF_GAUSSIAN functor. That way you could provide a default and let the user change it.
| */ | ||
| KOKKOS_INLINE_FUNCTION | ||
| void eval_basis_vector(const IntDeviceMatView& slice_length, const Coord& p, | ||
| void eval_basis_vector(const IntDeviceMatView& slice_length, const double* p, |
There was a problem hiding this comment.
Please use the Rank2View here. That's a mdspan with two indexes tagged for the memory space.
|
|
||
| // initilize U (orthogonal matrices) array | ||
| ScratchMatView U(team.team_scratch(1), row, row); | ||
| fill(-5.0, team, U); |
There was a problem hiding this comment.
Do you need to fill these arrays before use? It seems like this will be doing unnecessary work?
|
|
||
| if (team.team_rank() == 0) { | ||
| KokkosBatched::SerialSVD::invoke(KokkosBatched::SVD_USV_Tag(), matrix, U, | ||
| sigma, Vt, work, 1e-6); |
There was a problem hiding this comment.
Tolerance probably needs to be in the function interface since we will likely need to mess around with that depending on the data set.
| KokkosBatched::SerialGemm< | ||
| KokkosBatched::Trans::Transpose, KokkosBatched::Trans::NoTranspose, | ||
| KokkosBatched::Algo::Gemm::Unblocked>::invoke(1.0, Vt, temp_matrix, 0.0, | ||
| vsigmaInvUtMul); |
There was a problem hiding this comment.
I think Serial operation should be on the team thread 0. Otherwise, you are doing this operation on multiple threads.
| KokkosBlas::SerialGemv< | ||
| KokkosBlas::Trans::NoTranspose, | ||
| KokkosBlas::Algo::Gemv::Unblocked>::invoke(1.0, vsigmaInvUtMul, | ||
| rhs_values, 0.0, |
There was a problem hiding this comment.
same comment as above.
Previously we implemented a normal equation based solver using LU decomposition. The solver didn't work well for ill conditioned systems. The current implementation is based on requalrized SVD decomposition
| @@ -7,18 +7,82 @@ using namespace Omega_h; | |||
|
|
|||
There was a problem hiding this comment.
Don't know how we missed this on the first review, but this needs to be removed. @Fuad-HH already reported errors from this.
The problem with this is you are forcing all downstream code to bring all of Omega_h into it's namespace and that will create naming conflicts. For example with parallel_for
|
Looks like there are some build failures. Possibly missing header files? |
4 participants
This pull request introduces an SVD-based solver to improve robustness in mls interpolation routines. The solver has been tested across a variety of scenarios and performs reliably in most typical cases.
This offers a more stable alternative to existing LU solver that improves numerical stability in extreme cases.