Complete version of the aderrors library

aderrors is a fortran 2008 library for error analysis of MC data with
the following features:

- Exact linear error propagation, even in iterative algorithms
  (i.e. error propagation in fit parameters).
- Handles data from any number of ensembles (i.e. simulations with
  different parameters).
- Support for replicas (i.e. several runs with the same simulation
  parameters).
- Standalone portable implementation without any external
  dependencies.
- Fast computation of autocorrelation functions with the FFT
  package (see http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html).
- Exact determination of gradients and Hessians of arbitrary
  functions.

README.md

+
+# aderrors - Error analysis of Monte Carlo data with Automatic Differentiation
+
+`aderrors` is a fortran implementation of the $`\Gamma`$-method for
+analysis of Monte Carlo data. It uses Automatic Differentiation to
+perform **exact linear error propagation**. It preforms the computation
+of gradients and Hessians of arbitrary functions, which allows a
+robust and **exact error propagation even in iterative algorithms**. 
+
+- [Features](#features)
+- [Examples](#examples)
+  - [Simple analysis of MC data](#simple-analysis-of-mc-data)
+  - [A complete example](#a-complete-example)
+  - [A calculator with uncertainties](#a-calculator-with-uncertainties)
+- [Installation](#installation)
+  - [Requirements](#requirements)
+  - [Instructions](#instructions)
+  - [Using the library](#using-the-library)
+- [Full documentation](#full-documentation)
+- [How to cite](#how-to-cite)
+
+## Features
+
+- **Exact** linear error propagation, even in iterative algorithms
+  (i.e. error propagation in fit parameters).
+- Handles data from **any number of ensembles** (i.e. simulations with
+  different parameters).
+- Support for **replicas** (i.e. several runs with the same simulation
+  parameters). 
+- Standalone **portable** implementation without any external
+  dependencies.
+- **Fast** computation of autocorrelation functions with the [FFT
+  package](http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html) (included
+  in the distribution).
+- **Exact** determination of gradients and Hessians of arbitrary
+  functions. 
+
+## Examples 
+
+This is just a small collection of examples. Note that the basic data type is
+`uwreal`, that is able to handle MC histories and data with
+errors. The error is determined by calling the method `uwerr`
+on the data type. More examples can be found in the `test` directory
+of the distribution. They are explained in the documentation `doc/aderrors.pdf`.
+
+All these examples use the `module simulator` (available in
+`test/simulator.f90`) to generate autocorrelated data with
+autocorrelation function $`\Gamma(t) = \sum_k \lambda_k
+e^{-|t|/\tau_k}`$. 
+
+### Simple analysis of MC data
+
+This code performs a simple error analysis of the data with the default
+parameters ($`S_{\tau}=4`$ for automatic windowing, no tail added to the
+autocorrelation function). The result of the analysis with the
+$`\Gamma`$-method is compared with the exact values of the error and
+$`\tau_{\rm int}`$ returned by the `module simulator`. This example is
+included in the distribution in the file `test/simple.f90`.
+
+```fortran
+program simple
+
+  use ISO_FORTRAN_ENV, Only : error_unit, output_unit
+  use numtypes
+  use aderrors
+  use simulator
+  
+  implicit none
+
+  integer, parameter  :: nd = 20000
+  type (uwreal)       :: x
+  real (kind=DP)      :: data_x(nd), err, ti, tau(5), lam(5)
+
+
+  ! Fill arrays data_x(:) with autocorrelated
+  ! data from the module simulator.
+  tau = (/1.0_DP, 3.0_DP, 4.0_DP, 5.0_DP, 7.34_DP/)
+  lam = (/1.00_DP, 0.87_DP, 1.23_DP, 0.56_DP, 0.87_DP/)
+
+  call gen_series(data_x,  err, ti, tau, lam, 0.3_DP)
+
+  ! Load data_x(:) measurements in variable x. Use
+  ! default settings (Stau=4, texp=0, 1 replica)
+  x = data_x
+  
+  ! Perform error analysis (optimal window)
+  call x%uwerr() 
+
+  ! Print results
+  write(*,'(1A,1I6,1A)')'** Measurements: ', nd, ' **   '
+  write(*,100)'  - Gamma-method:  ', x%value(),  " +/- ", x%error(), '( tauint: ', &
+       x%taui(1), " +/- ", x%dtaui(1), ')'
+  write(*,100)'  - Exact:         ', 0.3_DP,     " +/- ", err, '( tauint: ', ti, " )"
+
+100 FORMAT((2X,1A,1F8.5,1A,1F7.5,5X,1A,1F0.2,1A,1F7.5,1A))
+  
+  
+  stop
+end program simple
+```
+
+Running this code gives as output
+```
+** Measurements:  20000 **   
+    - Gamma-method:   0.29883 +/- 0.04271     ( tauint: 4.13 +/- 0.48344)
+    - Exact:          0.30000 +/- 0.04074     ( tauint: 3.82 )
+```
+
+### A complete example
+
+This example is included in the distribution in the file
+`test/complete.f90`. The analysis is performed
+for a complicated non-linear function of the two primary
+observables. It also uses replica for the MC ensemble labeled 1. 
+
+```fortran
+program complete
+
+  use ISO_FORTRAN_ENV, Only : error_unit, output_unit
+  use numtypes
+  use constants
+  use aderrors
+  use simulator
+  
+  implicit none
+
+  integer, parameter  :: nd = 5000, nrep=4
+  type (uwreal)       :: x, y, z
+  integer             :: iflog, ivrep(nrep)=(/1000,30,3070,900/), i, is, ie
+  real (kind=DP)      :: data_x(nd), data_y(nd/2), err, ti, texp
+  real (kind=DP)      :: tau(4), &
+       lam_x(4)=(/1.0_DP, 0.70_DP, 0.40_DP, 0.40_DP/), &
+       lam_y(4)=(/2.3_DP, 0.40_DP, 0.20_DP, 0.90_DP/)
+  character (len=200) :: flog='history_z.log'
+
+
+  ! Fill arrays data_x(:) with autocorrelated
+  ! data from the module simulator. Use nrep replica
+  tau = (/1.0_DP, 3.0_DP, 12.0_DP, 75.0_DP/)
+  texp = maxval(tau)
+
+  is = 1
+  do i = 1, nrep
+     ie = is + ivrep(i) - 1
+     call gen_series(data_x(is:ie),  err, ti, tau, lam_x, 0.3_DP)
+     is = ie + 1
+  end do
+
+  ! Fill data_y(:) with different values of tau also using
+  ! module simulator
+  forall (i=1:4) tau(i) = real(2*i,kind=DP)
+  call gen_series(data_y, err, ti, tau, lam_y, 1.3_DP)
+
+
+  ! Load data_x(:) measurements in variable x.
+  ! Set replica vector, exponential autocorrelation time
+  ! and ensemble ID.
+  x = data_x
+  call x%set_id(1)
+  call x%set_replica(ivrep)
+  call x%set_texp(texp)
+  
+  ! Load data_y(:) measurements in variable y
+  y = data_y
+  call y%set_id(2)
+
+  ! Exact, transparent error propagation
+  z = sin(x)/(cos(y) + 1.0_DP)
+
+  ! Attach tail in ensemble with ID 1 when signal in the
+  ! normalized auto-correlation function equals its error
+  call z%set_dsig(1.0_DP,1) 
+  ! Set Stau=3 for automatic window in ensemble with ID 2
+  call z%set_stau(3.0_DP,2)
+  ! Perform error analysis (tails, optimal window,...)
+  call z%uwerr() 
+
+  ! Print results and output details to flog  
+  write(*,'(1A,1F8.5,1A,1F8.5)')'** Observable z:   ', z%value(),  " +/- ", z%error()
+  do i = 1, z%neid()
+     write(*,'(3X,1A,1I3,3X,1F5.2,"%")',advance="no")'Contribution to error from ensemble ID', &
+          z%eid(i), 100.0_DP*z%error_src(i)
+     write(*,'(2X,1A,1F0.4,1A,1F8.4,1A)')'(tau int: ', z%taui(i), " +/- ", z%dtaui(i), ")"
+  end do
+
+  open(newunit=iflog, file=trim(flog))
+  call z%print_hist(iflog)
+  close(iflog)
+  
+  stop
+end program complete
+```
+
+Running this code gives as output
+```
+** Observable z:    0.24426 +/-  0.05374
+   Contribution to error from ensemble ID  1   83.93%  (tau int: 5.8333 +/-   2.0772)
+   Contribution to error from ensemble ID  2   16.07%  (tau int: 2.5724 +/-   0.5268)
+```
+
+
+The file `history_z.log` contains the MC histories, normalized
+autocorrelation functions and $`\tau_{\rm int}`$. The format is
+a simple text file: it can be processed by shell scripts
+(see `tools/plot/plot_hist.sh`) to produce the following graphics. 
+
+<img src="plots/ID1-hist.png" alt="ID1-hist" width="800"/>
+</br>
+<img src="plots/ID1-auto.png" alt="ID1-hist" width="400"/>
+<img src="plots/ID1-taui.png" alt="ID1-hist" width="400"/>
+
+### A calculator with uncertainties
+
+The module also handles simple data with errors (i.e. **not** a MC
+history). This is extremely useful since in many occasions we have to
+use data from statistically independent sources where only the central
+value and the error is available. 
+
+With the module `aderrors` these kind of data in handled
+transparently, just as if it were another ensemble (see `test/calculator.f90`).
+
+```fortran
+program calculator
+
+  use aderrors
+
+  implicit none
+
+  type (uwreal) :: x, y, z, t
+  
+  x = (/1.223_8, 0.012_8/) ! x = 1.223 +/- 0.012
+  y = sin(2.0_8*x)
+  z = 1.0_8 + 2.0_8 * sin(x)*cos(x)
+  t = z - y
+  call t%uwerr()
+  write(*,'(1A,1F18.16,1A,1F18.16)')'Exactly one: ', t%value(), " +/- ", t%error() ! 1.0 +/- 0.0
+  
+  stop
+end program calculator
+```
+
+Running this code gives as output
+```
+Exactly one: 1.0000000000000000 +/- 0.0000000000000000
+```
+
+## Installation
+
+### Requirements
+
+The code is strict `fortran 2008` compliant. Any compiler that
+supports this standard can be used. The code has been tested with
+`gfortran 6.X`, `gfortran 7.X`, `gfortran 8.X`, `intel fortran 17`,
+`intel fortran 18`. 
+
+Note that `gfortran 4.X` and `gfortran 5.X` do not support
+`submodules` (part of the `fortran 2008` standard). This code will not
+work whith these versions. 
+
+### Instructions
+
+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
+   be placed in the `test` directory.
+1. If preferred, move the contents of the `include` and `lib`
+   directories somewhere else.
+
+### Using the library
+
+1. Compile your programs with `-I <dir>/include`.
+1. Link your programs with the `-L <dir>/lib` and `-laderr` options. 
+
+## Full documentation
+
+Look into the `doc/aderrors.pdf` file.
+
+## How to cite
+
+If you use this package for a scientific publication, please cite the
+original work ...

build/Makefile

+
+
+INC    = ./include
+DOBJ   = ./obj
+
+FC     = gfortran6 
+FOPT   = -J$(INC) -O3 --fast-math  -march=native -ftree-vectorize \
+         -funroll-loops -std=f2008 
+
+#
+# For intel fortran try
+#
+#FC     = ifort 
+#FOPT   = -module $(INC) -O2 -stand f08 
+
+
+FFLAGS = $(FOPT)
+
+MISC     = numtypes.o constants.o time.o nonnum.o bdio.o fourier.o
+LINKED   = linked_list.o linked_list_basic.o
+ADERRORS = aderrors.o aderrors_arith.o aderrors_cf.o aderrors_func.o \
+	aderrors_io.o aderrors_cov.o aderrors_op.o aderrors_diff.o \
+	aderrors_darith.o 
+
+
+OBJWP  = $(MISC) $(RANDOM) $(LINKED) $(ADERRORS)
+OBJ    = $(foreach ob,$(OBJWP),$(DOBJ)/$(ob))
+
+VPATH=../src/misc:../src/aderrors:../src/linked_list
+.SUFFIXES: .f90 .o .x
+.PHONY: clean, all, lib, test
+
+$(DOBJ)/%.o: %.f90
+	$(FC) $(FFLAGS) -c $< -o $@ -I$(INC)
+
+%.x: %.f90 $(OBJ)
+	$(FC) $(FFLAGS) ../test/simulator.f90 $< -o$@ -L./lib/ -I$(INC)  -laderr
+
+all:  lib
+
+test: lib ../test/calculator.x ../test/uwerr.x ../test/simple.x \
+      ../test/newton.x ../test/derived.x ../test/multi.x \
+      ../test/hessian.x ../test/complete.x ../test/newton_grad.x
+
+lib: $(OBJ)
+	ar rcs ./lib/libaderr.a $(DOBJ)/*.o
+#	rm *.o
+
+clean:
+	rm obj/*.o 
+	rm include/*.mod 
+	rm include/*.smod 
+	rm lib/libaderr.a

build/include/.gitignore

+# Ignore everything in this directory
+*
+# Except this file
+!.gitignore

build/lib/.gitignore

+# Ignore everything in this directory
+*
+# Except this file
+!.gitignore

build/obj/.gitignore

+# Ignore everything in this directory
+*
+# Except this file
+!.gitignore

doc/aderrors.bbl

+\providecommand{\href}[2]{#2}\begingroup\raggedright\begin{thebibliography}{1}
+
+\bibitem{aderrors-mod}
+{Ramos, Alberto}, ``{\tt aderrors}: Error analysis of monte carlo data with
+  automatic differentiation,'' 2018.
+\newblock \url{https://gitlab.ift.uam-csic.es/alberto/aderrors}.
+
+\bibitem{fft-package}
+{Ooura, Takuya}, ``{FFT} package,'' 1996.
+\newblock \url{http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html}.
+
+\bibitem{wiki:ad}
+{Wikipedia contributors}, ``Automatic differentiation --- {W}ikipedia{,} the
+  free encyclopedia,'' 2018.
+\newblock \url{https://en.wikipedia.org/wiki/Automatic_differentiation}.
+  [Online; accessed 22-June-2018].
+
+\bibitem{web:bdio}
+{Tomasz Korzec and Hubert Simma}, ``binary data i/o,'' 2014.
+\newblock \url{http://www.bdio.org/}.
+
+\end{thebibliography}\endgroup

info/DATE

+Wed Sep  5 02:43:08 CEST 2018

src/aderrors/aderrors_cf.f90

+!
+! MODULE uwerrors_cf
+!
+! "THE BEER-WARE LICENSE":
+! Alberto Ramos wrote this file. As long as you retain this 
+! notice you can do whatever you want with this stuff. If we meet some 
+! day, and you think this stuff is worth it, you can buy me a beer in 
+! return. <alberto.ramos@cern.ch>
+!
+! file:    uwerrors_cf.f90
+! created: Thu May 10 17:55:14 2018
+!                               
+
+submodule (aderrors) aderrors_cf
+
+  use constants
+  use fourier
+  
+  implicit none
+
+  type ws_alloc
+     real (kind=DP), allocatable     :: cf(:), w(:), cdata(:)
+     integer, allocatable :: ip(:)
+     integer :: nt = -1
+  end type ws_alloc
+  
+  type (ws_alloc), save :: ws
+  
+contains
+
+! ****************************************
+! *
+  subroutine init_ws(n1, nf, n)
+! *
+! ****************************************
+    integer, intent (in)  :: n1
+    integer, intent (out) :: nf, n
+    
+    nf = int(log(real(n1+1))/lge2_dp)+2
+    n = 2**(nf+1)
+    if (n1.gt.ws%nt) then
+       if (allocated(ws%cf))    deallocate(ws%cf)
+       if (allocated(ws%cdata)) deallocate(ws%cdata)
+       if (allocated(ws%ip))    deallocate(ws%ip)
+       if (allocated(ws%w))     deallocate(ws%w)
+
+       allocate( &
+            ws%cdata(0:n-1), &
+            ws%cf(0:n1-1), &
+            ws%ip(0:3+int(sqrt(real(n)))), &
+            ws%w(0:int(n-1)) &
+            )
+
+       ws%nt = n1
+    end if
+
+    return
+  end subroutine init_ws
+    
+! ****************************************
+! *
+  subroutine init_cf(p)
+! *
+! ****************************************
+    class (uwreal), intent (inout) :: p
+
+    integer :: n
+    
+    n = max(sum(p%tmax),1)
+    if ( (allocated(p%gamm)).and.(n.eq.size(p%gamm)).and. &
+         (allocated(p%drho)).and.(n.eq.size(p%drho)) ) return
+
+    if (allocated(p%gamm)) deallocate(p%gamm)
+    if (allocated(p%drho)) deallocate(p%drho)
+    allocate(p%gamm(0:n-1), p%drho(0:n-1))
+
+    return
+  end subroutine init_cf
+
+! ****************************************
+! *
+  module subroutine uwerr(p)
+! *
+! ****************************************
+    class (uwreal), intent (inout) :: p
+
+    integer :: n, irst, irnd
+
+    call init_cf(p)
+    call autocorrfunc_obs(p)
+    call wopt_texp(p)
+    call tauint(p)
+   
+    p%var  = 0.0_DP
+    p%derr = 0.0_DP
+    do n = 1, p%nid
+       if (p%nd(n).eq.1) then
+          irst     = 1 + sum(p%nd(1:n-1)) 
+          p%var(n) = p%data(irst)**2
+       else
+          irst = sum(p%tmax(1:n-1))
+          irnd = irst + p%iw(n)
+          
+          p%var(n) = p%gamm(irst)/real(p%nd(n),kind=dp) + &
+               2.0_dp * sum(p%gamm(irst+1:irnd))/real(p%nd(n),kind=dp)
+          if (p%texp(n) > 0.0_dp) p%var(n) = p%var(n) + &
+               2.0_dp*p%texp(n)*p%gamm(irnd)/real(p%nd(n),kind=dp)
+
+          p%derr = p%derr + &
+               p%var(n)*(real(p%iw(n),kind=DP)-0.5_DP) / real(p%nd(n),kind=DP)
+       end if
+    end do
+
+    p%err  = sqrt(sum(p%var))
+    p%derr = sqrt(p%derr)
+
+    return
+  end subroutine uwerr
+  
+! ****************************************
+! *
+  subroutine autocorrfunc_obs(p)
+! *
+! ****************************************
+    class (uwreal), intent (inout) :: p
+
+    integer        :: ist, ind, j, i, k, n, irst, irnd, icnd, iadd, nrcnt, irmax
+    real (kind=dp) :: cont, dbias
+
+    p%gamm = 0.0_DP
+    n = 0
+    do n = 1, p%nid
+       if (p%nd(n).eq.1) cycle
+       call p%get_cf_offset(irst,irnd,n)
+       irmax = irnd
+       do j = 1, p%nrep(n)
+          call p%get_offset(ist,ind,n,j)
+
+          call autocf(p%data(ist:ind))
+          irnd = min(irst + (ind-ist)/2, irmax)
+          icnd = irnd - irst
+          p%gamm(irst:irnd) = p%gamm(irst:irnd) + ws%cf(0:icnd)
+       end do
+    end do
+    
+    irst = 0
+    do n = 1, p%nid
+       if (p%nd(n).eq.1) cycle
+       call p%get_vrep_offset(ist,ind,n)
+       do j = 0, p%tmax(n)-1
+          nrcnt = count(p%ivrep(ist:ind)/2 .ge. j)
+          p%gamm(j+irst) = p%gamm(j+irst) / real(p%nd(n) - nrcnt*j,kind=dp)
+       end do
+       irst = irst + p%tmax(n)
+    end do
+
+    call wopt(p)
+    irst = 0
+    do n = 1, p%nid
+       if (p%nd(n).eq.1) cycle
+       irnd = irst + p%tmax(n) - 1
+       dbias = p%gamm(irst) + 2.0_DP*sum(p%gamm(irst+1:irst+p%iw(n)))
+       p%gamm(irst:irnd) = p%gamm(irst:irnd) + dbias/real(p%nd(n),kind=DP)
+       irst = irst + p%tmax(n)
+    end do
+    
+    irst = 0
+    do n = 1, p%nid
+       if (p%nd(n).eq.1) cycle
+       irnd = irst + p%tmax(n) - 1
+       p%drho(irst:irnd) = 0.0_dp
+
+       if (p%gamm(irst) == 0.0_dp) then
+          p%drho(irst:irnd) = 0.0_dp
+       else
+          iadd = p%iw(n)
+          do i = irst, irnd
+             ist  = max(1,i-irst-iadd) + irst
+             ind  = i + iadd
+
+             do k = ist, ind
+                if (k.lt.p%tmax(n)+irst) then
+                   cont = -2.0_dp*p%gamm(k)*p%gamm(i)/p%gamm(irst)**2
+                else
+                   cont = 0.0_DP
+                end if
+                
+                if ((i+k-2*irst) < p%tmax(n)) then
+                   cont = cont + p%gamm(i+k-irst)/p%gamm(irst)
+                end if
+                if (abs(i-k) < p%tmax(n)) then
+                   cont = cont + p%gamm(abs(i-k)+irst)/p%gamm(irst)
+                end if
+                
+                p%drho(i) = p%drho(i) + cont**2
+             end do
+             p%drho(i) = sqrt(p%drho(i)/real(p%nd(n),kind=dp))
+          end do
+       end if
+
+       irst = irst + p%tmax(n)
+    end do
+
+    return
+  end subroutine autocorrfunc_obs
+
+! ****************************************
+! *
+  subroutine autocf(data)
+! *
+! ****************************************
+
+    real (kind=dp), intent (in) :: data(:)
+
+    integer :: i, nt, nfourier, n
+    
+    nt = size(data)
+    call init_ws(nt, nfourier, n)
+    ws%cdata = 0.0_dp
+    do i = 0, nt-1
+       ws%cdata(2*i) = data(i+1)
+    end do
+
+    ws%ip(0) = 0
+    call cdft(n, 1, ws%cdata, ws%ip, ws%w)
+    do i = 0, n/2-1
+       ws%cdata(2*i)   = ws%cdata(2*i)**2 + ws%cdata(2*i+1)**2
+       ws%cdata(2*i+1) = 0.0_dp
+    end do
+    
+    call cdft(n, -1, ws%cdata, ws%ip, ws%w)
+    ws%cdata(0:n-1) = 2.0_dp*ws%cdata(0:n-1)/real(n,kind=dp)
+    
+    do i = 0, nt-1
+       ws%cf(i) = real(ws%cdata(2*i),kind=DP)
+    end do
+    
+    return
+  end subroutine autocf
+
+! ********************************  
+! *
+  subroutine wopt(p)
+! *
+! ********************************
+    type (uwreal), intent (inout) :: p
+    
+    integer :: i=0, n, ist, ind
+    Real (kind=DP) :: tiw, tau, gw, rw
+
+    do n = 1, p%nid
+       if (p%nd(n).eq.1) then
+          p%iw(n) = 0
+          cycle
+       end if
+       call p%get_cf_offset(ist,ind,n)
+       
+       tiw = 0.5_DP
+       if (p%gamm(ist)==0.0_DP) then
+          p%iw(n) = 1
+       else
+          do i = ist+1, ind
+             rw  = real(i-ist,kind=DP)
+             tiw = tiw + p%gamm(i)/p%gamm(ist)
+             if (tiw < 0.5_DP) then
+                tau = 0.01_DP
+             else
+                tau = p%s(n)/log((2.0_DP*tiw+1.0_DP)/(2.0_DP*tiw-1.0_DP))
+             end if
+             gw  = exp(-rw/tau) - tau/sqrt(rw*real(p%nd(n),kind=DP))
+
+             if (gw < 0.0_DP) exit
+          end do
+          p%iw(n) = i-ist
+       end if
+    end do
+
+    return
+  end subroutine wopt
+  
+! ********************************  
+! *
+  subroutine wopt_texp(p)
+! *
+! ********************************
+    type (uwreal), intent (inout) :: p
+    
+    integer :: i=0, n, ist, ind
+    Real (kind=DP) :: tiw
+
+    do n = 1, p%nid
+       if (p%nd(n).eq.1) then
+          p%iw(n) = 0
+          cycle
+       end if
+       call p%get_cf_offset(ist,ind,n)
+       
+       tiw = 0.5_DP
+       if (p%texp(n) > 0.0_DP) then
+          do i = ist+1, ind
+             if (p%dsig(n)*p%drho(i) > p%gamm(i)/p%gamm(ist)) exit
+          end do
+          p%iw(n) = i-ist
+       end if
+    end do
+
+    return