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