Differenze

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

--- wikipaom2017:051:030:000 [2017/03/31 11:54] – [prodotto matrice - matrice , flessibile] ebertocchi
+++ wikipaom2017: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>