elim.Gauss.nz <- function(A, dbg=FALSE) {
# caso di elementi tutti diversi da zero 
# (prudenziale! -> solo i pivot devono essere non nulli)
    if(dim(A)[1] != dim(A)[2]) {
    cat(sprintf("La matrice deve essere quadrata\n"))
      return(NULL)
    }
    if (any(A ==0))  {
        cat(sprintf("Questo metodo funziona solo con tutti elementi !=0\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)        
       for (i in k:nr) {
          A[i,] = A[i,] -  A[k-1,] *  A[i,k-1]/A[k-1,k-1]
       }
       if(dbg) print(A)
    }
    
    return(A)
}