Hitta flaskhals i Matlab?

Jag sitter och skriver på ett exjobb som går ut på att lösa lite differentialekvationer numeriskt i Matlab. I korthet löser jag ett linjärt ekvationssystem i ett antal celler som jag itererar över. Koden är inte optimerad på prestanda än så länge och kommer eventuellt inte bli på grund av tidsbrist men när jag satt och tittade lite på vad som händer när den räknar så blev jag förvånad. 100 celler går på under en sekund. 1000 celler tar ca 30 sekunder och 10000 celler orkade jag inte vänta på så gav upp efter tio minuter. När jag kollar på processorbelastningen så tar Matlab ca 20 %. Koden är enkeltrådad och jag kan rätt tydligt se vilken kärna det är som får jobba med den. RAM är inte heller något problem för det går bara åt en ca 700 MB för tillfället. Det kan bli ett problem senare men nu verkar Matlab kunna hålla alla variabler i RAM.

Så vad är det som tar så extremt lång tid? Vad väntar datorn på? Är det att läsa data från RAM till cache som tar lång tid eftersom själva beräkningen inte är tung? För det verkar ju inte som om det är RAM som krånglar då jag har 16 GB tillgängligt och Matlab borde använda allt den behöver.

Vore även trevligt om det fanns ett generellt verktyg för har en liknande fundering kring ett spel där spelet säger att prestandan begränsas av min processor men enhetshanterar säger att jag peakar på 80 % och oftast ligger jag runt 50 % så bevisligen finns det massa prestanda kvar att hämta.

Det som tar mycket tid verkar vara lösningen av det linjära ekvationssystemet som för tillfället bara har tre variabler och det kommer inte bli så mycket större. Sen är det Gauss-Seidel som görs över alla cellerna.
Jag vet redan att den bara använder en kärna för algoritmen går inte att köra seriellt på ett smidigt sätt. Varje iteration beror på lösningen av föregående cell och så vidare. Som sagt så kan jag se vilken kärna som används men den går inte över 20 %.

Som sagt Gauss-Seidel och tror inte jag kan byta hur som helst för det är ett CFD-problem som kräver en iterativ lösning. Eller tekniskt sett hade det gått att lösa som ett linjärt ekvationssystem men det skulle haverera för stora matriser senare.

Kommersiell CFD-mjukvara kör samma algoritm och klarar av att belasta processorn väldigt mycket mer. Men de kanske kör lite annan voodoo för att göra det mer effektivt.

Att den inte löste det för extremt många celler kan bero på numeriska problem som ställer till det. Jag analyserade inte lösningen ordentligt men det är möjligt att det krånglar med toleranser för konvergens i något steg. Eller att matlab bråkar i vanlig ordning. Men ja tiden ökar rätt mycket vid större körningar.

Nej inte riktigt för att den ensamma tråden bara belastar en kärna 20 %. Då är det ju inte uppenbart att fler trådar löser det.

Egentligen är det inte tidsåtgången som oroar mig utan det faktum att den inte belastar processorn alls.

Problemet är en massbalans i en dimension med konvektion, dispersion, reaktion och än så länge enkla källtermer. Programmet är helt egenskrivet och räknar fram en koncentrationsprofil baserat på massbalansen. I varje cell har jag olika storlekar på de partiklar som reagerar så i varje cell ska det inte bara lösas en massbalans med transport in och ut ur cellen för varje partikel storlek utan jag ska även lösa en balans för hur storleken förändras. Just förändringen av storlek löses med ett linjärt ekvationssystem som just nu har tre obekanta och kommer aldrig ha fler än 15. Den biten går fort. Det som tar tid är att lösa massbalansen för alla celler.

Metoden går ju ut på att jag gissar en koncentration för alla celler och sedan använder jag den gissningen och min massbalans för att räkna fram en bättre koncentration och sedan littererar man som en galning tills man får en lösning.

Den är helt egenbyggd eftersom det är det exjobbet går ut på.

Ja det är så jag måste göra för att lösa problemet. Det är ju en iterativ metod så lösningen av det 101a talet beror av värdet på det 100e. Tidigare uträknade tal sparas redan i en vektor eftersom algoritmen går igenom alla tal från början till slut och räknar fram ett nytt tal baserat på talen bredvid. Så ett nytt värde för talet i fås av i-1 och i+1. Så löser man från i upp till antalet celler och sedan börjar man om från ett och gör ett svep till fram tills dess man har hittat en konvergerande lösning. Så jag sparar de senaste lösningarna i en vektor och inget annat.

Jag använder rätt få av Matlabs funktioner för tillfället. Har kört profileraren för att få fram tiderna men den gav inte så mycket information. Ska ta och köra en omgång till i morgon och se om jag kan lura ut lite mer från den.

Ja det är ju det som gör mig så fundersam. Det är inte datorn i signaturen som används utan det är en fyrkärnig i7 med HT så har 8 trådar att leka med. Har för mig att det är en Nahelem men kommer inte ihåg exakt modell. Enligt enhetshanteraren så tar Matlab ca 20 % CPU och tittar jag på varje enskild logisk kärna så ser jag tydligt vilken det är som faktiskt arbetar och den ligger väldigt lågt. Så något annat sätter stopp.

