Preliminares

A solução numérica de sistemas lineares de equações é uma etapa indispensável em muitas áreas da matemática aplicada, engenharia e indústria de TI, seja trabalhando com gráficos, calculando a aerodinâmica de um avião ou otimizando a logística. A agora "máquina" na moda também não teria progredido muito sem ela. Além disso, a solução de sistemas lineares, como regra, consome a maior porcentagem de todos os custos de computação. É claro que esse componente crítico requer otimização da velocidade máxima.

Muitas vezes, trabalha com os chamados matrizes esparsas - aquelas cujos zeros são ordens de magnitude maiores que outros valores. Isso, por exemplo, é inevitável se você estiver lidando com equações diferenciais parciais ou em outros processos nos quais os elementos que surgem nas relações lineares definidoras são conectados apenas a "vizinhos". Aqui está um exemplo possível de uma matriz esparsa para a equação de Poisson unidimensional conhecida na física clássica $$$- \ phi ^ {''} = f$ em uma grade uniforme (sim, embora não haja muitos zeros nela, mas ao moer a grade, eles serão saudáveis!):

$$

$$\ begin {pmatrix} 2 & -1 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & -1 & 2 \ end {pmatrix}$$

As matrizes opostas a elas - aquelas nas quais elas não prestam atenção ao número de zeros e levam em consideração todos os componentes, sem exceção - são chamadas de densas.

Matrizes esparsas geralmente são apresentadas, por várias razões, em um formato de coluna compactada - CSC (Compressed Sparse Column). Esta descrição usa duas matrizes inteiras e um ponto flutuante. Deixe a matriz ter tudo $$$nnz_A$ diferente de zero e $$$N$ colunas. Os elementos da matriz são listados em colunas da esquerda para a direita, sem exceções. Primeira matriz $$$i_A$ comprimentos $$$nnz_A$ contém números de linhas de componentes de matriz diferentes de zero. Segunda matriz $$$j_A$ comprimentos $$$N + 1$ contém ... não, não os números das colunas, porque a matriz seria escrita no formato de coordenadas (coordenadas) ou ternário (trigêmeo). E a segunda matriz contém os números de série dos componentes da matriz com os quais as colunas começam, incluindo uma coluna fictícia adicional no final. Finalmente, a terceira matriz $$$v_A$ comprimentos $$$nnz_A$ já contém os próprios componentes, na ordem correspondente à primeira matriz. Aqui, por exemplo, sob a suposição de que a numeração de linhas e colunas de matrizes começa do zero, para uma matriz específica

$$

$$A = \ begin {pmatrix} 0 & 1 & 0 & 4 & -1 \\ 7 & 0 & 2.1 & 0 & 3 \\ 0 & 0 & 0 & 10 & 0 \ end {pmatrix}$$

matrizes serão $$$i_A = \ {1, 0, 1, 0, 2, 0, 1 \}$ , $$$j_A = \ {0, 1, 2, 3, 5, 7 \}$ , $$$v_A = \ {7, 1, 2,1, 4, 10, -1, 3 \}$ .

Os métodos para resolver sistemas lineares são divididos em duas grandes classes - direta e iterativa. As linhas retas são baseadas na possibilidade de representar a matriz como um produto de duas matrizes mais simples, para então dividir a solução em dois estágios simples. A iteração usa as propriedades únicas dos espaços lineares e trabalha no antigo como o método mundial de aproximação sucessiva de incógnitas à solução desejada e, no processo de convergência em si, as matrizes são usadas, em regra, apenas para multiplicá-las por vetores. Os métodos iterativos são muito mais baratos que os diretos, mas funcionam lentamente em matrizes mal condicionadas. Quando a confiabilidade do concreto armado é importante - use linhas retas. Aqui estou eu e quero tocar um pouco.

