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