http://www.mathworks.com/matlabcentral/newsreader/view_thread/292643
MATLAB Central Newsreader  floating point control settings
Feed for thread: floating point control settings
enus
©19942014 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Wed, 29 Sep 2010 02:22:06 +0000
floating point control settings
http://www.mathworks.com/matlabcentral/newsreader/view_thread/292643#783426
John
I have a complicated numerical procedure written in C that has been compiled into a library. When I call one function from a short, standalone C program, I get one answer. When I call the same compiled code from matlab (through a mex interface), I get a slightly different answer. The difference is very small, out in the last digit of accuracy, but it causes a later problem that is serious. (I have to take the result of this calculation and plug it into a function with restricted domain; the slight inaccuracy causes that function to blow up, e.g. sqrt(1.0e14).)<br>
<br>
I've spent several hours working on this, and the only guess have as to what is causing this is that matlab sets the floating point control word differently from the standalone program, and this setting makes operations round slightly differently. <br>
<br>
Questions:<br>
(1) Does anyone know how matlab sets the FP control word?<br>
(2) is there any way to find the value of the control word from matlab?<br>
(3) is there any way to find the value of the control word from within C. I am using MinGW, so gcc compiler on Windows.

Wed, 29 Sep 2010 03:15:31 +0000
Re: floating point control settings
http://www.mathworks.com/matlabcentral/newsreader/view_thread/292643#783438
Roger Stafford
"John" <jpnolan.nospam@mathworks.com> wrote in message <i7u7se$c4d$1@fred.mathworks.com>...<br>
> I have a complicated numerical procedure written in C that has been compiled into a library. When I call one function from a short, standalone C program, I get one answer. When I call the same compiled code from matlab (through a mex interface), I get a slightly different answer. The difference is very small, out in the last digit of accuracy, but it causes a later problem that is serious. (I have to take the result of this calculation and plug it into a function with restricted domain; the slight inaccuracy causes that function to blow up, e.g. sqrt(1.0e14).)<br>
> <br>
> I've spent several hours working on this, and the only guess have as to what is causing this is that matlab sets the floating point control word differently from the standalone program, and this setting makes operations round slightly differently. <br>
> <br>
> Questions:<br>
> (1) Does anyone know how matlab sets the FP control word?<br>
> (2) is there any way to find the value of the control word from matlab?<br>
> (3) is there any way to find the value of the control word from within C. I am using MinGW, so gcc compiler on Windows.<br>
         <br>
In my opinion you will not find the cause of your differing computation results in differing floating point rounding procedures. The great majority of computing systems' floating point algorithms adhere to the IEEE 754 Floating Point Standard of 1985, and though four different kinds of rounding are provided in that standard, the majority of all rounding is done in the very preciselydefined "roundtonearest" mode. The net effect of this is that rounding per se on your C system is very likely identical to the rounding done on a matlab system.<br>
<br>
The place to look for different results lies in the functions used to carry out mathematical operations which go beyond simple binary operations. To give a simple example, if we have three numbers, a, b, and c, you can get different results from (a+b)+c than from a+(b+c) in spite of the associative law of addition even though a consistent rounding method is used. Consequently you are all the more likely to get different answers as the result of two long sequences of operations which are supposedly mathematically identical but involve different orders of numerical operations. Also one sine function for example may use a different numerical approximation method than another. One iterative procedure may proceed a few steps further than another. The algorithms used to do the various kinds of computations in complex systems are probably as varied in their immediate detail as are the <br>
programmers that devise them, and all such variation is capable of producing different results out in the area of least significant bits. <br>
<br>
If you have written procedures that will "blow up" as a result of differences in these least bits, it would be advisable to do some further work on them to make them more robust against such diversity. As a rule you will not be able to use other computing systems if you do not do this.<br>
<br>
Roger Stafford

