So I think I managed to program the block diagonal function, but something is not working out:

```
matrix block_diag_matrix(matrix matrix_diag, int K, int J) {
matrix[K,K] O;
matrix[K,K] A[J,J];
matrix[J*K,J*K] rtn;
O = rep_matrix(0, K, K);
for (s in 1:J) {
for (t in 1:J) {
if (t==s){
A[s,t]= matrix_diag;
}
else{
A[s,t]= O;
}
}
}
for (i in 1:J){
for (j in 1:J){
for (n in 1:K){
for (m in 1:K){
rtn[i*m,j*n] = A[i, j, m, n];
}
}
}
}
return rtn;
}
```

I then code the model like this:

```
Beta ~ multi_normal(mu, Omega);
```

And I get this error:

```
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: multi_normal_lpdf: Covariance matrix is not symmetric. Covariance matrix[1,127] = nan, but Covariance matrix[127,1] = nan (in 'model2de871122f6a_Block_Diagonal_Work_in_progress' at line 155)
```

I do not understand what is wrong though.

Edit: I do understand what is wrong, the indices are not correct. However I do not understand how to correct them.