[Megoldva: -ffloat-store] Egyenletrendszer megoldás, i686 vs x86-64

 ( sebist | 2008. június 19., csütörtök - 14:21 )

A lenti program megoldana egy lineáris egyenletrendszert. A megoldás: 1 2 3 4. A kód olyan, amilyen. A probléma az, hogy i686-on rossz lesz az eredmény. Nevezetesen: 1.07407 2 3 3.48148 . x86_64-en és sparc-on jól futott/fordult le. Érdekesség: ha kihagyom az oszlopcserés részt(HIBA innen-HIBA eddig között), akkor jó lesz i686-on is. Jo akkor is, ha 7 helyett 6 es 16 helyett 15 van. Mellékelem a kódot, ha valakinek van egy kis ideje, akkor futtassa le, és írja meg, és hogy működött-e. Írja le hozzászólásban, hogy milyen környezetben használta. Bármilyen ötletet szívesen veszek.

Köszönöm a segítséget:
István

main.cpp

#include<iostream>
#include "gauss.h"

int main()
{
    mat A={{1.0,1.0,1.0,1.0},
           {1.0,3.0,1.0,4.0},
           {7.0,1.0,1.0,1.0},
           {3.0,3.0,1.0,3.0}};
    vect b={10.0,26.0,16.0,24.0};
    vect x;
    gauss(&A,&b,&x);
    for (int i=0;i<N;++i)  std::cout << x[i] <<" ";
    std::cout << std::endl;
    return 0;
}

gauss.h:

#ifndef GAUSS_H
#define GAUSS_H

#include<iostream>
#include<cmath>

#define N 4

typedef double vect[N];
typedef double mat[N][N];


int gauss(mat *A, vect *b, vect *x)
{
    
    int i,j,k;
    int xn[N];
    int mx,my;
    double d;
        
    for(i=0;i<N;i++) xn[i]=i;
    for(i=0;i<N;i++)
    {
        
        // foelem kivalasztas
        mx=i;
        my=i;
        for(j=i;j<N;j++)
        {
            for(k=i;k<N;k++)
            {
                if(fabs((*A)[j][k])>fabs((*A)[my][mx]))
                {
                    my=j;
                    mx=k;                    
                }
            }
        }

        if( my != i) // sorcsere
        {
            for(j=i;j<N;j++)
            {
                d=(*A)[i][j];
                (*A)[i][j]=(*A)[my][j];
                (*A)[my][j]=d; 
            }
            d=(*b)[i];
            (*b)[i]=(*b)[my];
            (*b)[my]=d;        
        } 

/* HIBA innen */
        
       if( mx != i ) // oszlopcsere
        {
            for(j=i;j<N;j++)
            {
                d=(*A)[j][i];
                (*A)[j][i]=(*A)[j][mx];
                (*A)[j][mx]=d;    
            }
            k=xn[i];
            xn[i]=xn[mx];
            xn[mx]=k;
        }
/* HIBA Eddig  */
        if( (*A)[i][i] == 0.0)
        {
            std::cerr << "Gauss Error" << std::endl;
            return 10;
        }
        d=(*A)[i][i];
        
        for(j=i;j<N;j++) (*A)[i][j]=(*A)[i][j]/d;
        (*b)[i]=(*b)[i]/d;
    
        if( i < N-1)
        {
            for(k=i+1;k<N;k++)
            {
                d=(*A)[k][i];
                for(j=i;j<N;j++)
                {
                    (*A)[k][j]=(*A)[k][j]-d*(*A)[i][j]/(*A)[i][i];
                }   
                (*b)[k]=(*b)[k]-d*((*b)[i]);        
            }
        }
    }
    
    for(i=N-2;i>0;i--)
    {
        for(j=N-1;j>i;j--)
        {
            (*b)[i]=(*b)[i]-(*A)[j][i]*(*b)[j];
        }    
    }
    
    for(i=N-2;i>=0;i--)
    {
        d=0.0;
        for(j=N-1;j>i;j--)
        {
            d+=(*A)[i][j]*(*b)[j];
        }    
        (*b)[i]-=d;
    }
    
    for(i=0;i<N;i++)
    {
        (*x)[xn[i]]=(*b)[i];
    }
    return 0;

}

#endif

Hozzászólás megjelenítési lehetőségek

A választott hozzászólás megjelenítési mód a „Beállítás” gombbal rögzíthető.

i686-on

majki@fa ~ $ g++ test01.cpp -o test01
majki@fa ~ $ ./test01 
1.07407 2 3 3.48148 

majki@fa ~ $ g++ test01.cpp -o test01 -msse2 -mfpmath=sse
majki@fa ~ $ ./test01 
1 2 3 4 

érdekes hp-ux 11.11 alatt jó

-bash-3.2$ g++ gauss.cpp
-bash-3.2$ ./a.out
1234

Core2Duo T7100, 2.5G, Ubuntu 8.04, 2.6.24

$ ./gauss
1.07407 2 3 3.48148

$ valgrind ./gauss 2>/dev/null
1 2 3 4

Tehát valami memóriával kapcsolatos probléma áll fent, de a valgrind nem taláta meg (más címet ad vissza mint a sima malloc). Talán az indexeknél van gond?

A gcc 4.1.2 és a 4.3.0 esetén is hibás az eredmény, csak a valgrind "segítségével" jó.

Nálam is ugyanez.
g++ (GCC) 4.2.3 (Ubuntu 4.2.3-2ubuntu7)

Érdekes, hogy a valgrind megoldja a problémát.
Valószínűleg az oszlopcsere kicímez a tömbböl, de ha amúgy is 0 van ott valamiért, akkor nem okoz problémát.

c:\Documents and Settings\saxus\Dokumentumok\Visual Studio 2005\Projects\gauss\debug>gauss.exe
1 2 3 4

