class Array
   def add(v)
      v.each_index{|i| self[i] += v[i]}
   end
   def sub(v)
      v.each_index{|i| self[i] -= v[i]}
   end
   def scalar(v)
      sum = 0.0
      v.each_index{|i| sum += self[i]*v[i]}
      return sum
   end
end

class Matrice

   attr_accessor :m, :n

   @@icall=0
   # usato per i numeri random
   def myrand()
      if @@icall==0
         @@icall=1
         srand(93467)
      end
      return rand
   end

   # costruttore
   def initialize (ndim)
      @n = ndim
      @m = Array.new(@n*@n,0.0)
   end

   def add(c)
      m.add(c.m)
   end

   def sub(c)
      m.sub(c.m)
   end

   def copia(a)
      (0...@n).each do |i|
         (0...@n).each do |j|
            m[i*@n+j] = a.get(i,j)
         end
      end
   end

   def riempiRandom(i)
   srand(i)
   # riempie un vettore con numeri random
      @m.each_index {|j| @m[j] = rand}
   end

   def get(i, j)
      return  @m[i*@n+j]
   end

   # scrittura di un elemento di un vettore
   def set(i, j, x)
      @m[i*@n+j]=x
   end

   # lettura della dimensione del vettore
   def getOrder
      return @n
   end

#  calcolo l'inversa di una matrice
   def inversa

      b = Vettore.new(@n)
      x = Vettore.new(@n)
      ainv = Matrice.new(@n)
      a = Matrice.new(@n)

      for col in 0...@n
#         b.v.each {|v| v=0.0}
         for j in 0...@n
            b.set(j,0.0)
         end
         b.set(col,1.0)
         ainv.stampa
         a.copia(self)
         risolviSistema(a,b,x)
         for j in 0...@n
#            ainv.@m[col+j*n] = x.get(j)
            ainv.set(j,col,x.get(j))
         end
      end
      ainv.stampa

      return ainv

   end

   def postoPivot(i)

      # Calcolo la riga in cui si trova l'elemento pivot
      # sulla colonna c e le righe tra c e n comprese

      res = (@m[i*n+i]).abs # elemento diagonale
      k=i
      for j in (i+1)...n
         x=(@m[j*n+i]).abs
         if(x > res )
            res=x
            k=j
         end
      end
      return k
   end

   def scambiaRighe(b, i, j)
   # scambia due righe della matrice e i corrispondenti termini noti

      for k in 0...(b.getOrder)
         @m[i*n+k], @m[j*n+k] = @m[j*n+k], @m[i*n+k]
      end
      b.v[i], b.v[j] = b.v[j], b.v[i]

   end

   def stampa
      print "\n"
      (0...@n).each do |j|
         (0...@n).each do |k|
            printf("%8.4f   ", @m[j*@n+k])
         end
         print "\n"
      end
   end

   def determinante
      # Calcola il determinante di una matrice in forma triangolare
      # e di cui sia nota la parita' d degli scambi di riga

      # eliminazione gaussiana
      (0...@n).each do |i| # per ogni riga
         (i+1...@n).each do |j|  # parto dall riga successiva
            c= @m[j*@n+i]/@m[i*@n+i]  # calcolo il rapporto
            # sottraggo un multiplo della riga i alla riga j
            for k in i...@n
               @m[j*@n+k] -= c*@m[i*@n+k]
            end
         end
      end

      # prodotto degli elementi diagonali
      prod = 1.0
      (0...@n).each do |i|
         prod *= @m[i*@n+i]
      end
      return prod

   end

   public   :postoPivot, :scambiaRighe, :inversa, :copia
end #fine classe Matrice


class Vettore
   attr_accessor :v, :n

   @@icall=0
   # usato per i numeri random
   def myrand()
      if @@icall==0
         @@icall=1
         srand(93467)
      end
      return rand
   end

   # costruttore
   def initialize (ndim)
      @n = ndim
      @v = Array.new(@n,0.0)
   end

   def add(c)
      @v.add(c.v)
   end

   def sub(c)
      @v.sub(c.v)
   end

   def copia(b)
      @v.each_index {|i| @v[i] = b.get(i)}
   end

   def scalar(c)
      return @v.scalar(c.v)
   end

   def riempiRandom(i)
   srand(i)
   # riempie un vettore con numeri random
      @v.each_index {|j| @v[j] = rand}
   end

   def get(i)
      return  @v[i]
   end

   # scrittura di un elemento di un vettore
   def set(i, x)
      @v[i]=x
   end

   # lettura della dimensione del vettore
   def getOrder
      return @v.size
   end

   def stampa
      print "\n"
      @v.each_index {|j| printf("%8.4f   \n", @v[j])}
      print "\n"
   end

