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

Groups > comp.lang.c > #138834 > unrolled thread

storing and retrieving a double may show a difference?

Back to article view | Back to comp.lang.c

Contents

  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

Page 1 of 3  [1] 2 3  Next page →

#138834 — storing and retrieving a double may show a difference?

All values are doubles.
		
        x = *(A + c * width_A + c);
	diagonal_A[c] = x;
        
When diagonal is stored and tested for equality immediately, it tests equal:

c 1, x 103997435021.915, diagonal 103997435021.915
diff 0

But when the diagonal was stored much earlier, an after some blocks of codes, x and diagonal will fail the equality test. 
for (c = 0; c < width_A; ++c) {
		x = *(A + c * width_A + c);
		//diagonal_A[c] = x;
		if (x != diagonal_A[c]) {
			printf("c %d, x %.15g != diagonal %.15g\n", c, x, diagonal_A[c]);
			printf("diff %.15g\n", fabs(x - diagonal_A[c]));
			assert(0);
		}

c 1, x 103997435021.915 != diagonal 103997435022.047
diff 0.131515502929688

Why a slight difference?

[toc] | [next] | [standalone]

#138837

W dniu piątek, 21 grudnia 2018 08:25:55 UTC+1 użytkownik chanr...@gmail.com napisał:
> All values are doubles.
> 		
>         x = *(A + c * width_A + c);
> 	diagonal_A[c] = x;
>         
> When diagonal is stored and tested for equality immediately, it tests equal:
> 
> c 1, x 103997435021.915, diagonal 103997435021.915
> diff 0
> 
> But when the diagonal was stored much earlier, an after some blocks of codes, x and diagonal will fail the equality test. 
> for (c = 0; c < width_A; ++c) {
> 		x = *(A + c * width_A + c);
> 		//diagonal_A[c] = x;
> 		if (x != diagonal_A[c]) {
> 			printf("c %d, x %.15g != diagonal %.15g\n", c, x, diagonal_A[c]);
> 			printf("diff %.15g\n", fabs(x - diagonal_A[c]));
> 			assert(0);
> 		}
> 
> c 1, x 103997435021.915 != diagonal 103997435022.047
> diff 0.131515502929688
> 
> Why a slight difference?

storing and retriving a double may show the difference as far as i know

this is becouse x86 fpu coprocessor registers are 80-bits when you store 
to ram it is stored as 64 bits then if you retrive it you may less precise value in register then

hovever more modern compilers make floating point arithmetic on simd unit
not on fpu as far as i know and as far as
i remember you should not see that difference there (and those more modern means something like AD 2005 and newer,
sorta like that though i dont know exact date when they generally switched)

[toc] | [prev] | [next] | [standalone]

#138860

On Friday, December 21, 2018 at 4:56:33 PM UTC+8, fir wrote:
> W dniu piątek, 21 grudnia 2018 08:25:55 UTC+1 użytkownik > 
> storing and retriving a double may show the difference as far as i know
> 
> this is becouse x86 fpu coprocessor registers are 80-bits when you store 
> to ram it is stored as 64 bits then if you retrive it you may less precise value in register then
> 
> hovever more modern compilers make floating point arithmetic on simd unit
> not on fpu as far as i know and as far as
> i remember you should not see that difference there (and those more modern means something like AD 2005 and newer,
> sorta like that though i dont know exact date when they generally switched)

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.


double gauss() {
	int i, j, k;
	double x, y, diagonal[256];
// augmented matrix A ;  variables x eq
// initialize A for system of equations for gaussian elimination/ MLR.
// must loop r nested within loop c else last row (r == num_var+1) not initialized.
	for (i = 0; i < width_A; ++i) {
		// c  == eq
		for (j = i; j < height_A; ++j) {
//r == variables
			x = 0;
			for (k = 0; k < num_obsv; ++k) {
				x += *(obsv + j * num_obsv + k) * *(obsv + i * num_obsv + k);
			}
			*(A + j * width_A + i) = x;
			if (i == j){
				diagonal[i] = x;
			}
			if (j > i && j < width_A) {
				// symmetry about diagonal; must; && needed, else assign overwrite earlier rows.
				*(A + i * width_A + j) = *(A + j * width_A + i);
			}
		}
	}
//print_A("A before gauss");

//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
		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);
			}
			assert(*(A + i * width_A + k) < 1.0e-300);
		}
	}

