[C/C++] Gaussian en Gauss-Jordan eliminatie performance

[Verwijderd]

Ik moet voor een project een systeem van lineaire vergelijkingen met 5 onbekenden oplossen. Geschreven in matrices levert dit systeem een vergelijking op in de vorm A * x = B, waarbij A een 5x5 matrix en B een 5x1 matrix is, de onbekenden staan in de 5x1 matrix x. Een dergelijk systeem kan o.a. opgelost worden met Gaussian of Gauss-Jordan eliminatie.

Matrix A wordt opgestelt uit 5 paar data samples uit een AD converter. De samples worden vooraf gecontroleerd op afhankelijkheid. Als de samples niet afhankelijk zijn, en de matrix vergelijking dus oplosbaar is, wordt de matrix opgesteld en de berekening uitgevoerd door een TMS320C6713 floating-point DSP van TI. De uitkomst van de matrix moet berekend zijn voordat de volgende samples door de AD converter aangeleverd worden en de rekentijd voor het oplossen van de matrix moet tot een minimum beperkt worden.

Volgens verschillende bronnen op internet (o.a. PlanetMath) zou een implementatie van gaussian eliminatie sneller zijn. Om dit te testen heb ik 2 functies geschreven om de executietijd te kunnen vergelijken. Bij beide implementaties is SamplesMatrix een 5x6 array, waarbij matrix A in de eerste 5 kolommen staat en matrix B in de laatste kolom. Na afloop van de functie staan in de laatste kolom de antwoorden van de 5 onbekenden.

C:

void PerformGaussian (float (*SamplesMatrix)[6])
{
    float Temp;
    
    // Walk through columns
    for (INT8 i = 0; i < 5; i++)
    {
        // Walk through rows beneath working row number
        for (INT8 j = i + 1; j < 5; j++)
        {
            // Walk through current column and those on the right it 
            for (INT8 k = 5; k >= i; k--)
            {
                SamplesMatrix[j][k] -= SamplesMatrix[i][k] *
                    (SamplesMatrix[j][i] / SamplesMatrix[i][i]);
            }
        }
    }
    
    // Back substitution
    for (INT8 j = 4; j >= 0; j--)
    {
        Temp = 0;
        for (INT8 k = j + 1; k < 5; k++)
        {
            Temp += SamplesMatrix[j][k] * SamplesMatrix[k][5];
        }
        SamplesMatrix[j][5] = (SamplesMatrix[j][5] - Temp) / SamplesMatrix[j][j];
    }
}

C:

void PerformGaussJordan (float (*SamplesMatrix)[6])
{
    float Calculator;
    
    // Walk through the rows of the matrix
    for (INT8 i = 0; i < 5; i++)
    {
        // Store working point value
        Calculator = SamplesMatrix[i][i];
        
        // Walk through columns in current row
        for (INT8 j = 0; j < 6; j++)
        {
            // Divide sample value by working point value.
            // Results in a 1 in the working point. 
            SamplesMatrix[i][j] /= Calculator;
        }
        
        // Walk through the rows of the matrix...
        for (INT8 j = 0; j < 5; j++)
        {
            // ...except the current row
            if (j != i)
            {
                // Store value in working row working point column
                Calculator = SamplesMatrix[j][i];

                // Walk through columns in working row 
                for (INT8 k = 0; k < 6; k++)
                {
                    SamplesMatrix[j][k] -= Calculator * SamplesMatrix[i][k];
                }
            }
        }
    }
}

Tegen verwachting in is de exectutietijd van de PerformGaussJordan functie korter (11843 klokcyles op de c6713) dan PerformGaussian functie (13717 klokcycles op de c6713). Waarschijnlijk zullen de geneste for-loops bij de Gaussian eliminatie voor veel vertraging zorgen, echter ik kan zelf geen oplossing bedenken waar deze for-loops, of een andere vorm van loop, niet nodig zijn.

Ziet een van jullie waar de performence van de Gaussian eliminatie verbeterd kan worden? Deze implementatie zou namelijk sneller moeten zijn. Als jullie nog verbeterpunten zien in de PerformGaussJordan functie of misschien een heel ander algoritme weten dan hoor ik dit ook erg graag, de huidige 11843 klokcycles is namelijk nog veel te veel!

