c f77 -o jacobi jacobi.f ~/public_html/f/m/lib/Matrice/source/Matrice100.f
program jacobi_lin
implicitnone
real m(100, 100), X(100), P(100), rags
integer righe, exc
c l'iterazione si ferma
c dopo 'ferma' volte che viene trovato lo stesso vettore
integer ferma
ferma = 5
righe = 6
write(6,*) 'Benvenuti in jacobi:'
write(6,*) 'risoluzione di sistemi lineari'
call assumiMAtrice100(m, righe, righe+1)
write(6,*) '\napplicazione di jacobi alla matrice:'
call stampaMatrice100(m, righe, righe+1)
write(6,*) '---------------------'
call jacobi(m, righe, X, ferma, rags, exc)
if(exc .eq. 0) then
c qui e' andato tutto bene:
write(6,*)'\nil vettore dei risultati e\':'
call stampaMatrice100(x, righe ,1)
write(6,*) '\nvettore M*X = '
call prodottoMatrice100(righe,righe,1,m,x,p)
call stampaMatrice100(p,righe,1)
else
if(exc .eq. 1) then
write(6,*) 'jacobi stava per tentare una divisione per'
write(6,*) 'ed e\' stato terminato'
endif
if(exc .eq. 2) then
write(6,*) 'il criterio del raggio spettrale mi dice che'
write(6,*) 'non e\' possibile applicare jacobi alla '
write(6,*) 'alla nostra matrice'
write(6,*) 'il raggio spettrale della matrice di'
1 ,' convergenza e\'', rags
endif
c qui c'e' stato qualche problema:
endif
stop
end
c m = matrice di dimensione lr * lr+1
c exc = parametro per errori
c exc = 1 si stava per tentare una duivisione per zero
c exc = 2 il raggio spettrale in modulo e' <1 => non converge
subroutine jacobi(m, lr, X, ferma, rags, exc)
implicitnone
integer lr, exc, i, j, cont,r, ecc
real m(100, 100),D(100,100),EpF(100,100), conv(100,100), rags, X(100)
real oldX(100), somm1, somm2, ultUg
c dopo che sono state trovate 'ferma' soluzioni uguali si ferma
integer ferma
exc = 0
ultUg = 0
ferma = 4
c riempio D:
do i = 1, lr
D(i, i) = 1/m(i,i)
enddo
write(6,*) '\nmatrice D invertita = '
call stampaMAtrice100(d,lr,lr)
write(6,*) ''
write(6,*)
c riempio E + F
do i = 1,lr
do j = i+1,lr
EpF(i,j) = -M(i,j)
enddo
do j =1, i - 1
EpF(i,j) = -M(i,j)
enddo
EpF(i,i) = 0
enddo
write(6,*) 'matrice E + F = '
call stampaMAtrice100(EpF,lr,lr)
write(6,*)
call prodottoMatrice100(lr,lr,lr,D,EpF,CONV)
write(6,*) '\nla matrice di convergenza (D invertita)*(E+F)
1 e\':'
call stampaMatrice100(conv,lr,lr)
write(6,*)'\nne calcolo il raggio spettrale col metodo delle potenze:'
call raggioSpettrale(conv, lr, rags, ecc)
if(ecc .ne. 0) then
if(ecc .eq. 1) then
write(6,*) 'il metodo delle potenze si e\' arrestato per'
write(6,*) '\tper aver fatto una divisione per zero'
exc = 1
return
endif
return
endif
if (rags .lt. 0) then
rags = -rags
endif
if(rags .ge. 1) then
exc = 2
return
endif
write(6,*) '\nil raggio spettrale della matrice di convergenza e\' ='
write(6,*) '\t', rags, ' < 1, quindi posso applicare'
write(6,*) '\t---Jacobi---'
Call initOldX(m, lr, oldX)
write(6,*) '\noldX iniziale:'
call stampaMatrice100(oldX,lr,1)
c ora si applica jacobi sul serio:
cont = 1
dowhile(cont .eq. 1)
do i = 1,lr
somm1=0.
do j = 1, lr
if(i .ne. j) then
somm1 = somm1 + (m(i,j)*oldx(j))/m(i,i)
endif
enddo
x(i) = m(i, lr+1)/m(i,i) - somm1
enddo
write(6,*) '\nvettore risultati x corrente:'
call stampaMatrice100(X,lr,1)
call uguali(X, oldX, lr, r)
if(r .eq. 1) then
ultug = ultug + 1
endif
if (r .ne. 1) then
ultug = 0
endif
if(ultUg .gt. ferma) then
write(6,*) '\nHo trovato ', ferma, ' volte lo stesso risultato:'
write(6,*) '\tposso fermarmi'
return
endif
do i = 1,lr
oldx(i) = X(i)
ENDDO
enddo
return
end
subroutine uguali(X1,X2,lr,r)
integer lr,r,i
real x1(lr), x2(lr)
do i=1,lr
if(x1(i) .ne. x2(i))then
r = 0
return
endif
enddo
r = 1
return
end
subroutine initOldX(m, lr, oldX)
implicitnone
integer lr, i, j
real m(lr,lr+1),oldX(lr), dmu(lr,lr), K(lr)
c inizializzo D :
do i = 1, lr
do j = 1,lr
dmu(i,j) = 0
enddo
if(i .eq. j) then
if(m (i,j) .ne. 0) then
dmu(i,j) = 1/m(i,j)
else
dmu(i,j) = 100
endif
endif
enddo
do i = 1, lr
K(i) =1
oldX(i) = i
c m(i, lr+1)
enddo
c call prodottoMAtrice100(lr, lr, 1, dmu, K, oldX)
return
end
c m = matrice (INPUT)
c lr = dimensione matrice (INPUT)
c r = raggio spettrale (OUTPUT)
c exc = eccezioni
c 0 = tutto bene
c 1 = si e' tenatata una divisione per zero
subroutine raggioSpettrale(m, lr, r, exc)
implicit none
integer i,lr,exc, ecc
real v(lr), toll,m(lr,lr), ri, scarto,r
toll = 0.01
exc =0
do i=1,lr
v(i) = 0.2 + i
enddo
call quipotenz(M,v,lr,toll,ri,scarto,ecc)
if(exc .ne. 0) then
if(ecc .eq. 1) then
exc = 1
return
endif
return
endif
r = ri
return
end
c calcola col metodo delle potenze il raggio spettrale di modulo massimo
c M = matrice in entrata
c v = vettore in entrata di iterazione
c n = dimensione della matrice
c toll = tolleranza superata la quale si ferma
c ri = raggio spettrale di modulo massimo (OUTPUT)
c scarto = ultimo scarto (OUTPUT)
c exc: parametro per le eccezioni:
c exc = 1 si stava per tentare una divisione per zero
subroutine quipotenz(M,v,n,toll,ri, scarto, exc)
implicit none
integer n,i,j,exc
real M(n,n),v(n), scarto, toll, ri, quinormax, w(n)
real norma, rip1
exc = 0
scarto = 2.0 * toll
ri = quinormax(n,v)
i = 0
dowhile(scarto .gt. toll)
i = i+1
call prodottoMatrice100(n,n,1,M,v,w)
norma = quinormax(n,w)
rip1= norma
if(rip1 .eq. 0) then
exc = 1
return
endif
scarto = abs(1-ri/rip1)
do j=1,n
if(norma .eq. 0) then
exc = 1
return
endif
v(j) = w(j)/norma
enddo
ri = rip1
write(6,*) 'ri = ', ri
enddo
return
end
real function quinormax(n,v)
c calcola la norma massima del vettore v
c parametri di scambio:n,v in input; normax in outpu
implicit none
integer n,i
real v(*)
c inizio calcolo
quinormax = abs(v(1))
doi=2,n
if(quinormax .lt. abs(v(i))) then
quinormax = abs (v(i))
endif
enddo
return
end