subroutine snfast22(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 does not assume that array t(i) is sorted in ascending
c*** order, and also allows for possible duplicate times t(i)
c*** that 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: hpsort_indx,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 ts(NDIM),err2s(NDIM),xs(NDIM),ys(NDIM)
DOUBLE PRECISION tt(NDIM),ee(NDIM),xx(NDIM),yy(NDIM)
INTEGER imap(NDIM+1),indx(NDIM)
DOUBLE PRECISION sum1,sum2,sum3,ZERO,ONE,ldetc,ldetc1,ldetm
PARAMETER (ZERO=0.0d0,ONE=1.0d0)
SAVE indx,tt,ee,err2s,nn,imap,ldetm
c
if(NDIM.lt.n)then
stop 'snfast22: NDIM too small'
endif
c
if(mode.eq.0)then
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***
c*** ***Note*** This sort takes of order n*log(n) steps in
c*** general, so if the times are already sorted, skip this
c*** step or use the routine snfast21.f instead. Also, if
c*** the times are unsorted, but in some simple way, it may
c*** be better to replace the general hpsort_indx routine
c*** with one that takes advantage of those simplifications.
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*** 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)=ts(1)
m=1
imap(1)=1
do i=2,n
if(ts(i)-tlimit.lt.tt(ii))then
m=m+1
else
imap(ii)=i-m
ii=ii+1
if(ii.gt.NDIM)then
stop 'snfast22: NDIM too small'
endif
tt(ii)=ts(i)
m=1
endif
enddo
nn=ii
if(nn.gt.NDIM)then
stop 'snfast22: NDIM too small'
endif
imap(nn)=n-m+1
imap(nn+1)=n+1
tt(nn)=ts(n)
c
ldetm=ZERO
do ii=1,nn
i1=imap(ii)
i2=imap(ii+1)-1
if(i1.eq.i2)then
ee(ii)=err2s(i1)
else
sum1=ZERO
sum2=ZERO
do i=i1,i2
sum1=sum1+ONE/err2s(i)
sum2=sum2+log(err2s(i))
enddo
ee(ii)=ONE/sum1
ldetm=ldetm+log(sum1)+sum2
endif
enddo
endif
c
do i=1,n
ys(i)=y(indx(i))
enddo
c
do ii=1,nn
i1=imap(ii)
i2=imap(ii+1)-1
if(i1.eq.i2)then
yy(ii)=ys(i1)
else
sum3=ZERO
do i=i1,i2
sum3=sum3+ys(i)/err2s(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
xs(i1)=xx(ii)
else
do i=i1,i2
xs(i)=(xx(ii)*ee(ii)+ys(i)-yy(ii))/err2s(i)
enddo
endif
enddo
c
do i=1,n
x(indx(i))=xs(i)
enddo
c
return
end