subroutine snfast21(n,t,err2,var,tcorr,ldetc,x,y,mode)
c
c***    Solves the linear system (S+N)*x=y, where S is a (signal)
c***    autocorrelation matrix with elements of the form,
c***
c***        S(i,j)=var*exp(-abs(t(i)-t(j))/tcorr),  i,j=1,...,n
c***
c***    where t(i), i=1,...,n, are an given array of "times", and
c***    where N is a diagonal (noise) matrix with elements
c***
c***		N(i,j)=delta(i,j)*err2(i)
c***
c***	Given a right hand side vector, y(i), i=1,...,n, the unknown
c***	vector x(i), i=1,...,n, is found.  Also calculated is ldetc,
c***	the logarithm of the determinant of C=S+N.  
c***
c***    The integer parameter "mode" must be set to 0 for the first call
c***    to the routine for a given set of input parameters n,t,err2,var,tcorr.
c***    Subsequent calls with these same parameters, but with different
c***    RHSs y(i), can be done with mode.ne.0, which avoids unnecessary
c***    computations, increasing the speed.
c***
c***	The fast method of Rybicki and Press (1995, Phys. Rev. Lett., 
c***	74, 1060; cf.eq. [17]) is used.  
c***
c***    This routine assumes that array t(i) is sorted in ascending 
c***    order, but allows for possible duplicate times t(i) that 
c***	are equal to within a fraction TINY of the correlation time 
c***	tcorr.  This is done by transforming to a new equivalent problem
c***    (with smaller n) in which duplicate times are eliminated.  
c***	The resulting problem is then solved using the routine snfast11.  
c***	For maximum speed one can use an alternative routine based on 
c***	the following table:
c***
c***                     times sorted       times unsorted
c***                   _____________________________________
c***                  |                 |                   |
c***    no duplicate  |    snfast11     |    snfast12       |
c***      times       |                 |                   |
c***                  |     O(n)        |   O(n log n)      |
c***                  |                 |                   |
c***                   _____________________________________
c***                  |                 |                   |
c***    duplicate     |  snfast21       |    snfast22       |
c***      times       |                 |                   |
c***                  |     O(n)        |   O(n log n)      |
c***                  |                 |                   |
c***                   _____________________________________
c***
c***
c***	The routine snfast22 is the most general, and thus could be used
c***	for all problems.  However, choosing the least complicated routine
c***	suitable for a given problem (i.e., the one with the most "1"s in
c***	its name) will minimize CPU time and storage requirements.  These
c***	are all "fast" routines, but the timing requirements of snfast11
c***	and snfast21 are O(n), while those of snfast12 and snfast22 are
c***	only O(n log n) due to the sorting.  If the times are unsorted,
c***	but in some simple way that can be sorted in less than O(n log n),
c***	one should consider the possibility of doing that simple sort
c***	first, then using snfast11 or snfast21.
c***
c***	G. Rybicki.  Version 24 April 1998.
c***
c***	USES: snfast11
c***	USES: trimult,tridiag,tridet   (in package subs.f)
c
	INTEGER i,ii,i1,i2,m,n,nn,NDIM,mode
	DOUBLE PRECISION t(*),err2(*),x(*),y(*),var,tcorr,TINY,tlimit
	PARAMETER (NDIM=10000)	!alter to suit problem
	PARAMETER (TINY=1.0d-7)	!alter to suit problem
	DOUBLE PRECISION tt(NDIM),ee(NDIM),xx(NDIM),yy(NDIM)
	INTEGER imap(NDIM+1)
	DOUBLE PRECISION sum1,sum2,sum3,ZERO,ONE,ldetc,ldetc1,ldetm
	PARAMETER (ZERO=0.0d0,ONE=1.0d0)
	SAVE tt,ee,nn,imap,ldetm
c
	if(mode.eq.0)then
c
c***		Create new array tt(ii) that omits duplicate times
c***		(to within fraction TINY of the correlation time).
c***		Also create the array imap(ii), which defines the
c***		relationship between the original and modified problems.
c
		tlimit=TINY*tcorr
c	
		ii=1
		tt(1)=t(1)
		m=1
		imap(1)=1
		do i=2,n
			if(t(i).lt.t(i-1))then
				stop 'snfast21: times not sorted'
			elseif(t(i)-tlimit.lt.tt(ii))then
				m=m+1
			else
				imap(ii)=i-m
				ii=ii+1
				if(ii.gt.NDIM)then
					stop 'snfast21: NDIM too small'
				endif
				tt(ii)=t(i)
				m=1
			endif
		enddo
		nn=ii
		if(nn.gt.NDIM)then
			stop 'snfast21: NDIM too small'
		endif
		imap(nn)=n-m+1
		imap(nn+1)=n+1
		tt(nn)=t(n)
c
		ldetm=ZERO
		do ii=1,nn
			i1=imap(ii)
			i2=imap(ii+1)-1
			if(i1.eq.i2)then
				ee(ii)=err2(i1)
			else
				sum1=ZERO
				sum2=ZERO
				do i=i1,i2
					sum1=sum1+ONE/err2(i)
					sum2=sum2+log(err2(i))
				enddo
				ee(ii)=ONE/sum1
				ldetm=ldetm+log(sum1)+sum2
			endif
		enddo
	endif
c
	do ii=1,nn
		i1=imap(ii)
		i2=imap(ii+1)-1
		if(i1.eq.i2)then
			yy(ii)=y(i1)
		else
			sum3=ZERO
			do i=i1,i2
				sum3=sum3+y(i)/err2(i)
			enddo
			yy(ii)=sum3*ee(ii)
		endif
	enddo
c
	call snfast11(nn,tt,ee,var,tcorr,ldetc1,xx,yy,mode)
	ldetc=ldetc1+ldetm
c
	do ii=1,nn
		i1=imap(ii)
		i2=imap(ii+1)-1
		if(i1.eq.i2)then
			x(i1)=xx(ii)
		else
			do i=i1,i2
				x(i)=(xx(ii)*ee(ii)+y(i)-yy(ii))/err2(i)
			enddo
		endif
	enddo
c
	return
	end