====== Routine ausiliarie ======
===== Matrici e vettori =====
==== Calcolo di determinante per matrici 2x2 o 3x3 ====
c234567---------------------------------------------------------------72
c===================== calcolo di determinante per matrici 2x2 o 3x3 ===
c calcolo il determinante di una matrice di ordine 2 o 3
c a : matrice della quale calcolare il determinante
c n : ordine della matrice
c det : determinante calcolato
c restituisce 'not a number' se l'ordine è diverso da 2 o 3.
c
subroutine detsm(a,n,det)
implicit none
integer n
double precision a(n,n),det
if (n.eq.1) then
det=a(1,1)
else if (n.eq.2) then
det=a(1,1)*a(2,2)-a(1,2)*a(2,1)
else if (n.eq.3) then
det= a(1,1)*a(2,2)*a(3,3)
$ +a(1,2)*a(2,3)*a(3,1)
$ +a(1,3)*a(2,1)*a(3,2)
$ -a(1,3)*a(2,2)*a(3,1)
$ -a(1,2)*a(2,1)*a(3,3)
$ -a(1,1)*a(2,3)*a(3,2)
else
write(0,*) 'ERR: detsm called with invalid n=',n
det=0.d0
det=det/det
end if
return
end subroutine
==== Calcolo dell'inversa di una matrice 2x2 o 3x3 ====
c234567---------------------------------------------------------------72
c============================== calcolo l'inversa di una matrice 2x2 ===
c calcolo la matrice inversa di una matrice di ordine 2 o 3
c a : matrice della quale calcolare il determinante
c n : ordine della matrice
c inva : matrice inversa calcolata
c restituisce una matrice 'not a number' se l'ordine è diverso da 2.
c
subroutine invsm(a,n,inva)
implicit none
integer n
double precision a(n,n),inva(n,n),det
c calcolo il determinante
call detsm(a,n,det)
if (n.eq.1) then
inva(1,1)= 1.d0/a(1,1)
else if (n.eq.2) then
inva(1,1)= a(2,2)/det
inva(1,2)=-a(1,2)/det
inva(2,1)=-a(2,1)/det
inva(2,2)= a(1,1)/det
c
c TODO: inversa matrice 3x3!
c
else
write(0,*) 'ERR: invsm called with invalid n=',n
det=0.d0
det=det/det
end if
return
end subroutine
==== Assemblaggio matrice su matrice ====
c234567---------------------------------------------------------------72
c======================== mappa, scala e assembla matrice su matrice ===
c
c INPUT:
c . sca : fattore scalare, double precision
c . a : matrice ma x na, double precision
c . b : matrice mb x nb, double precisione
c . imap : vettore di mappatura delle righe di a su b
c . jmap : vettore di mappatura delle colonne di a su b
c
c OUTPUT :
c
c b : modificata nei termini imap,jmap
c con b( imap, jmap) = b(imap,jmap) + sca *a
c e b(!imap,*),b(*,!jmap) non modificati
c
c i termini di imap fuori range 1 <= imap(i) <= mb verranno trascurati
c i termini di jmap fuori range 1 <= jmap(j) <= nb verranno trascurati
c
subroutine assmbl(sca,a,ma,na,b,mb,nb,imap,jmap)
implicit none
integer ma,na,mb,nb,imap(ma),jmap(na), i,j
double precision a(ma,na), b(mb,nb),sca
do i=1,ma
if (imap(i) .ge. 1 .and. imap(i) .le. mb) then
do j=1,na
if (jmap(j) .ge. 1 .and. imap(i) .le. nb) then
b(imap(i),jmap(j)) = b(imap(i),jmap(j)) + sca*a(i,j)
end if
enddo
end if
enddo
return
end
==== prodotto matrice - matrice , flessibile ====
c234567---------------------------------------------------------------72
c=========================== prodotto matrice - matrice , flessibile ===
c può essere usato con vettori
c
c ab = a . b
c
c a : matrice m x l double precision
c b : matrice l x n double precision
c ab: matrice m x n double precision
c
subroutine dot(m,l,n,a,b,ab)
implicit none
integer m,n,l,i,j,k
double precision a(m,l),b(l,n),ab(m,n)
c scorro sugli elementi (i,j) della matrice prodotto
do i=1,m
do j=1,n
c azzero il termine ab(i,j)
ab(i,j)=0.0d0
c accumulo i contributi di sommatoria interna
do k=1,l
ab(i,j)=ab(i,j) + a(i,k)*b(k,j)
enddo
enddo
enddo
return
end
==== prodotto matrice trasposta - matrice , flessibile ====
c234567---------------------------------------------------------------72
c================= prodotto matrice trasposta - matrice , flessibile ===
c può essere usato con vettori
c
c ab = a^T . b
c
c a : matrice m x l double precision
c b : matrice l x n double precision
c ab: matrice l x n double precision
c
subroutine tdot(m,l,n,a,b,ab)
implicit none
integer m,n,l,i,j,k
double precision a(m,l),b(l,n),ab(m,n)
c scorro sugli elementi (i,j) della matrice prodotto
do i=1,l
do j=1,n
c azzero il termine ab(i,j)
ab(i,j)=0.0d0
c accumulo i contributi di sommatoria interna
do k=1,m
ab(i,j)=ab(i,j) + a(k,i)*b(k,j)
enddo
enddo
enddo
return
end
==== azzeramento matrice/vettore ====
c234567---------------------------------------------------------------72
c======================================= azzeramento matrice/vettore ===
c passare un valore n pari al numero di celle totali allocate,
c ad es. 12 nel caso di una matrice 3x4
c
subroutine dzero(a,n)
integer n,i
double precision a(n)
do i=1,n
a(i)=0.0d0
enddo
end subroutine