#####
rZ_SR_rl<-function(S,Ranks,Z)
{
  n<-dim(Z)[1] ; p<-dim(Z)[2]
  for(j in sample(1:p))
  {
    Sjc<- S[j,-j]%*%solve(S[-j,-j])
    sdj<- sqrt( S[j,j] -S[j,-j]%*%solve(S[-j,-j])%*%S[-j,j] )
    muj<- Z[,-j]%*%t(Sjc)
    for(r in sample( unique(Ranks[  !is.na(Ranks[,j])  ,j]) ) )
    {
      ir<- (1:n)[Ranks[,j]==r & !is.na(Ranks[,j])]
      lb<-suppressWarnings(max( Z[ Ranks[,j]<r,j],na.rm=TRUE ))
      ub<-suppressWarnings(min( Z[ Ranks[,j]>r,j],na.rm=TRUE ))
      z<-qnorm(runif(length(ir), pnorm(lb,muj[ir],sdj),pnorm(ub,muj[ir],sdj)),
                muj[ir],sdj)
      z[z==Inf]<-max( lb,z[z!= Inf] )
      z[z==-Inf]<-min( ub,z[z!= -Inf] )
      Z[ir,j]<-z
    }
    ir<-(1:n)[is.na(Ranks[,j])]
    Z[ir,j]<-rnorm(length(ir),muj[ir],sdj)
  }

  Z                               }
#####