Update README.md. Corrected covariance computation.

- README.md contains the reference [arXiv:1809.01289]
- Corrected bug in te computations of the covariance for the case
  where there is no Monte Carlo data.

README.md

@@ -262,8 +262,8 @@ work whith these versions.
 1. Download or clone the repository.
 1. Edit the `Makefile` in the `build` directory. Change the compiler
    command/options (variables `FC` and `FOPT`). 
-1. Compile the library with `gmake`.
-1. Optionally build/run the test codes with `gmake test`. Executabes will
+1. Compile the library with `GNU make`.
+1. Optionally build/run the test codes with `make test`. Executabes will
    be placed in the `test` directory.
 1. If preferred, move the contents of the `include` and `lib`
    directories somewhere else.
@@ -280,4 +280,8 @@ Look into the `doc/aderrors.pdf` file.
 ## How to cite
 
 If you use this package for a scientific publication, please cite the
-original work ...
+original work:
+
+"Automatic differentiation for error analysis of Monte Carlo data"
+Alberto Ramos. 
+[arXiv:1809.01289](https://arxiv.org/abs/1809.01289)

src/aderrors/aderrors.f90

@@ -838,8 +838,8 @@ contains
     
     if (imode.eq.UW_CF_MODE) then
        call p%init(1, [1], [1], [1])
-       p%mean = 0.0_DP
-       p%data = d
+       p%mean = d
+       p%data = 0.0_DP
        p%ivrep(1) = p%nd(1)
     else
        call uwfree(p)

src/aderrors/aderrors_cov.f90

@@ -45,19 +45,22 @@ contains
     type (uwreal),  intent (in)    :: p(:)
     real (kind=DP), intent (inout) :: mat(:,:)
 
-    integer :: n, io1, io2, irst, irnd
+    integer :: n, io1, io2, irst, irnd, is1, is2, ie1, ie2
     
     call autocorrfunc_cov(p)
     call wopt()
     
-    
     mat = 0.0_DP
     do io1 = 1, wsid%no
        do io2 = io1, wsid%no
           do n = 1, wsid%nid
              if (wsid%nd(n).eq.1) then
-                irst     = 1 + sum(wsid %nd(1:n-1)) 
-!                mat(io1,io2) = wscf%data(irst)**2
+                if ( (wsid%map(n,io1).ne.0).and.(wsid%map(n,io2).ne.0) ) then
+                   call p(io1)%get_offset(is1,ie1,wsid%map(n,io1),0)
+                   call p(io2)%get_offset(is2,ie2,wsid%map(n,io2),0)
+                   mat(io1,io2) = mat(io1,io2) + &
+                        p(io1)%data(is1)*p(io2)%data(is2)
+                end if
              else
                 irst = sum(wsid%tmax(1:n-1))
                 irnd = irst + wsid%iw(n)
@@ -290,6 +293,7 @@ contains
 
     maxrl = sum(wsid%tmax(1:nid))
     call init_wscf(0, no, maxrl, n1, n2)
+    if (maxrl.eq.0) return
 
     wscf%gamm = 0.0_DP
     do i = 1, nid