//print_A("A after gauss");
// get soln B
	for (i = num_var; i >= 0; --i) {
		// i == col
		B[i] = *(A + (num_var + 1) * width_A + i);
		for (j = i + 1; j <= num_var; ++j) {
			B[i] -= B[j] * *(A + j * width_A + i);
		}
	}

	/* SSE = ∑(y - y-predicted)^2
	 * SSY = ∑ (y - y-mean)^2
	 *
	 * R2 = 1 - SSE / SSY
	 */

// get predicted y obsv
	x = 0;
	for (i = 0; i < num_obsv; ++i) {
		y = B[0];
		for (j = 0; j < num_var; ++j) {
			y += B[j + 1] * *(obsv + (j + 1) * num_obsv + i);
		}
		x += pow(obsv_y[i] - y, 2); // SSE
	}

	x = x / SSY; // 1 - R2
	R2 = 1 - x;
	//assert(R2 >= 0 && R2 <= 1.0);
	if (R2 < 0 || R2 > 1.0) {
		printf("R2 %g\n", R2);
		return HUGE_VAL;
	}

	return x;
}

[toc] | [prev] | [next] | [standalone]

#138863

On 21/12/2018 15:02, chanrasjid@gmail.com wrote:

> // 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.

float has about 23 bits of precision, double 53 and x87 fpu has 64.

The values in your first post matched in about their first 40 bits.

How many calculations approx are performed to arrive the diagonal value?

