subroutine snfast11(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 order, c*** and that there are no duplicate times, since this makes the matrix c*** S singular. If either or both of these two assumptions are violated, 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: trimult,tridiag,tridet (in package subs.f) c INTEGER i,n,n1,NDIM,mode PARAMETER (NDIM=10000) !alter to suit problem DOUBLE PRECISION t(n),err2(n),var,tcorr,ldetc,x(n),y(n) DOUBLE PRECISION r(NDIM),e(NDIM),z(NDIM),ldetmat,ldets,ldetm DOUBLE PRECISION a(NDIM),b(NDIM),c(NDIM),aa(NDIM),bb(NDIM),cc(NDIM) DOUBLE PRECISION varinv,ZERO,ONE,dum PARAMETER (ZERO=0.d0,ONE=1.d0) SAVE a,b,c,aa,bb,cc,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 'snfast11: NDIM too small' endif n1=n-1 c c*** Check if times are properly sorted. c do i=1,n1 if(t(i+1).lt.t(i))then stop 'snfast11: Array t(i) not properly sorted' endif enddo varinv=ONE/var 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(-(t(i+1)-t(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)*err2(1) cc(1)=c(1)*err2(1) do i=2,n1 aa(i)=a(i)*err2(i) bb(i)=ONE+b(i)*err2(i) cc(i)=c(i)*err2(i) enddo aa(n)=a(n)*err2(n) bb(n)=ONE+b(n)*err2(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. c call tridiag(aa,bb,cc,z,y,n) call trimult(a,b,c,z,x,n) ldetc=ldetmat return end