Strumenti Utente

Strumenti Sito


wikipaom2017:051:030:000

Differenze

Queste sono le differenze tra la revisione selezionata e la versione attuale della pagina.

Link a questa pagina di confronto

Entrambe le parti precedenti la revisioneRevisione precedente
Prossima revisione
Revisione precedente
wikipaom2017:051:030:000 [2017/03/31 11:53] ebertocchiwikipaom2017:051:030:000 [2017/03/31 11:55] (versione attuale) ebertocchi
Linea 1: Linea 1:
 +====== Routine ausiliarie ======
 +===== Matrici e vettori =====
 +==== Calcolo di determinante per matrici 2x2 o 3x3 ====
 +<code>
 +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
 +</code>
  
 +==== Calcolo dell'inversa di una matrice 2x2 o 3x3 ====
 +
 +<code>
 +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
 +</code>
 +
 +==== Assemblaggio matrice su matrice ====
 +
 +<code>
 +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
 +</code>
 +
 +==== prodotto matrice - matrice , flessibile ====
 +
 +<code>
 +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
 +
 +</code>
 +==== prodotto matrice trasposta - matrice , flessibile ====
 +
 +<code>
 +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
 +</code>
 +==== azzeramento matrice/vettore ====
 +
 +
 +<code>
 +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
 +</code>