[Infinitive]

Je hebt een stelsel voor máár 5 variabelen. Dan kan je denk ik weinig zeggen over het algoritme zelf. Je kan hooguit gaan mierenneuken over de implementatie

Trouwens, een stelsel van 5 variabelen kan je waarschijnlijk nog wel met de hand uitwerken en implementeren.

[ Voor 27% gewijzigd door Infinitive op 16-08-2004 09:52 ]

[T.T.]

Inderdaad, je moet dergelijke algoritmes niet vergelijken op martrices van deze orde. Probeer eens matrices van 100 keer zo groot. En dan liefst nog verschillende matrices zoals een dense en een sparse matrix. De uitkomsten moet je dan vergelijken met de typen matrices die jij moet oplossen.

Als je altijd dergelijk kleine matrices wilt oplossen, dan zou ik me helemaal geen zorgen maken over de ene of andere implementatie.

[ Voor 19% gewijzigd door T.T. op 16-08-2004 10:03 ]

[Verwijderd]

Bedankt voor jullie reactie Infinitive & T.T. Ik had er even geen rekening mee gehouden dat ik idd met een kleine matrix aan het rekenen ben en dat de performance voor grote matrices misschien heel anders ligt. Echter het is idd de bedoeling om altijd met dergelijk kleine matrices te rekenen en het betreft altijd een dense matrix.

Het oplossen van deze matrix vergelijking maakt onderdeel uit van een digitaal filter welke fouten in de ingangssignalen moet berekenen en corrigeren aan de hand van de samples van de AD converter. De samples worden aangeleverd met 16kHz en het heeft de voorkeur om voor elk nieuwe sample de correctiefactoren te bereken. Dit geeft dus iets meer dan 60us om de matrix op te stellen, door te rekenen en de ingangswaarden te compenseren. Deze 60us komt neer op 12000 klokcycles terwijl ik al 11843 cycles nodig heb om alleen de matrix door te rekenen.

Als er door te mierenneuken op de implementatie performence winst te halen is hoor ik dit dus graag.

[Zoijar]

Misschien dat unrollen (met een template meta program) wel voor snelheids winst kan zorgen als het om ordes als 5x5 gaat.

[klinz]

Zoijar schreef op 16 augustus 2004 @ 12:32:
Misschien dat unrollen (met een template meta program) wel voor snelheids winst kan zorgen als het om ordes als 5x5 gaat.

Met unrollen moet je altijd een beetje uitkijken. Dit is niet per definitie sneller (vanwege codecache), even benchmarken dus voordat je daar op overstapt.

[Zoijar]

Er staat toch ook "misschien"

[Confusion]

Verwijderd schreef op 16 augustus 2004 @ 11:57:
Het oplossen van deze matrix vergelijking maakt onderdeel uit van een digitaal filter welke fouten in de ingangssignalen moet berekenen en corrigeren aan de hand van de samples van de AD converter. De samples worden aangeleverd met 16kHz en het heeft de voorkeur om voor elk nieuwe sample de correctiefactoren te bereken.

Dat heeft natuurlijk de voorkeur, maar is het ook echt noodzakelijk? Ik kan me voldoende situaties voorstellen waarin slechts elke vierde sample corrigeren geen significante vergroting van de fout oplevert (ten opzichte van elk sample corrigeren), maar dat ligt heel erg aan de vorm van het signaal. Ik zou dat zeker even narekenen en toetsen aan de eisen.

Verder kan je veel optimaliseren als je iets weer over de te verwachten vorm van de matrices, maar in de meest algemene vorm kan je slechts optimaliseren door slim gebruik van datastructuren en operaties of eventueel een andere taal (assembler?). Over het laatste weet ik niets; over het eerste zou ik nog iets zinnigs kunnen zeggen.

[klinz]

Zoijar schreef op 16 augustus 2004 @ 12:39:
Er staat toch ook "misschien"

Er stond geen argument bij, daarom.

[Macros]

Mierenneuk optimalisaties:

