Concorso Energia in gioco IT

Calcolo numerico

Programmazione in Fortran 77.

Metodi di soluzione dei sistemi lineari.

In questa rubrica presento una completa serie di algoritmi implementati in Digital Visual Fortran ( totalmente compatibile con il FORTRAN 77 ).

Nel 2005 questi programmi sono serviti agli studenti dei corsi serali per studenti lavoratori per superare l'esame universitario "Calcolo numerico e programmazione" presso la facoltà di ingegneria di Padova ( corso della professoressa Michela Redivo Zaglia per allievi aereospaziali).

I programmi sono testati e funzionanti ma essi non assicurano il superamento della prova d'esame senza un approfondito studio teorico nel testo "Calcolo numerico" dell'autrice Redivo Zaglia.

Attenzione: Nell'anno in cui questi programmi sono stati creati sono stati prodotti in sede d'esame (quindi vanno ben memorizzati nella loro essenza e non solo a memoria) usando il compilatore in linea del sistema operativo Linux installato in aula Taliercio della facoltà di ingegneria in via Belzoni 7 (PD). E' indispensabile al fine di superare la prova d'esame, sapere compilare in ambiente Linux ed effettuare la visualizzazione degli eventuali grafici con i tools di default presenti in quel sistema operativo.

I programmi name.f posso essere tagliati ed incollati nell'editor di programmazione per un corretto funzionamento. Ove sia necessario bisognerà far trovare nel direttorio di compilazione i file di testo contenenti l'insieme di input "formattato" ovvero con righe e colonne al posto in cui l'algoritmo le aspetta. Nella stessa directory troveremo il file testuale di output.

Marco Gottardo.

Metodo di jacobi

jacobi.f

C====================================================================
C Programma per in calcolo della soluzione di un sistema con il
C METODO di Jacobi
C
C Va preparato un file di input che contiene :
C la matrice A +
C sulla riga successiva il vettore termine noto b +
C sulla riga successiva il vettore iniziale x0
C
C====================================================================
PROGRAM jacob
IMPLICIT NONE
C====================
C Parte dichiarativa
C====================
INTEGER n, maxdim, maxiter, k, kmax, i
PARAMETER (maxdim=100, maxiter=100)
REAL*8 A(1:maxdim,1:maxdim), b(1:maxdim), x(1:maxdim),
* iterate(1:maxdim)
REAL*8 rn(0:maxiter)
REAL*8 toll
CHARACTER*20 nomefilein, nomefileout
C====================
C Parte esecutiva
C====================
C
C Ingresso dati
C
print *, ' METODO Jacobi'
print *, ' Inserisci la dimensione n'
read *, n
print *, ' Inserisci la tolleranza'
read *, toll
print *, ' Inserisci il numero massimo di iterazioni'
read *, kmax
C
print *, 'Dammi il nome del file di ingresso'
read*, nomefilein
print *, 'Dammi il nome del file in uscita'
read*, nomefileout
open (unit=9,file=nomefilein)
open (unit=10,file=nomefileout)
print *, ' Lettura della matrice A per righe'
CALL leggimat (A, n, n, maxdim, 9)
print *, ' Lettura del vettore b su di una riga'
CALL leggimat (b, 1, n, 1, 9)
print *, ' Lettura del vettore x0 su di una riga'
CALL leggimat (x, 1, n, 1, 9)
C
C Uscita dati
C
print *
print *, ' METODO Jacobi'
print *
print *, ' Dimensione = ', n
print *, ' Tolleranza = ', toll
print *, ' Numero massimo di iterazioni = ', kmax
print *, ' Matrice A'
CALL scrivimat (A, n, n, maxdim, 6)
print *
print *, ' Vettore b'
CALL scrivimat (b, 1, n, 1, 6)
print *, ' Vettore x0'
CALL scrivimat (x, 1, n, 1, 6)
C
write(10,*)
write(10,*) ' METODO Jacobi'
write(10,*)
write(10,*) ' Dimensione = ', n
write(10,*) ' Tolleranza = ', toll
write(10,*) ' Numero massimo di iterazioni = ', kmax
write(10,*) ' Matrice A'
CALL scrivimat (A, n, n, maxdim, 10)
write(10,*)
write(10,*) ' Vettore b'
CALL scrivimat (b, 1, n, 1, 10)
write(10,*) ' Vettore x0'
CALL scrivimat (x, 1, n, 1, 10)
write(10,*)
C
C Calcolo della soluzione
C
CALL Jacobi (A, b, x, n, maxdim, toll, iterate, kmax, k )
C
C Uscita risultati
C
print *,'********** Risultati ************'
print *
print *, ' Numero di iterazioni eseguite = ', k-1
print *, ' Vettore soluzione x '
CALL scrivimat (x, 1, n, 1, 6)
C
write(10,*)
write(10,*)'********** Risultati ************'
write(10,*)
write(10,*) ' Numero di iterazioni eseguite = ', k-1
write(10,*) ' Vettore soluzione x '
CALL scrivimat (x, 1, n, 1, 10)
write(10,*)
write(10,*)' n. iterata',' valore test-tolleranza'
do i = 0, k-1
write(10,*) i,iterate(i)
end do
close(9)
close(10)
stop
end

