# NOT RUN {
m = 3
h = 1
J = Jn(n = m+1+h)
Q = Qm(m = m, nu = 0)
# using K
K = Km(J = J[1:(m+1),1:(m+1)], Q = Q)
Kkron0 = Kkronm(K = K, h = h)
# using m and nu
Kkron = Kkronm(m = m, nu = 0, h = h)
# using kronecker product from R
K = Km(J = J[1:(m+1),1:(m+1)], Q = Q)
Kh = rbind(matrix(0, nrow = h, ncol = m+1+h),
cbind(matrix(0, nrow = m+1, ncol = h), K))
KkronR = K %x% Kh
# using the definition K* = (Jm %x% J)'(Q %x% Q)(Jm %x% J)
J_m = J[1:(m+1),1:(m+1)]
J_h = J[(h+1):(m+1+h),1:(m+1+h)]
KkronD = t(J_m %x% J_h)%*%(Q %x% Q)%*%(J_m %x% J_h)
# comparing the results
sum(abs(Kkron0 - Kkron))
sum(abs(Kkron0 - KkronR))
sum(abs(Kkron0 - KkronD)) # difference due to rounding error
# }
# NOT RUN {
# Function Kkronm is computationaly faster than a pure implementation in R:
m = 100
h = 1
J = Jn(n = m+1)
Q = Qm(m = m, nu = 0)
# using Kkronm
t1 = proc.time()
Kkron = Kkronm(m = m, nu = 0, h = 1)
t2 = proc.time()
# elapsed time:
t2-t1
# Pure R implementation:
K = Km(J = J, Q = Q)
Kh = rbind(matrix(0, nrow = h, ncol = m+1+h),
cbind(matrix(0, nrow = m+1, ncol = h), K))
t3 = proc.time()
KkronR = K %x% Kh
t4 = proc.time()
# elapsed time
t4-t3
# }
Run the code above in your browser using DataLab