gammaprop=function(n){
        
    prop=rep(0,n)
    prop[1]=1/sqrt(pi)
    prop[2]=sqrt(pi)/2
    
    m=floor(n/2)
    
    for(i in 1:(m-1)){
        temp=2*i
        prop[temp+1]=prop[temp-1]*(temp)/(temp-1)
        prop[temp+2]=prop[temp]*(temp+1)/(temp)
    }

# if n is odd, input prop[n].
    if(m!=n/2){
        prop[n]=prop[n-2]*(n-1)/(n-2)
    }
    
    prop    
}

rngseed <- function(seed)
{
        seed <- as.integer(seed)
        if(seed <= 0)
                stop("Seed must be a positive integer")
        tmp <- .Fortran("rngs",
                seed)
        invisible()
}

#####Routines start here######
gibbsneteff=function(txx,tyx,tyy,p,nsize,gamma.ini,alpha,bi,bii,bmati,bmatii,nbi,nbii,bvarsq,
set,nset,setsort,nsort,nop,niter,nburnin,gammaout){

rngseed(floor(runif(1)*7483647))

set=as.integer(set)
nset=as.integer(nset)
setsort=as.integer(setsort)
nsort=as.integer(nsort)
nop=as.integer(nop)
modelsize=nop

start=as.integer(1)
fin=as.integer(nop)

gammain=as.integer(gamma.ini)

niter=as.integer(niter)
nburnin=as.integer(nburnin)
#p=ncol(xmat)
p=as.integer(p)
#nsize=nrow(xmat)
nsize=as.integer(nsize)
#nu=as.integer(nu)
diag(txx)=diag(txx)+(1/(bvarsq^2))
#tyx=t(y)%*%xmat
#tyx=as.vector(tyx)
#tyy=sum(y*y)



gmrgcnt=rep(0,p)
gmrgcnt=as.integer(gmrgcnt)


detanoi=0.0
detai=0.0
tyxanoitxy=0.0
tyxaitxy=0.0
bf=0.0

chgcount=rep(0,p)
chgcount=as.integer(chgcount)

prgpost=rep(0.0,p)

dveci=rep(0.0,nsort)
dvec=rep(0.0,nsort)
#obtain the first lmat and dvec
lmat=matrix(0.0,nrow=nsort,ncol=nsort)

settemp=setsort[1:nop]
mattemp=solve(txx[settemp,settemp])
cmat=t(chol(mattemp))
dvec[1:nop]=diag(cmat)

for(i in 1:nop){
    cmat[,i]=cmat[,i]/dvec[i]
}

lmat[1:nop,1:nop]=cmat

dvec[1:nop]=dvec[1:nop]^2

out=.Fortran("gibbseffnetf",niter=niter,nburnin=nburnin,gammain=gammain,n=nsize,p=p,txx=txx,tyx=tyx,tyy=tyy,
chgcount=chgcount,lmat=lmat,dvec=dvec,dveci=dveci,set=set,nset=nset,setsort=setsort,nsort=nsort,nop=nop,
start=start,fin=fin,detanoi=detanoi,detai=detai,tyxanoitxy=tyxanoitxy,tyxaitxy=tyxaitxy,bvarsq=bvarsq,
bf=bf,prgpost=prgpost,alpha=alpha,bi=bi,bii=bii,bmati=bmati,bmatii=bmatii,nbi=nbi,nbii=nbii,modelsize=modelsize,
gmrgcnt=gmrgcnt,gammaname=gammaout)
} 

gibbsneteff1=function(txx,tyx,tyy,p,nsize,gamma.ini,alpha,bi,bii,bmati,bmatii,nbi,nbii,bvarsq,
set,nset,setsort,nsort,nop,niter,nburnin,gammaout){

rngseed(floor(runif(1)*7483647))

set=as.integer(set)
nset=as.integer(nset)
setsort=as.integer(setsort)
nsort=as.integer(nsort)
nop=as.integer(nop)
modelsize=nop
indices=rep(0,nsort)
indices=as.integer(indices)

start=as.integer(1)
fin=as.integer(nop)

gammain=as.integer(gamma.ini)

niter=as.integer(niter)
nburnin=as.integer(nburnin)
#p=ncol(xmat)
p=as.integer(p)
#nsize=nrow(xmat)
nsize=as.integer(nsize)
#nu=as.integer(nu)
diag(txx)=diag(txx)+(1/(bvarsq^2))
#tyx=t(y)%*%xmat
#tyx=as.vector(tyx)
#tyy=sum(y*y)



gmrgcnt=rep(0,p)
gmrgcnt=as.integer(gmrgcnt)


detanoi=0.0
detai=0.0
tyxanoitxy=0.0
tyxaitxy=0.0
bf=0.0

chgcount=rep(0,p)
chgcount=as.integer(chgcount)

prgpost=rep(0.0,p)

dveci=rep(0.0,nsort)
dvec=rep(0.0,nsort)
#obtain the first lmat and dvec
lmat=matrix(0.0,nrow=nsort,ncol=nsort)

settemp=setsort[1:nop]
mattemp=solve(txx[settemp,settemp])
cmat=t(chol(mattemp))
dvec[1:nop]=diag(cmat)

for(i in 1:nop){
    cmat[,i]=cmat[,i]/dvec[i]
}

lmat[1:nop,1:nop]=cmat

dvec[1:nop]=dvec[1:nop]^2

out=.Fortran("gibbseffnet1f",niter=niter,nburnin=nburnin,gammain=gammain,n=nsize,p=p,txx=txx,tyx=tyx,tyy=tyy,
chgcount=chgcount,lmat=lmat,dvec=dvec,dveci=dveci,set=set,nset=nset,setsort=setsort,nsort=nsort,nop=nop,
start=start,fin=fin,detanoi=detanoi,detai=detai,tyxanoitxy=tyxanoitxy,tyxaitxy=tyxaitxy,bvarsq=bvarsq,
bf=bf,prgpost=prgpost,alpha=alpha,bi=bi,bii=bii,bmati=bmati,bmatii=bmatii,nbi=nbi,nbii=nbii,indices=indices,modelsize=modelsize,
gmrgcnt=gmrgcnt,gammaname=gammaout)
}