Wed, 29 Sep 2010 03:49:05 +0000
Re: floating point control settings
http://www.mathworks.com/matlabcentral/newsreader/view_thread/292643#783442
John
"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i7ub0j$3tm$1@fred.mathworks.com>...<br>
> "John" <jpnolan.nospam@mathworks.com> wrote in message <i7u7se$c4d$1@fred.mathworks.com>...<br>
> > I have a complicated numerical procedure written in C that has been compiled into a library. When I call one function from a short, standalone C program, I get one answer. When I call the same compiled code from matlab (through a mex interface), I get a slightly different answer. The difference is very small, out in the last digit of accuracy, but it causes a later problem that is serious. (I have to take the result of this calculation and plug it into a function with restricted domain; the slight inaccuracy causes that function to blow up, e.g. sqrt(1.0e14).)<br>
> > <br>
> > I've spent several hours working on this, and the only guess have as to what is causing this is that matlab sets the floating point control word differently from the standalone program, and this setting makes operations round slightly differently. <br>
> > <br>
> > Questions:<br>
> > (1) Does anyone know how matlab sets the FP control word?<br>
> > (2) is there any way to find the value of the control word from matlab?<br>
> > (3) is there any way to find the value of the control word from within C. I am using MinGW, so gcc compiler on Windows.<br>
>          <br>
> In my opinion you will not find the cause of your differing computation results in differing floating point rounding procedures. The great majority of computing systems' floating point algorithms adhere to the IEEE 754 Floating Point Standard of 1985, and though four different kinds of rounding are provided in that standard, the majority of all rounding is done in the very preciselydefined "roundtonearest" mode. The net effect of this is that rounding per se on your C system is very likely identical to the rounding done on a matlab system.<br>
> <br>
> The place to look for different results lies in the functions used to carry out mathematical operations which go beyond simple binary operations. To give a simple example, if we have three numbers, a, b, and c, you can get different results from (a+b)+c than from a+(b+c) in spite of the associative law of addition even though a consistent rounding method is used. Consequently you are all the more likely to get different answers as the result of two long sequences of operations which are supposedly mathematically identical but involve different orders of numerical operations. Also one sine function for example may use a different numerical approximation method than another. One iterative procedure may proceed a few steps further than another. The algorithms used to do the various kinds of computations in complex systems are probably as varied in their immediate detail as are the <br>
> programmers that devise them, and all such variation is capable of producing different results out in the area of least significant bits. <br>
> <br>
> If you have written procedures that will "blow up" as a result of differences in these least bits, it would be advisable to do some further work on them to make them more robust against such diversity. As a rule you will not be able to use other computing systems if you do not do this.<br>
> <br>
> Roger Stafford<br>
<br>
I will repeat: exactly the same binary code is operating on exactly the same input. There are not different functions, or different order of data elements. The only plausible reason I can see is different FP control settings from the calling routine (matlab vs standalone C program). I did design my code to be robust, and built in a threshold to avoid domain problems. My problem here is that I set the threshold after lots of testing under C, only to find out that this does not work when called from matlab. This makes me wary of running routines under matlab.

Wed, 29 Sep 2010 05:03:04 +0000
Re: floating point control settings
http://www.mathworks.com/matlabcentral/newsreader/view_thread/292643#783449
Mark Shore
> > <br>
> > If you have written procedures that will "blow up" as a result of differences in these least bits, it would be advisable to do some further work on them to make them more robust against such diversity. As a rule you will not be able to use other computing systems if you do not do this.<br>
> > <br>
> > Roger Stafford<br>
> <br>
> I will repeat: exactly the same binary code is operating on exactly the same input. There are not different functions, or different order of data elements. The only plausible reason I can see is different FP control settings from the calling routine (matlab vs standalone C program). I did design my code to be robust, and built in a threshold to avoid domain problems. My problem here is that I set the threshold after lots of testing under C, only to find out that this does not work when called from matlab. This makes me wary of running routines under matlab.<br>
<br>
To repeat what Roger said, if your algorithm fails due to the LSB of a floating point number, it is NOT robust. A link recommended from time to time by John D'Errico is <a href="http://docs.sun.com/source/8063568/ncg_goldberg.html">http://docs.sun.com/source/8063568/ncg_goldberg.html</a> or else look at Higham's Accuracy and Stability of Numerical Algorithms.

Wed, 29 Sep 2010 06:01:23 +0000
Re: floating point control settings
http://www.mathworks.com/matlabcentral/newsreader/view_thread/292643#783455
Roger Stafford
"John" <jpnolan.nospam@mathworks.com> wrote in message <i7ucvh$9mb$1@fred.mathworks.com>...<br>
> I will repeat: exactly the same binary code is operating on exactly the same input. There are not different functions, or different order of data elements. The only plausible reason I can see is different FP control settings from the calling routine (matlab vs standalone C program). I did design my code to be robust, and built in a threshold to avoid domain problems. My problem here is that I set the threshold after lots of testing under C, only to find out that this does not work when called from matlab. This makes me wary of running routines under matlab.<br>
        <br>
