module draw
!	!DIR$ NOOPTIMIZE
	use basicmodule
	implicit none

	contains
	
    subroutine draw_muU
		integer	:: ir,i
		real	:: m,prec,s2,v(1),u(0:q)

		prec=muprec_prior
		m=0
		do ir=1,n
			i=clubid(ir)
			s2=om2*(1-rhos(ir)**2)*kappa2(ir)
			u=Us(:,ir)-rhos(ir)*CFs(:,i)
			m=m+sum(Siu(:,0,cindu(ir))*u)/s2
			prec=prec+Siu(0,0,cindu(ir))/s2
		enddo
		call rnnoa(v)
		muU=m/prec+v(1)/sqrt(prec)
	end subroutine
	
	subroutine draw_clubid
		real	:: u(0:q),v(0:q),s2,p(nclubs),unif(1)
		integer	:: ir,i
!$omp parallel do private(i,u,v,s2,p,unif) 
		do ir=1,n
			v=Us(:,ir)
			v(0)=v(0)-muU
			s2=(1-rhos(ir)**2)*kappa2(ir)*om2
			do i=1,nclubs
				u=v-rhos(ir)*CFs(:,i)
				p(i)=-.5*sum(matmul(Siu(:,:,cindu(ir)),u)*u)/s2
			enddo
			p=exp(p-maxval(p))
			do i=2,nclubs
				p(i)=p(i)+p(i-1)
			enddo
			p=p/p(nclubs)
			call rnun(unif)
			clubid(ir)=count(unif(1)>p)+1
		enddo
	end subroutine
	
	subroutine draw_CoCid
		real	:: u(0:q),v(0:q),s2,p(nCoCs),unif(1)
		integer	:: i,k

!$omp parallel do private(k,u,v,s2,p,unif) 
		do i=1,nclubs
			v=CFs(:,i)
			s2=(1-CoCrhos(i)**2)*clubs2(i)*om2
			do k=1,nCoCs
				u=v-CoCrhos(i)*CoCFs(:,k)
				p(k)=-.5*sum(matmul(Siu(:,:,cindCF(i)),u)*u)/s2
			enddo
			p=exp(p-maxval(p))
			do k=2,nCoCs
				p(k)=p(k)+p(k-1)
			enddo
			p=p/p(nCoCs)
			call rnun(unif)
			CoCid(i)=count(unif(1)>p)+1
		enddo
	end subroutine

	subroutine draw_rho
		integer	:: ir,i,l
		real	:: p(nrhogrid),unif(1),s2
		real	:: u(0:q)
		
!$omp parallel do private(i,l,u,s2,p,unif) 
		do ir=1,n
			i=clubid(ir)
			do l=1,nrhogrid
				u=Us(:,ir)-rhogrid(l)*CFs(:,i)
				u(0)=u(0)-muU
				s2=(1-rhogrid(l)**2)*kappa2(ir)*om2
				p(l)=-.5*sum(matmul(Siu(:,:,cindu(ir)),u)*u)/s2-0.5*(q+1)*log(s2)+log(rhodist(l))
			enddo
			p=exp(p-maxval(p))
			do l=2,nrhogrid
				p(l)=p(l)+p(l-1)
			enddo
			p=p/p(nrhogrid)
			call rnun(unif)
			rhos(ir)=rhogrid(count(unif(1)>p)+1)
		enddo
	end subroutine

	subroutine draw_CoCrho
		integer	:: i,l,k
		real	:: p(nrhogrid),unif(1),s2
		real	:: u(0:q)
		
!$omp parallel do private(k,l,u,s2,p,unif) 	
		do i=1,nclubs
			k=CoCid(i)
			do l=1,nrhogrid
				u=CFs(:,i)-rhogrid(l)*CoCFs(:,k)
				s2=(1-rhogrid(l)**2)*clubs2(i)*om2
				p(l)=-.5*sum(matmul(Siu(:,:,cindCF(i)),u)*u)/s2-0.5*(q+1)*log(s2)+log(CoCrhodist(l))
			enddo
			p=exp(p-maxval(p))
			do l=2,nrhogrid
				p(l)=p(l)+p(l-1)
			enddo
			p=p/p(nrhogrid)
			call rnun(unif)
			CoCrhos(i)=rhogrid(count(unif(1)>p)+1)
		enddo
	end subroutine
	
	subroutine draw_ind(nmeas,measind,meas,measdist,meass2)
		real	:: meas(0:,:),measdist(:),meass2(:)
		integer	:: measind(:),nmeas
		
		real	:: p(nth),unif(1),u(0:q),s2
		integer	:: i,j,jtrans,j0
		
