Thanks again for these suggestions. After a lot of checking (unfortunately I cannot debug directly…), I found out that the problem was caused by the mismatch between the format of the reals and integers in R when calling the Fortran subroutines.
Luckily I have been able to get the full subroutine running (for now) on my Mac and on a Windows machine. I still get an error though when checking it on a Linux machine. I build the subroutine again in small steps and a problem seems to appear in the following matrix multiplication step:
meanMat(1:Njs(g1),1:P) = matmul(Xgroups(g1,1:Njs(g1),1:K), &
BDraws(g1,1:K,1:P))
Here Njs(g1)=20, g1=1, K=1, P=3. Even though this does not give problems on a Mac and Windows machine, I could imagine that this may cause a problem because K is 1. Based on what is said here (Solved: about matmul() - Intel Community), the dimensionality would not be a problem as 1:K would remain to be an array. But I also checked it using the reshape function:
meanMat(1:Njs(g1),1:P) = matmul(reshape(Xgroups(g1,1:Njs(g1),1:K),(/Njs(g1),K/)), &
reshape(BDraws(g1,1:K,1:P),(/K,P/)))
but I still get an error. Would there be other (less error prone) ways to do this matrix multiplication? Of course, it could also be that the error is caused somewhere else but for some reason this line kills the check (and then the provided information would be insufficient). Thanks again.