John, in your "complicated numerical procedure written in C" do you ever do anything other than add, subtract, multiply, or divide just two numbers together in one step? Of course you do! The minute you multiply two matrices together the results in the least bits will depend on just how your C compiler was designed to handle matrix multiplication. If you call on a trigonometric function your results will depend on the kinds of approximations your C library uses in its computation. I very seriously doubt the statement you made: "exactly the same binary code is operating on exactly the same input." It is almost certainly incorrect in the sense that we are talking about here.<br>
<br>
Roger Stafford

Wed, 29 Sep 2010 06:12:07 +0000
Re: floating point control settings
http://www.mathworks.com/matlabcentral/newsreader/view_thread/292643#783459
Bruno Luong
"John" <jpnolan.nospam@mathworks.com> wrote in message <i7u7se$c4d$1@fred.mathworks.com>...<br>
> I have a complicated numerical procedure written in C that has been compiled into a library. When I call one function from a short, standalone C program, I get one answer. When I call the same compiled code from matlab (through a mex interface), I get a slightly different answer. The difference is very small, out in the last digit of accuracy, but it causes a later problem that is serious. (I have to take the result of this calculation and plug it into a function with restricted domain; the slight inaccuracy causes that function to blow up, e.g. sqrt(1.0e14).)<br>
> <br>
> I've spent several hours working on this, and the only guess have as to what is causing this is that matlab sets the floating point control word differently from the standalone program, and this setting makes operations round slightly differently. <br>
> <br>
> Questions:<br>
> (1) Does anyone know how matlab sets the FP control word?<br>
> (2) is there any way to find the value of the control word from matlab?<br>
> (3) is there any way to find the value of the control word from within C. I am using MinGW, so gcc compiler on Windows.<br>
<br>
May be you can try this code that I use for VS<br>
<a href="http://www.mathworks.com/matlabcentral/newsreader/view_thread/274099#719196">http://www.mathworks.com/matlabcentral/newsreader/view_thread/274099#719196</a><br>
<br>
To know Matlab FPCW just query it within a "clean" MEX file.<br>
<br>
FPCW = _controlfp(0,0);<br>
<br>
Bruno

Wed, 29 Sep 2010 06:24:28 +0000
Re: floating point control settings
http://www.mathworks.com/matlabcentral/newsreader/view_thread/292643#783460
Bruno Luong
"John" <jpnolan.nospam@mathworks.com> wrote in message <i7ucvh$9mb$1@fred.mathworks.com>...<br>
<br>
> <br>
> I will repeat: exactly the same binary code is operating on exactly the same input. There are not different functions, or different order of data elements. The only plausible reason I can see is different FP control settings from the calling routine (matlab vs standalone C program). I did design my code to be robust, and built in a threshold to avoid domain problems. My problem here is that I set the threshold after lots of testing under C, only to find out that this does not work when called from matlab. This makes me wary of running routines under matlab.<br>
<br>
Something puzzle me, it your algorithm requires some particular setting of FPCW, why you are bother with Matlab FPCW? Just set it to the evalue the code is designed for.<br>
<br>
In any case I must agree with others have said. An algorithm that can blow due to the uncertainty of the 53th bit does not seem at all robust for any realistic application.<br>
<br>
Bruno