!$omp parallel do private(u,s2,j0,j,p,unif)
		do i=1,nmeas
			u=meas(:,i)
			s2=om2*meass2(i)
			do j=1,nth
				p(j)=-.5*sum(matmul(Siu(:,:,j),u)*u)/s2+log(measdist(j)/sqrt(detSu(j)))
			enddo
			p=exp(p-maxval(p))
			do j=2,nth
				p(j)=p(j-1)+p(j)
			enddo
			p=p/p(nth)	
			call rnun(unif)
			measind(i)=count(unif(1)>p)+1
		enddo
		
	end subroutine

	subroutine draw_cindu
		real	:: meas(0:q,n)
		integer	:: ir,i

		do ir=1,n
			i=clubid(ir)
			meas(:,ir)=Us(:,ir)-rhos(ir)*CFs(:,i)
			meas(0,ir)=meas(0,ir)-muU
		enddo
		call draw_ind(n,cindu,meas,thdist,kappa2*(1-rhos**2))
	end subroutine

	subroutine draw_cindCF
		real	:: meas(0:q,nclubs)
		integer	:: i,k
		do i=1,nclubs
			k=CoCid(i)
			meas(:,i)=CFs(:,i)-CoCrhos(i)*CoCFs(:,k)
		enddo
		call draw_ind(nclubs,cindCF,meas,CFthdist,clubs2*(1-CoCrhos**2))
	end subroutine

	subroutine draw_cindCoC
		call draw_ind(nCoCs,cindCoC,CoCFs,CoCthdist,CoCs2)
	end subroutine

	subroutine draw_kappa2grid
		real	:: s2,u(0:q),p(ns2grid),usu,unif(1)
		integer	:: ir,i,l
		real	:: pbase(ns2grid)

		do l=1,ns2grid
			pbase(l)=-.5*(q+1)*log(s2grid(l))+log(s2dist(l))
		enddo
!$omp parallel do private(i,u,s2,usu,p,unif,l)
		do ir=1,n
			i=clubid(ir)
			u=Us(:,ir)-rhos(ir)*CFs(:,i)
			u(0)=u(0)-muU
			s2=om2*(1-rhos(ir)**2)
			usu=sum(matmul(Siu(:,:,cindu(ir)),u)*u)
			do l=1,ns2grid
				p(l)=-.5*usu/(s2*s2grid(l))
			enddo
			p=p+pbase
			p=exp(p-maxval(p))
			do l=2,ns2grid
				p(l)=p(l)+p(l-1)
			enddo
			p=p/p(ns2grid)
			call rnun(unif)
			kappa2(ir)=s2grid(count(unif(1)>p)+1)
		enddo
	end subroutine

	subroutine draw_clubs2grid
		real	:: s2,u(0:q),p(ns2grid),usu,unif(1)
		integer	:: i,k,l
		real	:: pbase(ns2grid)

		do l=1,ns2grid
			pbase(l)=-.5*(q+1)*log(s2grid(l))+log(CFs2dist(l))
		enddo
!$omp parallel do private(k,u,s2,usu,p,unif,l)
		do i=1,nclubs
			k=CoCid(i)
			u=CFs(:,i)-CoCrhos(i)*CoCFs(:,k)
			s2=om2*(1-CoCrhos(i)**2)
			usu=sum(matmul(Siu(:,:,cindCF(i)),u)*u)
			do l=1,ns2grid
				p(l)=-.5*usu/(s2*s2grid(l))
			enddo
			p=p+pbase
			p=exp(p-maxval(p))
			do l=2,ns2grid
				p(l)=p(l)+p(l-1)
			enddo
			p=p/p(ns2grid)
			call rnun(unif)
			clubs2(i)=s2grid(count(unif(1)>p)+1)
		enddo
	end subroutine

	subroutine draw_CoCs2grid
		real	:: s2,u(0:q),p(ns2grid),usu,unif(1)
		integer	:: k,l
		real	:: pbase(ns2grid)

		do l=1,ns2grid
			pbase(l)=-.5*(q+1)*log(s2grid(l))+log(COCs2dist(l))
		enddo
