anchor.diff:  ! DAXPY(N,DA,DX,INCX,DY,INCY) Adding vectors
         EXTERNAL DCOPY
   ! DCOPY(N,DX,INCX,DY,INCY) copies x to y
  -      EXTERNAL DDOT
  +! NB DDOT from ATLAS, compiled with gfortran 4.2 on Ubuntu Gutsy,
  +! was producing bad output- hence the manual dot product.
         
   !     x <- (x-mu)      
         call DAXPY(n, -1.0D0, mu, 1, x, 1)
  
positive.code:       PARAMETER (infinity = 1.7976931348623157d308)      
        DOUBLE PRECISION PI
        PARAMETER (PI=3.141592653589793238462643d0) 
        DOUBLE PRECISION twopi_N, log_detC
  
        EXTERNAL DPOTRF
  ! DPOTRF( UPLO, N, A, LDA, INFO ) Cholesky factorization
        EXTERNAL DPOTRS
  ! DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO ) Solves triangular system
        EXTERNAL DAXPY
  ! DAXPY(N,DA,DX,INCX,DY,INCY) Adding vectors
        EXTERNAL DCOPY
  ! DCOPY(N,DX,INCX,DY,INCY) copies x to y
        EXTERNAL DDOT
        
  !     C <- cholesky(C)      
        call DPOTRF( 'L', n, C, n, info )
  
negative.code: !     mu <- tau * x
        call DSYMV('L',n,1.0D0,tau,n,x,1,0.0D0,mu,1)
  
        gd = 0.0D0
        do i=1,n
            gd=gd+x(i)*mu(i)
        end do
  
  !     like <- -.5 dot(x,mu) (.5 (x-mu) C^{-1} (x-mu)^T)
        like = -0.5D0 * gd
  
  !      Cholesky factorize tau for the determinant.      
         call DPOTRF( 'L', n, tau, n, info )
  