C====================================================================
C Subroutine per la risoluzione di un sistema con il
C METODO Jacobi
C====================================================================
SUBROUTINE Jacobi
* (A, b, x, n, maxdim, toll, iterate, kmax, k )
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed n colonne
C b vettore termine noto con n componenti
C n dimensione della matrice A
C maxdim numero massimo di righe di A (come da dimensionamento)
C toll tolleranza desiderata
C kmax numero massimo di iterazioni
C
C Variabili INGRESSO/USCITA:
C x vettore iniziale x0 in ingresso e
C vettore soluzione finale in uscita (n componenti)
C
C Variabili USCITA:
C k numero di iterazioni effettuate
C iterate vettore contenente i valori della norma
C vettore differenza
C
C Altre variabili
C test norma del vettore differenza
C
C====================
C Parte dichiarativa
C====================
INTEGER n, maxdim, kmax, k
INTEGER i, j
REAL*8 A(1:maxdim,*), b(1:*), x(1:*)
REAL*8 xb(1:maxdim), r(1:maxdim),xn(1:maxdim)
REAL*8 toll, test, iterate(1:maxdim)
REAL*8 norm2
C====================
C Parte esecutiva
C====================
k=0
test = toll+1
do while (test.ge.toll.and.k.lt.kmax)
do i = 1, n
xb(i) = 0.0d0
do j = 1, i-1
xb(i) = xb(i) + A(i,j)*xb(j)
end do
do j = i+1, n
xb(i) = xb(i) + A(i,j)*x(j)
end do
xb(i) = (b(i)-xb(i))/A(i,i)
enddo
c test
do i = 1, n
xn(i)=xb(i) - x(i)
enddo
test=norm2(xn,n)
iterate(k)=test
c assegnazione alle componenti di x
c dei valori del vettore appoggio xb
do i = 1, n
x(i) = xb(i)
enddo
k=k+1
enddo
return
end

C====================================================================
C Function per il calcolo della norma euclidea di un vettore
C====================================================================
REAL*8 FUNCTION norm2 (x, n)
IMPLICIT NONE
C
C Variabili INGRESSO:
C x vettore di dimensione n
C n dimensione di x
C
C====================
C Parte dichiarativa
C====================
INTEGER n
INTEGER i
REAL*8 x(1:*)
C====================
C Parte esecutiva
C====================
norm2 = 0.0d0
do i = 1, n
norm2 = norm2 + x(i)**2
end do
norm2 = dsqrt(norm2)
return
end

C====================================================================
C Subroutine per il calcolo del prodotto matrice - vettore
C====================================================================
SUBROUTINE matxvet (A, x, n, m, maxdim, y)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed m colonne
C x vettore m righe
C n numero di righe di A e numero di righe di y
C m numero di colonne di A e numero di righe di x
C maxdim numero massimo di righe di A (come da dimensionamento)
C
C Variabili USCITA:
C y vettore risultante di n righe. y=A*x
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim
INTEGER i, j
REAL*8 A(1:maxdim,*), x(1:*), y(1:*)
C====================
C Parte esecutiva
C====================
do i = 1, n
y(i) = 0.0d0
do j = 1, m
y(i) = y(i) + A(i,j) * x(j)
end do
end do
return
end

C====================================================================
C Subroutine per la LETTURA di una matrice a n righe ed m colonne
C====================================================================
SUBROUTINE leggimat (A, n, m, maxdim, io)
IMPLICIT NONE
C
C Variabili INGRESSO:
C n numero di righe di A
C m numero di colonne di A
C maxdim numero massimo di righe di A (come da dimensionamento)
C io numero dell'unita' da cui deve leggere la matrice
C
C Variabili USCITA:
C A matrice n righe ed m colonne
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim, io
INTEGER i, j
REAL*8 A(1:maxdim,*)
C====================
C Parte esecutiva
C====================
do i = 1, n
read (io,*) (A(i,j), j=1,m)
end do
return
end

C====================================================================
C Subroutine per la SCRITTURA di una matrice a n righe ed m colonne
C====================================================================
SUBROUTINE scrivimat (A, n, m, maxdim, io)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed m colonne
C n numero di righe di A
C m numero di colonne di A
C maxdim numero massimo di righe di A (come da dimensionamento)
C io numero dell'unita' da cui deve scrivere la matrice
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim, io
INTEGER i, j
REAL*8 A(1:maxdim,*)
C====================
C Parte esecutiva
C====================
do i = 1, n
write (io,*) (A(i,j), j=1,m)
end do
return
end

Programma per in calcolo della soluzione di un sistema con il metodo di GAUSS CON pivoting

Versione con lettura della matrice AUMENTATA (A|b).

gaussP.f