end # fine classe vettore

# moltiplica una matrice per un vettore
def mulMV(a, b)

   n = a.getOrder

   x = Vettore.new(n)

   (0...n).each do |i|
      sum = 0.0
      (0...n).each do |j|
         sum += a.get(i,j)*b.get(j)
      end
      x.set(i,sum)
   end

   return x
end


def eliminazioneGaussiana(a, b)

   # eliminazione gaussiana
   n = a.getOrder

   (0...n).each do |i| # per ogni riga
      (i+1...n).each do |j|  # parto dall riga successiva
         c = a.get(j,i)/a.get(i,i)  # calcolo il rapporto
         # sottraggo un multiplo della riga i alla riga j
         (i...n).each {|k| a.set(j,k,a.get(j,k)-c*a.get(i,k))}
         b.set(j,b.get(j)-c*b.get(i))
      end
   end

end


def Pivoting(a, b)

   n = a.getOrder

   # cerco l'elemento massimo nelle varie colonne
   (0...n).each do |i|
      j = a.postoPivot(i)
      # scambio la riga del pivot con quella iniziale
      if(i!=j)
         a.scambiaRighe(b,i,j)   #scambio le righe e i termini noti
      end
   end
end


def backSubstitution(a, b, x)

   n = a.getOrder()

   i = n-1;
   while i>=0
      temp = b.get(i)
      for j in (i+1)...n
         temp -= a.get(i,j)*x.get(j)
      end
      temp /= a.get(i,i)
      x.set(i,temp)
      i -= 1
   end
end

def risolviSistema(a, b, x)

   ap = Matrice.new(a.getOrder)
   bp = Vettore.new(b.getOrder)

   ap.copia(a)
   bp = b

   #
   # Risolve il sistema a*x=b
   # x contiene la soluzione
   #
   Pivoting(ap,bp)
   #
   eliminazioneGaussiana(ap,bp)
   backSubstitution(ap,bp,x)
end

def verificaSoluzione(a, b, x)
   print "Verifico la soluzione\n"
   mulMV(a,x).stampa
   b.stampa
end

# moltiplica due matrici
def mulMM(a1,a2)

   n = a1.getOrder

   a = Matrice.new(n)

   (0...n).each do |j|
      (0...n).each do |k|
         temp = 0.0
         (0...n).each do |l|
            temp += a1.get(j,l)*a2.get(l,k)
         end
         a.set(j,k,temp)
      end
   end
   return a
end

n=5   # ordine della matrice

# definisco gli oggetti che mi servono
a = Matrice.new(n)
ainv = Matrice.new(n)
a_copia = Matrice.new(n)

# riempio la matrice con numeri random
a.riempiRandom(60623)
a.stampa
# copia di sicurezza di 'a'
a_copia.copia(a)
print "\nCopia di a\n"
a_copia.stampa

# stesse operazioni con i vettori
b = Vettore.new(n)
p "Vettore nullo di #{n} elementi\n"
b.stampa
b.riempiRandom 8276
p "Vettore random di #{n} elementi\n"
b.stampa
b_copia = Vettore.new(n)
b_copia.copia(b)
p "Vettore copia del precedente\n"
b_copia.stampa
b.add(b_copia)
p "Vettore doppio del precedente\n"
b.stampa

# moltiplica la matrice a per il vettore b
# e stampa il risultato
mulMV(a,b).stampa

print "Prodotto scalare\n"
print b.scalar b

x = Vettore.new(n)

print "\nCopia di a 3a\n"
a_copia.stampa

a.copia(a_copia)
b.copia(b_copia)
# la risoluzione del sistema distrugge a e b
risolviSistema(a,b,x)
print "\nCopia di a 3b\n"
a_copia.stampa

print "\n\n Soluzione: \n"
x.stampa

a.copia(a_copia)
b.copia(b_copia)
# recupera la matrice originale
verificaSoluzione(a,b,x)

# Calcolo il determinante
print "\nDeterminante: \n"
print a.determinante()
print "\n"
print "A dopo Determinante"
a.stampa

a.copia(a_copia)
# calcolo l'inversa
print "\nOra calcolo l'inversa\n"
a.stampa
ainv.copia(a.inversa)
print "\nMatrice inversa\n"
ainv.stampa
print "\nA dopo Inversa: \n"
a.stampa
# verifico l'esattezza della matrice inversa
# a contiene l'inversa, a_copia la matrice originale
print "\nVerifica inversa\n\n"
mulMM(ainv,a).stampa

print "\nA dopo verifica inversa: \n"
a.stampa