integer lda, n, ipiv(n), ifail
real a(lda, n)
Il problema è quello dell'eliminazione di Gauss per una matrice quadrata A, secondo la formula P*A=L*U con P matrice delle permutazioni, L matrice triangolare inferiore, U matrice triangolare superiore.
Descrizione dell'algoritmo
La difficoltà maggiore da affrontare è quella di ridurre l'amplificazione degli errori di round-off, che si propagano con l'applicazione del metodo nel caso generale. La strategia adottata è quella del pivoting parziale, secondo la quale, ad ogni passo si controlla se l'elemento scelto come pivot della k-esima riga è il maggiore di tutti i coefficienti ad esso sottostanti, in caso negativo si scambiano le righe della matrice, in modo tale che nella k-esima riga venga a trovarsi quella con elemento maggiore in modulo, tra tutti quelli ad esso sottostanti. Si procede poi con la fattorizzazione di Gauss sulla "nuova" matrice, equivalente alla precedente. Il pivoting viene simulato, in realtà, non sarà effettuato nessuno scambio di righe, l'applicazione del pivoting alla matrice avviene virtualmente, utilizzando un vettore di interi (ipiv) nel quale si memorizza l'ordine effettivo delle righe. L'algoritmo lavora per colonne.
Raccomandazioni sull'uso
La maggiore raccomandazione sta nel fatto che la matrice da fattorizzare deve essere non singolare, per applicare l'eliminazione gaussiana. Nel caso in cui la matrice da fattorizzare sia malcondizionata, si suggerisce l'uso della subroutine Lups, che fornisce una soluzione più accurata. Per problemi legati alla allocazione in memoria di array bidimensionali, in fortran la dimensione lda (leading dimension) deve essere uguale alla dimensione di riga della matrice dichiarata nel programma chiamante [Appendice D]. Per ridurre la complessità di spazio la subroutine memorizza nel triangolo superiore 1 della matrice di ingresso la matrice triangolare U, in quello inferiore la matrice triangolare L, la matrice di ingresso è, quindi, persa. Come ulteriore output dalla fattorizzazione LU, sfruttando le proprietà matematiche del determinante e il teorema di Binet, si ha:
Si può ottenere, quindi, il determinante della matrice A come risultato indiretto della fattorizzazione LU, moltiplicando i termini della diagonale principale della matrice triangolare U, e comunque con una complessità di O(n3/3), contro una complessità di O(n!) della formula classica. Il vettore ipiv è impiegato per permutare le righe della matrice in uscita, in modo da "ricostruire" la matrice fattorizzata.
a - matrice di reali, matrice da fattorizzare
lda - intero, leading dimension
n - intero, ordine della matrice
output:
a - matrice di reali, matrice fattorizzata
ipiv - vettore di interi, vettore delle permutazioni
ifail - intero, indicatore d'errore
ifail = 3 la dimensione inserita è non positiva
ifail = 4 la matrice inserita è singolare
la complessità di tempo è somma di due contributi separati, quello del pivoting e quello dell'eliminazione. Il primo è pari a O(n2 ),il secondo è pari a O(n3).
Program lupd
c Programma dimostrativo della subroutine lup
external lup
integer dim
parameter (dim=6)
real a(dim,dim)
real*4 temp
integer ipiv(dim),n,i,j,ifail
character*12 fin,fout,foutLU
character sc
print*
print*,'*** PROGRAMMA DIMOSTRATIVO DELLA SUBROUTINE LUP ***'
print*
c inizializzo il generatore di numeri casuali
call seed(RND$TIMESEED)
write(*,10)'Ordine della matrice (max ',dim,')= '
read(*,*)n
print*,'Crea una nuova matrice ' n'
print*,'Carica una matrice da file: c'
read(*,*)sc
if (sc .eq. `c')
print*,'Nome del file di input= '
read(*,'(A)')fin
open(11,FILE=fin,STATUS='OLD')
do 4 i=1,n
do 3 j=1,n
read(11,20)a(i,j)
3 continue
4 continue
endif
goto 100
c creazione di una matrice n*n di valori generati casualmente
do 2 i=1,n
do 1 j=1,n
call random(temp)
a(i,j) = temp * i * j
1 continue
2 continue
print*,'Nome del file di output per la matrice da fattorizzare= '
read(*,'(A)')fout
open(12,FILE=fout)
write(12,20)((a(i,j),j=1,n),i=1,n)
close(12)
100
c richiamo la subroutine lup
call lup(a,dim,n,ipiv,ifail)
if (ifail .eq. 3) stop 'Errore! dimensione non valida'
if (ifail .eq. 4) stop 'Errore! matrice singolare'
print*,'Nome del file di output per la matrice fattorizzata= '
read(*,'(A)')foutLU
open(12,FILE=foutLU)
write(12,20)((a(ipiv(i),j),j=1,n),i=1,n)
write(12,*)
write(12,30)(ipiv(i),i=1,n)
close(12)
print*,'Matrice fattorizzata= '
write(*,20)((a(ipiv(i),j),j=1,n),i=1,n)
10 format (1X,A,I2,A)
20 format (1X,E14.7)
30 format (1X,I3)