> double gauss() {


(I tried to run this code but it's missing types and values for width_A, 
height_A, num_obsv, obsv, A, B, num_var, obsv_y, SSY and R2.)

[toc] | [prev] | [next] | [standalone]

#138871

On 12/21/2018 10:33 AM, Bart wrote:
> On 21/12/2018 15:02, chanrasjid@gmail.com wrote:
> 
>> // 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.
> 
> float has about 23 bits of precision, double 53 and x87 fpu has 64.
> 
> The values in your first post matched in about their first 40 bits.
> 
> How many calculations approx are performed to arrive the diagonal value?

Note, the x87 CPU core has a setting that allows a lesser number of bits
to be computed even if the loaded value is 64-bits or 80-bits.  Make sure
those bits are set to allow the full 53 or 64 bits of computation.

>> double gauss() {
> 
> (I tried to run this code but it's missing types and values for width_A, 
> height_A, num_obsv, obsv, A, B, num_var, obsv_y, SSY and R2.)

It's the PC bits:

     Page 12
     http://home.agh.edu.pl/~amrozek/x87.pdf

     00b -- 24-bit precision
     01b -- Reserved
     10b -- 53-bit precision
     11b -- 64-bit precision

-- 
Rick C. Hodgin

[toc] | [prev] | [next] | [standalone]

#138884

On Saturday, December 22, 2018 at 12:15:16 AM UTC+8, Rick C. Hodgin wrote:
> On 12/21/2018 10:33 AM, Bart wrote:
> > On 21/12/2018 15:02, chanrasjid@gmail.com wrote:
> > 
> >> // 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.
> > 
> > float has about 23 bits of precision, double 53 and x87 fpu has 64.
> > 
> > The values in your first post matched in about their first 40 bits.
> > 
> > How many calculations approx are performed to arrive the diagonal value?
> 
> Note, the x87 CPU core has a setting that allows a lesser number of bits
> to be computed even if the loaded value is 64-bits or 80-bits.  Make sure
> those bits are set to allow the full 53 or 64 bits of computation.
> 
> >> double gauss() {
> > 
> > (I tried to run this code but it's missing types and values for width_A, 
> > height_A, num_obsv, obsv, A, B, num_var, obsv_y, SSY and R2.)
> 
> It's the PC bits:
> 
>      Page 12
>      http://home.agh.edu.pl/~amrozek/x87.pdf
> 
>      00b -- 24-bit precision
>      01b -- Reserved
>      10b -- 53-bit precision
>      11b -- 64-bit precision
> 
> -- 
> Rick C. Hodgin

It would be an issue for me as I just rely on what Arch Linux gives me. But I can accept such minor differences.

But can I set full bits in gcc compiler flags?  

[toc] | [prev] | [next] | [standalone]

#138885

W dniu piątek, 21 grudnia 2018 18:25:58 UTC+1 użytkownik chanr...@gmail.com napisał:
> On Saturday, December 22, 2018 at 12:15:16 AM UTC+8, Rick C. Hodgin wrote:
> > On 12/21/2018 10:33 AM, Bart wrote:
> > > On 21/12/2018 15:02, chanrasjid@gmail.com wrote:
> > > 
> > >> // 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.
> > > 
> > > float has about 23 bits of precision, double 53 and x87 fpu has 64.
> > > 
> > > The values in your first post matched in about their first 40 bits.
> > > 
> > > How many calculations approx are performed to arrive the diagonal value?
> > 
> > Note, the x87 CPU core has a setting that allows a lesser number of bits
> > to be computed even if the loaded value is 64-bits or 80-bits.  Make sure
> > those bits are set to allow the full 53 or 64 bits of computation.
> > 
> > >> double gauss() {
> > > 
> > > (I tried to run this code but it's missing types and values for width_A, 
> > > height_A, num_obsv, obsv, A, B, num_var, obsv_y, SSY and R2.)
> > 
> > It's the PC bits:
> > 
> >      Page 12
> >      http://home.agh.edu.pl/~amrozek/x87.pdf
> > 
> >      00b -- 24-bit precision
> >      01b -- Reserved
> >      10b -- 53-bit precision
> >      11b -- 64-bit precision
> > 
> > -- 
> > Rick C. Hodgin
> 
> It would be an issue for me as I just rely on what Arch Linux gives me. But I can accept such minor differences.
> 
> But can I set full bits in gcc compiler flags?

some fpu assembly command stes that, its avaliable in compilers by some library calls, but note this is about fpu, as i asit like from aboput 2005 compilers generate not fpu code but scalard simd one, as far as i remember simd has its own control word you also probably could set..but note you can probably bring
small problems on your head tinkering with that as i heard some libreries do 
not check it and are written in asumptions that this control world(s)
are liek default (thats what i heard not checked it to much)

[toc] | [prev] | [next] | [standalone]

#138886

On 12/21/2018 12:05 PM, chanrasjid@gmail.com wrote:
> On Saturday, December 22, 2018 at 12:15:16 AM UTC+8, Rick C. Hodgin wrote:
>> It's the PC bits:
>>
>>       Page 12
>>       http://home.agh.edu.pl/~amrozek/x87.pdf
>>
>>       00b -- 24-bit precision
>>       01b -- Reserved
>>       10b -- 53-bit precision
>>       11b -- 64-bit precision
> 
> It would be an issue for me as I just rely on what Arch Linux gives me. But I can accept such minor differences.

The differences relate to how accurate the computations are.  If you
use 24-bit precision, you can have losses at 2^(-23) bits, which means
every value beyond the approximately 7th significant digit will be
trash.  If you use 53 bits, then every value beyond the approximately
15th digit will be trash.  And using 64-bit precision will get you
approximately 17-18 significant digits.

> But can I set full bits in gcc compiler flags?

I don't typically use GCC, but IIRC there are intrinsics that allow you
to set the values.  I couldn't find them just now on a quick Google
search.

See "x87 FP Assists" on how these bits are set:

 
https://software.intel.com/en-us/articles/x87-and-sse-floating-point-assists-in-ia-32-flush-to-zero-ftz-and-denormals-are-zero-daz

Here's the assembly code to set those bits in the PC field.  In GCC
you'll have to mangle it to get it to work with the way GCC works.
I'm not an expert in that mangling, so I can't help you there, but
these are the instructions:

     unsigned short x;

     _asm
     {
         // Load whatever it already is
         fstcw   x

         // Mask off all the PC bits
         and     x,0fcffh
         // Note:  Right now it's set to 00b --24-bit precision

         // Use one of these to set bit-level of precision
//      or      x,00200h  // Set to 53-bit precision
//      or      x,00400h  // Set to 64-bit precision

         // Store it back into the FPU
         fldcw   x
     }

Note:  This won't affect direct loads and stores.  It won't be converted.
But if you do any calculation then it will affect it.  To get exact IEEE-
754 rounding when storing values there are some other things you must do.
GCC has a host of flags which follow IEEE-754 rules, or bypass them for
speed.  Depending on your need, you may need to explore them.

-- 
Rick C. Hodgin

[toc] | [prev] | [next] | [standalone]

#138889

chanrasjid@gmail.com writes:

> On Saturday, December 22, 2018 at 12:15:16 AM UTC+8, Rick C. Hodgin wrote:
<snip>
>> Note, the x87 CPU core has a setting that allows a lesser number of bits
>> to be computed even if the loaded value is 64-bits or 80-bits.  Make sure
>> those bits are set to allow the full 53 or 64 bits of computation.
<snip>
> It would be an issue for me as I just rely on what Arch Linux gives
> me. But I can accept such minor differences.

It's not to do with floating point accuracy.  It can't be.  It's just a
programming error.  By the way, the difference is not always small.  By
choosing different data you can get huge differences.

> But can I set full bits in gcc compiler flags?

No need.

-- 
Ben.

[toc] | [prev] | [next] | [standalone]

#138876

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.

I hate to say it, but you are very unlikely to get sound advice from fir.

I've marked the line where it happens:

> double gauss() {
> 	int i, j, k;
> 	double x, y, diagonal[256];
> // augmented matrix A ;  variables x eq
> // initialize A for system of equations for gaussian elimination/ MLR.
> // must loop r nested within loop c else last row (r == num_var+1) not initialized.
> 	for (i = 0; i < width_A; ++i) {
> 		// c  == eq
> 		for (j = i; j < height_A; ++j) {
> //r == variables
> 			x = 0;
> 			for (k = 0; k < num_obsv; ++k) {
> 				x += *(obsv + j * num_obsv + k) * *(obsv + i * num_obsv + k);
> 			}
> 			*(A + j * width_A + i) = x;
> 			if (i == j){
> 				diagonal[i] = x;
> 			}
> 			if (j > i && j < width_A) {
> 				// symmetry about diagonal; must; && needed, else assign overwrite earlier rows.
> 				*(A + i * width_A + j) = *(A + j * width_A + i);
> 			}
> 		}
> 	}
> //print_A("A before gauss");
>
> //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
> 		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.

> 			}
> 			assert(*(A + i * width_A + k) < 1.0e-300);
> 		}
> 	}

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

-- 
Ben.

[toc] | [prev] | [next] | [standalone]

#138878

W dniu piątek, 21 grudnia 2018 18:11:09 UTC+1 użytkownik Ben Bacarisse napisał:
> 
> I hate to say it, but you are very unlikely to get sound advice from fir.
> 
such crappy coments on me never worked, dont work now, and im dare to say will never work (hovever, never say never)(some could say)

[toc] | [prev] | [next] | [standalone]

#138879

W dniu piątek, 21 grudnia 2018 18:18:04 UTC+1 użytkownik fir napisał:
> W dniu piątek, 21 grudnia 2018 18:11:09 UTC+1 użytkownik Ben Bacarisse napisał:
> > 
> > I hate to say it, but you are very unlikely to get sound advice from fir.
> > 
> such crappy coments on me never worked, dont work now, and im dare to say will never work (hovever, never say never)(some could say)

note my answer was an answer to question if storing and reading back may change result (and its interesting that in long double case this is simply the case afair) - i was not answering on bug source as i got no time to read into it

[toc] | [prev] | [next] | [standalone]

#138883

W dniu piątek, 21 grudnia 2018 18:18:04 UTC+1 użytkownik fir napisał:
> W dniu piątek, 21 grudnia 2018 18:11:09 UTC+1 użytkownik Ben Bacarisse napisał:
> > 
> > I hate to say it, but you are very unlikely to get sound advice from fir.
> > 
>(hovever, never say never)

the problem may be im recently feel older and with age it is expected to only get worse... i know its theoreticaly reversable but it needs training to be in good logical form and i dont train to much sadly ;c still hoping one day i will get back

[toc] | [prev] | [next] | [standalone]

#138899

On Saturday, December 22, 2018 at 1:11:09 AM UTC+8, 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.
> 
> I hate to say it, but you are very unlikely to get sound advice from fir.
> 
> I've marked the line where it happens:
> 
> > double gauss() {
> > 	int i, j, k;
> > 	double x, y, diagonal[256];
> > // augmented matrix A ;  variables x eq
> > // initialize A for system of equations for gaussian elimination/ MLR.
> > // must loop r nested within loop c else last row (r == num_var+1) not initialized.
> > 	for (i = 0; i < width_A; ++i) {
> > 		// c  == eq
> > 		for (j = i; j < height_A; ++j) {
> > //r == variables
> > 			x = 0;
> > 			for (k = 0; k < num_obsv; ++k) {
> > 				x += *(obsv + j * num_obsv + k) * *(obsv + i * num_obsv + k);
> > 			}
> > 			*(A + j * width_A + i) = x;
> > 			if (i == j){
> > 				diagonal[i] = x;
> > 			}
> > 			if (j > i && j < width_A) {
> > 				// symmetry about diagonal; must; && needed, else assign overwrite earlier rows.
> > 				*(A + i * width_A + j) = *(A + j * width_A + i);
> > 			}
> > 		}
> > 	}
> > //print_A("A before gauss");
> >
> > //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
> > 		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.
> 
> > 			}
> > 			assert(*(A + i * width_A + k) < 1.0e-300);
> > 		}
> > 	}
> 
> <snip rest of function>
> > 	return x;
> > }
> 
> -- 
> Ben.