- Ik zou minimaal de innerloops unrollen als dat mogelijk is. Als je denkt dat dat niet mogelijk is gebruik dan een switch. Die forloop indices worden dan niet meer bijgehouden wat aardig wat cycles scheelt. Ook verminderd dat de branches. Kijk wel uit met je cache, teveel unrollen kan je code cache over laten lopen want een enorme penalty oplevert.
- Vaak is vergelijken met 0 sneller dan met een constant of variable. Het is dan ook een klein beetje sneller om van n naar 0 te lopen dan andersom.

Een beginnetje:

C:

void PerformGaussJordan (float (*SamplesMatrix)[6]) 
{ 
    float Calculator; 
     
    // Walk through the rows of the matrix 
    // --i is sneller (in brakke compilers) dan i-- , 
    // omdat bij i-- de huidige waarde bewaard moet blijven
    for (INT8 i = 4; i >= 0; --i) 
    { 
        // Store working point value 
        Calculator = SamplesMatrix[i][i]; 
         
        // Divide sample value by working point value. 
        // Results in a 1 in the working point.  
        SamplesMatrix[i][0] /= Calculator; 
        SamplesMatrix[i][1] /= Calculator; 
        SamplesMatrix[i][2] /= Calculator; 
        SamplesMatrix[i][3] /= Calculator; 
        SamplesMatrix[i][4] /= Calculator; 
        SamplesMatrix[i][5] /= Calculator; 
         
        // Walk through the rows of the matrix... 
        for (INT8 j = 4; j >= 0; --j) 
        { 
            // ...except the current row 
            if (j != i) 
            { 
                // Store value in working row working point column 
                Calculator = SamplesMatrix[j][i]; 

                // Walk through columns in working row  

                SamplesMatrix[j][0] -= Calculator * SamplesMatrix[i][0]; 
                SamplesMatrix[j][1] -= Calculator * SamplesMatrix[i][1]; 
                SamplesMatrix[j][2] -= Calculator * SamplesMatrix[i][2]; 
                SamplesMatrix[j][3] -= Calculator * SamplesMatrix[i][3]; 
                SamplesMatrix[j][4] -= Calculator * SamplesMatrix[i][4]; 
                SamplesMatrix[j][5] -= Calculator * SamplesMatrix[i][5]; 
            } 
        } 
    } 
}

[Verwijderd]

Macros schreef op 16 augustus 2004 @ 13:03:
Mierenneuk optimalisaties:

- Ik zou minimaal de innerloops unrollen als dat mogelijk is. Als je denkt dat dat niet mogelijk is gebruik dan een switch. Die forloop indices worden dan niet meer bijgehouden wat aardig wat cycles scheelt. Ook verminderd dat de branches. Kijk wel uit met je cache, teveel unrollen kan je code cache over laten lopen want een enorme penalty oplevert.
- Vaak is vergelijken met 0 sneller dan met een constant of variable. Het is dan ook een klein beetje sneller om van n naar 0 te lopen dan andersom.

Ik heb de code even gesimuleerd. Het tellen van n naar 0 alsmede het veranderen van i-- naar --i had geen effect. Nu kan ik de huidige compiler van TI ook alles behalve brak noemen, maar toch bedankt voor de tip.

Het unrollen van van de binnenste for-loops had wel degelijk effect. De simulator van TI geeft aan dat de code gepost door Marcos nog slechts 8575 klok cycles nodig heeft, verder worden er nog geen problemen met de cache aangegeven. Ik kan de code echter nog niet emuleren/benchmarken op de DSP zelf, maar het ziet er veelbelovend uit.

Nooit geweten dat mierenneuken je een performance winst van 27% op kan leveren Mucho Gracias!

[.oisyn]

