P1 <- data.frame(split(rnorm(5*16, 1, 1), 1:16))
M1 <- matrix(unlist(P1[1,]), nrow = 4)
M1[upper.tri(M1)] <- t(M1)[upper.tri(M1)]

P2 <- data.frame(split(rnorm(5*16, 1, 1), 1:16))
M2 <- matrix(unlist(P2[1,]), nrow = 4)
M2[upper.tri(M2)] <- t(M2)[upper.tri(M2)]

B <- rnorm(4, 0, 1)

R1 <- M1 %*% B
R2 <- M2 %*% B
cor(R1, R2)

com_rsk_p2 <- function(m1, m2, n = 2){
  nitt <- length(m1[,1])
  k <- sqrt(length(m1))
  B <- split(rnorm(n*k, 0, 1), 1:n)

  rv_cor <- split(rep(NA, times = n*nitt), 1:nitt)

  for(i in 1:nitt){
    R1 <- sapply(B, function(x) x %*% matrix(unlist(m1[i,]), ncol = k))
    R2 <- sapply(B, function(x) x %*% matrix(unlist(m2[i,]), ncol = k))

    rv_cor[[i]] <- diag(matrix(mapply(cor, list(R1), list(R2)), ncol = n))
  }

  return(t(data.frame(rv_cor)))
}

> out <- com_rsk_p2(P1, P2)
> out
         [,1]       [,2]
X1  0.7622732  0.8156658
X2  0.4414054  0.4266134
X3  0.4388098 -0.1248999
X4  0.5438046  0.7723585
X5 -0.5833943 -0.5294521