Add routines for saving/loading arrays in npy format

doc/specs/stdlib_io.md

@@ -36,7 +36,7 @@ program demo_loadtxt
     use stdlib_io, only: loadtxt
     implicit none
     real, allocatable :: x(:,:)
-    call loadtxt('example.dat', x) 
+    call loadtxt('example.dat', x)
 end program demo_loadtxt
 ```
 
@@ -128,6 +128,90 @@ program demo_savetxt
     use stdlib_io, only: savetxt
     implicit none
     real :: x(3,2) = 1
-    call savetxt('example.dat', x) 
+    call savetxt('example.dat', x)
+end program demo_savetxt
+```
+
+
+## `load_npy`
+
+### Status
+
+Experimental
+
+### Description
+
+Loads am `array` from a npy formatted binary file.
+
+### Syntax
+
+`call [[stdlib_io_npy(module):load_npy(interface)]](filename, array[, iostat][, iomsg])`
+
+### Arguments
+
+`filename`: Shall be  a character expression containing the file name from which to load the `array`.
+
+`array`: Shall be an allocatable array of any rank of type `real`, `complex` or `integer`.
+
+`iostat`: Default integer, contains status of loading to file, zero in case of success.
+          Optional argument, in case not present the program will halt for non-zero status.
+
+`iomsg`: Deferred length character value, contains error message in case `iostat` is non-zero.
+         Optional argument, error message will be dropped if not present.
+
+### Return value
+
+Returns an allocated `array` with the content of `filename` in case of success.
+
+### Example
+
+```fortran
+program demo_loadtxt
+    use stdlib_io_npy, only: load_npy
+    implicit none
+    real, allocatable :: x(:,:)
+    call loadtxt('example.npy', x)
+end program demo_loadtxt
+```
+
+
+## `save_npy`
+
+### Status
+
+Experimental
+
+### Description
+
+Saves an `array` into a npy formatted binary file.
+
+### Syntax
+
+`call [[stdlib_io_npy(module):save_npy(interface)]](filename, array[, iostat][, iomsg])`
+
+### Arguments
+
+`filename`: Shall be  a character expression containing the name of the file that will contain the `array`.
+
+`array`: Shall be an array of any rank of type `real`, `complex` or `integer`.
+
+`iostat`: Default integer, contains status of saving to file, zero in case of success.
+          Optional argument, in case not present the program will halt for non-zero status.
+
+`iomsg`: Deferred length character value, contains error message in case `iostat` is non-zero.
+         Optional argument, error message will be dropped if not present.
+
+### Output
+
+Provides a text file called `filename` that contains the rank-2 `array`.
+
+### Example
+
+```fortran
+program demo_savetxt
+    use stdlib_io_npy, only: save_npy
+    implicit none
+    real :: x(3,2) = 1
+    call save_npy('example.npy', x)
 end program demo_savetxt
 ```

src/CMakeLists.txt