Thanks. Yes! The diagonal elements of A starts changing when it starts to subtract one equation with the multiple of another - the initial diagonals of A are all changed. 

It is an IQ question. My brain just missed out thinking the diagonals are not affected.   

[toc] | [prev] | [next] | [standalone]

#138918

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:
- 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.
****/
>> 		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);
>> 		}
>> 	}
> 
> <snip rest of function>
>> 	return x;
>> }
> 

[toc] | [prev] | [next] | [standalone]

#138919

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.

[toc] | [prev] | [next] | [standalone]

#138984

On 12/22/2018 10:08 PM, Ben Bacarisse wrote:
> Manfred <noname@add.invalid> writes:
> 
>> On 12/21/2018 6:10 PM, Ben Bacarisse wrote:
<snip>
>>> 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.
I don't fully agree. I find asserts in floating point math very useful, 
but they must be used in the proper way, which is (the proper way) 
somewhat counter-intuitive in that you must assume there are no /exact/ 
values, which is odd in a binary world.
It is true that picking the exact boundary is quite hard, but still some 
reasonable estimate of rounding errors is possible, and with such, 
asserts can help spot errors - you can't get 100% error detection (but 
no test does), but you can get a significant rate.

   Much better to do a full error analysis or to
> collect error data as you go and present that as part of the output.
This can be an additional debug strategy, yes.

   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)
