-
Notifications
You must be signed in to change notification settings - Fork 10
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Kahan 2x2 determinant computation reduces numerical imprecision. #239
base: develop
Are you sure you want to change the base?
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If we want to use fma
we should put it in LvArray::math
and wrap it such that it uses std::fma
on host and CUDA's fma
on device.
Also this modification makes calculating the determinant significantly more costly. I'm not sure you can answer this question but maybe the kind of matrices that pop up in the unit tests aren't representative of the normal use case.
src/fixedSizeSquareMatrixOpsImpl.hpp
Outdated
return matrix[ 0 ][ 0 ] * matrix[ 1 ][ 1 ] - matrix[ 0 ][ 1 ] * matrix[ 1 ][ 0 ]; | ||
|
||
// Aliases for matrix [[a, b],[c, d]] for improved readability | ||
auto const & a = matrix[0][0]; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think matrix[ 0 ][ 0 ]
is more readable.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
OK, I will change this.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Also this modification makes calculating the determinant significantly more costly. I'm not sure you can answer this question but maybe the kind of matrices that pop up in the unit tests aren't representative of the normal use case.
It is 3 flop vs 2 flop. Look at the second variant I listed in https://godbolt.org/z/4bGdjYbd5
I don't think this is a big deal...for 3d it may be a bit more costly...but not much. I don't think we will often run into this sort of round off problem, but then again, if it doesn't cost much to "fix", then I kind of like this approach.
src/fixedSizeSquareMatrixOpsImpl.hpp
Outdated
|
||
// Using the more precise Kahan method to compute the 2x2 determinant. | ||
auto const w = b * c; | ||
auto const e = fma( -b, c, w ); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is the fma
required to obtain adequate precision? I'd prefer to do auto const e = -b * c + w
and let the compiler do it's thing.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sure, I'll check too.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
https://godbolt.org/z/ao9GMhPsq
I think you need the fma
I could spot 3x3 dets in the code (FEM stuff), but only 2x2 in the unit tests. EDIT: no 2x2 dets in the code, only in tests. |
I've built an integer dummy fallback in order to prevent integers to be cast to double (I had issues).
Related to issue GEOS-DEV/GEOS#1440: some tests
testTensorOpsInverseTwoArgs
andtestTensorOpsInverseOneArg
are failing due to large numerical imprecision.Using the Kahan method to compute 2x2 determinants brings more precision.
3x3 case still untouched.
Warning, I sometimes get compilation warnings like
which I find surprising since this
fma
function is supposed to bedevice
friendly.fma
also seems pretty standard so one should not need any fallback.So I need to double-check.
Any opinion?