Digamos que para uma matriz quadrada $$$A$ decomposição do formulário está disponível para nós $$$A = LU$ onde $$$L$ e $$$U$ - respectivamente, as matrizes triangular inferior e triangular superior, respectivamente. O primeiro significa que ele tem um zeros acima da diagonal, o segundo significa que está abaixo da diagonal. Como exatamente conseguimos essa decomposição - não estamos interessados ​​agora. Aqui está um exemplo simples de decomposição:

$$

$$\ begin {pmatrix} 1 & -1 & -1 \\ 2 & - 1 & -0,5 \\ 4 & -2 & -1,5 \ end {pmatrix} = \ begin {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ begin {pmatrix} 1 & -1 & -1 \\ 0 & 1 & 1,5 \\ 0 & 0 & -0,5 \ end {pmatrix}$$

Como, neste caso, resolver o sistema de equações $$$Ax = f$ por exemplo, com o lado direito $$$f = \ begin {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix}$ ? O primeiro estágio é um movimento direto (resolução para frente = substituição para frente). Nós denotamos $$$y: = Ux$ e trabalhar com o sistema $$$Ly = f$ . Porque $$$L$ - triangular inferior, sucessivamente no ciclo encontramos todos os componentes $$$y$ de cima para baixo:

$$

$$\ begin {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ begin {pmatrix} y_1 \\ y_2 \\ y_3 \ end {pmatrix} = \ begin {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix} \ implica y = \ begin {pmatrix} 4 \\ -6 \\ -1 \ end {pmatrix}$$

A idéia central é que, ao encontrar $$$i$ componentes vetoriais $$$y$ é multiplicado por uma coluna com o mesmo número de matriz $$$L$ que é subtraído do lado direito. A própria matriz parece entrar em colapso da esquerda para a direita, diminuindo de tamanho à medida que mais e mais componentes vetoriais são encontrados $$$y$ . Esse processo é chamado eliminação de coluna.

O segundo estágio é um movimento para trás (resolução para trás = substituição para trás). Ter um vetor encontrado $$$y$ nós decidimos $$$Ux = y$ . Aqui já analisamos os componentes de baixo para cima, mas a ideia permanece a mesma: $$$i$ a coluna é multiplicada pelo componente encontrado $$$x_i$ e se move para a direita e a matriz entra em colapso da direita para a esquerda. Todo o processo é ilustrado para a matriz mencionada no exemplo na figura abaixo:

$$

$$\ small \ begin {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ begin {pmatrix} y_1 \\ y_2 \\ y_3 \ end {pmatrix} = \ begin {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix} \ implica \ begin {pmatrix} 1 & 0 \\ 2 & 1 \ end {pmatrix} \ begin {pmatrix} y_2 \\ y_3 \ end {pmatrix } = \ begin {pmatrix} -6 \\ -13 \ end {pmatrix} \ implica \ begin {pmatrix} 1 \ end {pmatrix} \ begin {pmatrix} y_3 \ end {pmatrix} = \ begin {pmatrix} -1 \ end {pmatrix}$$

$$

$$\ small \ begin {pmatrix} 1 & -1 & -1 \\ 0 & 1 & 1.5 \\ 0 & 0 & -0.5 \ end {pmatrix} \ begin {pmatrix} x_1 \\ x_2 \\ x_3 \ end { pmatrix} = \ begin {pmatrix} 4 \\ -6 \\ -1 \ end {pmatrix} \ Rightarrow \ begin {pmatrix} 1 & -1 \\ 0 & 1 \ end {pmatrix} \ begin {pmatrix} x_1 \ \ x_2 \ end {pmatrix} = \ begin {pmatrix} 6 \\ -9 \ end {pmatrix} \ implica \ begin {pmatrix} 1 \ end {pmatrix} \ begin {pmatrix} x_1 \ end {pmatrix} = \ begin {pmatrix} -3 \ end {pmatrix}$$

e nossa decisão será $$$x = \ begin {pmatrix} -3 \\ - 9 \\ 2 \ end {pmatrix}$ .

Se a matriz é densa, ou seja, é completamente representada na forma de algum tipo de matriz, unidimensional ou bidimensional, e o acesso a um elemento específico nela ocorre ao longo do tempo $$$O (1)$ , então um procedimento de solução semelhante com a decomposição existente não é difícil e é fácil de codificar, portanto, não perderemos tempo com isso. E se a matriz for escassa? Aqui, em princípio, também não é difícil. Aqui está o código C ++ para o avanço em que a solução $$$x$ Ele foi escrito para o local no lado direito, sem verificar a correção dos dados de entrada (os arrays CSC correspondem à matriz triangular inferior):

Algoritmo 1:

void forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x) { size_t j, p, n = x.size(); for (j = 0; j < n; j++) //    { x[j] /= vA[jA[j]]; for (p = jA[j]+1; p < jA[j+1]; p++) x[iA[p]] -= vA[p] * x[j] ; } } 

Todas as discussões posteriores dizem respeito apenas ao curso para frente para resolver o sistema triangular inferior $$$Lx = f$ .

Gravata

Mas e se o lado direito, ou seja, vetor à direita do sinal de igual $$$Lx = f$ Ele próprio possui um grande número de zeros? Então faz sentido pular os cálculos associados às posições zero. As alterações no código neste caso são mínimas:

Algoritmo 2:

 void fake_sparse_forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x) { size_t j, p, n = x.size(); for (j = 0; j < n; j++) //    { if(x[j]) { x[j] /= vA[jA [j]]; for (p = jA[j]+1; p < jA[j+1]; p++) x[iA[p]] -= vA[p] * x[j]; } } } 

A única coisa que adicionamos é a if , cujo objetivo é reduzir o número de operações aritméticas ao número real. Se, por exemplo, todo o lado direito consistir em zeros, você não precisará contar nada: a solução será o lado direito.

Tudo parece ótimo e, é claro, funcionará, mas aqui, após um longo período de preliminares, o problema é finalmente visível - o desempenho assintoticamente baixo desse solucionador no caso de grandes sistemas. E tudo por causa de um loop for . Qual é o problema? Mesmo que a condição dentro do if se torne extremamente rara, não há como fugir do próprio ciclo, e isso cria a complexidade do algoritmo $$$O (N)$ onde $$$N$ É o tamanho da matriz. Não importa como os loops sejam otimizados pelos compiladores modernos, essa complexidade se fará sentir com grandes $$$N$ . Especialmente ofensivo quando todo o vetor $$$f$ consiste inteiramente de zeros, porque, como dissemos, não faça nada neste caso! Nem uma única operação aritmética! Que diabos $$$O (N)$ ação?

Bem, digamos. Mesmo assim, por que você não pode simplesmente suportar o modo ocioso, porque cálculos reais com números reais, ou seja, aqueles que se enquadram if ainda serão poucos? Mas o fato é que esse procedimento de fluxo direto com o lado direito rarefeito é frequentemente usado em ciclos externos e subjacente à decomposição de Cholesky $$$A = LL ^ T$ e decomposição LU à esquerda. Sim, sim, uma dessas expansões, sem a capacidade de fazer o que todos esses movimentos diretos e reversos na solução de sistemas lineares, perdem todo o significado prático.

Teorema Se a matriz for simétrica positiva definida (SPD), ela poderá ser representada como $$$A = LL ^ T$ a única maneira de $$$L$ - matriz triangular inferior com elementos positivos na diagonal.

Para matrizes SPD altamente esparsas, é usada a decomposição de Cholesky. Representando esquematicamente a decomposição em forma de bloco de matriz

$$

$$\ begin {pmatrix} \ mathbf {L} _ {11} e \ mathbf {0} \\ \ mathbf {l} _ {12} ^ T & l_ {22} \ end {pmatrix} \ begin {pmatrix} \ mathbf {L} _ {11} ^ T & \ mathbf {l} _ {12} \\ \ mathbf {0} & l_ {22} \ end {pmatrix} = \ begin {pmatrix} \ mathbf {A} _ { 11} & \ mathbf {a} _ {12} \\ \ mathbf {a} _ {12} ^ T & a_ {22} \ end {pmatrix},$$

todo o processo de fatoração pode ser logicamente dividido em apenas três etapas.

Algoritmo 3:

1. método superior de decomposição de Cholesky de menor dimensão $$$\ mathbf {L} _ {11} \ mathbf {L} _ {11} ^ T = \ mathbf {A} _ {11}$ (recursão!)
2. sistema triangular inferior com lado direito escasso $$$\ mathbf {L} _ {11} \ mathbf {l} _ {12} = \ mathbf {a} _ {12}$ ( aqui está! )
3. cálculo $$$l_ {22} = \ sqrt {a_ {22} - \ mathbf {l} _ {12} ^ T \ mathbf {l} _ {12}}$ (operação trivial)

Na prática, isso é implementado de maneira que as etapas 3 e 2 sejam executadas em um ciclo grande e nesta ordem. Assim, a matriz é executada na diagonal de cima para baixo, aumentando a matriz $$$L$ linha por linha em cada iteração do loop.

Se, em um tal agloritmo, a complexidade avançada na etapa 2 for $$$O (N)$ onde $$$N$ É o tamanho da matriz triangular inferior $$$\ mathbf {L} _ {11}$ em uma iteração arbitrária de um ciclo grande, a complexidade de toda a decomposição será pelo menos $$$O (N ^ 2)$ ! Ah, eu não gostaria disso!

Estado da arte

Muitos algoritmos são de alguma forma baseados na imitação de ações humanas na resolução de problemas. Se você der a uma pessoa um sistema linear triangular inferior com o lado direito, no qual apenas 1-2 são diferentes de zero, é claro que ele primeiro executa o vetor do lado direito com os olhos de cima para baixo (esse ciclo maldito de complexidade $$$O (N)$ ! ) para encontrar esses nonzeros. Depois, ele os usará apenas, sem perder tempo com os componentes zero, porque o último não afetará a solução: não faz sentido dividir zeros nos componentes diagonais da matriz, além de mover a coluna multiplicada por zero para a direita. Este é o algoritmo acima 2. Não há milagres. Mas e se uma pessoa receber imediatamente o número de componentes diferentes de zero de algumas outras fontes? Por exemplo, se o lado direito for uma coluna de alguma outra matriz, como é o caso da decomposição de Cholesky, teremos acesso instantâneo a seus componentes diferentes de zero, se solicitado em sequência:

 //     j,     : for (size_t p = jA[j]; p < jA[j+1]; p++) //    vA[p]     iA[p]   j 

A complexidade desse acesso é $$$O (nnz_j)$ onde $$$nnz_j$ É o número de componentes diferentes de zero na coluna j. Graças a Deus pelo formato CSC! Como você pode ver, ele não é usado apenas para economizar memória.

Nesse caso, podemos captar a essência do que está acontecendo ao resolver pelo método do curso direto $$$Lx = f$ para matrizes esparsas $$$L$ lado direito $$$f$ . Prenda a respiração: tomamos um componente diferente de zero $$$f_i$ no lado direito, encontramos a variável correspondente $$$x_i$ dividindo por $$$L_ {ii}$ e, multiplicando a i-ésima coluna por essa variável encontrada, introduzimos não-zeros adicionais no lado direito, subtraindo essa coluna do lado direito! Este processo está bem descrito na linguagem dos ... gráficos. Além disso, orientado e não cíclico.

Para uma matriz triangular inferior que não possui zeros na diagonal, definimos um gráfico conectado. Assumimos que a numeração das linhas e colunas começa do zero.
Definição de Um gráfico de conectividade para uma matriz triangular inferior $$$L$ o tamanho $$$N$ , que não possui zeros na diagonal, é chamado de conjunto de nós $$$V = \ {0, ..., N-1 \}$ e arestas orientadas $$$E = \ {(j, i) | L_ {ij} \ ne 0 \}$ .

Aqui, por exemplo, é a aparência do gráfico de conexão para uma matriz triangular inferior.

em que os números simplesmente correspondem ao número ordinal do elemento diagonal e os pontos indicam componentes diferentes de zero abaixo da diagonal:



Definição de Gráfico orientado para acessibilidade $$$G$ em vários índices $$$W \ subconjunto V$ chamar tantos índices $$$R_G (W) \ subconjunto V$ que em qualquer índice $$$z \ em R_G (W)$ você pode passar passando pelo gráfico de algum índice $$$w (z) \ em W$ .

Exemplo: para um gráfico de uma imagem $$$R_G (\ {0, 4 \}) = \ {0, 4, 5, 6 \}$ . É claro que sempre mantém $$$W \ subconjunto R_G (W)$ .

Se representarmos cada nó do gráfico como o número da coluna da matriz que o gerou, os nós vizinhos para os quais suas bordas são direcionadas corresponderão aos números de linha de componentes diferentes de zero nesta coluna.

Vamos $$$nnz (x) = \ {j | x_j \ ne 0 \}$ denota o conjunto de índices correspondentes a posições diferentes de zero no vetor x.

Teorema de Hilbert (não, não por cujo nome os espaços são nomeados) $$$nnz (x)$ onde $$$x$ existe um vetor de solução de um sistema triangular inferior esparso $$$Lx = f$ com lado direito escasso $$$f$ coincide com $$$R_G (nnz (f))$ .

Além disso: no teorema, não levamos em conta a possibilidade improvável de que os números diferentes de zero no lado direito, ao destruir as colunas, possam ser reduzidos a zero, por exemplo, 3 - 3 = 0. Esse efeito é chamado de cancelamento numérico. Levar em conta esses zeros que ocorrem espontaneamente é uma perda de tempo, e eles são percebidos como todos os outros números em posições diferentes de zero.

Um método eficaz de conduzir um movimento direto com um dado lado direito escasso, sob a suposição de que temos acesso direto a seus componentes diferentes de zero por índices, é o seguinte.

1. Percorremos o gráfico usando o método "pesquisa em profundidade", iniciando sequencialmente a partir do índice $$$i \ in nnz (f)$ cada componente diferente de zero do lado direito. Gravando nós encontrados em uma matriz $$$R_G (nnz (f))$ enquanto fazemos isso na ordem em que retornamos ao gráfico! Por analogia com o exército de invasores: ocupamos o país sem crueldade, mas quando eles começaram a nos levar de volta, nós, furiosos, destruímos tudo em seu caminho.
   Vale ressaltar que não é absolutamente necessário que a lista de índices $$$nnz (f)$ foi classificado por ascensão quando foi alimentado com a entrada do algoritmo "busca pela profundidade da primeira pesquisa". Você pode começar em qualquer ordem no aparelho $$$nnz (f)$ . Sequência diferente pertencente ao conjunto $$$nnz (f)$ índices não afeta a solução final, como veremos no exemplo.
   
   Esta etapa não requer nenhum conhecimento de uma matriz real. $$$v_A$ ! Bem-vindo ao mundo da análise simbólica ao trabalhar com solucionadores esparsos diretos!
   
2. Prosseguimos para a solução em si, tendo à disposição uma matriz $$$R_G (nnz (f))$ da etapa anterior. Destruímos as colunas na ordem inversa do registro dessa matriz . Voila!

Exemplo

Considere um exemplo no qual as duas etapas são demonstradas em detalhes. Suponha que tenhamos uma matriz 12 x 12 da seguinte forma:

O gráfico de conectividade correspondente tem o formato:

Deixe o diferente de zero na parte direita estar apenas nas posições 3 e 5, ou seja, $$$nnz (f) = \ {3, 5 \}$ . Vamos percorrer o gráfico, começando com esses índices na ordem escrita. Em seguida, o método de "pesquisa profunda" terá a seguinte aparência. Primeiro vamos visitar os índices em ordem $$$3 \ a 8 \ a 11$ , não esquecendo de marcar os índices como visitados. Na figura abaixo, eles estão sombreados em azul:

Ao retornar, inserimos índices em nossa matriz de índices de componentes de solução diferentes de zero $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} 3 \}$ . Em seguida, tente executar $$$5 \ to8 \ to ...$ , mas encontramos um nó 8 já marcado, por isso não tocamos nessa rota e vamos para a ramificação $$$5 \ to9 \ to ...$ . O resultado dessa execução será $$$5 \ to9 \ to10$ . Não podemos visitar o nó 11, pois ele já está rotulado. Então, volte e complemente a matriz $$$nnz (x)$ nova captura marcada em verde: $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} 3, \ color {green} {10}, \ color {green} 9, \ color {verde} 5 \}$ . E aqui está a foto:

