trace <- function(X) {
## Function to calculate the trace of
#  a square matrix.

  dd = dim(X)
  if (dd[1]!=dd[2]) {
   stop("The trace is only defined for square matrices.")
  }

sum(diag(X))
}