hawkan skrev:En grej som kan testas, om det ändå ska testas, är att byta alla "float" mot "double". Substitute all eller vad det kan heta.
Då vet du skillnad både i resultat och beräkningstid. Bara en idé, men jag hade gjort det. Jag hade för övrigt aldrig
gjort linear algebra med float, vad en sån åsikt nu kan vara värd.
Den verkar inte göra någon skillnad för q-matrisen.
Det blir ju lite krångligt att förstå det där, med alla *-tecken överallt.
Du kan ju lägga in alla uttryck som liknar *((ptr_r + i*n)+k) i temporära variabler, och köra ovanstående beräkning på de temporära variablerna. Då kan du titta på och titta på varje sådant värde, och kontrollera att det stämmer, eller ifall någon av de där i,j,k,q,r, o.s.v är inknappat på rätt ställe.
Jag använder nästan aldrig variabler som är endast en bokstav, eftersom att det är så svårt att förstå koden och hitta buggar. Men matematisker gillar ju det.
Exakt! Som du ser så har jag bara skrivit av MATLAB-kod. Men jag tror att jag kanske måste hitta annan QR-kod för C.
Jag har ju använt mig av modifierad Gram-Schmidt, men kanske Householder är bättre. Den är mer avancerad.
Det blir kanske problem med dynamisk allokering om du en dag ska använda 100x100-matriser istället för 6x6, och inte tänkt på att dom inte får plats i minnet. Men jag tror inte heller att det är minnesallokeringen som kommer vara problemet.
bearing skrev:Det blir kanske problem med dynamisk allokering om du en dag ska använda 100x100-matriser istället för 6x6, och inte tänkt på att dom inte får plats i minnet. Men jag tror inte heller att det är minnesallokeringen som kommer vara problemet.
Jag kommer alltid att använda 36x36 i mitt nuvarande program. Det är alltid samma matriser jag ska beräkna hela tiden.
Det jag gör för att verifiera om det fungerar eller inte är:
1. Beräkna i MATLAB
2. Skriv min C kod
3. Verifiera med valgrind i PC
4. Flytta över C-koden till STM32
Both algorithms compute a QR factorisation of a matrix A, that is, A=QR. Due to rounding errors, we have A=Q~R~+E, where Q~ is no longer orthogonal (but close to in some sense). The residual E is not usually the problem, for both Householder and Gram-Schmidt (even the classic one), ∥E∥≤cϵ∥A∥, where c depends only on the dimension of A, ϵ is the machine precision. Here, ∥⋅∥ is some "rounding errors friendly" norm, e.g., any p-norm (usually 1, 2, or ∞).
The main difference is in the accuracy of the orthogonal factor. The Householder orthogonalisation gives Q~ which is almost orthogonal in the sense that ∥I−Q~TQ~∥≤cϵ. In the modified Gram-Schmidt, this "loss of orthogonality" is proportional to the condition number of A, that is, MGS yields ∥I−Q~TQ~∥≤cϵκ2(A) (note that the Q-factor from the classical Gram-Schmidt satisfies a similar bound with the condition number of A squared).
The difference between Householder and MGS is, roughly speaking, due to the fact that Householder computes the Q-factor as a product of accurate Householder reflections while MGS directly orthogonalises the columns of A.
By the way, both Householder and MGS compute an accurate R-factor R~ in the sense of backward error. For both, there is an (exact) orthogonal matrix Q^ and a residual matrix E^ (of course different for each) such that ∥E^∥≤cϵ∥A∥ and A=Q^R~+E^.
Den verkar inte göra någon skillnad för q-matrisen.
Det betyder ju att produkten till höger är 0. Och det kan vara ett tecken på en bugg, att den läser från fel ställe i minnet, där värdena är noll typ. Så som jag sa, sätt varje faktor i temporära variabler och kolla att de får rätt värden, och inte alltid blir noll.
Eller så är det helt enkelt så att alltid en av de två faktorerna i produkten är noll. Och då ger ju ovanstående inte någon skillnad för Q. Men det antar jag inte är fallet.
R ändras ju på andra rader, så det är väl inte så konstigt.
Om R ändras på den där raden, ja då är något väldigt konstigt. Men det tror jag inte.
Hur som helst verkar det ju som att du inte är så intresserad av att söka efter ev. buggar, utan ska nu börja på ett nytt kapitel, trots många återstående frågetecken. Vore det inte klokt att ta reda på vad som händer i din kod, och varför den inte gör som den ska, innan du byter? Den lärdomen kommer nog att komma till användning senare när den nya algoritmen inte fungerar.
Efter jag fixar till qr.c filen så funderar jag på att skriva om projektet så att man bara allokerar minnet en gång och därmed får man återanvända matriserna.
Det var ingen bugg i qr.c, utan det var helt enkelt C's egna sätt att få fram decimaler. Jag har jämfört med Octave och det stämmer.
Det jag gjorde är att om jag har matrisen H så utfär jag bara denna ekvation: \(H = QR\)
Där efter så multiplicerar jag ihop dem och då fick jag tillbaka mitt H. Bra! \(Q*R\)
Jag testade även med SVD. \(H = USV^{T}\)
Sedan satt jag ihop \(USV^{T}\) och fick mitt H, men det var rätt mycket fel i det H:et, vilket är normalt då det är SVD.
Så slutsatsen kan jag säga att jag hittar inga fel. Lite sämre avrundningar i C och C verkar vara dåligare på flyttal än vad GNU Octave är. Men annars så fungerar det.
Jag säger så här. Mitt C-bibliotek är fullt godkänt för matrisberäkningar, men nackdelen med biblioteket är att den har svårt att hantera små små små tal typ e-7 och nedåt. Vad det kan bero på tror jag det har med själva C programmeringsspråk att göra. Om C++ skulle göra bättre ifrån sig, vet jag inte. Jag tvivlar lite på det faktiskt. Rätta mig om jag har fel.
Annars funderar jag på att skriva om strukturen på biblioteket så att man bara allokerar matriser EN gång. Där efter återanvänder man dem. Vad tror ni om det? Tror ni det skulle vara bra om jag anpassar det för c89 standard? Eller tjänar jag något på det?
Nu blev jag sugen på lite debatt här. Motivera VARFÖR jag inte ska använda C++ i STM32? Alla motiveringar accepteras.
C ska inte ha några problem att hantera flyttal inom de specificerade ramarna för
den aktuella flyttals typen (float, double, whatever). e-7 är ju inte speciellt
litet för i princip alla flyttals typer.
Jag tror mer på att det är något temporärt värde inne i formlerna som sticker iväg
uppåt eller nedåt.
Så det borde inre ha någonting med C som språk i sig. Kanske med någon implementering
av flyttalsoperationer i ett libc eller liknande, men om rutinerna använder en FPU så borde
de bara skicka samma FPU instruktioner oavsett programmeringsspråk.
Så kan det vara. Men med tanke på att jag delar upp H = QR och sedan Q*R -> H och jag gör samma sak med SVD där jag sätter H = USV^T och sedan sätter jag U*S*V^T -> H, med lite avrundingar.
Jag kör samma algoritmer i Octave och dem visar samma. Skillnaden mellan i C och Octave är att precisionen är bättre i Octave än i C. Jag mätte upp just nu. Om precisionen är e-10 i C för SVD så är precisionen e-18 i Octave för just samma tal. Vi alla vet att e-10 är riktigt litet, men det kan göra påverkningar ändå.
Vad tror du Sojjan? Kan det vara så att om jag initialiserar matriserna bara en gång och sedan återanvänder jag dem. Då kan det bli bättre resultat?