Os nós verdes enlameados 8 e 11 são os que queríamos visitar durante a segunda corrida, mas não pudemos, porque já visitou durante o primeiro.

Então a matriz $$$nnz (x) = R_G (nnz (f))$ formado. Passamos para o segundo passo. Um passo simples: percorremos a matriz $$$nnz (x)$ na ordem inversa (da direita para a esquerda), localizando os componentes correspondentes do vetor de solução $$$x$ divisão pelos componentes diagonais da matriz $$$L$ e movendo as colunas para a direita. Outros componentes $$$x$ como eram zeros, continuavam assim. Esquematicamente, isso é mostrado abaixo, onde a seta inferior indica a ordem na qual as colunas são destruídas:

Nota: nesta ordem de destruição de colunas, o número 3 será encontrado depois dos números 5, 9 e 10! As colunas não são destruídas em ordem crescente, o que seria um erro para as densas, ou seja, matrizes incompletas. Mas para matrizes esparsas, isso é comum. Como pode ser visto na estrutura da matriz diferente de zero neste exemplo, o uso das colunas 5, 9 e 10 na coluna 3 não distorcerá os componentes na resposta $$$x_5$ , $$$x_9$ , $$$x_ {10}$ de forma alguma, eles têm c $$$x_3$ apenas "sem cruzamentos". Nosso método levou isso em consideração! Da mesma forma, o uso da coluna 8 após a coluna 9 não estragará o componente $$$x_9$ . Ótimo algoritmo, não é?