Wed, 29 Sep 2010 06:26:05 +0000
Re: floating point control settings
http://www.mathworks.com/matlabcentral/newsreader/view_thread/292643#783461
James Tursa
"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i7uknj$s3i$1@fred.mathworks.com>...<br>
> "John" <jpnolan.nospam@mathworks.com> wrote in message <i7ucvh$9mb$1@fred.mathworks.com>...<br>
> > I will repeat: exactly the same binary code is operating on exactly the same input. There are not different functions, or different order of data elements. The only plausible reason I can see is different FP control settings from the calling routine (matlab vs standalone C program). I did design my code to be robust, and built in a threshold to avoid domain problems. My problem here is that I set the threshold after lots of testing under C, only to find out that this does not work when called from matlab. This makes me wary of running routines under matlab.<br>
>         <br>
> John, in your "complicated numerical procedure written in C" do you ever do anything other than add, subtract, multiply, or divide just two numbers together in one step? Of course you do! The minute you multiply two matrices together the results in the least bits will depend on just how your C compiler was designed to handle matrix multiplication. If you call on a trigonometric function your results will depend on the kinds of approximations your C library uses in its computation. I very seriously doubt the statement you made: "exactly the same binary code is operating on exactly the same input." It is almost certainly incorrect in the sense that we are talking about here.<br>
> <br>
> Roger Stafford<br>
<br>
To my understanding of OP's first post, he has a compiled library ... only one file. He calls functions in that file from C and gets one answer, and then calls the same functions from the same compiled library with the same inputs in MATLAB and gets a different answer. As far as the library is concerned, it doesn't matter what the compiler did ... it's the same binary file being used by both C and MATLAB ... so where is the difference coming from? My first question would be how has OP verified that the inputs are indeed exactly the same? Are they read from a file or what?<br>
<br>
On a side issue, your statement "... do you ever do anything other than add, subtract, multiply, or divide just two numbers together in one step?" is unnecessarily restrictive. Even if you never did anything other than that in your source code, that would not necessarily stop the compiler from saving temporary results in an extended bit precision register and then using *that* downstream instead of a saved lower precision result. So something as simple as A = B + C followed downstream by D = A + E could still give a different answer on one system vs another if the compiler used an extended register temporary saved value for A in the 2nd expression rather than getting A from memory.<br>
<br>
James Tursa

Wed, 29 Sep 2010 13:37:04 +0000
Re: floating point control settings
http://www.mathworks.com/matlabcentral/newsreader/view_thread/292643#783583
John
"James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <i7um5t$o6$1@fred.mathworks.com>...<br>
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i7uknj$s3i$1@fred.mathworks.com>...<br>
> > "John" <jpnolan.nospam@mathworks.com> wrote in message <i7ucvh$9mb$1@fred.mathworks.com>...<br>
> > > I will repeat: exactly the same binary code is operating on exactly the same input. There are not different functions, or different order of data elements. The only plausible reason I can see is different FP control settings from the calling routine (matlab vs standalone C program). I did design my code to be robust, and built in a threshold to avoid domain problems. My problem here is that I set the threshold after lots of testing under C, only to find out that this does not work when called from matlab. This makes me wary of running routines under matlab.<br>
> >         <br>
> > John, in your "complicated numerical procedure written in C" do you ever do anything other than add, subtract, multiply, or divide just two numbers together in one step? Of course you do! The minute you multiply two matrices together the results in the least bits will depend on just how your C compiler was designed to handle matrix multiplication. If you call on a trigonometric function your results will depend on the kinds of approximations your C library uses in its computation. I very seriously doubt the statement you made: "exactly the same binary code is operating on exactly the same input." It is almost certainly incorrect in the sense that we are talking about here.<br>
> > <br>
> > Roger Stafford<br>
> <br>
> To my understanding of OP's first post, he has a compiled library ... only one file. He calls functions in that file from C and gets one answer, and then calls the same functions from the same compiled library with the same inputs in MATLAB and gets a different answer. As far as the library is concerned, it doesn't matter what the compiler did ... it's the same binary file being used by both C and MATLAB ... so where is the difference coming from? My first question would be how has OP verified that the inputs are indeed exactly the same? Are they read from a file or what?<br>
> <br>
> On a side issue, your statement "... do you ever do anything other than add, subtract, multiply, or divide just two numbers together in one step?" is unnecessarily restrictive. Even if you never did anything other than that in your source code, that would not necessarily stop the compiler from saving temporary results in an extended bit precision register and then using *that* downstream instead of a saved lower precision result. So something as simple as A = B + C followed downstream by D = A + E could still give a different answer on one system vs another if the compiler used an extended register temporary saved value for A in the 2nd expression rather than getting A from memory.<br>
> <br>
> James Tursa<br>
<br>
<br>
Thank you to Bruno for the information on the function _controlfp( ), which comes from float.h. You can read the floating point control word via<br>
unsigned int FPCW;<br>
FPCW = _controlfp(0,0);<br>
printf("FPCW=%012x \n",FPCW);<br>
This does confirm that the issues is rounding. When I run the code from a standalone C (using MinGW flavor of gcc) program I get <br>
FPCW=00000008001f<br>
When I run the same code under matlab using mex, I get<br>
FPCW=00000009001f<br>
<br>
If I manually set the FPCW inside the C program to the the matlab setting, my function gives a NaN as in matlab; if i manually set the FPCW to the C setting inside my mex file, I get the correct value out of the matlab form of the function. <br>
<br>
A few comments on the thread:<br>
(1) yes, I did verify that all inputs are exactly the same<br>
(2) the standard thresholds for determining the smallest double that you can add to 1 and get a result bigger than 1 is DBL_EPSILON in float.h. That value is the same, regardless of the FPCW setting. So, you cannot rely on DBL_EPSILON to help in this kind of problem.<br>
(3) The calculations I am doing are complicated, with thousands of lines of code computing quantities whose values and domains of definition depend on the input parameters, and then doing further calculations based on those results. I designed my code to test for these issues, and built into the code a tolerance that prevented problems. I set that tolerance after a lot of work evaluating the code from C. I would not be surprised if that tolerance was hardware dependent, or compiler dependent, or compiler flag dependent. I was surprised that the same binary code on the same machine gave different answers.<br>
(4) I know I am pushing the limits of double precision calculations on these calculations. I want to get maximum precision with the fastest possible execution time. I can solve the problem by increasing my threshold, but I lose accuracy (for the same compiled code on the same machine) by doing this. And I know I can put in checks in the latter phase of my calculations to do domain checking, but that slows things down. Indeed, I've spent a lot of time in the first phase to set up things precisely so that I didn't have to do constant checking in the second phase. <br>
<br>
I know there is no free lunch. I was just surprised at the way this problem occurred.