>>

<snip>
>> /****
>> 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.
You are technically right about this, but my point is about practical 
use of == and != .
In practice, since FP numbers are typically (at least all that need be 
verified) the result of some expression, and they are affected by 
rounding errors, it makes no sense to use == to compare the result of 
two expressions that /in theory/ do compare equal, but in finite 
precision we know they don't.
We know we must use something like (fabs(x-v) < something*DBL_EPSILON), 
and in this perspective the equivalence of positive zero and negative 
zero is not really relevant.
My point about using == as an identity test above is substantially about 
using it to compare a number with itself.
In the example at hand, this means verifying e.g. that an index in an 
array of FP numbers points to the expected element, or that some value 
is not modified between the point it is stored and the point it is 
retrieved (which is the problem you spotted).
The identity test, as a further remark, is guaranteed to work only on 
architectures that do not have the 80/64 bit issue, or something else 
that does not preserve the binary representation of a number between 
storing and retrieving.

> 
> 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!
Yes, I did not explain it, I only recalled the conclusion. Done now, in 
case anyone is interested.

> 
>> ****/
>>>> 		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.

Maybe we see a different purpose for the assert. I think it is 
legitimate to use it in places like this, where it could spot some 
logical errors.
OTOH testing for /accuracy/ of the result, like "is this algorithm 
suitable for the kind of matrices it is used with?" is a different 
matter, and you are right that a simple assert() can't do, while the 
alternative you indicate is more appropriate.

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