C====================================================================
C Programma per in calcolo della soluzione di un sistema
C con il metodo di GAUSS CON pivoting
C Versione con lettura della matrice AUMENTATA (A|b)
C====================================================================
C
C Ab matrice n righe ed (n+1) colonne composta dalla matrice
C A (n righe ed n colonne) e dal vettore termine noto b
C (ultima colonna di A)
C x vettore soluzione
C
PROGRAM gaussnp
IMPLICIT NONE
C====================
C Parte dichiarativa
C====================
INTEGER n, maxdim
INTEGER i
PARAMETER (maxdim=100)
REAL*8 Ab(1:maxdim,1:maxdim+1), L(1:maxdim,1:maxdim)
REAL*8 LU(1:maxdim,1:maxdim)
REAL*8 y(1:maxdim), x(1:maxdim)
CHARACTER*20 nomefilein, nomefileout
C====================
C Parte esecutiva
C====================
C
C Ingresso dati
C
C
print *, 'Dammi il nome del file di ingresso'
read*, nomefilein
print *, 'Dammi il nome del file in uscita'
read*, nomefileout
open (unit=9,file=nomefilein)
open (unit=10,file=nomefileout)
print *, ' Inserisci la dimensione n'
read *, n
CALL leggimat (Ab, n, n+1, maxdim, 9)
C
C Uscita dati
C
write(10,*)
write(10,*) ' Matrice aumentata Ab'
CALL scrivimat (Ab, n, n+1, maxdim, 10)
write(10,*)
C
C Fattorizzazione di Gauss
C
CALL gaussp (Ab, n, maxdim, L)
C
C Costruzione del vettore y
C
do i = 1, n
y(i) = Ab(i,n+1)
end do
C
C Calcolo della soluzione
C
CALL trisup (Ab, y, n, maxdim, x)
C
C Calcolo del prodotto LU
C
CALL matxmat (L, Ab, n, n, n, maxdim, LU)
C
C Uscita risultati
C
write(10,*)
write(10,*) ' Matrice L'
CALL scrivimat (L, n, n, maxdim, 10)
write(10,*)
write(10,*) ' Matrice U'
CALL scrivimat (Ab, n, n, maxdim, 10)
write(10,*)
write(10,*) ' Matrice L*U'
CALL scrivimat (LU, n, n, maxdim, 10)
write(10,*)
write(10,*) ' Vettore y'
CALL scrivimat (y, 1, n, 1, 10)
write(10,*)
write(10,*) ' Vettore soluzione x '
CALL scrivimat (x, 1, n, 1, 10)
stop
end

C====================================================================
C Subroutine per la LETTURA di una matrice a n righe ed m colonne
C====================================================================
SUBROUTINE leggimat (A, n, m, maxdim, io)
IMPLICIT NONE
C
C Variabili INGRESSO:
C n numero di righe di A
C m numero di colonne di A
C maxdim numero massimo di righe di A (come da dimensionamento)
C io numero dell'unita' da cui deve leggere la matrice
C
C Variabili USCITA:
C A matrice n righe ed m colonne
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim, io
INTEGER i, j
REAL*8 A(1:maxdim,*)
C====================
C Parte esecutiva
C====================
do i = 1, n
read (io,*) (A(i,j), j=1,m)
end do
return
end

C====================================================================
C Subroutine per la SCRITTURA di una matrice a n righe ed m colonne
C====================================================================
SUBROUTINE scrivimat (A, n, m, maxdim, io)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed m colonne
C n numero di righe di A
C m numero di colonne di A
C maxdim numero massimo di righe di A (come da dimensionamento)
C io numero dell'unita' da cui deve scrivere la matrice
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim, io
INTEGER i, j
REAL*8 A(1:maxdim,*)
C====================
C Parte esecutiva
C====================
do i = 1, n
write (io,*) (A(i,j), j=1,m)
end do
return
end

C====================================================================
C Subroutine per la risoluzione di un sistema triangolare superiore
C====================================================================
SUBROUTINE trisup (U, b, n, maxdim, x)
IMPLICIT NONE
C
C Variabili INGRESSO:
C U matrice n righe ed n colonne, triangolare superiore
C b vettore termine noto con n componenti
C n dimensione della matrice U
C maxdim numero massimo di righe di U (come da dimensionamento)
C
C Variabili USCITA:
C x vettore soluzione con n componenti
C
C====================
C Parte dichiarativa
C====================
INTEGER n, maxdim
INTEGER i, j
REAL*8 U(1:maxdim,*), b(1:*), x(1:*)
C====================
C Parte esecutiva
C====================
do i = n, 1, -1
x(i) = b(i)
do j = i+1, n
x(i) = x(i) - U(i,j) * x(j)
end do
x(i) = x(i)/U(i,i)
end do
return
end

C====================================================================
C Subroutine per il calcolo del prodotto matrice - matrice
C====================================================================
SUBROUTINE matxmat (A, B, n, m, r, maxdim, C)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed m colonne
C B matrice m righe ed r colonne
C n numero di righe di A e numero di righe di C
C m numero di colonne di A e numero di righe di B
C r numero di colonne di B e numero di colonne di C
C maxdim numero massimo di righe di A, B, C (come da dimensionamento)
C
C Variabili USCITA:
C C matrice prodotto n righe ed r colonne
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, r, maxdim
INTEGER i, j, k
REAL*8 A(1:maxdim,*), B(1:maxdim,*), C(1:maxdim,*)
REAL*8 acc
C====================
C Parte esecutiva
C====================
do i = 1, n
do j = 1, r
acc = 0.0d0
do k = 1, m
acc = acc + A(i,k)*B(k,j)
end do
C(i,j) = acc
end do
end do
return
end

C====================================================================
C Subroutine per la risoluzione di un sistema lineare con il
C Metodo di Gauss CON pivoting
C====================================================================
SUBROUTINE gaussp (Ab, n, maxdim, L)
IMPLICIT NONE
C
C Variabili INGRESSO:
C n dimensione della matrice U
C maxdim numero massimo di righe di U (come da dimensionamento)
C
C Variabili INGRESSO/USCITA:
C Ab matrice n righe ed n+1 colonne
C in ingresso deve contenere la matrice aumentata
C in uscita assume la struttura (U|y)
C
C Variabili USCITA:
C L matrice n righe ed n colonne