Se percorrermos o gráfico primeiro a partir do índice 5 e depois do índice 3, nossa matriz será $$$nnz (x): = \ {\ color {blue} {11}, \ color {blue} 8, \ color {blue} {10}, \ color {blue} {9}, \ color {blue} 5, \ color {green} 3 \}$ . Destruir as colunas na ordem inversa dessa matriz fornecerá exatamente a mesma solução que no primeiro caso!

A complexidade de todas as nossas operações é dimensionada pelo número de operações aritméticas reais e pelo número de componentes diferentes de zero no lado direito, mas não pelo tamanho da matriz! Atingimos nosso objetivo.

Crítica

NO ENTANTO! Um leitor crítico notará que a suposição inicial de que temos “acesso direto aos componentes diferentes de zero do lado direito pelo índice” já significa que uma vez antes de executarmos o lado direito de cima para baixo para encontrar esses mesmos índices e organizá-los na matriz $$$nnz (f)$ isto é, já gastaram $$$O (N)$ ação! Além disso, o próprio gráfico exige que alocemos primeiro a memória do tamanho máximo possível (em outro lugar, é necessário anotar os índices encontrados pesquisando em profundidade!) Para não sofrer realocação adicional à medida que a matriz cresce $$$nnz (x)$ , e isso também requer $$$O (N)$ operações! Por que, então, eles dizem, todo esse barulho?