!$omp parallel do private(k,u,s2,usu,p,unif,l)
		do k=1,nCoCs
			u=CoCFs(:,k)
			s2=om2
			usu=sum(matmul(Siu(:,:,cindCoC(k)),u)*u)
			do l=1,ns2grid
				p(l)=-.5*usu/(s2*s2grid(l))
			enddo
			p=p+pbase
			p=exp(p-maxval(p))
			do l=2,ns2grid
				p(l)=p(l)+p(l-1)
			enddo
			p=p/p(ns2grid)
			call rnun(unif)
			CoCs2(k)=s2grid(count(unif(1)>p)+1)
		enddo
	end subroutine
	
	
	subroutine draw_CF
		real	:: u(0:q),s2,Sigis(0:q,0:q,nclubs),ms(0:q,nclubs)
		integer	:: ir, i, k
		
!$omp parallel do private(k) 	
		do i=1,nclubs
			k=CoCid(i)
			Sigis(:,:,i)=Siu(:,:,cindCF(i))/((1-CoCrhos(i)**2)*om2*clubs2(i))
			ms(:,i)=matmul(Sigis(:,:,i),CoCrhos(i)*CoCFs(:,k))
		enddo
!$omp parallel do private(i,u,s2) reduction(+:Sigis,ms)
		do ir=1,n
			i=clubid(ir)
			u=Us(:,ir)
			u(0)=u(0)-muU
			s2=om2*(1-rhos(ir)**2)*kappa2(ir)
			Sigis(:,:,i)=Sigis(:,:,i)+rhos(ir)**2*SiU(:,:,cindu(ir))/s2
			ms(:,i)=ms(:,i)+rhos(ir)*matmul(SiU(:,:,cindu(ir)),u)/s2
		enddo
!$omp parallel do 
		do i=1,nclubs
			call linds(Sigis(:,:,i),Sigis(:,:,i))
			call rnnoa(CFs(:,i))
			CFs(:,i)=matmul(choleski(Sigis(:,:,i)),CFs(:,i))+matmul(Sigis(:,:,i),ms(:,i))
		enddo
	end subroutine			

	subroutine draw_CoCF
		real	:: u(0:q),s2,Sigis(0:q,0:q,nCoCs),ms(0:q,nCoCs)
		integer	:: i, k
		do k=1,nCoCs
			Sigis(:,:,k)=Siu(:,:,cindCoC(k))/(om2*CoCs2(k))
			ms(:,k)=0
		enddo
!$omp parallel do private(k,u,s2) reduction(+:Sigis,ms)
		do i=1,nclubs
			k=CoCid(i)
			u=CFs(:,i)
			s2=om2*(1-CoCrhos(i)**2)*clubs2(i)
			Sigis(:,:,k)=Sigis(:,:,k)+CoCrhos(i)**2*SiU(:,:,cindCF(i))/s2
			ms(:,k)=ms(:,k)+CoCrhos(i)*matmul(SiU(:,:,cindCF(i)),u)/s2
		enddo
!$omp parallel do 
		do k=1,nCoCs
			call linds(Sigis(:,:,k),Sigis(:,:,k))
			call rnnoa(CoCFs(:,k))
			CoCFs(:,k)=matmul(choleski(Sigis(:,:,k)),CoCFs(:,k))+matmul(Sigis(:,:,k),ms(:,k))
		enddo
	end subroutine			

	subroutine draw_sigma2_F
		real	:: ssum
		real	:: snu
		real	:: v(1),u(0:q),p(ns2fm),usu,unif(1)
		integer	:: l
		
		usu=sum(fm*matmul(Sifm(:,:,cindf),fm))
		do l=1,ns2fm
			p(l)=-.5*usu/s2fmtab(l)-.5*(q+1)*log(s2fmtab(l))
		enddo
		p=exp(p-maxval(p))
		do l=1,ns2fm
			p(l)=p(l)*merge(l,ns2fm+1-l,2*l<=ns2fm)
		enddo
		do l=2,ns2fm
			p(l)=p(l)+p(l-1)
		enddo
		p=p/p(ns2fm)
		call rnun(unif)
		sigma2_F(1)=s2fmtab(count(unif(1)>p)+1)
		
		u=f-fm
		u(0:1)=u(0:1)-mt_F
		ssum=sigma2_Faprior+sum(u*matmul(Sifa(:,:),u))
		snu=sigma2_Fanuprior+q+1
		call rnchi(snu,v)
		sigma2_F(2)=ssum/v(1)
	end subroutine
	
	subroutine draw_cindf
		real	:: p(nf),unif(1),u(0:q)
		integer	:: i
		
		u=fm