C
C====================
C Parte dichiarativa
C====================
INTEGER n, maxdim
INTEGER k, i, r, j
REAL*8 Ab(1:maxdim,1:*), L(1:maxdim,1:*)
REAL*8 tmp
C====================
C Parte esecutiva
C====================
do k = 1, n-1
tmp = dabs(Ab(k,k))
r = k
do i = k+1, n
if (dabs(Ab(i,k)).gt.tmp) then
tmp = dabs(Ab(i,k))
r = i
print*, 'Scambio riga',r,' con riga', k
end if
end do
if (r.ne.k) then
do i = 1, n+1
tmp = Ab(k,i)
Ab(k,i) = Ab(r,i)
Ab(r,i) = tmp
end do
end if
do i = k+1, n
Ab(i,k) = Ab(i,k)/Ab(k,k)
do j = k+1, n+1
Ab(i,j) = Ab(i,j) - Ab(i,k)*Ab(k,j)
end do
end do
end do
do i = 1, n
L(i,i) =1.0d0
do j = i+1, n
L(i,j) = 0.0d0
L(j,i) = Ab(j,i)
Ab(j,i) = 0.0d0
end do
end do
return
end

Programma per in calcolo della soluzione di un sistema con il METODO GAUSS-SEIDEL.

Gauss-Seidel.f

C====================================================================
C Programma per in calcolo della soluzione di un sistema con il
C METODO GAUSS-SEIDEL
C====================================================================
PROGRAM gaussseidel
IMPLICIT NONE
C====================
C Parte dichiarativa
C====================
INTEGER n, maxdim, maxiter, k, kmax, i
PARAMETER (maxdim=100, maxiter=100)
REAL*8 A(1:maxdim,1:maxdim), b(1:maxdim), x(1:maxdim)
REAL*8 rn(0:maxiter)
REAL*8 toll
CHARACTER*20 nomefilein, nomefileout
C====================
C Parte esecutiva
C====================
C
C Ingresso dati
C
print *, ' METODO GAUSS-SEIDEL'
print *, ' Inserisci la dimensione n'
read *, n
print *, ' Inserisci la tolleranza'
read *, toll
print *, ' Inserisci il numero massimo di iterazioni'
read *, kmax
C
print *, 'Dammi il nome del file di ingresso'
read*, nomefilein
print *, 'Dammi il nome del file in uscita'
read*, nomefileout
open (unit=9,file=nomefilein)
open (unit=10,file=nomefileout)
print *, ' Lettura della matrice A per righe'
CALL leggimat (A, n, n, maxdim, 9)
print *, ' Lettura del vettore b su di una riga'
CALL leggimat (b, 1, n, 1, 9)
print *, ' Lettura del vettore x0 su di una riga'
CALL leggimat (x, 1, n, 1, 9)
C
C Uscita dati
C
print *
print *, ' METODO GAUSS-SEIDEL'
print *
print *, ' Dimensione = ', n
print *, ' Tolleranza = ', toll
print *, ' Numero massimo di iterazioni = ', kmax
print *, ' Matrice A'
CALL scrivimat (A, n, n, maxdim, 6)
print *
print *, ' Vettore b'
CALL scrivimat (b, 1, n, 1, 6)
print *, ' Vettore x0'
CALL scrivimat (x, 1, n, 1, 6)
C
write(10,*)
write(10,*) ' METODO GAUSS-SEIDEL'
write(10,*)
write(10,*) ' Dimensione = ', n
write(10,*) ' Tolleranza = ', toll
write(10,*) ' Numero massimo di iterazioni = ', kmax
write(10,*) ' Matrice A'
CALL scrivimat (A, n, n, maxdim, 10)
write(10,*)
write(10,*) ' Vettore b'
CALL scrivimat (b, 1, n, 1, 10)
write(10,*) ' Vettore x0'
CALL scrivimat (x, 1, n, 1, 10)
write(10,*)
C
C Calcolo della soluzione
C
CALL GSEIDEL (A, b, x, n, maxdim, toll, kmax, k, rn)
C
C Uscita risultati
C
print *
print *, ' Tolleranza * norm2(b) = ', toll
print *, ' Numero di iterazioni eseguite = ', k
print *, ' Vettore soluzione x '
CALL scrivimat (x, 1, n, 1, 6)
print *, ' Norma residuo finale = ', rn(k)
C
write(10,*)
write(10,*) ' Tolleranza * norm2(b) = ', toll
write(10,*) ' Numero di iterazioni eseguite = ', k
write(10,*) ' Vettore soluzione x '
CALL scrivimat (x, 1, n, 1, 10)
write(10,*) ' Norma residuo finale = ', rn(k)
write(10,*)
do i = 0, k
write(10,'(i5,f7.2,e25.16)') i, log10(rn(i)), rn(i)
end do
close(9)
close(10)
stop
end