[toc] | [prev] | [next] | [standalone]

#139177

Manfred <noname@add.invalid> writes:

> On 12/22/2018 10:08 PM, Ben Bacarisse wrote:
>
>> Manfred <noname@add.invalid> writes:
>>
>>> On 12/21/2018 6:10 PM, Ben Bacarisse wrote:
>
> [snip] [some blank lines added]
>
>>> 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.
>
> You are technically right about this, but my point is about
> practical use of == and != .
>
> In practice, since FP numbers are typically (at least all that need
> be verified) the result of some expression, and they are affected by
> rounding errors, it makes no sense to use == to compare the result
> of two expressions that /in theory/ do compare equal, but in finite
> precision we know they don't.
>
> We know we must use something like (fabs(x-v) <
> something*DBL_EPSILON), and in this perspective the equivalence of
> positive zero and negative zero is not really relevant.
>
> My point about using == as an identity test above is substantially
> about using it to compare a number with itself.
>
> In the example at hand, this means verifying e.g. that an index in
> an array of FP numbers points to the expected element, or that some
> value is not modified between the point it is stored and the point
> it is retrieved (which is the problem you spotted).
>
> The identity test, as a further remark, is guaranteed to work only
> on architectures that do not have the 80/64 bit issue, or something
> else that does not preserve the binary representation of a number
> between storing and retrieving.

In ISO C, assigning a floating-point value to a 64-bit double
must convert any intermediate value (eg, of greater range or
precision) to a 64-bit value, and that 64-bit value must
subsequently be retrieved upon a subsequent access of the
assigned-to object.  That rule applies regardless of whatever
quirks the architecture might have.

Note that the same requirement applies for values that are
explicitly converted (that is, cast) to a particular type.
Thus for example if we have

    double a = ... ;
    double b = ... ;

    if(  (double)(a*b) == 1  )  ...

the result of the multiplication 'a*b' may have additional range
or precision (or both) over the type 'double', but the cast to
'(double)' requires any greater range or precision be discarded
(that is, the value is converted in the same way as if it were
being assigned to a 'double' object) before comparing the result.

[toc] | [prev] | [next] | [standalone]

#138932

On Sunday, December 23, 2018 at 4:41:15 AM UTC+8, Manfred wrote:
> 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:
> - 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.
> ****/
> >> 		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);
> >> 		}
> >> 	}
> > 
> > <snip rest of function>
> >> 	return x;
> >> }
> >

I am aware there is (or should be) a proper way to deal with comparison of floating values or testing for ==0 , etc! But I am not a professional programmer! I think I won't pass to get a job as a programmer - unless the most elementary of project. I use 1.0e-300 to catch logic error.

In my original post, I was thinking along the reasoning as fir; that somehow the double get passed through various memory and the floating representation changes give the difference. If the original error was huge - not just the last few digits - I probably would have found the bug Ben pointed out.  

[toc] | [prev] | [next] | [standalone]

#138877