Visual Studio 2005, XP SP3, 32 bit, full alap C++ Console Application projekt, debug mód

----

saxus@arrakis ~/gauss                                                           
$ uname -a
FreeBSD arrakis.muportal.hu 6.3-RELEASE FreeBSD 6.3-RELEASE #0: Wed Jan 16 04:18:52 UTC 2008     
root@dessler.cse.buffalo.edu:/usr/obj/usr/src/sys/GENERIC  i386
$ g++ --version
g++ (GCC) 3.4.6 [FreeBSD] 20060305
Copyright (C) 2006 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

$ g++ gauss.cpp -o gauss
$ ./gauss
1 2 3 4

"Hiba az ön készülékében van." Sztem GCC-ddel van valami.

Ez memóriakezelési hiba. Foglalj helyet malloc- vagy calloc al. csak a végén legyen free is.

typedef int klaszterdef [307300][2];

klaszterdef *klaszter;
klaszter = (klaszterdef *) calloc(1, sizeof (klaszterdef));
....
....
....
free ();

Nekem is mátrixműveletekre írt progim van, de más a feladata. Nrm írtam, át a te változóidra, gondolom nem lesz nehéz átírnod.

Mar a hibasnak jelolt resz elott hiba van. (3. alkalom)

a:

 1 0.142857 0.142857 0.142857 
 0 1 0.222222 0.740741 
-0 0 -1.13679e-16 0.666667 
 0 0 0.666667 0.222222 
+0 0 -1.11022e-16 0.666667 
 
-b:2.28571 6.14815 1.33333 2.44444 
+b:2.28571 6.14815 2.44444 1.33333

szerk:
if( i < N-1) {} -en belul keletkezhet az elteres.

szerk2:

(*A)[k][j]=(*A)[k][j]-d*(*A)[i][j]/(*A)[i][i];
std::cout << i <<" "<< j <<" "<< k <<" ered: " << (*A)[k][j]  << std::endl ;

-1 2 3 ered: -1.13679e-16
+1 2 3 ered: -1.11022e-16

Köszönöm mindenki segítségét !

Ezzel a sorral van gond:
if(fabs((*A)[j][k])>fabs((*A)[my][mx]))

Ha beszúrjuk elé pl:

double jk = fabs( (*A)[j][k] );
double yx = fabs( (*A)[my][mx] );
if( fabs( jk - yx ) < EPSILON ) continue;

,ahol EPSILON egy kis szám, akkor már jó lesz.

Üdv:
István

Érdekes hiba. Ilyet még nem láttam.

A hiba akkor jön elő, amikor a két összehasonlított oldal értéke egyaránt 2/3 (0,6666667). Valószínűleg egy kerekítési hiba miatt (talán az átviteljelző bittől is függ a kerekítés) a két oldal nem egyezik meg és pont a bal oldal több. A valgrind mint nyomkövető biztosan törölhetett egy bitet, amitől más lett a kerekítés.

Aha így logikus, ahogy írod. Nekem a kolegám futott ilyen hibába (most mondta), amikor ez egyik függvényillesztőkódomat írta át delphibe, c++-ból. Érdekes hiba egyébként.

máshol van a hiba ..
oliver@lenny:~/develop/gaus$ diff o1 o2
1,2c1,2
< i -1209655308
< j -1208255264
---
> i -1210089484
> j -1208689440
72,74d71
< if(fabs((*A)[j][k])>fabs((*A)[my][mx]))
< ___mx 3
< ___my 3
83c80
< 1.07407 2 3 3.48148
---
> 1 2 3 4


debian gnu/linux @ linux-2.6.22.24-op1 | patch
info

use:
-mfpmath=sse

Ugyan az az erdmeny, mint x86_64 -en.

387 80 bittes prceizitasatol nem varnek ekkora elterest, foleg nem a rossz iranyba :)
BUG szagot erzek.
Lassuk, mit mond a dissassembler. :)

A 80 bit neha 64, ugyanis az FPU regiszter barmikor kikerulhet a memoriaba, es akkor 64 bitre lesz kerekitve. Visszatolteskor ismet 80 bit lesz, de ez ugye nem ugyanaz a 80 bit, mint az eredeti. Ha jol emlexem, van egy -ffloat-store opcio, azzal is lehet jatszani.

Kedves efel!

>A 80 bit neha 64, ugyanis az FPU regiszter barmikor
> kikerulhet a memoriaba, es akkor 64 bitre lesz
> kerekitve.

Visszakérdezek: tudnád ezt bizonyítani is?

Ezt azért kérdem, mert az i80x86-os processzorok assemblerei tudomásom szerint tartalmaznak egy ún. 10 byte-os (80 bit, ,,ten bytes'') adatterület lefoglaló operátort, amelyet ,,direkt'' a lebegőpontos számítások, azaz a FPU-val való kommunikációhoz lehet használni.

Önmagában a 80 bit --> 64 bit --> 80 bit érdekes megoldás, de akkor -- azon felül hogy IEE szabvány is (IEE 754) -- az adatvesztések miatt mi értelme lenne 80 bitnek???

G.
============================================
"Share what you know. Learn what you don't."

see also:
long double @x86

GCC docs (az alahuzas tolem)

* On 68000 and x86 systems, for instance, you can get paradoxical results if you test the precise values of floating point numbers. For example, you can find that a floating point value which is not a NaN is not equal to itself. This results from the fact that the floating point registers hold a few more bits of precision than fit in a double in memory. Compiled code moves values between memory and floating point registers at its convenience, and moving them into memory truncates them.

You can partially avoid this problem by using the `-ffloat-store' option (see section 3.7 Options That Control Optimization).

Rendben, meggyőztél.

G.
============================================
"Share what you know. Learn what you don't."