C====================================================================
C Subroutine per la risoluzione di un sistema con il
C METODO GAUSS-SEIDEL
C====================================================================
SUBROUTINE GSEIDEL
* (A, b, x, n, maxdim, toll, kmax, k, resnorm)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed n colonne
C b vettore termine noto con n componenti
C n dimensione della matrice A
C maxdim numero massimo di righe di A (come da dimensionamento)
C toll tolleranza desiderata
C kmax numero massimo di iterazioni
C
C Variabili INGRESSO/USCITA:
C x vettore iniziale x0 in ingresso e
C vettore soluzione finale in uscita (n componenti)
C
C Variabili USCITA:
C k numero di iterazioni effettuate
C resnorm vettore norma del residuo - dimensione (0:k)
C
C====================
C Parte dichiarativa
C====================
INTEGER n, maxdim, kmax, k
INTEGER i, j
REAL*8 A(1:maxdim,*), b(1:*), x(1:*)
REAL*8 xb(1:maxdim), r(1:maxdim)
REAL*8 resnorm(0:*)
REAL*8 toll
REAL*8 test, norm2
C====================
C Parte esecutiva
C====================
k=0
CALL residuo (A, b, x, n, maxdim, r)
test = norm2(r,n)
toll = toll * norm2(b,n)
do while (test.ge.toll.and.k.lt.kmax)
do i = 1, n
xb(i) = 0.0d0
do j = 1, i-1
xb(i) = xb(i) + A(i,j)*xb(j)
end do
do j = i+1, n
xb(i) = xb(i) + A(i,j)*x(j)
end do
xb(i) = (b(i)-xb(i))/A(i,i)
end do
do i = 1, n
x(i) = xb(i)
end do
resnorm(k) = test
CALL residuo (A, b, x, n, maxdim, r)
test = norm2(r,n)
k=k+1
enddo
resnorm(k) = test
return
end

C====================================================================
C Subroutine per il calcolo del vettore residuo
C====================================================================
SUBROUTINE residuo (A, b, x, n, maxdim, r)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice di dimensione n
C b vettore termine noto di dimensione n
C x vettore iterata di dimensione n
C n dimensione di A, di b, di x e di r
C maxdim numero massimo di righe di A (come da dimensionamento)
C
C Variabili USCITA:
C r vettore residuo di dimensione n
C
C====================
C Parte dichiarativa
C====================
INTEGER n, maxdim
INTEGER i
REAL*8 A(1:maxdim,*), b(1:*), x(1:*), r(1:*)
C====================
C Parte esecutiva
C====================
C
C Calcolo di A*x
C
CALL matxvet (A, x, n, n, maxdim, r)
C
C Calcolo del vettore residuo
C
do i = 1, n
r(i) = b(i) - r(i)
end do
return
end

C====================================================================
C Function per il calcolo della norma euclidea di un vettore
C====================================================================
REAL*8 FUNCTION norm2 (x, n)
IMPLICIT NONE
C
C Variabili INGRESSO:
C x vettore di dimensione n
C n dimensione di x
C
C====================
C Parte dichiarativa
C====================
INTEGER n
INTEGER i
REAL*8 x(1:*)
C====================
C Parte esecutiva
C====================
norm2 = 0.0d0
do i = 1, n
norm2 = norm2 + x(i)**2
end do
norm2 = dsqrt(norm2)
return
end

C====================================================================
C Subroutine per il calcolo del prodotto matrice - vettore
C====================================================================
SUBROUTINE matxvet (A, x, n, m, maxdim, y)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed m colonne
C x vettore m righe
C n numero di righe di A e numero di righe di y
C m numero di colonne di A e numero di righe di x
C maxdim numero massimo di righe di A (come da dimensionamento)
C
C Variabili USCITA:
C y vettore risultante di n righe. y=A*x
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim
INTEGER i, j
REAL*8 A(1:maxdim,*), x(1:*), y(1:*)
C====================
C Parte esecutiva
C====================
do i = 1, n
y(i) = 0.0d0
do j = 1, m
y(i) = y(i) + A(i,j) * x(j)
end do
end do
return
end

C====================================================================
C Subroutine per la LETTURA di una matrice a n righe ed m colonne
C====================================================================
SUBROUTINE leggimat (A, n, m, maxdim, io)
IMPLICIT NONE
C
C Variabili INGRESSO:
C n numero di righe di A
C m numero di colonne di A
C maxdim numero massimo di righe di A (come da dimensionamento)
C io numero dell'unita' da cui deve leggere la matrice
C
C Variabili USCITA:
C A matrice n righe ed m colonne
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim, io
INTEGER i, j
REAL*8 A(1:maxdim,*)
C====================
C Parte esecutiva
C====================
do i = 1, n
read (io,*) (A(i,j), j=1,m)
end do
return
end

C====================================================================
C Subroutine per la SCRITTURA di una matrice a n righe ed m colonne
C====================================================================
SUBROUTINE scrivimat (A, n, m, maxdim, io)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed m colonne
C n numero di righe di A
C m numero di colonne di A
C maxdim numero massimo di righe di A (come da dimensionamento)
C io numero dell'unita' da cui deve scrivere la matrice
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim, io
INTEGER i, j
REAL*8 A(1:maxdim,*)
C====================
C Parte esecutiva
C====================
do i = 1, n
write (io,*) (A(i,j), j=1,m)
end do
return
end

Programma per in calcolo della soluzione di un sistema con il METODO JACOBI.

jacobi.f