Mas, de fato, para uma solução única de um sistema triangular inferior esparso com o lado direito esparso, originalmente definido como um vetor denso, não faz sentido desperdiçar tempo do desenvolvedor em todos os algoritmos mencionados acima. Eles podem ser ainda mais lentos que o método da testa apresentado pelo algoritmo 2 acima. Mas, como já mencionado, este aparelho é indispensável para a fatoração de Cholesky, portanto você não deve jogar tomates em mim. De fato, antes de executar o algoritmo 3, toda a memória necessária de tamanho máximo é alocada imediatamente e requer $$$O (N)$ hora. Em outro ciclo de coluna $$$A$ todos os dados são substituídos apenas em uma matriz de comprimento fixo $$$N$ , e apenas nas posições em que essa reescrita é relevante, graças ao acesso direto a elementos diferentes de zero. E é justamente por isso que a eficiência surge!

Implementação C ++

O próprio gráfico como uma estrutura de dados no código não é necessário para compilar. É suficiente usá-lo implicitamente ao trabalhar diretamente com a matriz. Toda a conectividade necessária será levada em consideração pelo algoritmo. Obviamente, a pesquisa de profundidade é convenientemente implementada usando recursão banal. Aqui, por exemplo, parece uma pesquisa de profundidade recursiva com base em um único índice inicial:

 /* j -   iA, jA -    ,    CSC top -     nnz(x);     0 result -   N,    nnz(x)   0  top-1  marked -    N;      */ void DepthFirstSearch(size_t j, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t& top, std::vector<size_t>& result, std::vector<int>& marked) { size_t p; marked[j] = 1; //   j   for (p = jA[j]; p < jA[j+1]; p++) //       j { if (!marked[iA[p]]) //  iA[p]   { DepthFirstSearch(iA[p], iA, jA, top, result, marked); //      iA[p] } } result[top++] = j; //  j   nnz(x) } 

Se você passar a variável top igual a zero na primeira chamada para DepthFirstSearch , após a conclusão de todo o procedimento recursivo, a variável top será igual ao número de índices encontrados na matriz de result . Trabalho de casa: escreva outra função que pegue uma matriz de índices de componentes diferentes de zero no lado direito e os passe sequencialmente para a função DepthFirstSearch . Sem isso, o algoritmo não está completo. Atenção: uma matriz de números reais $$$vA$ não passamos para a função, porque não é necessário no processo de análise simbólica.

Apesar de sua simplicidade, um defeito na implementação é óbvio: para sistemas grandes, o excesso de pilha está ao virar da esquina. Bem, então há uma opção baseada em um loop em vez de recursão. É mais difícil de entender, mas já adequado para todas as ocasiões:

 /* j -   N -   iA, jA -    ,    CSC top -     nnz(x) result -   N,     nnz(x)   top  ret-1 ,  ret -    marked -    N work_stack -     N */ size_t DepthFirstSearch(size_t j, size_t N, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t top, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack) { size_t p, head = N - 1; int done; result[N - 1] = j; //    while (head < N) { j = result[head]; //  j     if (!marked[j]) { marked[j] = 1; //   j   work_stack[head] = jA[j]; } done = 1; //    j      for (p = work_stack[head] + 1; p < jA[j+1]; p++) //    j { //   c   iA[p] if (marked[iA[p]]) continue; //  iA[p]  ,   work_stack[head] = p; //       j result[--head] = iA[p]; //       iA[p] done = 0; //   j    break; } if (done) //      j  { head++; //  j    result[top++] = j; //  j    } } return (top); } 

E é assim que o gerador de uma estrutura vetorial de solução diferente de zero já se parece $$$nnz (x)$ :

 /* iA, jA -   ,    CSC iF -        f nnzf -      f result -   N,    nnz(x)   0  ret-1 ,  ret -    marked -    N,    work_stack -     N */ size_t NnzX(const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<size_t>& iF, size_t nnzf, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack) { size_t p, N, top; N = jA.size() - 1; top = 0; for (p = 0; p < nnzf; p++) if (!marked[iF[p]]) //        top = DepthFirstSearch(iF[p], N, iA, jA, vA, top, result, marked, work_stack); for (p = 0; p < top; p++) marked[result[p]] = 0; //   return (top); } 

Finalmente, combinando tudo em um, escrevemos o solucionador triangular inferior para o caso do lado direito rarefeito:

 /* iA, jA, vA -  ,    CSC iF -        f nnzf -      f vF -       f result -   N,    nnz(x)   0  ret-1 ,  ret -    marked -    N,    work_stack -     N x -    N,    ;     */ size_t lower_triangular_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, const std::vector<size_t>& iF, size_t nnzf, const std::vector<double>& vF, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack, std::vector<double>& x) { size_t top, p, j; ptrdiff_t px; top = NnzX(iA, jA, iF, nnzf, result, marked, work_stack); for (p = 0; p < nnzf; p++) x[iF[p]] = vF[p]; //    for (px = top; px > -1; px--) //     { j = result[px]; // x(j)   x [j] /= vA[jA[j]]; //   x(j) for (p = jA[j]+1; p < jA[j+1]; p++) { x[iA[p]] -= vA[p]*x[j]; //  j-  } } return (top) ; } 

Vemos que nosso ciclo é executado apenas através dos índices da matriz $$$nnz (x)$ , não todo o conjunto $$ . Feito!

Existe uma implementação que não usa uma matriz de tags markedpara economizar memória. Em vez disso, um conjunto adicional de índices é usado.$$$V_1$ não se cruza com o conjunto $$$V = \{0, ..., N - 1\}$no qual existe um mapeamento individual por uma operação algébrica simples como um procedimento de marcação de nó. No entanto, em nossa era da memória barata, economizando-a em uma única matriz de comprimento$$$N$ parece completamente redundante.

Em conclusão

O processo de solução de um sistema esparso de equações lineares pelo método direto é dividido em três etapas:

1. Análise de caracteres
2. Fatoração numérica baseada nos dados da análise sivol
3. A solução dos sistemas triangulares obtidos com um lado denso do lado direito

O segundo passo, a fatoração numérica, é a parte que consome mais recursos e consome a maior parte (> 90%) do tempo estimado. O objetivo do primeiro passo é reduzir o alto custo do segundo. Um exemplo de análise simbólica foi apresentado neste post. No entanto, é o primeiro passo que requer o maior tempo de desenvolvimento e custos mentais máximos por parte do desenvolvedor. Uma boa análise simbólica requer conhecimento da teoria de gráficos e árvores e posse de um "instinto algorítmico". O segundo passo é incomparavelmente mais fácil de implementar.
Bem, o terceiro passo, tanto na implementação quanto no tempo estimado, é na maioria dos casos o mais despretensioso.

Uma boa introdução aos métodos diretos para sistemas esparsos está em Tim Davis , professor do Departamento de Ciência da Computação da Universidade A&M do Texas , intitulado "Métodos diretos para sistemas lineares esparsos ". O algoritmo descrito acima é descrito lá. Não estou familiarizado com Davis, embora eu próprio trabalhasse na mesma universidade (embora em uma faculdade diferente). Se não me engano, o próprio Davis participou do desenvolvimento de solucionadores (!). usado em Matlab Além disso, ele estava envolvido até mesmo para geradores de energia imagem no Street View do Google (OLS), a propósito, na mesma faculdade nenhum listado menos que o próprio Bjarne Stroustrup - .. o criador do C ++

Talvez um dia Escreverei algo mais sobre matrizes esparsas ou métodos numéricos em geral. Bom para todos!