PROGRAM vecgen
!
!  This program randomly generates an N-length vector and writes it to
!  a file, with N on the first line, and the rest of the array on
!  following line(s).
!
   IMPLICIT NONE
   INTEGER :: N, k, KK(1)
   DOUBLE PRECISION, ALLOCATABLE :: V(:), Z(:)
   EXTERNAL :: save_vec_to_file, read_vec_from_file
   INTEGER :: read_vec_from_file

   PRINT *, 'Enter N (> 0)'
   READ *, N
   IF (N .lt. 1) THEN
      PRINT*, 'N must be greater than 0, but you used:', N
      STOP
   END IF
   ALLOCATE(V(N))               ! allocate array
   CALL RANDOM_NUMBER(V)        ! randomly generate array
!
!  Write array to file, then read back into Z
!
   CALL save_vec_to_file(N, V)
   PRINT *, 'Results written to file ''vector.dat'''

   ALLOCATE(Z(N))               ! allocate new array
   k = read_vec_from_file(N, Z) ! should be duplicate of V
   PRINT *, 'sizediff = ', N-k
!
!  Do error check; can be slightly different due to rounding
!
   Z = abs(V - Z)                ! zero with no rounding
   kk = maxloc(Z)                ! biggest difference from zero
   k = kk(1)
   PRINT *, 'Maximum error found in location = ', k
   PRINT *, 'Error, Relative error = ', Z(k), Z(k) / abs(V(k))
   PRINT *, 'diffsum = ', sum(Z)
   DEALLOCATE(V)
   DEALLOCATE(Z)

END PROGRAM vecgen

SUBROUTINE save_vec_to_file(N, X)
!
!  This routine saves an N-length array to 'vector.dat'.  The first line
!  will hold the vector length, and the next line will have all elements of
!  X separated by spaces
!
   IMPLICIT NONE
   INTEGER, INTENT(IN) :: N
   DOUBLE PRECISION, INTENT(IN) :: X(*)

   OPEN(unit=10, file='vector.dat', status='replace', action='write')
   WRITE(10, *) N
   WRITE(10, *) X(1:N)
   CLOSE(10)

END SUBROUTINE save_vec_to_file

FUNCTION read_vec_from_file(N, X)
!
!  This routine reads an array from 'vector.dat'.  The first line
!  will hold the vector length, and the next line will have all elements of
!  X separated by spaces.  N is the length of X, and will be checked against
!  the file length, to ensure X is big enough to hold the file's array
!
   IMPLICIT NONE
   INTEGER read_vec_from_file
   INTEGER, INTENT(IN) :: N
   DOUBLE PRECISION, INTENT(OUT) :: X(*)
   INTEGER :: M

   OPEN(unit=10, file='vector.dat', status='old', action='read')
   READ(10, *) M
   IF (M > N) THEN
      PRINT*, 'FATAL: Need array >= ', M, ' in size, but only have', N, '!'
      STOP
   END IF
   READ(10, *) X(1:M)
   CLOSE(10)
   read_vec_from_file = M

END FUNCTION read_vec_from_file