układ równań liniowych - rozkład LU

Witam,
muszę napisać drobny program rozwiązujący układ równań liniowych metodą dokładną. Prowadzący polecił mi Rozkład LU, znalazłem zaimplementowany w 4-5 funkcji C# algorytm, rozwiązuje on małe macierze. problemem jest ze zmiennymi, których używa dany algorytm:

double[,] A = new double[n,n];
double[,] L = new double[n,n];
double[,] U = new double[n,n];

Od prowadzącego otrzymałem dane wejściowe w postaci:

wiersz, kolumna, wartość

[pusta linia - oddzielająca macierz wejściową od wektora B] 
wiersz wartość

wektor zawiera 10200 wierszy, zamiast macierzy wejściowej jest 3 kolumnowa tablica z koordynatami i wartością. jest to macierz rzadka - reszta pól to 0 (zera) i nie są brane pod uwagę.
konieczne jest więc przerobienie algorytmu aby poruszał się po "macierzy wejściowej" zgodnie z koordynatami danych wejściowych zamiast wymagać 3 zmiennych tablicowych - każdej z ponad milionem pól

no fajnie... i?

...
czy ktoś na tyle ogarnia, żeby przekształcić algorytm poruszający się po macierzy 10200x10200, tworząc przy okazji 2 nowe macierze (L i U) o tych samych rozmiarach.

Zapomnij o: double[,] A = new double[n,n];
Ty masz macierz rzadką - dużo zer, wiec taka alokacja jest po prostu marnowaniem pamięci i czasu obliczeń (dużo mnożeń przez zero).
Musisz przerobić algorytm tak (przykłady normalnych implementacji znajdziesz w internecie), by operował na macierzy rzadkiej, czyli na strukturze danych, która zapamiętuje w jakich wierszach występują jakie liczby w zadanej kolumnie.

właśnie umieszczam plik w liście tablic

List<double[]> A = new List<double[]>();
for (int i = 0; i < n; i++)
{
       string[] line = strAllLines[i].Split(new[] { '\t' });
       if (line.Length < 4) { break; }
       A.Add(new double[3] { Double.Parse(line[1]), Double.Parse(line[2]), Double.Parse(line[3]) }); //\t0\t0\t0,999999980558641 - pierwszy wiersz(line[0] to null - bo dane są takie a nie inne
}

teraz tylko pytanie jak się poruszać po macierzy zgodnie z algorytmem LU, jeżeli ten wymaga macierze nxn, będę przeszukiwał koordynaty w liście A, zgodnie z iteratorami w pętlach...

public void LUDecompozition(double[,] A, double[,] L, double[,] U, int n)
        {
            for (int i = 0; i < L.GetLength(0); i++)
            {
                L[i, i] = 1;
            }


            double tmp;
            for (int i = 0; i < n; i++)
            {
                for (int j = i; j < n; j++)
                {
                    tmp = 0;
                    for (int k = 0; k < i; k++)
                    {
                        tmp += L[i, k] * U[k, j];
                    }
                    U[i, j] = A[i, j] - tmp;
                }

                for (int j = i + 1; j < n; j++)
                {
                    tmp = 0;
                    for (int k = 0; k < i; k++)
                    {
                        tmp += L[j, k] * U[k, i];
                    }
                    L[j, i] = (A[j, i] - tmp) / U[i, i];
                }
            }
        }

Utwórz sobie obiekt z metodą lub indeksatorem (operator[]) aby tłumaczył tobie odniesienia jak w normalnej macierzy na wewnętrzne indeksy i tyle...

Musisz zaimplementować macierz rzadką tak, aby można było łatwo iterować po wierszach i kolumnach. W tym celu użyj np. SortedDictionary. Dwóch SortedDictionary, w których kluczem jest indeks elementu. Niech jeden słownik sortuje elementy najpierw względem kolumny a później względem wiersza, nim będziesz się posługiwać iterując po wierszach. Drugi niech sortuje odwrotnie - najpierw względem wiersza, później względem kolumny, posłużysz się nim iterując po kolumnach.

z skrajności w skrajność, zużycie RAMu spadło z 12GB+ na 25MB, ale Procesor wyskakuje z lapka...

public class SpCoor
    {
        public Point wsp { get; set; }
        public double val { get; set; }
    }

nie wiem czy najszybsza metoda...

//Lambda Expressions
 SpCoor resultL = L.Find( sc => (sc.wsp == new Point(i, k)));
SpCoor resultU = U.Find(sc => (sc.wsp == new Point(i, k)));
double l = resultL != null ? resultL.val : 0;
double u = resultU != null ? resultU.val : 0;
if (l != 0 && u != 0)
{
    tmp += l * u;
}

Krok pierwszy tworzysz rzadki wektor! Tak jak ci podpowiada adf88 :

