subroutine sfast2(n,t,var,tcorr,ldets,x,y,mode)
c
c***	Solves the linear system S*x=y, where S is a matrix with elements
c***
c***		S(i,j)=var*exp(-abs(t(i)-t(j))/tcorr),  i,j=1,...,n,
c***
c***    The array t(i), i=1,...,n, are given "times", no two of which can
c***	be identical, since S is then a singular matrix.
c***
c***    Given a right hand side vector, y(i), i=1,...,n, the unknown
c***    vector x(i), i=1,...,n, is found, using the fast inversion
c***	method of Rybicki and Press (1995, Phys. Rev. Lett., 74, 1060).
c***	with O(n) timings.  Also calculated is ldets, the logarithm of 
c***	the determinant of S.
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,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***	If the array t(i) is not given in ascending order,
c***
c***		t(1) < t(2) < ..... < t(n)
c***
c***	a new sorted array ts(i) is internally generated (necessary to apply
c***	the Rybicki and Press method).  Sorting takes O(n log n) steps,
c***	so it makes the overall method less "fast".  If the times
c***	are out of order, but in some simple way that allows a fast
c***	preliminary sorting by the user, it may be more efficient to do
c***	that sorting before using this routine.  In that case, one might
c***	also consider using the simpler routine sfast1.
c***
c***    G. Rybicki.  Version 24 April 1998.
c***
c***    USES: hpsort_indx,trimult   (in package subs.f)
c
	INTEGER i,n,n1,mode,NDIM
	PARAMETER (NDIM=1000)	!alter to suit problem
	DOUBLE PRECISION t(n),var,tcorr,ldets,x(n),y(n)
	DOUBLE PRECISION r(NDIM),e(NDIM),a(NDIM),b(NDIM),c(NDIM)
	INTEGER indx(NDIM)
	DOUBLE PRECISION ts(NDIM),xs(NDIM),ys(NDIM)
	DOUBLE PRECISION varinv,ldet,ZERO,ONE
	PARAMETER (ZERO=0.d0,ONE=1.d0)
	LOGICAL sorted
	SAVE a,b,c,ldet,indx
c
	if(mode.eq.0)then
		if(NDIM.lt.n)then
			stop 'sfast2: NDIM too small'
		endif
		n1=n-1
c
c***		Copy t(i) into ts(i), and set up index array indx.
c
                do i=1,n
                        ts(i)=t(i)
			indx(i)=i
                enddo
c
c***		Check if times are already sorted.
c
		sorted=.true.
		do i=1,n1
			if(ts(i+1).lt.ts(i))then
				sorted=.false.
			endif
		enddo
c
c***		If not, then sort array ts(i), creating the index array
c***		indx relating sorted and unsorted quantities.
c
		if(.not.sorted)then
	                call hpsort_indx(n,ts,indx)
		endif
	        ldet=log(var)
		do i=1,n1
			r(i)=exp(-(ts(i+1)-ts(i))/tcorr)
			e(i)=r(i)/(ONE-r(i)**2)
	                ldet=ldet+log(var*(ONE-r(i)**2))
		enddo
c
		varinv=ONE/var
		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))
	endif
c
c***    Use the quantities determined above to solve the linear system Sx=y
c***    for x, given y.  Note that the calculation involves simple
c***    tridiagonal matrices only when expressed in terms of the sorted
c***    versions xs and ys.
c
        do i=1,n
                ys(i)=y(indx(i))
        enddo
c
	call trimult(a,b,c,ys,xs,n)
c
        do i=1,n
                x(indx(i))=xs(i)
        enddo
	ldets=ldet
	return
	end