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;
    }
  }
  r <- 0;
  for (i in 1:n) {
    for (j in 1:n) {
      for (k in 1:n) {
        r <- r+M[i,k]*M[j,k];
      }
    }
  }
  r;
}