This is from Legarra & Misztal 2008

• Fortran code:

Double precision:: xpx(neq),y(ndata),e(ndata),X(ndata,neq), &
sol(neq),lambda,lhs,rhs,val
do i=1,neq
	xpx(i)=dot_product(X(:,i),X(:,i)) !form diagonal of X′X
enddo
lambda=vare/vara
e=y
do until convergence
	do i=1,neq
		!form lhs
		lhs=xpx(i)+lambda
		! form rhs with y corrected by other effects (formula 1)
		rhs=dot_product(X(:,i),e) +xpx(i) *sol(i)
		! do Gauss Seidel
		val=rhs/lhs
		! MCMC sample solution from its conditional (commented out here)
		! val=normal(rhs/lhs,1d0/lhs)
		! update e with current estimate (formula 2)
		e=e - X(:,i)*(val-sol(i))
		!update sol
		sol(i)=val
	enddo
enddo

• R code:

#get diagonal of X’X
for (i in 1:neq) {  
  xpx[i]=crossprod(X[,i],X[,i])
}
for (iter in 1:1000) {
  #Gauss Seidel
  for (i in 1:neq){
    lhs=xpx[i]+lambda
    rhs=crossprod(X[,i],e) + xpx[i]*ahat[i]
    val=rhs/lhs
    e = e - X[,i]*(val - ahat[i])
    ahat[i]=val
  }
}