subroutine sfast1(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 singular. This routine requires c*** that the given array t(i) be sorted in ascending order, c*** c*** t(1) < t(2) < ..... < t(n) c*** c*** If this is not the case, then the user must perform a preliminary c*** sort, or use the alternate routine sfast2. 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*** G. Rybicki. Version 24 April 1998. c*** c*** USES: trimult (in package subs.f) c*** 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) DOUBLE PRECISION varinv,ldet,ZERO,ONE PARAMETER (ZERO=0.d0,ONE=1.d0) SAVE a,b,c,ldet c if(mode.eq.0)then if(NDIM.lt.n)then stop 'sfast1: NDIM too small' endif n1=n-1 c c*** Check if times are properly sorted. c do i=1,n1 if(t(i+1).le.t(i))then stop 'sfast1: Array t(i) not properly sorted' endif enddo ldet=log(var) do i=1,n1 r(i)=exp(-(t(i+1)-t(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 call trimult(a,b,c,y,x,n) ldets=ldet c return end