Invers av triangulär matris, C algoritm

[Korken]

Godagens!

Är det någon här som sitter på kod för inversen av en triangulär matris?
Det skulle vara mycket skönt om man slapp skriva en och koncentrera sig på resten.
Man brukar nog ofta använda http://www.cs.berkeley.edu/~knight/knig ... poster.pdf

[Korken]

Det är värt att lägga till att jag har testat LINPACK med BLAS. Men när jag kompierar det så får jag dessa errors:
Edit: Detta kompileras i Ubuntu med arm-none-eabi-gcc från CodeBench för en STM32F405.

source/blas1_s.o: In function `snrm2':
/home/korken/Programming/KFly_STM32F4/source/blas1_s.c undefined reference to `sqrt'
source/blas1_s.o: In function `srotg':
/home/korken/Programming/KFly_STM32F4/source/blas1_s.c undefined reference to `sqrt'
source/linpack_s.o: In function `schdc':
/home/korken/Programming/KFly_STM32F4/source/linpack_s.c:249: undefined reference to `sqrt'
source/linpack_s.o: In function `schdd':
/home/korken/Programming/KFly_STM32F4/source/linpack_s.c:421: undefined reference to `sqrt'
/home/korken/Programming/KFly_STM32F4/source/linpack_s.c:431: undefined reference to `sqrt'
source/linpack_s.o:/home/korken/Programming/KFly_STM32F4/source/linpack_s.c:431: more undefined references to `sqrt' follow
/usr/local/CodeBench_1802_EABI/bin/../lib/gcc/arm-none-eabi/4.6.1/../../../../arm-none-eabi/lib/thumb2/libc.a(lib_a-exit.o): In function `exit':
exit.c:(.text+0x1c): undefined reference to `_exit'
/usr/local/CodeBench_1802_EABI/bin/../lib/gcc/arm-none-eabi/4.6.1/../../../../arm-none-eabi/lib/thumb2/libc.a(lib_a-writer.o): In function `_write_r':
writer.c:(.text+0x16): undefined reference to `_write'
/usr/local/CodeBench_1802_EABI/bin/../lib/gcc/arm-none-eabi/4.6.1/../../../../arm-none-eabi/lib/thumb2/libc.a(lib_a-closer.o): In function `_close_r':
closer.c:(.text+0x12): undefined reference to `_close'
/usr/local/CodeBench_1802_EABI/bin/../lib/gcc/arm-none-eabi/4.6.1/../../../../arm-none-eabi/lib/thumb2/libc.a(lib_a-fstatr.o): In function `_fstat_r':
fstatr.c:(.text+0x14): undefined reference to `_fstat'
/usr/local/CodeBench_1802_EABI/bin/../lib/gcc/arm-none-eabi/4.6.1/../../../../arm-none-eabi/lib/thumb2/libc.a(lib_a-isattyr.o): In function `_isatty_r':
isattyr.c:(.text+0x12): undefined reference to `_isatty'
/usr/local/CodeBench_1802_EABI/bin/../lib/gcc/arm-none-eabi/4.6.1/../../../../arm-none-eabi/lib/thumb2/libc.a(lib_a-lseekr.o): In function `_lseek_r':
lseekr.c:(.text+0x16): undefined reference to `_lseek'
/usr/local/CodeBench_1802_EABI/bin/../lib/gcc/arm-none-eabi/4.6.1/../../../../arm-none-eabi/lib/thumb2/libc.a(lib_a-readr.o): In function `_read_r':
readr.c:(.text+0x16): undefined reference to `_read'

De skulle annars vara perfekt, men lyckas inte bli av med dessa errors.

[ylle]

