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