The routines in this directory implement some of the fast inversion
methods of Rybicki and Press (1995, Phys. Rev. Lett., 74, 1060)
involving matrices of the exponential form

        S(i,j)=var*exp(-abs(t(i)-t(j))/tcorr),  i,j=1,...,n

where t(i), i=1,...,n, are an given array of "times".  Such matrices
appear as autocorrelation functions in "Ornstein-Uhlenbeck" random
processes, where "var" is the variance of the process and "tcorr"
is the correlation time.

The routines

	sfast1(n,t,var,tcorr,ldets,x,y,mode)
	sfast2(n,t,var,tcorr,ldets,x,y,mode)
      
find the solution of the linear system Sx=y, where the vector y(i),
i=1,...,n, is given and the vector x(i), i=1,...,n, is to be found.
Also found is dets, the logarithm of the determinant of S.  

The times t(i) must be distinct; otherwise the matrix S is singular.
The routine sfast1 requires the given array t(i) to be sorted in
ascending order,

               t(1) < t(2) < ... < t(n).  

The routine sfast2 does not require the array to satisfy this
requirement, since it does the necessary sorting internally.

The integer parameter "mode" must be set to 0 for the first call to the
routine for a given set of input parameters n,t,var,tcorr.
Subsequent calls with these same parameters, but with different RHSs
y(i), can be done with mode.ne.0, which avoids unnecessary
computations, increasing the speed.

The routines

	snfast11(n,t,err2,var,tcorr,ldetc,x,y,mode)
	snfast12(n,t,err2,var,tcorr,ldetc,x,y,mode)
	snfast21(n,t,err2,var,tcorr,ldetc,x,y,mode)
	snfast22(n,t,err2,var,tcorr,ldetc,x,y,mode)

solve the linear system (S+N)*x=y, where S is a (signal) autocorrelation
matrix with elements of the form,

        S(i,j)=var*exp(-abs(t(i)-t(j))/tcorr),  i,j=1,...,n

where t(i), i=1,...,n, are an given array of "times", and where N is a
diagonal (noise) matrix with elements

        N(i,j)=delta(i,j)*err2(i)

where the noise variances err2(i) are positive.  Given a right hand
side vector, y(i), i=1,...,n, the unknown vector x(i), i=1,...,n, is
found.  Also calculated is ldetc, the logarithm of the determinant of
C=S+N.  

The integer parameter "mode" must be set to 0 for the first call to the
routine for a given set of input parameters n,t,err2,var,tcorr.
Subsequent calls with these same parameters, but with different RHSs
y(i), can be done with mode.ne.0, which avoids unnecessary
computations, increasing the speed.

The routines differ in whether they assume that array t(i) of times is
sorted in ascending order, or whether there are any duplicate times.
Duplicate times do not make the matrix S+N singular, but a 
straightforward application of the method of Rybicki and Press
requires the inversion of S, which is singular for duplicate times.
Thus the case of duplicate times requires some extra manipulation of
the equations.

The choice of routine should be based on the following table:

                 times sorted       times unsorted
               _____________________________________
              |                 |                   |
no duplicate  |    snfast11     |    snfast12       |
  times       |                 |                   |
              |     O(n)        |   O(n log n)      |
              |                 |                   |
               _____________________________________
              |                 |                   |
duplicate     |  snfast21       |    snfast22       |
  times       |                 |                   |
              |     O(n)        |   O(n log n)      |
              |                 |                   |
               _____________________________________

The routine snfast22 is the most general, and thus could be used for
all problems.  However, choosing the least complicated routine
suitable for a given problem (i.e., the one with the most "1"s in its
name) will minimize CPU time and storage requirements.  These are all
"fast" routines, but the timing requirements of snfast11 and snfast21
are O(n), while those of snfast12 and snfast22 are only O(n log n) due
to the sorting.  If the times are unsorted, but in some simple way
that can be sorted in less than O(n log n), one should consider the
possibility of doing that simple sort first, then using snfast11 or
snfast21.