subroutine snfast12(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***	Here the array t(i) is not assumed to be sorted in ascending order,
c***    but that there are no duplicate times, since this makes the matrix
c***    S singular.  For other possibilities, it may be more efficient to
c***    use an alternative routine based on 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: hpsort_indx,trimult,tridiag,tridet  (in package subs.f)
c
	INTEGER i,n,n1,NDIM,mode
	PARAMETER (NDIM=7200)
	INTEGER indx(NDIM)
	DOUBLE PRECISION t(n),err2(n),var,tcorr,ldetc,x(n),y(n)
	DOUBLE PRECISION ts(NDIM),err2s(NDIM),xs(NDIM),ys(NDIM),zs(NDIM)
	DOUBLE PRECISION r(NDIM),e(NDIM),ldetmat,ldets,ldetm
	DOUBLE PRECISION a(NDIM),b(NDIM),c(NDIM),aa(NDIM),bb(NDIM),cc(NDIM)
	DOUBLE PRECISION exp,log
	DOUBLE PRECISION varinv,ZERO,ONE,dum
	PARAMETER (ZERO=0.d0,ONE=1.d0)
	SAVE a,b,c,aa,bb,cc,indx,ldetmat
c
	if(n.eq.1)then
		x(1)=y(1)/(var+err2(1))
		ldetc=log(var+err2(1))
		return
	endif
	if(mode.eq.0)then
		if(NDIM.lt.n)then
			stop 'snfast12: NDIM too small'
		endif
		n1=n-1
		varinv=ONE/var
c
c***		Create sorted versions of t and err2 (ts and err2s), using
c***		t as the key.  Index array indx can be similarly
c***		used for creating sorted version of other arrays.
c
		do i=1,n
			ts(i)=t(i)
		enddo
		call hpsort_indx(n,ts,indx)
		do i=1,n
			err2s(i)=err2(indx(i))
		enddo
c
c***		Calculate log-determinant of S using special formulas
c***		for Ornstein-Unlenbeck process.
c
		ldets=log(var)
		do i=1,n1
			r(i)=exp(-(ts(i+1)-ts(i))/tcorr)
			e(i)=r(i)/(ONE-r(i)**2)
			ldets=ldets+log(var*(ONE-r(i)**2))
		enddo
c
c***		Calculate the arrays a, b, and c, which define the 
c***		tridiagonal matrix inverse S^{-1}.
c
		b(1)= varinv*(ONE+r(1)*e(1))
		c(1)=-varinv*e(1)
		do i=2,n1
			a(i)=-varinv*e(i-1)
			b(i)= varinv*(ONE+r(i)*e(i)+r(i-1)*e(i-1))
			c(i)=-varinv*e(i)
		enddo
		a(n)=-varinv*e(n1)
		b(n)= varinv*(ONE+r(n1)*e(n1))
c
c***		Calcutate the arrays aa, bb, and cc, which define the
c***		tridiagonal matrix (1+N*S^{-1})
c
		bb(1)=ONE+b(1)*err2s(1)
		cc(1)=c(1)*err2s(1)
		do i=2,n1
			aa(i)=a(i)*err2s(i)
			bb(i)=ONE+b(i)*err2s(i)
			cc(i)=c(i)*err2s(i)
		enddo
		aa(n)=a(n)*err2s(n)
		bb(n)=ONE+b(n)*err2s(n)
c
c***		Calculate log-determinant of (1+N*S^{-1}).  Finally 
c***		get ldetmat, the log-determinant of the full correlation
c***		matrix C=S+N=(1+N*S^{-1})*S
c
		call tridet(aa,bb,cc,n,ldetm,dum)
		ldetmat=ldets+ldetm
c
	endif
c
c***	Use the quantities determined above to solve system (S+N)x=y
c***	for x, given y.  Note that the calculation involves only simple
c***	tridiagonal matrices when expressed in terms of the sorted
c***	versions xs and ys.
c
	do i=1,n
		ys(i)=y(indx(i))
	enddo
	call tridiag(aa,bb,cc,zs,ys,n)
	call trimult(a,b,c,zs,xs,n)
	do i=1,n
		x(indx(i))=xs(i)
	enddo
	ldetc=ldetmat
	return
	end