matrixmult <- function(n) {
  M <- array(0.,dim=c(n,n));
  for (i in 1:n) {
    for (j in 1:n) {
      M[i,j] <- i+j+0.5;
    }
  }
  sum(M %*% M); 
}