C====================================================================
C Programma per in calcolo della soluzione di un sistema con il
C METODO JACOBI
C====================================================================
PROGRAM jacobi
IMPLICIT NONE
C====================
C Parte dichiarativa
C====================
INTEGER n, maxdim, maxiter, k, kmax, i
PARAMETER (maxdim=100, maxiter=100)
REAL*8 A(1:maxdim,1:maxdim), b(1:maxdim), x(1:maxdim)
REAL*8 toll
CHARACTER*20 nomefilein, nomefileout
C====================
C Parte esecutiva
C====================
C
C Ingresso dati
C
print *, ' METODO JACOBI'
print *, ' Inserisci la dimensione n'
read *, n
print *, ' Inserisci la tolleranza'
read *, toll
print *, ' Inserisci il numero massimo di iterazioni'
read *, kmax
C
print *, 'Dammi il nome del file di ingresso'
read*, nomefilein
print *, 'Dammi il nome del file in uscita'
read*, nomefileout
open (unit=9,file=nomefilein)
open (unit=10,file=nomefileout)
print *, ' Lettura della matrice A per righe'
CALL leggimat (A, n, n, maxdim, 9)
print *, ' Lettura del vettore b su di una riga'
CALL leggimat (b, 1, n, 1, 9)
print *, ' Lettura del vettore x0 su di una riga'
CALL leggimat (x, 1, n, 1, 9)
C
C Uscita dati
C
print *
print *, ' METODO JACOBI'
print *
print *, ' Dimensione = ', n
print *, ' Tolleranza = ', toll
print *, ' Numero massimo di iterazioni = ', kmax
print *, ' Matrice A'
CALL scrivimat (A, n, n, maxdim, 6)
print *
print *, ' Vettore b'
CALL scrivimat (b, 1, n, 1, 6)
print *, ' Vettore x0'
CALL scrivimat (x, 1, n, 1, 6)
C
write(10,*)
write(10,*) ' METODO JACOBI'
write(10,*)
write(10,*) ' Dimensione = ', n
write(10,*) ' Tolleranza = ', toll
write(10,*) ' Numero massimo di iterazioni = ', kmax
write(10,*) ' Matrice A'
CALL scrivimat (A, n, n, maxdim, 10)
write(10,*)
write(10,*) ' Vettore b'
CALL scrivimat (b, 1, n, 1, 10)
write(10,*) ' Vettore x0'
CALL scrivimat (x, 1, n, 1, 10)
write(10,*)
C
C Calcolo della soluzione
C
CALL IACOBI (A, b, x, n, maxdim, toll, kmax, k)
C
C Uscita risultati
C
print *
print *, ' Tolleranza = ', toll
print *, ' Numero di iterazioni eseguite = ', k
print *, ' Vettore soluzione x '
CALL scrivimat (x, 1, n, 1, 6)
C
write(10,*)
write(10,*) ' Tolleranza = ', toll
write(10,*) ' Numero di iterazioni eseguite = ', k
write(10,*) ' Vettore soluzione x '
CALL scrivimat (x, 1, n, 1, 10)
write(10,*)
close(9)
close(10)
stop
end

C====================================================================
C Subroutine per la risoluzione di un sistema con il
C METODO JACOBI
C====================================================================
SUBROUTINE IACOBI
* (A, b, x, n, maxdim, toll, kmax, k)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed n colonne
C b vettore termine noto con n componenti
C n dimensione della matrice A
C maxdim numero massimo di righe di A (come da dimensionamento)
C toll tolleranza desiderata
C kmax numero massimo di iterazioni
C
C Variabili INGRESSO/USCITA:
C x vettore iniziale x0 in ingresso e
C vettore soluzione finale in uscita (n componenti)
C
C Variabili USCITA:
C k numero di iterazioni effettuate
C
C====================
C Parte dichiarativa
C====================
INTEGER n, maxdim, kmax, k
INTEGER i, j
REAL*8 A(1:maxdim,*), b(1:*), x(1:*)
REAL*8 xb(1:maxdim), r(1:maxdim)
REAL*8 toll
REAL*8 test, norm2
C====================
C Parte esecutiva
C====================
k=0
test = toll + 1.0d0
do while (test.ge.toll.and.k.lt.kmax)
do i = 1, n
xb(i) = 0.0d0
do j = 1, i-1
xb(i) = xb(i) + A(i,j)*xb(j)
end do
do j = i+1, n
xb(i) = xb(i) + A(i,j)*x(j)
end do
xb(i) = (b(i)-xb(i))/A(i,i)
end do
do i = 1, n
test = abs(xb(i)-x(i))
enddo
do i = 1, n
x(i) = xb(i)
end do
k=k+1
enddo
return
end

C====================================================================
C Subroutine per il calcolo del prodotto matrice - vettore
C====================================================================
SUBROUTINE matxvet (A, x, n, m, maxdim, y)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed m colonne
C x vettore m righe
C n numero di righe di A e numero di righe di y
C m numero di colonne di A e numero di righe di x
C maxdim numero massimo di righe di A (come da dimensionamento)
C
C Variabili USCITA:
C y vettore risultante di n righe. y=A*x
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim
INTEGER i, j
REAL*8 A(1:maxdim,*), x(1:*), y(1:*)
C====================
C Parte esecutiva
C====================
do i = 1, n
y(i) = 0.0d0
do j = 1, m
y(i) = y(i) + A(i,j) * x(j)
end do
end do
return
end