@@ -7,6 +7,9 @@ set(fppFiles
     stdlib_bitsets_64.fypp
     stdlib_bitsets_large.fypp
     stdlib_io.fypp
+    stdlib_io_npy.fypp
+    stdlib_io_npy_load.fypp
+    stdlib_io_npy_save.fypp
     stdlib_kinds.fypp
     stdlib_linalg.fypp
     stdlib_linalg_diag.fypp

src/Makefile.manual

@@ -4,6 +4,9 @@ SRCFYPP = \
         stdlib_bitsets_large.fypp \
         stdlib_bitsets.fypp \
         stdlib_io.fypp \
+        stdlib_io_npy.fypp \
+        stdlib_io_npy_load.fypp \
+        stdlib_io_npy_save.fypp \
         stdlib_kinds.fypp \
         stdlib_linalg.fypp \
         stdlib_linalg_diag.fypp \
@@ -87,6 +90,16 @@ stdlib_io.o: \
     stdlib_optval.o \
     stdlib_kinds.o \
     stdlib_ascii.o
+stdlib_io_npy.o: \
+    stdlib_kinds.o
+stdlib_io_npy_load.o: \
+    stdlib_io_npy.o \
+    stdlib_error.o \
+    stdlib_strings.o
+stdlib_io_npy_save.o: \
+    stdlib_io_npy.o \
+    stdlib_error.o \
+    stdlib_strings.o
 stdlib_linalg.o: \
     stdlib_kinds.o
 stdlib_linalg_diag.o: \

src/stdlib_io_npy.fypp

@@ -0,0 +1,120 @@
+! SPDX-Identifer: MIT
+
+#:include "common.fypp"
+#:set RANKS = range(1, MAXRANK + 1)
+#:set KINDS_TYPES = REAL_KINDS_TYPES + INT_KINDS_TYPES + CMPLX_KINDS_TYPES
+
+!> Description of the npy format taken from
+!> https://numpy.org/doc/stable/reference/generated/numpy.lib.format.html
+!>
+!>## Format Version 1.0
+!>
+!> The first 6 bytes are a magic string: exactly \x93NUMPY.
+!>
+!> The next 1 byte is an unsigned byte:
+!> the major version number of the file format, e.g. \x01.
+!>
+!> The next 1 byte is an unsigned byte:
+!> the minor version number of the file format, e.g. \x00.
+!> Note: the version of the file format is not tied to the version of the numpy package.
+!>
+!> The next 2 bytes form a little-endian unsigned short int:
+!> the length of the header data HEADER_LEN.
+!>
+!> The next HEADER_LEN bytes form the header data describing the array’s format.
+!> It is an ASCII string which contains a Python literal expression of a dictionary.
+!> It is terminated by a newline (\n) and padded with spaces (\x20) to make the total
+!> of len(magic string) + 2 + len(length) + HEADER_LEN be evenly divisible by 64 for
+!> alignment purposes.
+!>
+!> The dictionary contains three keys:
+!>
+!> - “descr”: dtype.descr
+!>   An object that can be passed as an argument to the numpy.dtype constructor
+!>   to create the array’s dtype.
+!>
+!> - “fortran_order”: bool
+!>   Whether the array data is Fortran-contiguous or not. Since Fortran-contiguous
+!>   arrays are a common form of non-C-contiguity, we allow them to be written directly
+!>   to disk for efficiency.
+!>
+!> - “shape”: tuple of int
+!>   The shape of the array.
+!>
+!> For repeatability and readability, the dictionary keys are sorted in alphabetic order.
+!> This is for convenience only. A writer SHOULD implement this if possible. A reader MUST
+!> NOT depend on this.
+!>
+!> Following the header comes the array data. If the dtype contains Python objects
+!> (i.e. dtype.hasobject is True), then the data is a Python pickle of the array.
+!> Otherwise the data is the contiguous (either C- or Fortran-, depending on fortran_order)
+!> bytes of the array. Consumers can figure out the number of bytes by multiplying the
+!> number of elements given by the shape (noting that shape=() means there is 1 element)
+!> by dtype.itemsize.
+!>
+!>## Format Version 2.0
+!>
+!> The version 1.0 format only allowed the array header to have a total size of 65535 bytes.
+!> This can be exceeded by structured arrays with a large number of columns.
+!> The version 2.0 format extends the header size to 4 GiB. numpy.save will automatically
+!> save in 2.0 format if the data requires it, else it will always use the more compatible
+!> 1.0 format.
+!>
+!> The description of the fourth element of the header therefore has become:
+!> “The next 4 bytes form a little-endian unsigned int: the length of the header data
+!> HEADER_LEN.”
+!>
+!>## Format Version 3.0
+!>
+!> This version replaces the ASCII string (which in practice was latin1) with a
+!> utf8-encoded string, so supports structured types with any unicode field names.
+module stdlib_io_npy
+    use stdlib_kinds, only : int8, int16, int32, int64, sp, dp, xdp, qp
+    implicit none
+    private
+
+    public :: save_npy, load_npy
+
+
+    !> Save multidimensional array in npy format
+    interface save_npy
+    #:for k1, t1 in KINDS_TYPES
+      #:for rank in RANKS
+        module subroutine save_npy_${t1[0]}$${k1}$_${rank}$(filename, array, iostat, iomsg)
+            character(len=*), intent(in) :: filename
+            ${t1}$, intent(in) :: array${ranksuffix(rank)}$
+            integer, intent(out), optional :: iostat
+            character(len=:), allocatable, intent(out), optional :: iomsg
+        end subroutine save_npy_${t1[0]}$${k1}$_${rank}$
+      #:endfor
+    #:endfor
+    end interface save_npy
+
+    !> Load multidimensional array in npy format
+    interface load_npy
+    #:for k1, t1 in KINDS_TYPES
+      #:for rank in RANKS
+        module subroutine load_npy_${t1[0]}$${k1}$_${rank}$(filename, array, iostat, iomsg)
+            character(len=*), intent(in) :: filename
+            ${t1}$, allocatable, intent(out) :: array${ranksuffix(rank)}$
+            integer, intent(out), optional :: iostat
+            character(len=:), allocatable, intent(out), optional :: iomsg
+        end subroutine load_npy_${t1[0]}$${k1}$_${rank}$
+      #:endfor
+    #:endfor
+    end interface load_npy
+
+
+    character(len=*), parameter :: nl = achar(10)
+
+    character(len=*), parameter :: &
+        type_iint8 = "<i1", type_iint16 = "<i2", type_iint32 = "<i4", type_iint64 = "<i8", &
+        type_rsp = "<f4", type_rdp = "<f8", type_rxdp = "<f10", type_rqp = "<f16", &
+        type_csp = "<c8", type_cdp = "<c16", type_cxdp = "<c20", type_cqp = "<c32"
+
+    character(len=*), parameter :: &
+        & magic_number = char(int(z"93")), &
+        & magic_string = "NUMPY"
+
+
+end module stdlib_io_npy