!		u(0:1)=u(0:1)-mt_F
!$omp parallel do 
		do i=1,nf
			p(i)=exp(-(0.5/sigma2_F(1))*sum(matmul(Sifm(:,:,i),u)*u))/sqrt(detSfm(i))
		enddo
		do i=2,nf
			p(i)=p(i-1)+p(i)
		enddo
		p=p/p(nf)
		call rnun(unif)
		cindf=count(unif(1)>p)+1
	end subroutine
	

	subroutine draw_om2
		real	:: ssum,snu,v(1),s2,u(0:q)
		integer	:: ir,i,k

		ssum=om2prior
		snu=om2nuprior
!$omp parallel do private(i,u,s2) reduction(+:ssum,snu)
		do ir=1,n
			i=clubid(ir)
			u=Us(:,ir)-rhos(ir)*CFs(:,i)
			u(0)=u(0)-muU
			s2=(1-rhos(ir)**2)*kappa2(ir)
			ssum=ssum+sum(matmul(Siu(:,:,cindu(ir)),u)*u)/s2
			snu=snu+q+1
		enddo
		
!$omp parallel do private(k,u,s2) reduction(+:ssum,snu)
		do i=1,nclubs			
			k=CoCid(i)
			u=CFs(:,i)-CoCrhos(i)*CoCFs(:,k)
			s2=clubs2(i)*(1-CoCrhos(i)**2)
			ssum=ssum+sum(matmul(Siu(:,:,cindCF(i)),u)*u)/s2
			snu=snu+q+1
		enddo
!$omp parallel do private(u,s2) reduction(+:ssum,snu)
		do k=1,nCoCs			
			u=CoCFs(:,k)
			s2=CoCs2(k)
			ssum=ssum+sum(matmul(Siu(:,:,cindCoC(k)),u)*u)/s2
			snu=snu+q+1
		enddo
		call rnchi(snu,v)
		om2=ssum/v(1)
	end subroutine
	
	subroutine draw_dist(nmeas,mid,dist)
		use rngam_int
		integer	:: nmeas,mid(:)
		real	:: dist(:)

		real	:: alpha(size(dist))
		integer	:: i
		
		alpha=alphatot_prior/size(dist)
		do i=1,nmeas
			alpha(mid(i))=alpha(mid(i))+1
		enddo
!$omp parallel do 
		do i=1,size(dist)
			call rngam(alpha(i),dist(i:i))
		enddo
		dist=dist/sum(dist)
	end subroutine

	subroutine draw_thdist
		call draw_dist(n,cindu,thdist)
	end subroutine

	subroutine draw_CFthdist
		call draw_dist(nclubs,cindCF,CFthdist)
	end subroutine

	subroutine draw_CoCthdist
		call draw_dist(nCoCs,cindCoC,CoCthdist)
	end subroutine

	subroutine draw_rhodist
		use rngam_int
		real	:: alpha(nrhogrid)
		integer	:: ir,j
				
		alpha=alphatot_prior/size(alpha)
		do ir=1,n
			j=nint((nrhogrid-1)*rhos(ir)/maxval(rhogrid))+1
			alpha(j)=alpha(j)+1
		enddo
!$omp parallel do 
		do j=1,nrhogrid
			call rngam(alpha(j),rhodist(j:j))
		enddo
		rhodist=rhodist/sum(rhodist)
	end subroutine
	
	subroutine draw_CoCrhodist
		use rngam_int
		real	:: alpha(nrhogrid)
		integer	:: i,j
				
		alpha=alphatot_prior/size(alpha)
		do i=1,nclubs
			j=nint((nrhogrid-1)*CoCrhos(i)/maxval(rhogrid))+1
			alpha(j)=alpha(j)+1
		enddo
