elim.Gauss.dnz <- function(A, dbg=FALSE) {
# caso di elementi tutti gli elementi della diagonale diversi da zero
# (anche questo controllo è prudenziale, in quanto se ci sono    
#  altri elementi nulli sulla diagonale, lo swap delle
#  righe potrebbe risolvere il problema)

    if(dim(A)[1] != dim(A)[2]) {
    cat(sprintf("La matrice deve essere quadrata\n"))
      return(NULL)
    }
    if (any(diag(A) ==0))  {
        cat(sprintf("Ci sono elementi nulli sulla diagonale nulli!\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)
}