C====================================================================
C Subroutine per la LETTURA di una matrice a n righe ed m colonne
C====================================================================
SUBROUTINE leggimat (A, n, m, maxdim, io)
IMPLICIT NONE
C
C Variabili INGRESSO:
C n numero di righe di A
C m numero di colonne di A
C maxdim numero massimo di righe di A (come da dimensionamento)
C io numero dell'unita' da cui deve leggere la matrice
C
C Variabili USCITA:
C A matrice n righe ed m colonne
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim, io
INTEGER i, j
REAL*8 A(1:maxdim,*)
C====================
C Parte esecutiva
C====================
do i = 1, n
read (io,*) (A(i,j), j=1,m)
end do
return
end

C====================================================================
C Subroutine per la SCRITTURA di una matrice a n righe ed m colonne
C====================================================================
SUBROUTINE scrivimat (A, n, m, maxdim, io)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed m colonne
C n numero di righe di A
C m numero di colonne di A
C maxdim numero massimo di righe di A (come da dimensionamento)
C io numero dell'unita' da cui deve scrivere la matrice
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim, io
INTEGER i, j
REAL*8 A(1:maxdim,*)
C====================
C Parte esecutiva
C====================
do i = 1, n
write (io,*) (A(i,j), j=1,m)
end do
return
end

Programma per in calcolo della soluzione di un sistema con il METODO SOR

sor.f

C====================================================================
C Programma per in calcolo della soluzione di un sistema con il
C METODO SOR
C====================================================================
PROGRAM mainsor
IMPLICIT NONE
C====================
C Parte dichiarativa
C====================
INTEGER n, maxdim, maxiter, k, kmax, i
PARAMETER (maxdim=100, maxiter=100)
REAL*8 A(1:maxdim,1:maxdim), b(1:maxdim), x(1:maxdim)
REAL*8 rn(0:maxiter)
REAL*8 toll, omega
CHARACTER*20 nomefilein, nomefileout
C====================
C Parte esecutiva
C====================
C
C Ingresso dati
C
print *, ' METODO SOR'
print *, ' Inserisci la dimensione n'
read *, n
print *, ' Inserisci la tolleranza'
read *, toll
print *, ' Inserisci il numero massimo di iterazioni'
read *, kmax
print *, ' Inserisci il parametro omega'
read *, omega
C
print *, 'Dammi il nome del file di ingresso'
read*, nomefilein
print *, 'Dammi il nome del file in uscita'
read*, nomefileout
open (unit=9,file=nomefilein)
open (unit=10,file=nomefileout)
print *, ' Lettura della matrice A per righe'
CALL leggimat (A, n, n, maxdim, 9)
print *, ' Lettura del vettore b su di una riga'
CALL leggimat (b, 1, n, 1, 9)
print *, ' Lettura del vettore x0 su di una riga'
CALL leggimat (x, 1, n, 1, 9)
C
C Uscita dati
C
print *
print *, ' METODO SOR'
print *
print *, ' Dimensione = ', n
print *, ' Tolleranza = ', toll
print *, ' Numero massimo di iterazioni = ', kmax
print *, ' Omega = ', omega
print *, ' Matrice A'
CALL scrivimat (A, n, n, maxdim, 6)
print *
print *, ' Vettore b'
CALL scrivimat (b, 1, n, 1, 6)
print *, ' Vettore x0'
CALL scrivimat (x, 1, n, 1, 6)
C
write(10,*)
write(10,*) ' METODO SOR'
write(10,*)
write(10,*) ' Dimensione = ', n
write(10,*) ' Tolleranza = ', toll
write(10,*) ' Numero massimo di iterazioni = ', kmax
write(10,*) ' Omega = ', omega
write(10,*) ' Matrice A'
CALL scrivimat (A, n, n, maxdim, 10)
write(10,*)
write(10,*) ' Vettore b'
CALL scrivimat (b, 1, n, 1, 10)
write(10,*) ' Vettore x0'
CALL scrivimat (x, 1, n, 1, 10)
write(10,*)
C
C Calcolo della soluzione
C
CALL SOR (A, b, x, n, maxdim, toll, kmax, omega, k, rn)
C
C Uscita risultati
C
print *
print *, ' Tolleranza * norm2(b) = ', toll
print *, ' Numero di iterazioni eseguite = ', k
print *, ' Vettore soluzione x '
CALL scrivimat (x, 1, n, 1, 6)
print *, ' Norma residuo finale = ', rn(k)
C
write(10,*)
write(10,*) ' Tolleranza * norm2(b) = ', toll
write(10,*) ' Numero di iterazioni eseguite = ', k
write(10,*) ' Vettore soluzione x '
CALL scrivimat (x, 1, n, 1, 10)
write(10,*) ' Norma residuo finale = ', rn(k)
write(10,*)
do i = 0, k
write(10,'(i5,f7.2,e25.16)') i, log10(rn(i)), rn(i)
end do
close(9)
close(10)
stop
end