[.edit: bliep, onderstaande reactie is gebaseerd op huidige moderne CPU's, ik had er even overheen gelezen dat het om een TI geval ging . Voor zover het niet relevant blijkt kun je het als niet gepost beschouwen ]

Macros schreef op 16 augustus 2004 @ 13:03:

- Ik zou minimaal de innerloops unrollen als dat mogelijk is.

Mwoa, op zich is dat iets wat een compiler voor je zou moeten doen als dat echt sneller is (en de compiler heeft daar over het algemeen wat meer verstand van). Ervanuitgaande dat er natuurlijk een degelijke compiler gebruikt wordt

- Vaak is vergelijken met 0 sneller dan met een constant of variable. Het is dan ook een klein beetje sneller om van n naar 0 te lopen dan andersom.

Dit is een mierenneuk optimalisatie die bovendien zelfs meer vertraging op kan leveren dan dat het een echte winst geeft. Ten eerste is het verschil tussen een vergelijking tegen 0 en tegen een andere waarde echt minimaal, in veel gevallen scheelt het slechts een enkele instructie, en vaak is het ook gewoon gelijk. Ten tweede is een loop omkeren terwijl je in die loop array access doet aan de hand van het tellertje echt evil: je kampt continu met cache misses bij reads, en je memory writes zijn niet meer sequential.

Wat ik overigens als (minimale, maar globale) optimalisatie zou willen aanhalen is dat er gewoon gebruik wordt gemaakt van voor de CPU native types. Een INT8 is natuurlijk mooi als je getallen toch niet buiten het bereik van -128..127 vallen, maar over het algemeen is de cpu/het instructieset veel beter geoptimaliseerd voor het gebruik van types die gelijk zijn aan de registergrootte... Over het algemeen betekent dat dus gewoon het gebruik van int

[ Voor 10% gewijzigd door .oisyn op 16-08-2004 14:19 ]

[Zoijar]

Hmm je kan proberen iets te vinden om die'if (j != i)' op regels 26 te omzeilen. Misschien toch alles unrollen. Ik weet niet precies hoe het op jouw architectuur zit, maar bij de moderne intel/amd's kan die conditional zorgen voor branch prediction failures, en dan moet al je cache geflushed worden. Aangezien de check nogal willekeurig is (steeds 4 keer wel, 1 keer niet, op een ander punt) denk ik niet dat dat goed voorspeld wordt.

[Macros]

Elke keer als je een mogelijke optimalisatie doorvoert moet je hem natuurlijk testen. Veel optimalisaties werkten vroeger wel, maar worden tegenwoordig door de compiler gedaan, bijvoorbeeld statische loop unrolling. Dynamische unrolling wordt hier niet gedaan, maar dat kunnen compilers denk ik niet.
Nog een mogelijke optimalisatie:
Ook is het sneller om pointers te gebruiken ipv. array indexering. Elke keer als ie een array indexeerd moet hij het begin van de (array addres) + (size van een element) * (index) doen. Met een pointer index hoef je alleen maar een (pointer) + (size van element). Dat scheelt 1 operatie. In hoeverre dat echt helpt moet je testen, grote kans dat de compiler zelf ook al en optimalisatie toepast.

[Soultaker]

Ik vind het wel logisch dat de implementatie van de Gauss-Jordan-eliminatie sneller is, aangezien die een complexiteit van O(N²) heeft, terwijl de implementatie van de Gauss-eliminatie een complexiteit van O(N³) heeft. Ik zou gevoelsmatig ook altijd voor Gauss-Jordan-eliminatie kiezen.

Overigens gaat deze implementatie van de Gauss-Jordan-eliminatie er van uit dat er op de diagonaal van de matrix geen nullen voorkomen (anders wordt Calculator nul en krijg je deling door nul). In het algemeen kan dat best voorkomen terwijl het systeem toch lineair onafhankelijk is en dan moet je eerst wat rijen verwisselen. Doe je die stap al eerder (bij de controle op lineaire onafhankelijkheid misschien) of is het toeval dat het steeds goed gaat?

[Zoijar]

Probeer deze code eens. Ik heb het niet kunnen testen, dus het zou kunnen dat ik iets over het hoofd heb gezien dat nog steeds zorgt voor runtime calls etc. Maar ik denk dat dit wel sneller zou moeten draaien. (misschien nog even inline toevoegen voor alle functies. Hoewel dat wel al automatisch moet gebeuren aangezien het class body definitions zijn)

C++:

typedef float** Matrix;

template <int dim, int i, int j, int k>
struct Cols2 {
    static void doit(Matrix m, float calc) {
        m[j][k] -= calc * m[i][k];

        // recurse to k+1
        Cols2<dim, i, j, k+1>::doit(m, calc);
    }
};

template <int dim, int i, int j>
struct Cols2<dim, i, j, dim> {
    static void doit(Matrix m, float calc) {
    }
};

////////
//
// general
template <int dim, int i, int j>
struct Rows2 {
    static void doit(Matrix m) {
        // if i!=j
        float calc = m[j][i];
        Cols2<dim+1, i, j, 0>::doit(m, calc);

        // recurse to j+1
        Rows2<dim, i, j+1>::doit(m);
    }
};

// when i==j
template <int dim, int i>
struct Rows2<dim, i, i> {
    static void doit(Matrix m) {
        // if i==j do nothing, and proceed to next
        Rows2<dim, i, i+1>::doit(m);
    }
};

template <int dim, int i>
struct Rows2<dim, i, dim> {
    static void doit(Matrix m) {
    }
};

////////
//
template <int dim, int i, int j>
struct Cols {
    static void doit(Matrix m, float calc) {
        m[i][j] /= calc;

        // recurse to j+1
        Cols<dim, i, j+1>::doit(m, calc);
    }
};

template <int dim, int i>
struct Cols<dim, i, dim> {
    static void doit(Matrix m, float calc) {
    }
};

////////
//
template <int dim, int i>
struct Rows {
    static void doit(Matrix m) {
        float calc = m[i][i];
        Cols<dim+1, i, 0>::doit(m, calc);

        Rows2<dim, i, 0>::doit(m);

        // recurse to i+1
        Rows<dim, i+1>::doit(m);
    }
};

template <int dim>
struct Rows<dim, dim> {
    static void doit(Matrix m) {
    }
};

////////
// Convenience function
template <int dim>
void PerformGaussJordan(Matrix m) {
    Rows<dim, 0>::doit(m);
}

int main(int argc, char* argv[]) {
    Matrix m; // load some matrix here to test, now it's uninitialized!

    // Perform the algorithm for a 5x5 matrix (6 columns)
    PerformGaussJordan<5>(m);
}

Ok, ik heb het getimed, en dit is niet sneller Waarschijnlijk dat de outerloop niet unrollen scheelt, het wordt nu te groot. Verder zit het meeste tijd bij het origineel in die deling, die kan je als 1/x vermenigvuldiging doen. Ik zal het even aanpassen en weer timen.

Goed... bij mij maakt het allemaal niets uit. VC7 doet namelijk idd al de unloop optimalisatie die Macros ook hiervoor al posten; en dus ook die deling opt. Maar aangezien jij wel snelheids verschil zag tussen dat en het origineel, blijkt dat jouw compiler dat waarschijnlijk niet doet. En dus zou het bij jou wel voor snelheids winst kunnen zorgen.

Nope. Zonder outerloop unrolling wordt het ongeveer even snel als macros oplossing (logisch, het is vrijwel hetzelfde). Dus ik denk dat je er niet meer snelheid uit kan halen qua code.

[ Voor 41% gewijzigd door Zoijar op 16-08-2004 17:14 ]

[Verwijderd]

Je moet er natuurlijk ook rekening mee houden dat een DSP architectuur heel anders is dan een 80x86 architectuur. Het optimaliseren verschilt dan ook tussen deze platformen. Vermenigvuldigingen en matrix bewerkingen zullen op een DSP veel sneller worden uitgevoerd dan op een PC.
Het zou dus kunnen dat het ene algoritme beter bij een DSP architectuur past en het andere beter bij een PC architectuur. Daarnaast is de performance van de algoritmes erg afhankelijk zijn de matrix dimensies.

[Verwijderd]

Soultaker schreef op 16 augustus 2004 @ 15:52:
Overigens gaat deze implementatie van de Gauss-Jordan-eliminatie er van uit dat er op de diagonaal van de matrix geen nullen voorkomen (anders wordt Calculator nul en krijg je deling door nul). In het algemeen kan dat best voorkomen terwijl het systeem toch lineair onafhankelijk is en dan moet je eerst wat rijen verwisselen. Doe je die stap al eerder (bij de controle op lineaire onafhankelijkheid misschien) of is het toeval dat het steeds goed gaat?

Ik controleer de data samples, voordat ik er een matrix mee maak, of ze 0 zijn. Als dit het geval is neem ik de sample niet mee in de berekening, aangezien een waarde 0 door de vorm van de ingangssignalen en door de AD converter niet mogelijk is. Ik weet dus vooraf zeker dat erg geen enkele 0 in de matrix voor komt.

Ik zal hier later nog iets van een foutmelding of interrupt aan koppelen want als ik overwachts wel een 0 binnen krijg is hoogst waarschijnlijk de AD converter defect.

[Infinitive]

Ik controleer de data samples, voordat ik er een matrix mee maak, of ze 0 zijn. Als dit het geval is neem ik de sample niet mee in de berekening, aangezien een waarde 0 door de vorm van de ingangssignalen en door de AD converter niet mogelijk is. Ik weet dus vooraf zeker dat erg geen enkele 0 in de matrix voor komt.

Dat geeft geen enkele garantie dat tijdens het elimineerproces geen 0 onstaat! Waarschijnlijk heb je door de onnauwkeurigheid van floats nog nooit een echte 0 gehad, zodat het niet op valt. En dat geeft meteen probleem twee: als je op een bijna-nul cel pivoteert, dan kan je verschrikkelijke afrondfouten onnauwkeurigheid krijgen, onafhankelijk van of hij nu wel of niet echt 0 is.

Als je aanneemt dat dat laatste niet zo vaak voorkomt en afrondfouten niet zo'n probleem zijn, dan kan je als oplossing nemen dat als de waarde in de kolom kleiner is dan een bepaalde waarde, dat je het dan als 0 aanziet. Zodra je een 0 tegenkomt, dan kan je er niet mee pivoteren. Je kan dan in dezelfde kolom zoeken of er lager nog een niet-0 cel is. Zo ja, dan kan je verder door de twee rijen om te wisselen. Zo niet, dan heeft je stelsel geen oplossing... en daar zul je mogelijk ook rekening mee moeten houden...

edit:
op afhankelijkheid controleer je al. Maar dat geeft dus nog geen garantie dat je geen 0 tegenkomt, alleen dat je tijdens proces altijd een andere rij kan vinden om mee op te wisselen wanneer je er een tegenkomt!

[ Voor 12% gewijzigd door Infinitive op 17-08-2004 15:22 ]

[Verwijderd]

Ik heb in de PerformGaussian functie ook enkele aanpassingen gedaan om te testen.

C:

void PerformGaussian (float (*SamplesMatrix)[6])
{
    float Temp;
    
    // Walk through columns
    for (INT8 i = 0; i < 5; i++)
    {
        // Walk through rows beneath working row number
        for (INT8 j = i + 1; j < 5; j++)
        {
            // Walk through current column and those on the right
            for (INT8 k = 5; k >= i; k--)
            {
                SamplesMatrix[j][k] -= SamplesMatrix[i][k] *
                    divsp(SamplesMatrix[j][i], SamplesMatrix[i][i]);
            }
        }
    }

    // Back substitution
    SamplesMatrix[4][5] = divsp((SamplesMatrix[4][5]), SamplesMatrix[4][4]);

    Temp = SamplesMatrix[3][4] * SamplesMatrix[4][5];
    SamplesMatrix[3][5] = divsp((SamplesMatrix[3][5] - Temp), SamplesMatrix[3][3]);

    Temp = SamplesMatrix[2][3] * SamplesMatrix[3][5]
             + SamplesMatrix[2][4] * SamplesMatrix[4][5];
    SamplesMatrix[2][5] = divsp((SamplesMatrix[2][5] - Temp), SamplesMatrix[2][2]);

    Temp = SamplesMatrix[1][2] * SamplesMatrix[2][5]
             + SamplesMatrix[1][3] * SamplesMatrix[3][5]
             + SamplesMatrix[1][4] * SamplesMatrix[4][5];
    SamplesMatrix[1][5] = divsp((SamplesMatrix[1][5] - Temp), SamplesMatrix[1][1]);

    Temp = SamplesMatrix[0][1] * SamplesMatrix[1][5]
             + SamplesMatrix[0][2] * SamplesMatrix[2][5]
             + SamplesMatrix[0][3] * SamplesMatrix[3][5]
             + SamplesMatrix[0][4] * SamplesMatrix[4][5];
    SamplesMatrix[0][5] = divsp((SamplesMatrix[0][5] - Temp), SamplesMatrix[0][0]);
}

Ik wist al dat er een behoorlijke bottleneck zat bij de delingen die gedaan moeten worden, zowel bij de Gaussian als de Gauss-Jordan implementatie. Ik heb nu alle delingen vervangen door de divsp functie uit de fastRTS library van TI. Dit is een library met implementaties voor rekenkundige operaties als div, sqrt, sin etc. speciaal geoptimaliseerd voor de 67xx series DSP. Dit gaat overigens enigsinds ten koste van de nauwkeurigheid, echter het had wel een major impact op de performence, vele malen groter dan ik verwacht had. Het aantal klok cycles voor de Gaussian methode was hierna nog slechts 5907.

Daarnaast heb ik het stukje van de backsubstitution unrolled dit leverde nog een extra winst op waardoor de Gaussian methode nu nog slechts 5312 cycles nodig heeft.

Ik heb de delingen bij de Gauss-Jordan methode ook vervangen door de divsp functie. De Gauss-Jordan methode heeft nu 6059 cycles nodig om de matrix door te rekenen. Voor de volledigheid ook nog even de code:

C:

void PerformGaussJordan (float (*SamplesMatrix)[6])  
{  
    float Calculator;  
      
    // Walk through the rows of the matrix  
    for (INT8 i = 0; i < 5; i++)  
    {  
        // Store working point value  
        Calculator = SamplesMatrix[i][i];  
          
        // Divide sample value by working point value.  
        // Results in a 1 in the working point.   
        SamplesMatrix[i][0] = divsp(SamplesMatrix[i][0], Calculator);  
        SamplesMatrix[i][1] = divsp(SamplesMatrix[i][1], Calculator);  
        SamplesMatrix[i][2] = divsp(SamplesMatrix[i][2], Calculator);  
        SamplesMatrix[i][3] = divsp(SamplesMatrix[i][3], Calculator);
        SamplesMatrix[i][4] = divsp(SamplesMatrix[i][4], Calculator);
        SamplesMatrix[i][5] = divsp(SamplesMatrix[i][5], Calculator);  
          
        // Walk through the rows of the matrix...  
        for (INT8 j = 0; j < 5; j++)  
        {  
            // ...except the current row  
            if (j != i)  
            {  
                // Store value in working row working point column  
                Calculator = SamplesMatrix[j][i];  

                // Walk through columns in working row   
                SamplesMatrix[j][0] -= Calculator * SamplesMatrix[i][0];  
                SamplesMatrix[j][1] -= Calculator * SamplesMatrix[i][1];  
                SamplesMatrix[j][2] -= Calculator * SamplesMatrix[i][2];  
                SamplesMatrix[j][3] -= Calculator * SamplesMatrix[i][3];  
                SamplesMatrix[j][4] -= Calculator * SamplesMatrix[i][4];  
                SamplesMatrix[j][5] -= Calculator * SamplesMatrix[i][5];  
            }  
        }  
    }  
}

[Soultaker]

Wat je trouwens ook kunt doen om inefficiënte delingen te vermijden, is eerst 1/pivot uitrekenen en daar vervolgens mee vermenigvuldigen, in plaats van steeds door de pivot te delen. Daarmee vervang je je delingen door vermenigvuldiging en voeg je een extra deling toe; netto zijn het meer operaties maar als je deling (veel) kostbaarder is dan je vermenigvuldiging (wat meestal het geval is) dan kan het waarschijnlijk uit.

[Apollo_Futurae]

Soultaker schreef op 16 augustus 2004 @ 15:52:
Ik vind het wel logisch dat de implementatie van de Gauss-Jordan-eliminatie sneller is, aangezien die een complexiteit van O(N²) heeft, terwijl de implementatie van de Gauss-eliminatie een complexiteit van O(N³) heeft.

Dat is niet waar. Gauss met substitutie en Gauss-Jordan zijn beide O(N³).
In flops gerekend is Gauss met substitutie ~ N³/3; Gauss-Jordan is 50% langzamer (N³/2).
_{[bron: Linear Algebra 3rd ed., Fraleigh and Beauregard, Addison-Wesley 1995]}

Het is niet erg ontopic, omdat de matrices zo klein zijn, maar ik wilde het toch even melden .

[Soultaker]

Apollo_Futurae schreef op 17 augustus 2004 @ 16:18:
Dat is niet waar. Gauss met substitutie en Gauss-Jordan zijn beide O(N³).
In flops gerekend is Gauss met substitutie ~ N³/3; Gauss-Jordan is 50% langzamer (N³/2).
_{[bron: Linear Algebra 3rd ed., Fraleigh and Beauregard, Addison-Wesley 1995]}

Ah, je hebt gelijk. Ik had het for-lusje op regel 29-31 over het hoofd gezien, vandaar dat ik op O(N²) uitkwam.

[Verwijderd]

Macros schreef op 16 augustus 2004 @ 15:01:
Ook is het sneller om pointers te gebruiken ipv. array indexering. Elke keer als ie een array indexeerd moet hij het begin van de (array addres) + (size van een element) * (index) doen. Met een pointer index hoef je alleen maar een (pointer) + (size van element). Dat scheelt 1 operatie. In hoeverre dat echt helpt moet je testen, grote kans dat de compiler zelf ook al en optimalisatie toepast.

Ik heb het even getest, maar ik merk geen enkel verschil tussen pointer notatie en array indexing. Misschien dat de compiler hier al iets mee doet?

[Soultaker]

Het hangt nogal van de processor en de compiler af. Ten eerste kan een goede compiler recht-toe-recht-aan for-lusjes zoals die hier gebruikt worden wel omschrijven naar een pointer waar steeds wat bij opgeteld wordt. Ten tweede kan in de x86-architectuur bijna elke operand worden geschreven in een vorm *(ptr+b*c) (met b een beperkt aantal 2-machten, maar toch); de berekening van het effectieve adres kos daarmee dus geen aparte instructies en levert vaak ook nauwelijks verschil in performance op.

Je werkt echter met een ander soort chip en specialistische compiler. Hoe het daarmee werkt is moeilijk te zeggen. Misschien is het sowieso verstandig om deze functie gewoon in assembly te schrijven, als de snelheid ervan zo kritiek is? Ik neem aan dat je een handleiding bij je processor(s) hebt waar ook implementatietips in staan.

[ Voor 3% gewijzigd door Soultaker op 17-08-2004 17:27 ]

[Macros]

Schrijf een genetisch algoritme dat je code optimaliseerd voor je gebruikte dsp

[Soultaker]

offtopic:
Dat is voor watjes die zelf niet kunnen coden.

[Zoijar]

offtopic:
burs

[ Voor 4% gewijzigd door Zoijar op 17-08-2004 20:53 ]

[Verwijderd]

Soultaker schreef op 17 augustus 2004 @ 17:26:
Misschien is het sowieso verstandig om deze functie gewoon in assembly te schrijven, als de snelheid ervan zo kritiek is? Ik neem aan dat je een handleiding bij je processor(s) hebt waar ook implementatietips in staan.

De DSP heeft een VLIW-architectuur met 2 datapaden en 4 rekenunits per pad. Om hier assembler voor te programmeren is een erg grote uitdaging, aangezien je alle 8 de rekenunits bezig moet zien te houden, maar ondertussen ook deadlock moet zien te voorkomen. Vandaar dat ik graag eerst een zo optimaal mogelijke implementatie in C code wil hebben. De assembler code die de compiler hieruit genereerd kan ik dan gebruiken als referentie.

Ik heb de tips van Zoijar en Soultaker, om 1/pivot te gebruiken, doorgevoerd in de PerformGaussJordan functie. Dit leverde nog wat extra winst op, waardoor deze functie nog maar 5191 cycles nodig heeft voor het oplossen van de matrix.

[ Voor 3% gewijzigd door Verwijderd op 18-08-2004 09:15 ]