Invers av triangulär matris, C algoritm

[mri]

Det går att ändra en rekursiv funktion till iterativ genom att bygga en egen "stack" som du pushar och popar data från, men den stacken kommer ju likväl att behöva minne, men det kunde allokeras och friges efter behov.

[Korken]

Tackar! Det ska jag göra!
Ska testa detta lite mer på en MCU också. Så jag vet mer exakt vad det blir för exekveringstid och stack djup som du säger.
Nu har jag räknar på antalet klockcykler som beräkningar tar + att ladda och spara variabler. Real world data är alltid bättre.

[thebolt]

En viktig fråga är om din matris har enhetsdiagonal (dvs alla diagonal-element är 1) eller är de olika?

Så vitt jag vet resulterar den mest effektiva metoden att beräkna inversen på en generell triangulär NxN matris i N divisioner och totalt ung. N^3/3 + 2N^2 operationer. Formuleringen och härledningen i sig är rekursiv men det är inga problem att implementera iterativt, med enda data som behövs är in och ut-matriserna (kan tom vara så att man kan göra det "in-place", men det måste undersökas i så fall).

I korthet:
u - in, v - ut
k = n..1
i = k-1..1

v_kk = u_kk^-1 = 1/u_kk
v_ik = -v_ii * sum(u_ij * v_jk, j = i+1 to k)

Om du gör den beräkningen på papper kommer du se att du alltid har datan du behöver redan beräknad.

[Korken]

Tackar! De är så min första algoritm fungerar jag tycker dock inte om att den kräver så många divisioner så jag försökte fixa det.
Slänger upp koden när jag kommit till jobbet så sett ni hur jag har gjort.

Edit: såg fel! Jag ska undersöka den algoritmen också!

[thebolt]

Jag räknade lite snabbt och insåg att det går använda den algoritmen jag visade "in-place" också, dvs in och utmatrisen i samma data-utrymme.

Om du är intresserad kan jag skriva ihop en enkel C-implementation av det (har en nu, men den beror på lite andra grejjer) när arbetsdagen är slut.

[Korken]

Det får du gärna göra!

Jag kan lägga till den i mitt bibliotek och jämföra med mina andra två implementationer, så har vi nog en ganska bra forumsbaserad matris-beräknings-grej.

[Korken]

Nu när jag kollar på min algoritm, är det inte så den här fungerar fast skriven på olika vis?
Man kan nog bli av med den andra divisionen nu när jag tänker efter som du säger.

Kod: Markera allt

public static void InvertUMatrix(float u[][], int n) {
	float sum;
	int i,j,k;
		
	for (j = n - 1; j > -1; j--) {
		u[j][j] = 1.0f/u[j][j];

		for (i = j - 1; i > -1; i--) {
			sum = 0.0f;
				
			for (k = j; k > i; k--)
				sum -= u[i][k]*u[k][j];
				
			u[i][j] = sum/u[i][i];
		}
	}
}

[Korken]

Jo, de va samma! De jag märkt är att min prestanda bov är min LU-faktorisering.
Går det att skriva denna på något bättre vis? Den drar nu dubbelt så mycket som en invers och har många divisioner.
Just nu fungerar den så att den praktiskt taget gör gausseliminering.

Kod: Markera allt

public static void LUFactorize(float[][] a, float[][] L, float[][] U) {
	int k, i, j, n;
	n = a.length;
		
	for(k = 0; k < n; k++)
	{
		L[k][k] = 1;

		for(i = k + 1; i < n; i++)
		{
			L[i][k]=a[i][k]/a[k][k];

			for(j = k + 1; j < n; j++) 
				a[i][j] = a[i][j]-L[i][k]*a[k][j];

		}
		for(j = k; j < n; j++)
			U[k][j] = a[k][j];
	}
}

[Korken]

Strunta i denna post! Jag hade tankevurpa!

Utökar med en extra fråga.
Hur skickar man bäst 2D array pekare till funktioner i C?

Jag gjorde såhär nu, men det fugnerar inte:
Borde jag bygga om det till en float *L och beräkna offseten manuellt i funktionen?

Kod: Markera allt

void InvertLMatrix(float **L, int n)
{
	float sum;
	int i,j,k;

	for (j = 0; j < n; j++)
	{
		for (i = j + 1; i < n; i++)
		{
			sum = 0.0f;

			for (k = j; k < i; k++)
				sum -= L[i][k]*L[k][j];

			L[i][j] = sum;
		}
	}
}

[Korken]

Nu så har jag testat det lite också!
Inversen av en generell 6x6 matris tar 5 648 klockcykler. Inte så bra som jag hade hoppats, med det är okej.
Inversen av en 12x12 matris tar 30 055 klockcykler.
Inversen av en 50x50 matris tar 1 606 259 klockcykler.
Testat på en STM32F405 i 168MHz med FPUn igång.
Där är några referenser.

Koden: https://github.com/korken89/KFly_STM32F ... _algebra.c

Kom gärna med förslag till förbättringar! Jag tror att det kan bli mycket snabbare!

[Andax]

Bläddrade lite i en gammal version av "Numerical recepies in C++" (second edition) och såg att kunde göra en snabbare än N^3 och med framförallt få divisioner med Strassens metod (kapitel 2.11 i min version av 'NR in C++'). Kanske kan vara en väg framåt?
Har du mycket programminne kan det ju vara en variant att hårdkoda en funktion för varje storlek av matris som du vill köra. Då kan man ju göra en hel del loop-unroll etc för att kvicka upp...

Förresten, hur många klockcykler tar respektive multiplikation, addition och division?

[Korken]

Tyvärr har jag inte den boken, men ska se om jag inte kan få tag på den.
Jag har också funderat på att göra flera för olika storlekar. Ska testa det imorgon för att se hur stor skillnad det gör.

När det kommer till klockcykler så ser det ut så här:
Add/Sub/Mul (float och int): 1 klockcykel
Div (float): 14 klockcykler
Div (int): 7 klockcykler

[Andax]

Se http://www.nr.com/

Det verkar som de har så att man kan läsa 30 sidor per månad fritt ur den senaste versionen.

Reflekterade inte över antal klockcykler du nämnde tidigare för 6x6, 12x12 och 50x50. Generellt låter det som många cykler för något som är O(N^3)

EDIT: Det verkar som NR in C från 1992 är helt fri att titta i: http://apps.nrbook.com/c/index.html

[swp]

Den äldre versionen (second edition) finns i pdf-format:

http://www.nrbook.com/a/bookcpdf.php

[Korken]

Plottar man mina resultat så ser de ut såhär:
Antal klockcykler som funktion av marisdimension: 12*N^3+40*N^2+340*N har jag för mig de va.