fam <- structure(list(ID = c(1L, 2L, 3L, 4L, 6L, 5L, 7L), dad = c(0L,
                                                                  0L, 1L, 1L, 1L, 3L, 5L), mum = c(0L, 0L, 0L, 2L, 4L, 4L, 6L),
                      GEN = c(1L, 1L, 2L, 2L, 3L, 3L, 4L)), class = "data.frame", row.names = c(NA,
                                                                                                -7L))

library(Matrix)

hom = function(fam) {
  t1 <- min(which.max(fam$dad > 0), which.max(fam$mum > 0))
  t2 <- max(fam[["ID"]])

  A<-Matrix(0, nrow=t2,ncol=t2, sparse=TRUE)
  diag(A) <- 2-0.5^(fam[["GEN"]]-1)

  for (t in t1:t2) {
    A[t,t]<- sum(c(A[t,t],     0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]))
    for(j in 1:(t-1))  {
      A[t,j]<- 0.5 * sum(c(A[j,fam[t,"dad"]],A[j,fam[t,"mum"]]))
      A <- tril(A) + t(tril(A, -1))
    }
  }
  A
}