C====================================================================
C Subroutine per la risoluzione di un sistema con il
C METODO SOR
C====================================================================
SUBROUTINE SOR
* (A, b, x, n, maxdim, toll, kmax, omega, k, resnorm)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed n colonne
C b vettore termine noto con n componenti
C n dimensione della matrice A
C maxdim numero massimo di righe di A (come da dimensionamento)
C toll tolleranza desiderata
C kmax numero massimo di iterazioni
C omega parametro di sovrarilassamento
C
C Variabili INGRESSO/USCITA:
C x vettore iniziale x0 in ingresso e
C vettore soluzione finale in uscita (n componenti)
C
C Variabili USCITA:
C k numero di iterazioni effettuate
C resnorm vettore norma del residuo - dimensione (0:k)
C
C====================
C Parte dichiarativa
C====================
INTEGER n, maxdim, kmax, k
INTEGER i, j
REAL*8 A(1:maxdim,*), b(1:*), x(1:*)
REAL*8 xb(1:maxdim), r(1:maxdim)
REAL*8 resnorm(0:*)
REAL*8 toll, omega
REAL*8 test, norm2
C====================
C Parte esecutiva
C====================
k=0
CALL residuo (A, b, x, n, maxdim, r)
test = norm2(r,n)
toll = toll * norm2(b,n)
do while (test.ge.toll.and.k.lt.kmax)
do i = 1, n
xb(i) = 0.0d0
do j = 1, i-1
xb(i) = xb(i) + A(i,j)*xb(j)
end do
do j = i+1, n
xb(i) = xb(i) + A(i,j)*x(j)
end do
xb(i) = (b(i)-xb(i))/A(i,i)
xb(i) = x(i) + omega*(xb(i)-x(i))
end do
do i = 1, n
x(i) = xb(i)
end do
resnorm(k) = test
CALL residuo (A, b, x, n, maxdim, r)
test = norm2(r,n)
k=k+1
enddo
resnorm(k) = test
return
end

C====================================================================
C Subroutine per il calcolo del vettore residuo
C====================================================================
SUBROUTINE residuo (A, b, x, n, maxdim, r)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice di dimensione n
C b vettore termine noto di dimensione n
C x vettore iterata di dimensione n
C n dimensione di A, di b, di x e di r
C maxdim numero massimo di righe di A (come da dimensionamento)
C
C Variabili USCITA:
C r vettore residuo di dimensione n
C
C====================
C Parte dichiarativa
C====================
INTEGER n, maxdim
INTEGER i
REAL*8 A(1:maxdim,*), b(1:*), x(1:*), r(1:*)
C====================
C Parte esecutiva
C====================
C
C Calcolo di A*x
C
CALL matxvet (A, x, n, n, maxdim, r)
C
C Calcolo del vettore residuo
C
do i = 1, n
r(i) = b(i) - r(i)
end do
return
end

C====================================================================
C Function per il calcolo della norma euclidea di un vettore
C====================================================================
REAL*8 FUNCTION norm2 (x, n)
IMPLICIT NONE
C
C Variabili INGRESSO:
C x vettore di dimensione n
C n dimensione di x
C
C====================
C Parte dichiarativa
C====================
INTEGER n
INTEGER i
REAL*8 x(1:*)
C====================
C Parte esecutiva
C====================
norm2 = 0.0d0
do i = 1, n
norm2 = norm2 + x(i)**2
end do
norm2 = dsqrt(norm2)
return
end

C====================================================================
C Subroutine per il calcolo del prodotto matrice - vettore
C====================================================================
SUBROUTINE matxvet (A, x, n, m, maxdim, y)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed m colonne
C x vettore m righe
C n numero di righe di A e numero di righe di y
C m numero di colonne di A e numero di righe di x
C maxdim numero massimo di righe di A (come da dimensionamento)
C
C Variabili USCITA:
C y vettore risultante di n righe. y=A*x
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim
INTEGER i, j
REAL*8 A(1:maxdim,*), x(1:*), y(1:*)
C====================
C Parte esecutiva
C====================
do i = 1, n
y(i) = 0.0d0
do j = 1, m
y(i) = y(i) + A(i,j) * x(j)
end do
end do
return
end

C====================================================================
C Subroutine per la LETTURA di una matrice a n righe ed m colonne
C====================================================================
SUBROUTINE leggimat (A, n, m, maxdim, io)
IMPLICIT NONE
C
C Variabili INGRESSO:
C n numero di righe di A
C m numero di colonne di A
C maxdim numero massimo di righe di A (come da dimensionamento)
C io numero dell'unita' da cui deve leggere la matrice
C
C Variabili USCITA:
C A matrice n righe ed m colonne
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim, io
INTEGER i, j
REAL*8 A(1:maxdim,*)
C====================
C Parte esecutiva
C====================
do i = 1, n
read (io,*) (A(i,j), j=1,m)
end do
return
end

C====================================================================
C Subroutine per la SCRITTURA di una matrice a n righe ed m colonne
C====================================================================
SUBROUTINE scrivimat (A, n, m, maxdim, io)
IMPLICIT NONE
C
C Variabili INGRESSO:
C A matrice n righe ed m colonne
C n numero di righe di A
C m numero di colonne di A
C maxdim numero massimo di righe di A (come da dimensionamento)
C io numero dell'unita' da cui deve scrivere la matrice
C
C====================
C Parte dichiarativa
C====================
INTEGER n, m, maxdim, io
INTEGER i, j
REAL*8 A(1:maxdim,*)
C====================
C Parte esecutiva
C====================
do i = 1, n
write (io,*) (A(i,j), j=1,m)
end do
return
end