Axel Ruhe
18 januari 2005
Tillämpade numeriska metoder II, 2D1250, våren 2005

LABORATION 1, Uppgift 2

Glesa Matriser, direkta metoder

En introduktion till glesa matriser ges i det utdelade materialet, Topics in Numerical Linear Algebra, kapitel 1.  

Matlabs
help/ MATLAB/Using Matlab/Mathematics/Sparse matrices
ger även en ganska instruktiv beskrivning. Några detaljer har jag utrett i Extra tips.

Den som verkligen vill veta mera om varför Matlabs rutiner ser ut just som de gör, kan läsa originaluppsatsen:
Gilbert, John R., Cleve Moler, and Robert Schreiber, "Sparse Matrices in MATLAB: Design and Implementation," SIAM J. Matrix Anal. Appl., Vol. 13, No. 1, January 1992, pp. 333-356.  gilbert92sparse.ps gilbert92sparse.pdf

Testmatriser hittar man i Matrix Market hos National Institute of Standards and Technology (NIST).

1. Enkelt exempel: Jämför omordningar!

Vi börjar med att  studera den matris man får vid finit differensapproximation av Laplace ekvation. Generera ett nät med G=numgrid(,)  och en matris med A=delsq(G)! Titta på den glesa matrisen med kommandot spy(A).

Ta lämplig storlek och använd  antingen en kvadrat, det blir enklast, eller ett L membran, lite mera komplicerat. Det är ett L membran man ser i Matlabs skrifter.
 
Studera nu hur olika omordningar av matrisen ser ut. Man kan räkna fram permutationsvektorer för bandviddsreduktion (RCM) och "minimum degree" (MMD) med anropen pr=symrcm(A) och pm=symmmd(A). Substrukturering, "Nested dissection" finns inte som rutin, det finns en särskild option i numgrid som ger den ordningen för en kvadrat. Titta på de omordnade matriserna med spy  och jämför med den bild jag visade på föreläsningen!

I grafen G ser man hur noderna räknas upp från början, se efter att det blir kolonnvis. Gör en bild av den omordnade grafen G för Reversed Cuthill Mc Kee. Det går att göra med rutinen v2g som beskrivs i extra tips. Observera att v2g inte kan anropas direkt med pm som parameter, vi måste ha den inversa permutationen, låt oss kalla den mp. Den fås enkelt genom anropet mp(pm)=iv, där iv är en vektor iv=1:n där n är ordningen på A.

När jag fick ordning på RCM så blev det precis som i föreläsningen. Se efter om Du lyckas med det!

Gör samma sak med MMD. Här liknar det inte alls det vi kommit fram till.

Pröva ett par större matriser. Kan man säga något om hur antalet påfyllda element ökar med ordningen n?

2. Allmänna glesa matriser

Jag har lagrat in ett par ganska stora matriser som Matlabfiler.

Välj några stycken från listan nedan!


1. Titta på matrisen och se hur mycket som fylls på vid Gausselimination för olika omordningar.
2. Hur mycket mindre blir det med RCM eller MMD jämfört med ingen omordning?
3. Är det någon märkbar skillnad mellan RCM och MMD?

Se upp, ser man sig inte för kan det bli alldeles för stort och fylla minnet!
Skriv inte ut resultatet av spy i redogörelsen, det kan bli förfärligt utrymmeskrävande.


OBS:
.m filer exekveras
.mat filer laddas in med kommandot load
.mtx filer är i Matrix  Market format och laddas in av rutinen A=mmread('filnamn') Om rutinen mmread inte finns kan den hämtas hem via Matrix Market .

Några testmatriser :

  1. fleigAB.mat styvhets och massmatris för ett egenvärdesproblem modellerat med finita elementmetoden med triangulära element. Skall likna de differensapproximationer vi hade i deluppgift 1 ganska mycket.
  2. Kör en 3 dimensionell Laplace  genererad med delsq3d.m t ex för 10*10*20 punkter.
  3. Ett kraftledningsnätverk för västra USA bcspwr
  4. Vi har en VLSI krets för en klocka i clock.mat Bara den ena matrisen G har tillräckligt många element för att vara intressant, jag tror varje element står för en resistans. En VLSI krets liknar i stora delar ett träd, det blir få slutna cykler. Därför blir matriserna mycket glesa och en algoritm som bygger på nivåstrukturer som RCM ger mycket liten påfyllnad.
  5. Leta gärna rätt på någon annan matris hos Matrix Market, det finns från alla möjliga tillämpningsområden!