# metodo di eliminazione Gaussiana

elim.Gauss <- function(A, dbg=FALSE) {
# (con gestione di trattamento di pivot nulli) 
    if(dim(A)[1] != dim(A)[2]) {
    cat(sprintf("La matrice deve essere quadrata\n"))
      return(NULL)
    }    
    nr = nc = dim(A)[1]
    if(nr == 1) {
        return(A)
    }

    if(dbg) print(A)
    for (k in 2:nr) {
        if(dbg) print(k)
        pivot = A[k-1,k-1]
        if(dbg) cat(sprintf("Pivot riga %d: %f\n", k-1, pivot))
        if (pivot == 0) {   # cerca la prima riga con candidato pivot != 0
          if(dbg) cat(sprintf("Pivot nullo in riga = %d\n", k-1))
          new.row = 0  
          for (i in k:nr) {
              if (A[i,k-1] != 0) {
                  new.row = i
                  if(dbg)  cat(sprintf("new.row %d \n", new.row ))
                  tmp = A[k-1,]
                  if(dbg) print(tmp)
                  A[k-1,] = A[i,]
                  A[i,] = tmp
                  if(dbg) print(A)
              }
              if(dbg) cat(sprintf("dopo scambio fra riga %d e riga %d\n", k-1, i))
              if(dbg) print(A)
              if(new.row != 0) break
          }
          if (new.row == 0) {
              cat(sprintf("Non è stato possibile scambiare due righe\n"))
              return(NULL)
          }
        }
        pivot = A[k-1,k-1]
        if(dbg) cat(sprintf("Pivot riga %d: %f\n", k-1, pivot))
        for (i in k:nr) {
            A[i,] = A[i,] -  A[k-1,] *  A[i,k-1] / pivot
        }        
        if(dbg) print(A)        
    }
    
    return(A)
}