Y<-matrix(c(110,100,110,100,100,110,110,100,100),nrow=9,ncol=1)
X<-matrix(c(1,1,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,1,1,1,1),nrow=9,ncol=3)

Z<-matrix(c(1,0,0,0,0,0,0,0,0,
            0,0,1,0,0,0,0,0,0,
            0,0,0,0,0,1,1,0,0,
            0,1,0,1,1,0,0,1,1),nrow=9,ncol=4)
R<-diag(1,nrow=9,ncol=9)
G<-diag(0.1,nrow=4,ncol=4)
A<-t(X)%*%solve(R)%*%X
B<-t(X)%*%solve(R)%*%Z
C<-t(Z)%*%solve(R)%*%X
D<-t(Z)%*%solve(R)%*%Z+solve(G)

E<-cbind(A,B)
F<-cbind(C,D)

M<-rbind(E,F)

t(X)%*%solve(R)%*%Y

##      [,1]
## [1,]  210
## [2,]  310
## [3,]  420

t(Z)%*%solve(R)%*%Y

##      [,1]
## [1,]  110
## [2,]  110
## [3,]  220
## [4,]  500

N<-rbind(t(X)%*%solve(R)%*%Y,t(Z)%*%solve(R)%*%Y)

solve(M,N)

##             [,1]
## [1,] 105.6386574
## [2,] 104.2757378
## [3,] 105.4584366
## [4,]   0.3964857
## [5,]   0.5203875
## [6,]   0.7569272
## [7,]  -1.6738004