W dniu piątek, 21 grudnia 2018 16:03:11 UTC+1 użytkownik chanr...@gmail.com napisał:
> On Friday, December 21, 2018 at 4:56:33 PM UTC+8, fir wrote:
> > W dniu piątek, 21 grudnia 2018 08:25:55 UTC+1 użytkownik > 
> > storing and retriving a double may show the difference as far as i know
> > 
> > this is becouse x86 fpu coprocessor registers are 80-bits when you store 
> > to ram it is stored as 64 bits then if you retrive it you may less precise value in register then
> > 
> > hovever more modern compilers make floating point arithmetic on simd unit
> > not on fpu as far as i know and as far as
> > i remember you should not see that difference there (and those more modern means something like AD 2005 and newer,
> > sorta like that though i dont know exact date when they generally switched)
> 
> 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.
> 

i answered on your question if storing and retrieving double may bring the difference (and note im not sure if in practive in c what you usually priint is this 64-bit stripped double, so im not sure if this in practice bring such differences for someone, usually you would use "long double" type probably to see that differences , not all compilers support that)

"The 53-bit significand precision gives from 15 to 17 significant decimal digits precision " - tour erreor is a bit bigger so it seems the reason my be some other maybe you just has normal bug...make sure that what you store is realy not what you read back maybe, as you might probably do some stupid mistake like mistaking index in array or something like that

simply debug it and trace down minimal surprising behaviour

> 
> double gauss() {
> 	int i, j, k;
> 	double x, y, diagonal[256];
> // augmented matrix A ;  variables x eq
> // initialize A for system of equations for gaussian elimination/ MLR.
> // must loop r nested within loop c else last row (r == num_var+1) not initialized.
> 	for (i = 0; i < width_A; ++i) {
> 		// c  == eq
> 		for (j = i; j < height_A; ++j) {
> //r == variables
> 			x = 0;
> 			for (k = 0; k < num_obsv; ++k) {
> 				x += *(obsv + j * num_obsv + k) * *(obsv + i * num_obsv + k);
> 			}
> 			*(A + j * width_A + i) = x;
> 			if (i == j){
> 				diagonal[i] = x;
> 			}
> 			if (j > i && j < width_A) {
> 				// symmetry about diagonal; must; && needed, else assign overwrite earlier rows.
> 				*(A + i * width_A + j) = *(A + j * width_A + i);
> 			}
> 		}
> 	}
> //print_A("A before gauss");
> 
> //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
> 		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);
> 			}
> 			assert(*(A + i * width_A + k) < 1.0e-300);
> 		}
> 	}
> 
> //print_A("A after gauss");
> // get soln B
> 	for (i = num_var; i >= 0; --i) {
> 		// i == col
> 		B[i] = *(A + (num_var + 1) * width_A + i);
> 		for (j = i + 1; j <= num_var; ++j) {
> 			B[i] -= B[j] * *(A + j * width_A + i);
> 		}
> 	}
> 
> 	/* SSE = ∑(y - y-predicted)^2
> 	 * SSY = ∑ (y - y-mean)^2
> 	 *
> 	 * R2 = 1 - SSE / SSY
> 	 */
> 
> // get predicted y obsv
> 	x = 0;
> 	for (i = 0; i < num_obsv; ++i) {
> 		y = B[0];
> 		for (j = 0; j < num_var; ++j) {
> 			y += B[j + 1] * *(obsv + (j + 1) * num_obsv + i);
> 		}
> 		x += pow(obsv_y[i] - y, 2); // SSE
> 	}
> 
> 	x = x / SSY; // 1 - R2
> 	R2 = 1 - x;
> 	//assert(R2 >= 0 && R2 <= 1.0);
> 	if (R2 < 0 || R2 > 1.0) {
> 		printf("R2 %g\n", R2);
> 		return HUGE_VAL;
> 	}
> 
> 	return x;
> }

[toc] | [prev] | [next] | [standalone]

Page 1 of 3  [1] 2 3  Next page →

Back to top | Article view | comp.lang.c

csiph-web