From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: floating point control settings
Date: Wed, 29 Sep 2010 03:49:05 +0000 (UTC)
Organization: Robust Analysis Inc
Lines: 21
Message-ID: <i7ucvh$9mb$>
References: <i7u7se$c4d$> <i7ub0j$3tm$>
Reply-To: <HIDDEN>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: 1285732145 9931 (29 Sep 2010 03:49:05 GMT)
NNTP-Posting-Date: Wed, 29 Sep 2010 03:49:05 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 122612
Xref: comp.soft-sys.matlab:674076

"Roger Stafford" <> wrote in message <i7ub0j$3tm$>...
> "John" <> wrote in message <i7u7se$c4d$>...
> > 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.0e-14).)
> > 
> > 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. 
> > 
> > Questions:
> > (1)  Does anyone know how matlab sets the FP control word?
> > (2)  is there any way to find the value of the control word from matlab?
> > (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.
> - - - - - - - - - -
>   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 precisely-defined "round-to-nearest" 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.
>   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 
> programmers that devise them, and all such variation is capable of producing different results out in the area of least significant bits. 
>   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.
> Roger Stafford

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.