Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.c > #138919

Re: storing and retrieving a double may show a difference?

Show all headers | View raw

Manfred <noname@add.invalid> writes:

> On 12/21/2018 6:10 PM, Ben Bacarisse wrote:
>> chanrasjid@gmail.com writes:
>> <snip>
>>> int width_A, height_A;// compute earlier
>>>
>>>
>>> // gaussian elimination for MLR, num_war,etc...
>>>
>>> i 1, x 158161154320.539 != diagonal 158161154320.547
>>> test: mlr-v1.c:561: gauss: Assertion `0' failed.
>>>
>>> I think the first answer abou fpu 80 bits and ram 64 bits may be the
>>> reason.
>>
>> No, it's a straight forward logic error.  You change A[1][1] after
>> calculating diagonal[1] so when you compare the two, they are no longer
>> the same.
> Besides, also by simple code inspection, there appears to be more to
> be corrected - error propagation is not correctly estimated, which
> affects the assertions. More precisely at least two ingredients of
> this are missing:

Yes, but I think the asserts are there only because an error was
spotted.  I'd go further -- it's rarely useful to use any floating point
comparisons in an assert because picking the boundary between OK and not
OK is almost impossible.  Much better to do a full error analysis or to
collect error data as you go and present that as part of the output.  Of
course, if this is some sort of real-time calculation for a game, say,
then you would probably not do either.

> - using DBL_EPSILON
> - considering how many operations are involved in the checked result
> that can affect accuracy (esp. additions and subtractions)
>
> I have marked two examples below
>
>>
>> I hate to say it, but you are very unlikely to get sound advice from fir.
>>
>> I've marked the line where it happens:
>>
> <snip>
>
>>> //reduce UT upper triangle to 0
>>> 	for (i = 0; i < width_A; ++i) {
>>> 		// div col i by leading non-zero element == diag element a-ii
>>> 		x = *(A + i * width_A + i);
>>> 		//diagonal
> /****
> Ben is right that there is a logic error, but still in floating point == and != tests require extra attention: in fact they are more correctly used as /identity/ tests rather than /compare equal/ tests.
> Here this is correct only for architectures that are guaranteed to
> preserve the FP binary representation.

I'm not entirely sure I know what you mean here.  In C, == and != /are/
equality tests and they should ignore any changes or differences in
representation.  For example, negative zero should equal positive zero,
and 1, say, in any bit-length should equal 1 in any other bit-length.

I think I know what you are getting at, but I don't think you've really
explained it, and I won't even try since these are complex matters and I
don't have anything like enough time!

> ****/
>>> 		if (i > 0 && x != diagonal[i]){
>>> 			printf("i %d, x %.15g != diagonal %.15g\n",i, x ,diagonal[i]);
>>> 			assert(0);
>>> 		}
>>> 		for (j = i; j < height_A; ++j) {
>>> 			*(A + j * width_A + i) /= x;
>>> 		}
>>> 		assert(*(A + i * width_A + i) == 1.0);
>>> 		// now diag elem = 1;
>>> 		for (k = i + 1; k < width_A; ++k) {
>>> 			x = *(A + i * width_A + k); // mult
>>> 			for (j = i; j < height_A; ++j) {
>>> 				*(A + j * width_A + k) -= x * *(A + j * width_A + i);
>>
>> **** Here you can have k == 1 and j == 1 and you change the element of A
>> that is on the diagonal.
>>
>>> 			}
> /*****************
> This is another example:
> - The constant 1.0e-300 is arbitrary. The test should be based on DBL_EPSILON an (some estimate of) the size of x.
> - Is this meant to be a test like == 0 or < 0 ? If it is the former, then a fabs() call is missing
> ******************/
>>> 			assert(*(A + i * width_A + k) < 1.0e-300);

I agree about the test, but I'd say get rid of the assert!  Either ignore
small non-zero values, do a more thorough analysis, or accumulate data
on the errors for later reporting.

Only the middle option is really robust.  There are a ill-conditioned
matrices that produce entirely unreliable results using this method, but
there are also well-understood tests for these.

>>> 		}
>>> 	}
>>
>> <snip rest of function>
>>> 	return x;
>>> }
>>
>

-- 
Ben.

Back to comp.lang.c | Previous | Next — Previous in thread | Next in thread | Find similar | Unroll thread

Thread

storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-20 23:25 -0800
  Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-21 00:56 -0800
    Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-21 07:02 -0800
      Re: storing and retrieving a double may show a difference? Bart <bc@freeuk.com> - 2018-12-21 15:33 +0000
        Re: storing and retrieving a double may show a difference? "Rick C. Hodgin" <rick.c.hodgin@gmail.com> - 2018-12-21 11:16 -0500
          Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-21 09:05 -0800
            Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-21 09:43 -0800
            Re: storing and retrieving a double may show a difference? "Rick C. Hodgin" <rick.c.hodgin@gmail.com> - 2018-12-21 12:57 -0500
            Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-21 20:06 +0000
      Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-21 17:10 +0000
        Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-21 09:17 -0800
          Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-21 09:20 -0800
          Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-21 09:24 -0800
        Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-21 21:46 -0800
        Re: storing and retrieving a double may show a difference? Manfred <noname@add.invalid> - 2018-12-22 21:41 +0100
          Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-22 21:08 +0000
            Re: storing and retrieving a double may show a difference? Manfred <noname@add.invalid> - 2018-12-23 20:34 +0100
              Re: storing and retrieving a double may show a difference? Tim Rentsch <txr@alumni.caltech.edu> - 2018-12-28 09:14 -0800
          Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-22 19:16 -0800
      Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-21 09:12 -0800
      Re: storing and retrieving a double may show a difference? James Kuyper <jameskuyper@alumni.caltech.edu> - 2018-12-21 09:23 -0800
        Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-21 20:03 +0000
          Re: storing and retrieving a double may show a difference? James Kuyper <jameskuyper@alumni.caltech.edu> - 2018-12-22 00:10 -0500
            Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-22 11:06 +0000
  Re: storing and retrieving a double may show a difference? mark.bluemel@gmail.com - 2018-12-21 01:00 -0800
    What type of problem is it? (Was: storing and retrieving a double may show a difference?) gazelle@shell.xmission.com (Kenny McCormack) - 2018-12-21 09:10 +0000
      Re: What type of problem is it? (Was: storing and retrieving a double may show a difference?) mark.bluemel@gmail.com - 2018-12-21 01:33 -0800
  Re: storing and retrieving a double may show a difference? Bart <bc@freeuk.com> - 2018-12-21 10:16 +0000
  Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-21 11:25 +0000
  Re: storing and retrieving a double may show a difference? James Kuyper <jameskuyper@alumni.caltech.edu> - 2018-12-21 08:06 -0500
    Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-21 07:13 -0800
      Re: storing and retrieving a double may show a difference? James Kuyper <jameskuyper@alumni.caltech.edu> - 2018-12-21 09:10 -0800
    Re: storing and retrieving a double may show a difference? Tim Rentsch <txr@alumni.caltech.edu> - 2018-12-28 08:16 -0800
  Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-21 08:09 -0800
  Re: storing and retrieving a double may show a difference? "Chris M. Thomasson" <invalid_chris_thomasson@invalid.invalid> - 2018-12-22 17:23 -0800
    Re: storing and retrieving a double may show a difference? "Chris M. Thomasson" <invalid_chris_thomasson@invalid.invalid> - 2018-12-22 17:58 -0800
  Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-22 19:30 -0800
    Re: storing and retrieving a double may show a difference? Richard Damon <Richard@Damon-Family.org> - 2018-12-22 22:48 -0500
      Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-22 20:10 -0800
        Re: storing and retrieving a double may show a difference? "Chris M. Thomasson" <invalid_chris_thomasson@invalid.invalid> - 2018-12-22 20:41 -0800
        Re: storing and retrieving a double may show a difference? Barry Schwarz <schwarzb@dqel.com> - 2018-12-22 22:02 -0800
          Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-23 02:43 -0800
            Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-23 12:38 +0000
              Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-23 05:12 -0800
                Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-23 05:32 -0800
                Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-23 14:55 +0000
            Re: storing and retrieving a double may show a difference? Richard Damon <Richard@Damon-Family.org> - 2018-12-23 07:43 -0500
    Re: storing and retrieving a double may show a difference? Manfred <noname@add.invalid> - 2018-12-23 20:42 +0100

csiph-web