!$omp parallel do 
		do j=1,nrhogrid
			call rngam(alpha(j),CoCrhodist(j:j))
		enddo
		CoCrhodist=CoCrhodist/sum(CoCrhodist)
	end subroutine

	subroutine draw_gens2dist(gens2dist,gens2)
		use rngam_int
		real	:: gens2dist(ns2grid), gens2(:)
		real	:: alpha(ns2grid)
		integer	:: j,ind
				
		alpha=alphatot_prior/size(alpha)
		do j=1,size(gens2)
			ind=count(gens2(j)>=s2grid)
			alpha(ind)=alpha(ind)+1
		enddo
		do ind=1,ns2grid
			call rngam(alpha(ind),gens2dist(ind:ind))
		enddo
		gens2dist=gens2dist/sum(gens2dist)
	end subroutine
	
	subroutine draw_s2dist
		call draw_gens2dist(s2dist,kappa2)
	end subroutine
	subroutine draw_CFs2dist
		call draw_gens2dist(CFs2dist,clubs2)
	end subroutine
	subroutine draw_CoCs2dist
		call draw_gens2dist(CoCs2dist,CoCs2)
	end subroutine
	
	subroutine draw_fcsts
		integer	:: ir,i,k
		real	:: fcstU(nh),u(0:q),s2
		
		call rnnoa(fcstf)
		fcstf=matmul(Xfcstf,mt_F)+matmul(mfcstfm(:,:,cindf),fm)+sqrt(sigma2_F(1))*matmul(cholfcstfm(:,:,cindf),fcstf)
		call rnnoa(fcstU)
		u=f-fm
		u(0:1)=u(0:1)-mt_F
		fcstf=fcstf+matmul(mfcstfa,u)+sqrt(sigma2_f(2))*matmul(cholfcstfa,fcstu)
!$omp parallel do private(fcstU,s2)
		do k=1,nCoCs
			call rnnoa(fcstU)
			s2=om2*CoCs2(k)
			fcstCoCFs(:,k)=matmul(mfcstu(:,:,cindCoC(k)),CoCFs(:,k))+sqrt(s2)*matmul(cholfcstu(:,:,cindCoC(k)),fcstU)
		enddo
!$omp parallel do private(k,fcstU,s2,u)
		do i=1,nclubs
			k=CoCid(i)
			call rnnoa(fcstU)
			s2=om2*(1-CoCrhos(i)**2)*clubs2(i)
			u=CFs(:,i)-CoCrhos(i)*CoCFs(:,k)
			fcstCFs(:,i)=CoCrhos(i)*fcstCoCFs(:,k)+matmul(mfcstu(:,:,cindCF(i)),u)+sqrt(s2)*matmul(cholfcstu(:,:,cindCF(i)),fcstU)
		enddo
		
!$omp parallel do private(i,fcstU,u,s2) 
		do ir=1,n
			i=clubid(ir)
			call rnnoa(fcstU)
			u=Us(:,ir)-rhos(ir)*CFs(:,i)
			u(0)=u(0)-muU
			s2=(1-rhos(ir)**2)*kappa2(ir)*om2
			fcstU=matmul(mfcstu(:,:,cindu(ir)),u)+sqrt(s2)*matmul(cholfcstu(:,:,cindu(ir)),fcstU)
			fcstU=fcstU+rhos(ir)*fcstCFs(:,i)
			fcstUs(:,ir)=fcstU
			fcsts(:,ir)=fcstf+fcstUs(:,ir)
		enddo
	end subroutine
	
	subroutine draw_F
		real	:: m(0:q), u(0:q), Sigi(0:q,0:q), s2, fhat(0:q)
		integer	:: ir,i
		m=0; fhat=0
		Sigi=sigma2_f(1)*Sfm(:,:,cindf)+sigma2_f(2)*Sfa
		call linds(Sigi,Sigi)