/home/korken/Programming/KFly_STM32F4/source/blas1_s.c undefined reference to `sqrt'

Är det så enkelt som?
#include <math.h>

[Norberg]

Gissar på att du missat länka mot libm (math library), läknas enklast med -lm flaggan till gcc

[Korken]

Båda är länkade, så det är okej där, ska se vad som inte vill sig med länkaren.

Jag tog och knegade ihop en egen funktion så det fixade sig.
Inversen av en 6x6 triangel-matris tar 35 multiplikationer och 15 divisioner.
Jag ska optimera lite mer, sedan lägger jag upp koden här för de som är sugna.

Edit: fel matris storlek va 13x13

[Korken]

Här står det hur man fixar kompieringsfelen: http://balau82.wordpress.com/2010/12/16 ... -programs/
srqt fugerade inte då det är double. jag måste föra bre float versionen.

Mer tester av algoritmen:
13x13 matriser tar 364 multiplikationer och 78 divisioner.

[Korken]

Okej, nu har jag räknat lite på min funktion.
Den tar 2.78 * n^2.975 antal klockcykler att exekvera, där n är matrisens dimension.

Är detta bra eller dåligt? Jag har inga referenser alls tyvärr.

Edit: Detta är för en generell matris!

[Korken]

Jag fick en ny idé, vad tror ni om detta?
Om man kollar på bilden i mitten från länken jag kollade på så kan man hela tiden dela sin matris till mindre och mindre trianglar tills man får en 1x1 matris.
Den kan man sedan invertera och börja vandra uppåt. Detta ger ett rekursivt programanrop. Kan man ha det i en MCU med tanke på hur mycket minne som kommer gå åt?
Om inte, hur kan man lättast översätta en rekursiv funktion som denna till en itterativ?

Detta get att man bara behöver dividera lika många gånger som matrisens dimension. Tex en 13x13 matris behöver 13 divisioner. Resten är bara Muls och Summor.

[Korken]

Okej, jag skrev ihop den rekursiva matrisinversen och såhär vart resultatet:
För små matriser < 50x50 (vid matriser som är ca 500x500 så är de ish lika snabba) så använder den rekursiva ca 30% färre klockcykler, men den är rekursiv.
Hur kan man avgöra om det är bra att använda den eller om man ska använda något annat?

[LHelge]

Rekursiva funktioner äter stackminne till frukost, för en stm32f4 är det nog inga problem men det kan nog vara bra att lägga en begränsning i koden och kommentera varför om det är något du tänker släppa fritt.

[Korken]

Jo jag täkter släppa det idag
Jag har en idé kvar och det är att göra den rekursiva funktionen iterativ. Dock så är det mycket svårare att dela upp triangel matrisen i mindre matriser när man kör iterativt.
Man kan då råka förstöra symmetrin och då fallerar algoritmen.

Finns det något sätt att veta hur mycket stack varje anrop tar? Då kan jag förutspå det och skriva in i dokumentationen.

[Andax]

Hur mycket stack som används är i princip returadress + backup av extra register som används av funktionen, dvs det beror på!

[Korken]

Sant, jag ska ta och kompilera och se vad disassemblyn säger.
Undertiden så behöver jag hjälp med en funktion.
När jag multiplicerar U_invers med L_invers så finns det massa 0or man kan hoppa över i matriserna.
Just nu gör jag det såhär (testar i Java först):

Kod: Markera allt

public static float[][] MultiplyULMatrix(float a[][], float b[][]) { 
	int n = a.length;
	float[][] resultant = new float[n][n];
	/* *
	 * 
	 * Illustration why we don't have to multiply and add every element:
	 * x is some number in the matrix, but the x:es are not necessarily the same.
	 * 
	 * | x x x x x |   | x 0 0 0 0 |
	 * | 0 x x x x |   | x x 0 0 0 |
	 * | 0 0 x x x | * | x x x 0 0 | = ... 
	 * | 0 0 0 x x |   | x x x x 0 |
	 * | 0 0 0 0 x |   | x x x x x |
	 * 
	 * Many zeroes that don't have to be taken into account when multiplying!
	 * 
	 * */
		
	for(int i = 0; i < n; i++) {
		for(int j = 0; j < n; j++) {
			for(int k = j; k < n; k++) { /* Now it will not multiply with half of the 0s */
				resultant[i][j] += a[i][k] * b[k][j];
			}
		}
	}
	   
	return resultant;
}

Problemet är att jag bara lyckas bli av med 0:orna i ena av de två matriserna. sätter jag k = i så struntar den i nollorna i den första matrisen och sätter jag k = j så struntar den i 0:orna i den andra matrisen.
Men jag lyckas inte få den att strunta i nollorna i båda. Kan någon vänlig själ hjälpa mig här? Jag har suttit fast så länge på detta nu.

[Korken]

Man kan lägga till detta innan den sista loopen:

Kod: Markera allt

if (j < i)
	k = i;
else
	k = j;

Men att använda if-satser är inte optimalt. Finns det någon bättre lösning?

Edit:
Jag undersökte också rekursionsdjupet. Det går lika djupt som matrisens storlek. Dvs är matrisen 50x50 så gåt rekursionen 50 anrop som djupast.
Detta borde betyda att minnesåtgången för rekursionen ökar linjärt med matrisens storlek.

[mri]

Printa ut (eller på nåt annat sätt ta reda på) en variabels adress som finns på stacken. Låt den rekursiva funktionen anropa sig själv och kolla hur mycket variabelns adress ändrade. Där har du stackåtgången per nivå.