Wed, 29 Sep 2010 18:44:33 +0000
Re: floating point control settings
http://www.mathworks.com/matlabcentral/newsreader/view_thread/292643#783730
James Tursa
"John" <jpnolan.nospam@mathworks.com> wrote in message <i7vfe0$pnc$1@fred.mathworks.com>...<br>
> <br>
> Thank you to Bruno for the information on the function _controlfp( ), which comes from float.h. You can read the floating point control word via<br>
> unsigned int FPCW;<br>
> FPCW = _controlfp(0,0);<br>
> printf("FPCW=%012x \n",FPCW);<br>
> This does confirm that the issues is rounding. When I run the code from a standalone C (using MinGW flavor of gcc) program I get <br>
> FPCW=00000008001f<br>
> When I run the same code under matlab using mex, I get<br>
> FPCW=00000009001f<br>
> <br>
> If I manually set the FPCW inside the C program to the the matlab setting, my function gives a NaN as in matlab; if i manually set the FPCW to the C setting inside my mex file, I get the correct value out of the matlab form of the function. <br>
<br>
So it looks like the Precision Control for intermediate values is different between the two. If I am reading the tables correctly (e.g., <a href="http://msdn.microsoft.com/enus/library/e9b52ceh.aspx">http://msdn.microsoft.com/enus/library/e9b52ceh.aspx</a>) the 8001f means intermediate values are extended precision 80bit with 64bit mantissa, whereas the 9001f means intermediate values are double precision 64bit with 53bit mantissa. I don't know the reason why MATLAB would do one method vs the other ... maybe related the the 3rd party libraries they use for BLAS & LAPACK, or maybe something else???<br>
<br>
James Tursa

Wed, 29 Sep 2010 18:52:01 +0000
Re: floating point control settings
http://www.mathworks.com/matlabcentral/newsreader/view_thread/292643#783734
Walter Roberson
On 100929 01:44 PM, James Tursa wrote:<br>
<br>
> So it looks like the Precision Control for intermediate values is<br>
> different between the two. If I am reading the tables correctly (e.g.,<br>
> <a href="http://msdn.microsoft.com/enus/library/e9b52ceh.aspx">http://msdn.microsoft.com/enus/library/e9b52ceh.aspx</a>) the 8001f means<br>
> intermediate values are extended precision 80bit with 64bit mantissa,<br>
> whereas the 9001f means intermediate values are double precision 64bit<br>
> with 53bit mantissa. I don't know the reason why MATLAB would do one<br>
> method vs the other ... maybe related the the 3rd party libraries they<br>
> use for BLAS & LAPACK, or maybe something else???<br>
<br>
IEEE 754 requires that 80/64 *not* be used internally for "double precision", <br>
as it changes the results and there is no rigorous definition of when <br>
intermediate results might be flushed to 64/53 and how that might effect the <br>
calculations.<br>
<br>
A further note along the theme of this topic: as it appears that MS Windows is <br>
being used, then Matlab's system_dependant('setprecision',53) could be used on <br>
the Matlab side to be sure of what one is getting. That call is valid only for <br>
MS Windows, but it would save mucking around with the floating point control <br>
structures.