Soluzione di sistemi

Quando si risolvono sistemi, come $A{\bf x}={\bf b}$, tutte le operazioni svolte su A vanno ripetute pari pari sul vettore di termini noti. Ad esempio, dentro il ciclo i precedente, andrebbe inserito

        /* Esegue combinazione lineare sul vettore dei termini noti */
        b[i] = b[i] - m*b[k];

In questo modo il sistema originale viene via via trasformato in sistemi equivalenti. Alla fine del procedimento il vettore b (uguale a ${\bf b}$ all'inizio) sarà stato trasformato in ${\bf c}=L^{-1}{\bf b}$. Per completare la soluzione del sistema si risolve infine il sistema triangolare $U{\bf x}={\bf c}$ con la back substitution descritta più avanti.

Come ulteriore esempio, consideriamo il calcolo dell'inversa. In questo caso B è una matrice, che va inizializzata all'identità. Ad ogni passo dovremo poi agire su un'intera riga di B, ad esempio con un altro ciclo j:

        /* Esegue combinazione lineare sulla matrice dei termini noti */
        for (j=0; j<n; j++)
            B[i][j] = B[i][j] - m*B[k][j];

Si può osservare che nel caso del calcolo dell'inversa la matrice B rimane sempre triangolare (durante questa fase di eliminazione), e quindi il precedente ciclo j si può fermare a k (compreso) anziché arrivare fino ad n. Alla fine del processo di eliminazione la matrice B contiene L^-1, e avremo ottenuto il sistema triangolare UX=L^-1.

Un'implementazione del calcolo dell'inversa tramite il metodo di Gauss si trova in inversa.c.