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

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

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

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

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

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