SortedDictionary<int, double>

Drugi krok łączysz wektory w macierz rzadką, np tak:

SortedDictionary<int, double>[] spsparseMatrix = new SortedDictionary<int, double>[n];

Zwróć uwagę, że elementy są posortowane, więc łatwo po nich iterować podczas mnożeń macierzy, to co ty zrobiłeś jest bez sensu bo musisz robić poszukiwania określonej komórki.

jest bez sense bo Find po 60k obiektów powtarzany a każdej iteracji to śmiechu warte,
muszę tylko rozkminić jak używać SortedDictionary, bo mając plik z danymi w formacie wiersz, kolumna, wartość, muszę jakoś wrzucić te dane odpowiednio do tej zmiennej..

                matA = new SortedDictionary<int, double>[n];
                for (int i = 0; i < strAllLines.Length - 1; i++)
                {
                    string[] line = strAllLines[i].Split(new[] { '\t' });
                    if (line.Length < 4) { break; }
                    SortedDictionary<int, double> sd = new SortedDictionary<int,double>();
                    
                    sd[int.Parse(line[2])] = Double.Parse(line[3]);
                    matA[int.Parse(line[1])] = sd;
                }

dobrze to zrozumiałem? bo coś mi się zdaje, że matA powinna być inaczej zadeklarowana niż :

SortedDictionary<int, double>[] matA;

Możesz zrobić coś na ten kształt:

public class SparseMatrix
{
    private SortedDictionary<uint, double>[] rows;
    private SortedDictionary<uint, double>[] cols;

    public SparseMatrix(uint n)
        : rows(new SortedDictionary<uint, double>[n])
        , cols(new SortedDictionary<uint, double>[n])
    {
    }

    public double this[uint row, uint col]
    {
        get
        { 
            double ret = 0;
            this.rows[row].TryGetValue(out col);
            return ret;
        }

        set
        { 
            if (IsZero(value)) { // tutaj sprawdzamy czy liczba jest bliska zeru
                this.rows[row].Remove(col) && this.cols[col].Remove(row);
            } else {
                this.rows[row][col] = value;
                this.cols[col][row] = value;
            }
        }
    }

    using RowEnumerator = SortedDictionary<uint, double>.Enumerator;
    using ColEnumerator = SortedDictionary<uint, double>.Enumerator;

    public RowEnumerator GetRowEnumerator(uint row) { return this.rows[row].GetEnumerator(); }
    public ColEnumerator GetColEnumerator(uint col) { return this.cols[col].GetEnumerator(); }
}

public void LUDecompozition(SparseMatrix A, SparseMatrix L, SparseMatrix U)
{
    for (int i = 0; i < L.GetLength(0); i++)
    {
        L[i, i] = 1;
    }

    for (int i = 0; i < n; i++)
    {
        for (int j = i; j < n; j++)
        {
            double tmp = 0;
            SparseMatrix.RowEnumerator l_row = L.GetRowEnumerator(i);
            SparseMatrix.ColEnumerator u_col = U.GetColEnumerator(j);
            bool valid = l_row.MoveNext() && u_col.MoveNext();

            while (valid) {
                if (l_row.Current.Key < u_col.Current.Key) {
                    valid &= l_row.MoveNext();
                } else if (l_row.Current.Key > u_col.Current.Key) {
                    valid &= u_col.MoveNext();
                } else {
                    tmp += l_row.Current.Value * u_col.Current.Value;
                }
            }

            U[i, j] = A[i, j] - tmp;
        }

        for (int j = i + 1; j < n; j++)
        {
            double tmp = 0;
            SparseMatrix.RowEnumerator u_row = U.GetRowEnumerator(i);
            SparseMatrix.ColEnumerator l_col = L.GetColEnumerator(j);
            bool valid = u_row.MoveNext() && l_col.MoveNext();

            while (valid) {
                if (u_row.Current.Key < l_col.Current.Key) {
                    valid &= u_row.MoveNext();
                } else if (u_row.Current.Key > l_col.Current.Key) {
                    valid &= l_col.MoveNext();
                } else {
                    tmp += u_row.Current.Value * l_col.Current.Value;
                }
            }

            L[j, i] = (A[j, i] - tmp) / U[i, i];
        }
    }
}

to pokazałeś, już sprawdzam. z góry już dziękuję za pomoc. zobaczymy jak pójdzie

    public SparseMatrix(uint n)
        : rows(new SortedDictionary<uint, double>[n])
        , cols(new SortedDictionary<uint, double>[n])
    {
    }

