function(datlist,alternative="greater",align="mean",
    Nsim=10000,PDF=F){
#================================================================
# This function needs as input a list of lists, say b lists.
# Each of the b lists should consist of two vectors x and y
# in that order, where y represents the treatment scores
# and x represents the control scores for the block 
# represented by this sublist.

# It is assumed that the random assigment of treatments 
# from block to block are independent.
# alternative="greater" means that we expect that the y-scores
# will tend to be larger than the x-scores under the alternative.
# Other values for alternative are "less" and "two.sided", with
# corresponding meanings.
#
# The argument align chooses the block alignment function. 
# align="mean" causes the subtraction of the block mean 
# from all scores in a block. Similar actions are taken when
# align="median" and "std.residual". In the latter case the 
# block mean is subtracted from all observations in a given 
# block and then these differences are divided by the block 
# Other alignment procedures could easily be added, see the first
# four lines in the code below.
#================================================================
if(align=="mean") AlignFun=function(z){z-mean(z)}
if(align=="median") AlignFun=function(z){z-median(z)}
if(align=="std.residual") AlignFun=function(z){(z-mean(z))/sqrt(var(z))}
if(align!="mean" & align != "median" & 
 align!= "std.residual") return("choose `mean', `median', or `std.residual' for block alignment")
if(alternative!="greater" & alternative!="less" & alternative!="two.sided"){
return("alternative must be `greater' or `less' or `two.sided'")
}
if(PDF==T) pdf(file="AlignedBlockedWilcoxon.pdf",width=11)
b=length(datlist)
zvec=NULL
mvec=NULL
nvec=NULL
Nvec=NULL
MM=1
for(i in 1:b){
   dat.i=datlist[[i]]
   x=dat.i[[1]]
   y=dat.i[[2]]
   nvec=c(nvec,length(y))
   mvec=c(mvec,length(x))
   Nvec=c(Nvec,length(x)+length(y))
   MM=MM*choose(Nvec[i],nvec[i])
}
mean.W=0
var.W=0
W.obs=0
if(MM<=Nsim){
  Wvec=rep(0,MM)
  for(i in 1:b){
     dat.i=datlist[[i]]
     x=dat.i[[1]]
     y=dat.i[[2]]
     z=c(x,y)
     z=AlignFun(z)
     zvec=c(zvec,z)
  }
  Rz=rank(zvec)
  Nz=0
  for(i in 1:b){
     mi=mvec[i]
     ni=nvec[i]
     Ni=mi+ni
     R=Rz[Nz+(1:Ni)]
     W.obs=W.obs+sum(R[mi+1:ni])
     dist.new=combn(R,ni,sum)
     if(i>1){
        Wvec=outer(Wvec,dist.new,"+")}else{
        Wvec=dist.new
     }
     mean.W=mean.W+ni*mean(R)
     var.W=var.W+var(R)*ni*mi/Ni
     Nz=Nz+Ni
  }
}else{
  # computing W.obs and aligned rank set Rz
  for(i in 1:b){
     dat.i=datlist[[i]]
     x=dat.i[[1]]
     y=dat.i[[2]]
     z=c(y,x)
     z=AlignFun(z)
     zvec=c(zvec,z)
  }
  Rz=rank(zvec)
  R.list=list()
  Nz=0
  for(i in 1:b){
        mi=mvec[i]
        ni=nvec[i]
        Ni=mi+ni
        R=Rz[Nz+1:Ni]
        mean.W=mean.W+ni*mean(R)   
        var.W=var.W+var(R)*ni*mi/Ni
        R.list[[i]]=R
        W.obs=W.obs+sum(R[1:ni])
        Nz=Nz+Ni
  }
  # simulation loop
  Wvec=rep(0,Nsim)
  for(j in 1:Nsim){
       WW=0
       for(i in 1:b){
           WW=WW+sum(sample(R.list[[i]],nvec[i]))
       }
       Wvec[j]=WW  
  }
}
sig.W=sqrt(var.W)
m=min(Wvec,mean.W-5*sig.W)
M=max(Wvec,mean.W+5*sig.W)
W.u=unique(sort(Wvec))
d.u=diff(W.u)
k=length(d.u)
br=c(W.u[1]-d.u[1]/2,W.u[1:k]+d.u/2,W.u[1+k]+d.u[k]/2)
if(MM<=Nsim){
out.hist=hist(Wvec,breaks=br,probability=T,xlim=c(m,M),col=c("grey","grey90"),
xlab="",main="Aligned Block Combined Wilcoxon Test")}else{
out.hist=hist(Wvec,nclass=151,probability=T,xlim=c(m,M),col=c("grey","grey90"),
xlab="",main="Aligned Block Combined Wilcoxon Test")
}
mtext(substitute(hat(W)[s]==sum(hat(W)[s]^"(i)", i==1,bx),list(bx=b)),1,3.5)
if(alternative=="greater" | alternative=="less"){
  abline(v=W.obs,col="red")}else{
  offset=abs(W.obs-mean.W)
  abline(v=mean.W+c(-offset,offset),col="red")
}
if(alternative=="less"){
  pval=mean(Wvec<=W.obs)
  pval.norm=pnorm((W.obs-mean.W)/sig.W)
  if(MM<=Nsim){
    text(W.obs-.05*sig.W,max(out.hist$density),"exact",adj=1)}else{
    text(W.obs-.05*sig.W,max(out.hist$density),"simulated",adj=1)
    text(W.obs-.05*sig.W,.70*max(out.hist$density),substitute(N[sim]==xNsim,list(xNsim=Nsim)),adj=1)
  }
  text(W.obs-.05*sig.W,.95*max(out.hist$density),paste("p-value =",round(pval,4)),adj=1)
  text(W.obs-.05*sig.W,.85*max(out.hist$density),paste("normal approximation"),adj=1)
  text(W.obs-.05*sig.W,.80*max(out.hist$density),paste("p-value =",round(pval.norm,4)),adj=1)
}
if(alternative=="greater"){
  pval=mean(Wvec>=W.obs)
  pval.norm=1-pnorm((W.obs-mean.W)/sig.W)
  if(MM<=Nsim){
    text(W.obs+.05*sig.W,max(out.hist$density),"exact",adj=0)}else{
    text(W.obs+.05*sig.W,max(out.hist$density),"simulated",adj=0)
    text(W.obs+.05*sig.W,.70*max(out.hist$density),substitute(N[sim]==xNsim,list(xNsim=Nsim)),adj=0)
  }
  text(W.obs+.05*sig.W,.95*max(out.hist$density),paste("p-value =",round(pval,4)),adj=0)
  text(W.obs+.05*sig.W,.85*max(out.hist$density),paste("normal approximation"),adj=0)
  text(W.obs+.05*sig.W,.80*max(out.hist$density),paste("p-value =",round(pval.norm,4)),adj=0)
}
if(alternative=="two.sided"){
  pval=mean(abs(Wvec-mean.W)>=abs(W.obs-mean.W))
  pval.norm=2*pnorm(-abs(W.obs-mean.W)/sig.W)
  if(MM<=Nsim){
    text(W.obs+.05*sig.W,max(out.hist$density),"exact",adj=0)}else{
    text(W.obs+.05*sig.W,max(out.hist$density),"simulated",adj=0)
    text(W.obs+.05*sig.W,.70*max(out.hist$density),substitute(N[sim]==xNsim,list(xNsim=Nsim)),adj=0)
  }
  text(W.obs+.05*sig.W,.95*max(out.hist$density),paste("2-sided p-value =",round(pval,4)),adj=0)
  text(W.obs+.05*sig.W,.85*max(out.hist$density),paste("normal approximation"),adj=0)
  text(W.obs+.05*sig.W,.80*max(out.hist$density),paste("2-sided p-value =",round(pval.norm,4)),adj=0)
}
x=seq(m,M,length.out=201)
y=dnorm(x,mean.W,sig.W)
lines(x,y,col="blue")

if(PDF==T) dev.off()
}