!$omp parallel do private(i,s2,u) reduction(+:m,Sigi,fhat)
		do ir=1,n
			i=clubid(ir)
			s2=(1-rhos(ir)**2)*kappa2(ir)*om2
			u=Xs(:,ir)-rhos(ir)*CFs(:,i)
			u(0)=u(0)-muU
			u(0:1)=u(0:1)-mt_F
			m=m+matmul(Siu(:,:,cindu(ir)),u)/s2
			Sigi=Sigi+Siu(:,:,cindu(ir))/s2
			if(fweights(ir)>0.0) then
				fhat=fhat+fweights(ir)*Xs(:,ir)
			endif			
		enddo
		Sigi=Sigi+Deltainv
		fhat(0:1)=fhat(0:1)-mt_F
		m=m+matmul(Deltainv,fhat-Ydummy)	
		call linds(Sigi,Sigi)
		call rnnoa(f)
		f=matmul(Sigi,m)+matmul(choleski(Sigi),f)
		f(0:1)=f(0:1)+mt_F
		do ir=1,n
			Us(:,ir)=Xs(:,ir)-F
		enddo
	end subroutine
	
	subroutine draw_mt_F
		real	:: m(2),prec(2,2),Sigi(0:q,0:q)
		Sigi=sigma2_f(1)*Sfm(:,:,cindf)+sigma2_f(2)*Sfa
		call linds(Sigi,Sigi)
		prec=muprec_prior*eye(2)
		prec=prec+Sigi(0:1,0:1)
		m=matmul(Sigi(0:1,:),f)
		call linds(prec,prec)
		call rnnoa(mt_F)
		mt_F=matmul(prec,m)+matmul(choleski(prec),mt_f)
	end subroutine
		
	subroutine draw_fm
		real	:: u(0:q),S(0:q,0:q),Sigi(0:q,0:q),mfm(0:q,0:q)
		S=sigma2_f(1)*Sfm(:,:,cindf)
		Sigi=S+sigma2_f(2)*Sfa
		call linds(Sigi,Sigi)
		u=f
		u(0:1)=u(0:1)-mt_F
		mfm=matmul(S,Sigi)
		call rnnoa(fm)
		fm=matmul(mfm,u)+matmul(choleski(S-matmul(mfm,S)),fm)
	end subroutine

	subroutine draw_X
		integer	:: ir
		real	:: m(0:q),s2
		integer	:: i
		real	:: fhat(0:q),sVs(0:q,0:q),e(0:q),u(0:q)
		
		fhat=0;sVs=0
!$omp parallel do private(i,s2,m,u) reduction(+:fhat,sVs)
		do ir=1,n
			i=clubid(ir)
			s2=(1-rhos(ir)**2)*kappa2(ir)*om2
			m=rhos(ir)*CFs(:,i)
			m(0)=m(0)+muU
			call rnnoa(u)
			u=sqrt(s2)*matmul(cholSu(:,:,cindu(ir)),u)+m
			Us(:,ir)=u-matmul(SuAA(:,:,cindu(ir),ir),u-Us(:,ir))
			if(fweights(ir)>0.0) then
				fhat=fhat+fweights(ir)*Us(:,ir)
				sVs=sVs+fweights(ir)**2*s2*SuAAS(:,:,cindu(ir),ir)
			endif
		enddo
		call rnnoa(e)
		e=Ydummy+matmul(choleski(Delta),e)
		fhat=matmul(invertpd(sVs+Delta),fhat-e)
		do ir=1,n
			if(fweights(ir)>0.0) then
				s2=(1-rhos(ir)**2)*kappa2(ir)*om2
				Us(:,ir)=Us(:,ir)-fweights(ir)*s2*matmul(SuAAS(:,:,cindu(ir),ir),fhat)
			endif
			Xs(:,ir)=Us(:,ir)+F
		enddo
	end subroutine

	subroutine draw_Y
		integer	:: ir
		real	:: m(0:q),s2
		integer	:: i
		real	:: fhat(0:q),sVs(0:q,0:q),e(0:q),u(0:q)
		
		fhat=0;sVs=0
!$omp parallel do private(i,s2,m,u) reduction(+:fhat,sVs)
		do ir=1,n
			i=clubid(ir)
			s2=(1-rhos(ir)**2)*kappa2(ir)*om2
			m=rhos(ir)*CFs(:,i)
			m(0)=m(0)+muU
			call rnnoa(u)
			u=sqrt(s2)*matmul(cholSu(:,:,cindu(ir)),u)+m
			Us(:,ir)=u-matmul(SuBB(:,:,cindu(ir),ir),u-Us(:,ir))
			if(fweights(ir)>0.0) then
				fhat=fhat+fweights(ir)*Us(:,ir)
				sVs=sVs+fweights(ir)**2*s2*SuBBS(:,:,cindu(ir),ir)
			endif
			Xs(:,ir)=Us(:,ir)+f
		enddo
		call rnnoa(e)
		Ydummy=matmul(choleski(Delta),e)+fhat
	end subroutine
end module