Det var cacheminnet som var min hypotes. Frågan är om det finns ett sätt att kolla om så är fallet för då är det nog inte mycket jag kan göra åt det.

Jag använder egenskrivna saker för att det är en del av exjobbet. Ingen utmaning att bara använda färdiga saker utan tanken är att jag ska lära mig lite på vägen också.

Egentligen är det inte viktigt att optimera för det har jag inte tid med i exjobbet. Det är en hård tidsplan som det är om jag ska göra allt som min examinator vill ha gjort. Hårdare optimering får han göra själv sen om jag lyckas få till arbetet tillräckligt bra för att det faktiskt ska användas senare. Tar körningen för lång tid får jag helt enkelt använda en grövre mesh. Hade varit grymt kul att pilla med optimering också för det har jag aldrig gjort innan men tyvärr har jag inte riktigt möjligheten för tillfället.

Nej inga matriser som ökar storlek. Allt är fördefinierat innan jag börjar räkna med det.

Tyvärr hinner jag inte lära mig C eller för den delen FORTRAN innan arbetet ska vara klart. Dessutom är det en modul som ska in i en större matlabmodel i slutändan så tror inte det är så populärt att porta den till ett annat språk.

Edit: Enligt profiler så tar raden:

global const mesh

hela 20 % av den totala körtiden. const och mesh är båda structures med konstanter respektive detaljer om min mesh som båda behövs för att köra funktionsfilen. Det är alltså en underfunktion som tar större delen av tiden. Jag gjorde en provkörning med 1000 celler och funktionen anropas då 230 000 gånger. Jag trodde inte att det skulle ta så lång tid att sätta upp globala variabler. Tyvärr känns det inte som om jag har något val då det är för mycket data för att det ska vara rimligt att skicka det som indata till funktionen. Vidare tar en rad som säger:

C(:,i) = balance_m(C,i)

Väldigt lång tid med en bra bit över hälften av körtiden. Matrisen C är redan ansatt till rätt storlek innan loopen körs så den ändrar inte storlek. Det är bara värdena inne i den som ändras. Det är för övrigt den funktionen som körs så många gånger och där global-kommandot tar så mycket tid.

Nej de deklareras i två andra funktioner. En som heter load_mesh och en som heter load_constants. Där ansätts massa varieabler i de två strukturerna const och mesh. Dessa variabler behövs som indata för att räkna funktionen och för att slippa att ha tre rader långt anrop på funktionen så är de satta som globala. I början av funktionen balance_m som körs 230000 gånger så måste jag ju säga att mesh och const är globala.

Har inte testat det faktiskt. Det var intressant läsning det där. Ska se om jag kan lägga lite tid på att göra om mina globala variabler till MATLAB handle object som det stod något om i en av länkarna. Får se hur mycket jobb det är bara för fördelen med global är att den är väldigt smidig att använda i kombination med strukturer. Jag har ju rätt olika typer av data och alla delar behöver inte ha tillgång till allt.

Jag har en for-loop som ligger under en while-loop. For-loopen tugger igenom alla kollumner i matrisen från ena sidan till andra en gång och sedan kollar while-loopen efter konvergens. Det är alla loopar som är nästlade. Har lite andra för att sätta upp lite matriser som indata men de körs enbart en gång.

Har du lite mer information om detta? Just nu itererar jag över kolumnerna så i varje loop så räknas en ny kolumn ut. Skulle det hjälpa om jag helt enkelt transponerade hela alltet?

A = rand(3,3)  % A = rand(n,m)
% A_ny = summa av alla grannar: (A_ny(1,1) = A(1,2) + A(2,1), A_ny(2,2) = A(1,2)+A(2,1)+A(1,3)+A(3,1))

rowX = [0 1 0; 1 0 1; 0 1 0];
% diag(ones(n-1,1),-1)+diag(ones(n-1,1),1)

rowY = [0 1 0; 1 0 1; 0 1 0];
% diag(ones(m-1,1),-1)+diag(ones(m-1,1),1) % för större matriser
while delta(A_ny) < eps  
   A_ny = rowX*A + A*rowY;
   M = A_ny - A; 
   A = f(M); 
end %end-while

C_new(2:end-1,2:end-1) = ( C(1:end-2, 2:end-1) + C(3:end, 2:end-1) + C(2:end-1, 1:end-2) +  C(2:end-1, 3:end) ) * 0.25;

Jag har inte satt mig in så mycket djupare i lösningen men det jag fått höra av min handlerade, som jobbat i många år med den här typen av problem, är att det är fullt möjligt att ställa upp mitt problem som ett linjärt ekvationssystem och lösa det rakt av men att det inte är praktiskt när antalet celler blir stort. Just nu kör jag i en dimension och använder mellan 10 och 100 celler för att det är lätt att få en överblick. När det sedan blir i tre dimensioner med miljontals celler så klarar Matlab tydligen inte av det längre. Så det finns en tanke med att köra iterativt förutom vissa pedagogiska anledningar.

Jag ska experimentera lite med att transponera och se vad som händer. Det finns egentligen ingen anledning att ställa upp det på ena eller andra sättet som det ser ut nu mer än att jag började. Det lär dock bli annorlunda senare i 3D då jag kommer behöva räkna med en matris i 4D. Hur vill Matlab ha det då?