nie mogę tego rozgryźć
...
ogólnie jest więcej błędów itd.

                matA = new SortedDictionary<int,SortedDictionary<int,double>>();

                for (int i = 0; i < strAllLines.Length - 1; i++)
                {
                    string[] line = strAllLines[i].Split(new[] { '\t' });
                    if (line.Length < 4) { break; }
                    SortedDictionary<int, double> sd = new SortedDictionary<int,double>();
                    sd.Add(int.Parse(line[2]), Double.Parse(line[3]));

                    try{
                        matA.Add(int.Parse(line[1]),sd);
                    }
                    catch{
                        SortedDictionary<int,double> sdTemp = new SortedDictionary<int, double>(); 
                        matA.TryGetValue(int.Parse(line[1]),out sdTemp);
                        sdTemp.Add(int.Parse(line[2]), Double.Parse(line[3]));
                    }
                }

nie wiem za chiny jak się poruszać po tych słownikach

using (var l_row = matL[i].GetEnumerator()); // to nie dziala
                    using (var u_col = matU[i].GetEnumerator()); // tez nie....
                    //zmienne matL matU matA to słowniki słowników...SortedDictionary<int,SortedDictionary<int,double>>
                    //SparseMatrix.RowEnumerator l_row = matL.GetEnumerator(i);// GetRowEnumerator(i);
                    //SparseMatrix.ColEnumerator u_col = matU.GetColEnumerator(j);
                    
                    bool valid = l_row.MoveNext() && u_col.MoveNext();

                    while (valid)
                    {
                        if (l_row.Current.Key < u_col.Current.Key)
                        {
                            valid &= l_row.MoveNext();
                        }
                        else if (l_row.Current.Key > u_col.Current.Key)
                        {
                            valid &= u_col.MoveNext();
                        }
                        else
                        {
                            tmp += l_row.Current.Value * u_col.Current.Value;
                        }
                    }

                    matU[i][j] = matA[i][j] - tmp;

Pamiętaj, że nie każdą macierz da się rozłożyć czystym rozkładem LU. Poczytaj o pivoted-LU.

thx za informacje, ale nie będę się aż tak zagłębiał. dostałem od prowadzącego macierz, którą da się wyliczyć, i podaruję sobie wszelkie "if'y". a tymczasem nadal mam problem z poruszaniem się po tych słownikach.. nie wiem czy nie będzie konieczne utworzenie po 2 słowniki dla każdej macierzy (wiersz-kolumna oraz kolumna-wiersz), w celu szybszego wyszukiwania - kosztem RAMu). nie wiem nawet za bardzo jak wyszukiwać po słowniku-słowników. może założę temat w dziale c# ( w samej tylko kwestii poruszania się po słownikach)

ps. ciekawy podpis:)

Poprawiłem kod tak, aby się kompilował.
http://4programmers.net/Pastebin/1725

naglik napisał(a):

nie wiem czy nie będzie konieczne utworzenie po 2 słowniki dla każdej macierzy
SparseMatrix ma już 2 słowniki (cols i rows). Są to tablice słowników (po jednym słowniku na każdy wiersz i na każdą kolumnę). Zazwyczaj wszystkie wiersze i kolumny będą posiadać wartości niezerowe (wystarczy, że diagonala jest cała niezerowa). Dlatego są to tablice słowników, a nie słowniki słowników jak to próbowałeś zrobić.

thx bardzo, poradziłem sobie z moim słownikiem słowników bez problemu :) ale masz rację - niepotrzebnie. jak uda mi się odpalić drugi thread, który będzie sprawdzał zużycie CPU, RAMu, i ogólnie czas wykonania to wezmę się za optymalizację... z tym, że jest to tylko projekt do oddania na ocenę - więc nie muszę tego jakoś specjalnie opracowywać...

naglik napisał(a):

jak uda mi się odpalić drugi thread, który będzie sprawdzał zużycie CPU, RAMu, i ogólnie czas wykonania to wezmę się za optymalizację...

Nie możesz tego zrobić Menedżerem Zadań?

miałbym wtedy za mało textbox'ów w aplikacji, i w ogóle byłaby za biedna:), zrobiłem sobie timer i bawię się w backgroundworker'ach. jedynie co muszę zrobić to w wymaganiach napisać: zalecane 2 rdzeniowy procesor, i 30MB ramu

/edit
macierz 100x100 - 80 sekund...
macierz 10200x10200 - 5389 sekund (89,8 minuty)

ps. pól z wartościami różnymi od 0 było 50539

        private void backgroundWorker2_DoWork(object sender, DoWorkEventArgs e)
        {
            lud = new LUDecompozition(matA, n, sender as BackgroundWorker, e as DoWorkEventArgs);
        }

gra i buczy
trochę zoptymalizowałem kod. odpalę i sprawdzę ile teraz potrwa rozwiązanie tych równań

ps2. nie wiem czy jest się po co chwalić tutaj kodem. nie mam problemu z udostępnieniem go. może komuś się przyda, albo ktoś go przejrzy sprawdzi i poprawi :)