Merge pull request #420 from ejovo13/lin_log_space

First implementation of real-valued linspace.

Author: jvdp1

doc/specs/stdlib_math.md

@@ -19,7 +19,7 @@ title: math
 
 #### Description
 
-Returns a value which lies in the given interval [`xmin`, `xmax`] (interval is `xmin` and `xmax` inclusive) and is closest to the input value `x`.  
+Returns a value which lies in the given interval [`xmin`, `xmax`] (interval is `xmin` and `xmax` inclusive) and is closest to the input value `x`.
 
 #### Syntax
 
@@ -35,8 +35,8 @@ Elemental function.
 
 #### Argument(s)
 
-`x`: scalar of either `integer` or `real` type. This argument is `intent(in)`.  
-`xmin`: scalar of either `integer` or `real` type. This argument is `intent(in)`.  
+`x`: scalar of either `integer` or `real` type. This argument is `intent(in)`.
+`xmin`: scalar of either `integer` or `real` type. This argument is `intent(in)`.
 `xmax`: scalar of either `integer` or `real` type, which must be greater than or equal to `xmin`. This argument is `intent(in)`.
 
 Note: All arguments must have same `type` and same `kind`.
@@ -90,3 +90,188 @@ program demo_clip_real
   ! clipped_value <- 3.02500010
 end program demo_clip_real
 ```
+
+### `linspace` - Create a linearly spaced rank one array
+
+#### Description
+
+Returns a linearly spaced rank 1 array from [`start`, `end`]. Optionally, you can specify the length of the returned array by passing `n`.
+
+#### Syntax
+
+`res = [[stdlib_math(module):linspace(interface)]] (start, end [, n])`
+
+#### Status
+
+Experimental
+
+#### Class
+
+Function.
+
+#### Argument(s)
+
+`start`: Shall be scalar of any numeric type or kind. This argument is `intent(in)`.
+`end`: Shall be the same `type` and `kind` as `start`. This argument is `intent(in)`.
+`n`: Shall be an integer specifying the length of the output. This argument is `optional` and `intent(in)`.
+
+#### Output value or Result value
+
+The output is a rank 1 array whose length is either 100 (default value) or `n`.
+
+If `n` == 1, return a rank 1 array whose only element is `end`.
+If `n` <= 0, return a rank 1 array with length 0.
+
+If `start`/`end` are `real` or `complex` types, the `result` will be of the same type and kind as `start`/`end`.
+If `start`/`end` are `integer` types, the `result` will default to a `real(dp)` array.
+
+#### Examples
+
+##### Example 1:
+
+Here inputs are of type `complex` and kind `dp`
+```fortran
+program demo_linspace_complex
+  use stdlib_math, only: linspace
+  use stdlib_kinds, only: dp
+  implicit none
+
+  complex(dp) :: start = complex(10.0_dp, 5.0_dp)
+  complex(dp) :: end = complex(-10.0_dp, 15.0_dp)
+
+  complex(dp) :: z(11)
+
+  z = linspace(start, end, 11)
+end program demo_linspace_complex
+```
+
+##### Example 2:
+
+Here inputs are of type `integer` and kind `int16`, with the result defaulting to `real(dp)`.
+```fortran
+program demo_linspace_int16
+  use stdlib_math, only: linspace
+  use stdlib_kinds, only: int16, dp
+  implicit none
+
+  integer(int16) :: start = 10_int16
+  integer(int16) :: end = 23_int16
+
+  real(dp) :: r(15)
+
+  r = linspace(start, end, 15)
+end program demo_linspace_int16
+```
+
+### `logspace` - Create a logarithmically spaced rank one array
+
+#### Description
+
+Returns a logarithmically spaced rank 1 array from [`base`^`start`, `base`^`end`]. The default size of the array is 50. Optionally, you can specify the length of the returned array by passing `n`. You can also specify the `base` used to compute the range (default 10).
+
+#### Syntax
+
+`res = [[stdlib_math(module):logspace(interface)]] (start, end [, n [, base]])`
+
+#### Status
+
+Experimental
+
+#### Class
+
+Function.
+
+#### Argument(s)
+
+`start`: Shall be a scalar of any numeric type. All kinds are supported for real and complex arguments. For integers, only the default kind is currently implemented. This argument is `intent(in)`.
+`end`: Shall be the same `type` and `kind` as `start`. This argument is `intent(in)`.
+`n`: Shall be an integer specifying the length of the output. This argument is `optional` and `intent(in)`.
+`base` : Shall be a scalar of any numeric type. All kinds are supported for real and complex arguments. For integers, only the default kind is currently implemented. This argument is `optional` and `intent(in)`.
+
+#### Output value or Result value
+
+The output is a rank 1 array whose length is either 50 (default value) or `n`.
+
+If `n` == 1, return a rank 1 array whose only element is `base`^`end`.
+If `n` <= 0, return a rank 1 array with length 0
+
+The `type` and `kind` of the output is dependent on the `type` and `kind` of the passed parameters.
+
+For function calls where the `base` is not specified: `logspace(start, end)`/`logspace(start, end, n)`, the `type` and `kind` of
+the output follows the same scheme as above for `linspace`.
+>If `start`/`end` are `real` or `complex` types, the `result` will be the same type and kind as `start`/`end`.
+>If `start`/`end` are integer types, the `result` will default to a `real(dp)` array.
+
+For function calls where the `base` is specified, the `type` and `kind` of the result is in accordance with the following table:
+
+| `start`/`end` | `n` | `base` | `output` |
+| ------------- | --- | ------ | -------- |
+| `real(KIND)` | `Integer` | `real(KIND)` | `real(KIND)` |
+| "          " | "       " | `complex(KIND)` | `complex(KIND)` |
+| "          " | "       " | `Integer` | `real(KIND)` |
+| `complex(KIND)` | "       " | `real(KIND)` | `complex(KIND)` |
+| "             " | "       " | `complex(KIND)` | `complex(KIND)` |
+| "             " | "       " | `Integer` | `complex(KIND)` |
+| `Integer` | "        " | `real(KIND)` | `real(KIND)` |
+| "              " | "        " | `complex(KIND)` | `complex(KIND)` |
+| "              " | "        " | `Integer` | `Integer` |
+
+#### Examples
+
+##### Example 1:
+
+Here inputs are of type `complex` and kind `dp`. `n` and `base` is not specified and thus default to 50 and 10, respectively.
+```fortran
+program demo_logspace_complex
+  use stdlib_math, only: logspace
+  use stdlib_kinds, only: dp
+  implicit none
+
+  complex(dp) :: start = (10.0_dp, 5.0_dp)
+  complex(dp) :: end = (-10.0_dp, 15.0_dp)
+
+  complex(dp) :: z(11) ! Complex values raised to complex powers results in complex values
+
+  z = logspace(start, end, 11)
+end program demo_logspace_complex
+```
+
+##### Example 2:
+
+Here inputs are of type `integer` and default kind. `base` is not specified and thus defaults to 10.
+```fortran
+program demo_logspace_int
+  use stdlib_math, only: logspace
+  use stdlib_kinds, only: dp
+  implicit none
+
+  integer :: start = 10
+  integer :: end = 23
+  integer :: n = 15
+
+  real(dp) :: r(n) ! Integer values raised to real powers results in real values
+
+  r = logspace(start, end, n)
+end program demo_logspace_int
+```
+
+##### Example 3:
+
+Here `start`/`end` are of type `real` and double precision. `base` is type `complex` and also double precision.
+```fortran
+program demo_logspace_rstart_cbase
+  use stdlib_math, only: logspace
+  use stdlib_kinds, only: dp
+  implicit none
+
+  real(dp) :: start = 0.0_dp
+  real(dp) :: end = 3.0_dp
+  integer :: n = 4
+  complex(dp) :: base = (0.0_dp, 1.0_dp)
+
+  complex(dp) :: z(n) ! complex values raised to real powers result in complex values
+
+  z = logspace(start, end, n, base)
+
+end program demo_logspace_rstart_cbase
+```

src/CMakeLists.txt

@@ -29,6 +29,8 @@ set(fppFiles
     stdlib_quadrature_simps.fypp
     stdlib_stats_distribution_PRNG.fypp
     stdlib_math.fypp
+    stdlib_math_linspace.fypp
+    stdlib_math_logspace.fypp
     stdlib_string_type.fypp
 )
 

src/Makefile.manual

@@ -25,6 +25,8 @@ SRCFYPP =\
         stdlib_stats_moment_scalar.fypp \
         stdlib_stats_var.fypp \
         stdlib_math.fypp \
+		stdlib_math_linspace.fypp \
+		stdlib_math_logspace.fypp \
         stdlib_stats_distribution_PRNG.fypp \
         stdlib_string_type.fypp
 
@@ -73,6 +75,7 @@ stdlib_error.o: stdlib_optval.o
 stdlib_specialfunctions.o: stdlib_kinds.o
 stdlib_specialfunctions_legendre.o: stdlib_kinds.o stdlib_specialfunctions.o
 stdlib_io.o: \
+	stdlib_ascii.o \
 	stdlib_error.o \
 	stdlib_optval.o \
 	stdlib_kinds.o
@@ -141,4 +144,8 @@ stdlib_strings.o: stdlib_ascii.o \
     stdlib_string_type.o \
 	stdlib_optval.o
 stdlib_math.o: stdlib_kinds.o
+stdlib_math_linspace.o: \
+	stdlib_math.o
+stdlib_math_logspace.o: \
+	stdlib_math_linspace.o
 stdlib_linalg_outer_product.o: stdlib_linalg.o

src/stdlib_math.fypp

@@ -1,21 +1,268 @@
 #:include "common.fypp"
 #:set IR_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
 
 module stdlib_math
     use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, qp
 
     implicit none
     private
-    public :: clip
+    public :: clip, linspace, logspace
+    public :: EULERS_NUMBER_SP, EULERS_NUMBER_DP, EULERS_NUMBER_QP
+    public :: DEFAULT_LINSPACE_LENGTH, DEFAULT_LOGSPACE_BASE, DEFAULT_LOGSPACE_LENGTH
+
+    integer, parameter :: DEFAULT_LINSPACE_LENGTH = 100
+    integer, parameter :: DEFAULT_LOGSPACE_LENGTH = 50
+    integer, parameter :: DEFAULT_LOGSPACE_BASE = 10
+
+    ! Useful constants for lnspace
+    real(sp), parameter :: EULERS_NUMBER_SP = exp(1.0_sp)
+    real(dp), parameter :: EULERS_NUMBER_DP = exp(1.0_dp)
+    real(qp), parameter :: EULERS_NUMBER_QP = exp(1.0_qp)
 
     interface clip
     #:for k1, t1 in IR_KINDS_TYPES
         module procedure clip_${k1}$
     #:endfor
     end interface clip
 
+    interface linspace
+    !! Version: Experimental
+    !!
+    !! Create rank 1 array of linearly spaced elements
+    !! If the number of elements is not specified, create an array with size 100. If n is a negative value,
+    !! return an array with size 0. If n = 1, return an array whose only element is end
+    !!([Specification](../page/specs/stdlib_math.html#linspace-create-a-linearly-spaced-rank-one-array))
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:set RName = rname("linspace_default", 1, t1, k1)
+      module function ${RName}$(start, end) result(res)
+        ${t1}$, intent(in) :: start
+        ${t1}$, intent(in) :: end
+
+        ${t1}$ :: res(DEFAULT_LINSPACE_LENGTH)
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:set RName = rname("linspace_n", 1, t1, k1)
+      module function ${RName}$(start, end, n) result(res)
+        ${t1}$, intent(in) :: start
+        ${t1}$, intent(in) :: end
+        integer, intent(in) :: n
+
+        ${t1}$ :: res(n)
+      end function ${RName}$
+    #:endfor
+
+
+    ! Add support for integer linspace
+    !!
+    !! When dealing with integers as the `start` and `end` parameters, the return type is always a `real(dp)`.
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("linspace_default", 1, t1, k1)
+      #! The interface for INT_KINDS_TYPES cannot be combined with RC_KINDS_TYPES
+      #! because the output for integer types is always a real with dp.
+      module function ${RName}$(start, end) result(res)
+        ${t1}$, intent(in) :: start
+        ${t1}$, intent(in) :: end
+
+        real(dp) :: res(DEFAULT_LINSPACE_LENGTH)
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("linspace_n", 1, t1, k1)
+      module function ${RName}$(start, end, n) result(res)
+        ${t1}$, intent(in) :: start
+        ${t1}$, intent(in) :: end
+        integer, intent(in) :: n
+
+        real(dp) :: res(n)
+      end function ${RName}$
+    #:endfor
+
+  end interface
+
+  interface logspace
+  !! Version: Experimental
+  !!
+  !! Create rank 1 array of logarithmically spaced elements from base**start to base**end.
+  !! If the number of elements is not specified, create an array with size 50. If n is a negative value,
+  !! return an array with size 0. If n = 1, return an array whose only element is base**end. If no base
+  !! is specified, logspace will default to using a base of 10
+  !!
+  !!([Specification](../page/specs/stdlib_math.html#logspace-create-a-logarithmically-spaced-rank-one-array))
+  #!=========================================================
+  #!=     logspace(start, end)                              =
+  #!=========================================================
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:set RName = rname("logspace", 1, t1, k1, "default")
+    module function ${RName}$(start, end) result(res)
+
+      ${t1}$, intent(in) :: start
+      ${t1}$, intent(in) :: end
+
+      ${t1}$ :: res(DEFAULT_LOGSPACE_LENGTH)
+
+    end function ${RName}$
+  #:endfor
+  #! Integer support
+  #:set RName = rname("logspace", 1, "integer(int32)", "int32", "default")
+    module function ${RName}$(start, end) result(res)
+
+      integer, intent(in) :: start
+      integer, intent(in) :: end
+
+      real(dp) :: res(DEFAULT_LOGSPACE_LENGTH)
+
+  end function ${RName}$
+
+  #!=========================================================
+  #!=     logspace(start, end, n)                           =
+  #!=========================================================
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:set RName = rname("logspace", 1, t1, k1, "n")
+    module function ${RName}$(start, end, n) result(res)
+      ${t1}$, intent(in) :: start
+      ${t1}$, intent(in) :: end
+      integer, intent(in) :: n
+
+      ${t1}$ :: res(n)
+    end function ${RName}$
+  #:endfor
+  #! Integer support
+  #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n")
+    module function ${RName}$(start, end, n) result(res)
+      integer, intent(in) :: start
+      integer, intent(in) :: end
+      integer, intent(in) :: n
+
+      real(dp) :: res(n)
+    end function ${RName}$
+
+  #!=========================================================
+  #!=     logspace(start, end, n, base)                     =
+  #!=========================================================
+  #! Need another function where base is not optional,
+  #! otherwise the compiler can not differentiate between
+  #! generic calls to logspace_n where a base is not present
+  #! ========================================================
+  #:for k1, t1 in REAL_KINDS_TYPES
+    ! Generate logarithmically spaced sequence from ${k1}$ base to the powers
+    ! of ${k1}$ start and end. [base^start, ... , base^end]
+    ! Different combinations of parameter types will lead to different result types.
+    ! Those combinations are indicated in the body of each function.
+    #:set RName = rname("logspace", 1, t1, k1, "n_rbase")
+    module function ${RName}$(start, end, n, base) result(res)
+      ${t1}$, intent(in) :: start
+      ${t1}$, intent(in) :: end
+      integer, intent(in) :: n
+      ${t1}$, intent(in) :: base
+      ! real(${k1}$) endpoints + real(${k1}$) base = real(${k1}$) result
+      ${t1}$ :: res(n)
+    end function ${RName}$
+
+    #:set RName = rname("logspace", 1, t1, k1, "n_cbase")
+    module function ${RName}$(start, end, n, base) result(res)
+      ${t1}$, intent(in) :: start
+      ${t1}$, intent(in) :: end
+      integer, intent(in) :: n
+      complex(${k1}$), intent(in) :: base
+      ! real(${k1}$) endpoints + complex(${k1}$) base = complex(${k1}$) result
+      ${t1}$ :: res(n)
+    end function ${RName}$
+
+    #:set RName = rname("logspace", 1, t1, k1, "n_ibase")
+    module function ${RName}$(start, end, n, base) result(res)
+      ${t1}$, intent(in) :: start
+      ${t1}$, intent(in) :: end
+      integer, intent(in) :: n
+      integer, intent(in) :: base
+      ! real(${k1}$) endpoints + integer base = real(${k1}$) result
+      ${t1}$ :: res(n)
+    end function ${RName}$
+  #:endfor
+  #! ========================================================
+  #! ========================================================
+  #:for k1, t1 in CMPLX_KINDS_TYPES
+    ! Generate logarithmically spaced sequence from ${k1}$ base to the powers
+    ! of ${k1}$ start and end. [base^start, ... , base^end]
+    ! Different combinations of parameter types will lead to different result types.
+    ! Those combinations are indicated in the body of each function.
+    #:set RName = rname("logspace", 1, t1, k1, "n_rbase")
+    module function ${RName}$(start, end, n, base) result(res)
+      ${t1}$, intent(in) :: start
+      ${t1}$, intent(in) :: end
+      integer, intent(in) :: n
+      real(${k1}$), intent(in) :: base
+      ! complex(${k1}$) endpoints + real(${k1}$) base = complex(${k1}$) result
+      ${t1}$ :: res(n)
+    end function ${RName}$
+
+    #:set RName = rname("logspace", 1, t1, k1, "n_cbase")
+    module function ${RName}$(start, end, n, base) result(res)
+      ${t1}$, intent(in) :: start
+      ${t1}$, intent(in) :: end
+      integer, intent(in) :: n
+      complex(${k1}$), intent(in) :: base
+      ! complex(${k1}$) endpoints + complex(${k1}$) base = complex(${k1}$) result
+      ${t1}$ :: res(n)
+    end function ${RName}$
+
+    #:set RName = rname("logspace", 1, t1, k1, "n_ibase")
+    module function ${RName}$(start, end, n, base) result(res)
+      ${t1}$, intent(in) :: start
+      ${t1}$, intent(in) :: end
+      integer, intent(in) :: n
+      integer, intent(in) :: base
+      ! complex(${k1}$) endpoints + integer base = complex(${k1}$) result
+      ${t1}$ :: res(n)
+    end function ${RName}$
+  #:endfor
+  #! ========================================================
+  #! ========================================================
+  #! Provide support for Integer start/endpoints
+    ! Generate logarithmically spaced sequence from ${k1}$ base to the powers
+    ! of ${k1}$ start and end. [base^start, ... , base^end]
+    ! Different combinations of parameter types will lead to different result types.
+    ! Those combinations are indicated in the body of each function.
+  #:for k1 in REAL_KINDS
+    #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n_r" + str(k1) + "base")
+    module function ${RName}$(start, end, n, base) result(res)
+      integer, intent(in) :: start
+      integer, intent(in) :: end
+      integer, intent(in) :: n
+      real(${k1}$), intent(in) :: base
+      ! integer endpoints + real(${k1}$) base = real(${k1}$) result
+      real(${k1}$) :: res(n)
+    end function ${RName}$
+
+    #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n_c" + str(k1) + "base")
+    module function ${RName}$(start, end, n, base) result(res)
+      integer, intent(in) :: start
+      integer, intent(in) :: end
+      integer, intent(in) :: n
+      complex(${k1}$), intent(in) :: base
+      ! integer endpoints + complex(${k1}$) base = complex(${k1}$) result
+      complex(${k1}$) :: res(n)
+    end function ${RName}$
+  #:endfor
+
+    #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n_ibase")
+    module function ${RName}$(start, end, n, base) result(res)
+      integer, intent(in) :: start
+      integer, intent(in) :: end
+      integer, intent(in) :: n
+      integer, intent(in) :: base
+      ! integer endpoints + integer base = integer result
+      integer :: res(n)
+    end function ${RName}$
+
+
+  end interface
+
 contains
-    
+
     #:for k1, t1 in IR_KINDS_TYPES
     elemental function clip_${k1}$(x, xmin, xmax) result(res)
         ${t1}$, intent(in) :: x
@@ -25,6 +272,6 @@ contains
 
         res = max(min(x, xmax), xmin)
     end function clip_${k1}$
-    
+
     #:endfor
 end module stdlib_math

src/stdlib_math_linspace.fypp

@@ -0,0 +1,97 @@
+#:include "common.fypp"
+submodule (stdlib_math) stdlib_math_linspace
+
+implicit none
+
+contains
+
+  #:for k1, t1 in REAL_KINDS_TYPES
+    #:set RName = rname("linspace_default", 1, t1, k1)
+    module function ${RName}$(start, end) result(res)
+      ${t1}$, intent(in) :: start
+      ${t1}$, intent(in) :: end
+
+      ${t1}$ :: res(DEFAULT_LINSPACE_LENGTH)
+
+      res = linspace(start, end, DEFAULT_LINSPACE_LENGTH)
+
+    end function ${RName}$
+  #:endfor
+
+  #:for k1, t1 in REAL_KINDS_TYPES
+    #:set RName = rname("linspace_n", 1, t1, k1)
+    module function ${RName}$(start, end, n) result(res)
+      ${t1}$, intent(in) :: start
+      ${t1}$, intent(in) :: end
+      integer, intent(in) :: n
+
+      ${t1}$ :: res(n)
+
+      integer :: i    ! Looping index
+      ${t1}$ :: interval ! Difference between adjacent elements
+
+
+      if(n <= 0) return ! If passed length is less than or equal to 0, return an empty (allocated with length 0) array
+      if(n == 1) then
+        res(1) = end
+        return
+      end if
+
+      interval = (end - start) / real((n - 1), ${k1}$)
+
+      res(1) = start
+      res(n) = end
+
+      do i = 2, n - 1
+
+        res(i) = real((i-1), ${k1}$) * interval + start
+
+      end do
+
+    end function ${RName}$
+  #:endfor
+
+
+    #:for k1, t1 in CMPLX_KINDS_TYPES
+      #:set RName = rname("linspace_default", 1, t1, k1)
+      module procedure ${RName}$
+
+        res = linspace(start, end, DEFAULT_LINSPACE_LENGTH)
+
+      end procedure ${RName}$
+    #:endfor
+
+    #:for k1, t1 in CMPLX_KINDS_TYPES
+      #:set RName = rname("linspace_n", 1, t1, k1)
+      module procedure ${RName}$
+
+        real(${k1}$) :: x(n) ! array of the real part of complex number
+        real(${k1}$) :: y(n) ! array of the imaginary part of the complex number
+
+        x = linspace(start%re, end%re, n)
+        y = linspace(start%im, end%im, n)
+
+        res = cmplx(x, y, kind=${k1}$)
+
+      end procedure ${RName}$
+    #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("linspace_default", 1, t1, k1)
+      module procedure ${RName}$
+
+        res = linspace(real(start, kind=dp), real(end, kind=dp), DEFAULT_LINSPACE_LENGTH)
+
+      end procedure ${RName}$
+    #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("linspace_n", 1, t1, k1)
+      module procedure ${RName}$
+
+        res = linspace(real(start, kind=dp), real(end, kind=dp), n)
+
+      end procedure ${RName}$
+    #:endfor
+
+end submodule

src/stdlib_math_logspace.fypp

@@ -0,0 +1,93 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+
+submodule (stdlib_math) stdlib_math_logspace
+
+implicit none
+
+contains
+
+  #!=========================================================
+  #!=     logspace(start, end)                              =
+  #!=========================================================
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:set RName = rname("logspace", 1, t1, k1, "default")
+    module procedure ${RName}$
+      res = logspace(start, end, DEFAULT_LOGSPACE_LENGTH, real(DEFAULT_LOGSPACE_BASE, ${k1}$))
+    end procedure
+  #:endfor
+  #! Integer support
+  #:set RName = rname("logspace", 1, "integer(int32)", "int32", "default")
+    module procedure ${RName}$
+      res = logspace(start, end, DEFAULT_LOGSPACE_LENGTH, DEFAULT_LOGSPACE_BASE)
+    end procedure
+
+  #!=========================================================
+  #!=     logspace(start, end, n)                           =
+  #!=========================================================
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:set RName = rname("logspace", 1, t1, k1, "n")
+    module procedure ${RName}$
+      res = logspace(start, end, n, real(DEFAULT_LOGSPACE_BASE, ${k1}$))
+    end procedure
+  #:endfor
+  #! Integer support
+  #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n")
+    module procedure ${RName}$
+      res = logspace(start, end, n, DEFAULT_LOGSPACE_BASE)
+    end procedure
+
+  #!=========================================================
+  #!=     logspace(start, end, n, base)                     =
+  #!=========================================================
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:set RName = rname("logspace", 1, t1, k1, "n_rbase")
+    module procedure ${RName}$
+      ${t1}$ :: exponents(n)
+      exponents = linspace(start, end, n)
+      res = base ** exponents
+    end procedure
+
+    #:set RName = rname("logspace", 1, t1, k1, "n_cbase")
+    module procedure ${RName}$
+      ${t1}$ :: exponents(n)
+      exponents = linspace(start, end, n)
+      res = base ** exponents
+    end procedure
+
+    #:set RName = rname("logspace", 1, t1, k1, "n_ibase")
+    module procedure ${RName}$
+      ${t1}$ :: exponents(n)
+      exponents = linspace(start, end, n)
+      res = base ** exponents
+    end procedure
+  #:endfor
+  #! Integer support:
+    ! Generate logarithmically spaced sequence from ${k1}$ base to the powers
+    ! of ${k1}$ start and end. [base^start, ... , base^end]
+    ! RName = ${RName}$
+  #:for k1 in REAL_KINDS
+    #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n_r" + str(k1) + "base")
+    module procedure ${RName}$
+      integer :: exponents(n)
+      exponents = linspace(start, end, n)
+      res = base ** exponents
+    end procedure
+
+    #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n_c" + str(k1) + "base")
+    module procedure ${RName}$
+      integer :: exponents(n)
+      exponents = linspace(start, end, n)
+      res = base ** exponents
+    end procedure
+  #:endfor
+
+    #:set RName = rname("logspace", 1, "integer(int32)", "int32", "n_ibase")
+    module procedure ${RName}$
+      integer :: exponents(n)
+      exponents = linspace(start, end, n)
+      res = base ** exponents
+    end procedure
+
+
+end submodule

src/tests/math/CMakeLists.txt

@@ -1 +1,3 @@
 ADDTEST(stdlib_math)
+ADDTEST(linspace)
+ADDTEST(logspace)

src/tests/math/Makefile.manual

@@ -1,4 +1,4 @@
-PROGS_SRC = test_stdlib_math.f90
+PROGS_SRC = test_stdlib_math.f90 test_linspace.f90 test_logspace.f90
 
 
 include ../Makefile.manual.test.mk

src/tests/math/test_linspace.f90

@@ -0,0 +1,383 @@
+program test_linspace
+    use stdlib_error, only: check
+    use stdlib_kinds, only: sp, dp, int8, int16
+    use stdlib_math, only: linspace, DEFAULT_LINSPACE_LENGTH
+
+    implicit none
+
+    integer :: iunit
+    logical :: warn = .false.
+    real(sp), parameter :: TOLERANCE_SP = 1000 * epsilon(1.0_sp)
+    real(dp), parameter :: TOLERANCE_DP = 1000 * epsilon(1.0_dp) ! Percentage of the range for which the actual gap must not exceed
+
+    ! Testing linspace.
+    !
+    ! For single and double precision, check if the beginning and end values are properly recorded
+    ! and make sure that the size of the result array is as expected.
+    !
+    ! This testing suite makes use of the a repeated section of code that will check to make
+    ! sure that every element is linearly spaced (i.e., call check(|array(i+1) - array(i)| < |expected_value| * TOLERANCE)).
+    ! I would convert this repeated code into a subroutine but that would require the implementation of a
+    ! generic procedure given that each linear space will have a different expected_value type and kind.
+
+    open(newunit=iunit, file="test_linspace_log.txt", status="unknown") ! Log the results of the functions
+
+    write(iunit,*) "Writing to unit #: ", iunit
+
+    call test_linspace_sp
+    call test_linspace_dp
+    call test_linspace_neg_index ! Make sure that when passed a negative index the result is an empty array
+    call test_linspace_cmplx
+    call test_linspace_cmplx_2
+    call test_linspace_cmplx_3
+    call test_linspace_cmplx_sp
+    call test_linspace_cmplx_sp_2
+    call test_linspace_int16
+    call test_linspace_int8
+
+    close(unit=iunit)
+
+contains
+
+    subroutine test_linspace_sp
+
+        integer, parameter :: n = 20
+        real(sp) :: start = 1.0_sp
+        real(sp) :: end = 10.0_sp
+        real(sp) :: expected_interval
+        real(sp) :: true_difference
+
+        integer :: i
+        real(sp), allocatable :: x(:)
+
+        x = linspace(start, end, n)
+
+        expected_interval =( end - start ) / real(( n - 1 ), sp)
+
+        call check(x(1) == start, msg="Initial value of array is not equal to the passed start parameter", warn=warn)
+        call check(x(n) == end, msg="Final array value is not equal to end parameter", warn=warn)
+        call check(size(x) == n, msg="Array not allocated to appropriate size", warn=warn)
+
+        print *, "Made it through first round of tests"
+
+        ! Due to roundoff error, it is possible that the jump from x(n-1) to x(n) is slightly different than the expected interval
+        do i = 1, n-1
+
+            true_difference = x(i + 1) - x(i)
+            call check(abs(true_difference - expected_interval) < abs(expected_interval) * TOLERANCE_SP)
+
+        end do
+
+    end subroutine
+
+    subroutine test_linspace_dp
+
+        real(dp) :: start = 1.0_dp
+        real(dp) :: end = 10.0_dp
+        integer, parameter :: n = DEFAULT_LINSPACE_LENGTH
+        real(dp) :: expected_interval
+        real(dp) :: true_difference
+
+        real(dp), allocatable :: x(:)
+        integer :: i
+
+        x = linspace(start, end)
+
+        expected_interval =( end - start ) / ( n - 1 )
+
+        call check(x(1) == start, msg="Initial value of array is not equal to the passed start parameter", warn=warn)
+        call check(x(n) == end, msg="Final array value is not equal to end parameter", warn=warn)
+        call check(size(x) == n, msg="Array not allocated to default size", warn=warn)
+
+        ! Due to roundoff error, it is possible that the jump from x(n-1) to x(n) is slightly different than the expected interval
+        do i = 1, n-1
+
+            true_difference = x(i + 1) - x(i)
+            call check(abs(true_difference - expected_interval) < abs(expected_interval) * TOLERANCE_DP)
+
+        end do
+
+    end subroutine
+
+    subroutine test_linspace_neg_index
+
+        real(dp) :: start = 1.0_dp
+        real(dp) :: end = 10.0_dp
+
+        real(dp), allocatable :: x(:)
+
+        x = linspace(start, end, -15)
+
+        call check(size(x) == 0, msg="Allocated array is not empty", warn=warn)
+
+    end subroutine
+
+    subroutine test_linspace_cmplx
+
+        complex(dp) :: start = (0.0_dp, 10.0_dp)
+        complex(dp) :: end = (1.0_dp, 0.0_dp)
+        complex(dp) :: expected_interval
+        integer, parameter :: n = 10
+
+        complex(dp), allocatable :: z(:)
+
+        integer :: i
+
+        z = linspace(start, end, n)
+
+        expected_interval =( end - start ) / ( n - 1 )
+
+        call check(z(1) == start, msg="Initial value of array is not equal to the passed start parameter", warn=warn)
+        call check(z(n) == end, msg="Final array value is not equal to end parameter", warn=warn)
+        call check(size(z) == n, msg="Array not allocated to correct size", warn=warn)
+
+        ! Due to roundoff error, it is possible that the jump from x(n-1) to x(n) is slightly different than the expected interval
+        do i = 1, n-1
+
+            call check(abs( ( z(i + 1) - z(i) ) - expected_interval) < abs(expected_interval) * TOLERANCE_DP)
+
+        end do
+
+        write(unit=iunit, fmt=*) "linspace((0.0_dp, 10.0_dp), (1.0_dp, 0.0_dp), 10): "
+        write(unit=iunit,fmt='(70("="))')
+        do i = 1, n
+            write(unit=iunit,fmt=*) z(i)
+        end do
+        write(iunit,*)
+        write(iunit,*)
+
+    end subroutine
+
+    subroutine test_linspace_cmplx_2
+
+        complex(dp) :: start = (10.0_dp, 10.0_dp)
+        complex(dp) :: end = (1.0_dp, 1.0_dp)
+        complex(dp) :: expected_interval
+
+        integer, parameter :: n = 5
+
+        complex(dp), allocatable :: z(:)
+
+        integer :: i
+
+        z = linspace(start, end, n)
+
+        expected_interval =( end - start ) / ( n - 1 )
+
+        call check(z(1) == start, msg="Initial value of array is not equal to the passed start parameter", warn=warn)
+        call check(z(n) == end, msg="Final array value is not equal to end parameter", warn=warn)
+        call check(size(z) == n, msg="Array not allocated to correct size", warn=warn)
+
+        ! Due to roundoff error, it is possible that the jump from x(n-1) to x(n) is slightly different than the expected interval
+        do i = 1, n-1
+
+            call check(abs( ( z(i + 1) - z(i) ) - expected_interval) < abs(expected_interval) * TOLERANCE_DP)
+
+        end do
+
+        write(unit=iunit, fmt=*) "linspace((10.0_dp, 10.0_dp), (1.0_dp, 1.0_dp), 5): "
+        write(unit=iunit,fmt='(70("="))')
+        do i = 1, n
+            write(unit=iunit,fmt=*) z(i)
+        end do
+        write(iunit,*)
+        write(iunit,*)
+
+    end subroutine
+
+    subroutine test_linspace_cmplx_3
+
+        complex(dp) :: start = (-5.0_dp, 100.0_dp)
+        complex(dp) :: end = (20.0_dp, 13.0_dp)
+        complex(dp) :: expected_interval
+
+        integer, parameter :: n = 20
+
+        complex(dp), allocatable :: z(:)
+
+        integer :: i
+
+        z = linspace(start, end, n)
+
+        expected_interval = ( end - start ) / ( n - 1 )
+
+        call check(z(1) == start, msg="Initial value of array is not equal to the passed start parameter", warn=warn)
+        call check(z(n) == end, msg="Final array value is not equal to end parameter", warn=warn)
+        call check(size(z) == n, msg="Array not allocated to correct size", warn=warn)
+
+        ! Due to roundoff error, it is possible that the jump from x(n-1) to x(n) is slightly different than the expected interval
+        do i = 1, n-1
+
+            call check(abs( ( z(i + 1) - z(i) ) - expected_interval) < abs(expected_interval) * TOLERANCE_DP)
+
+        end do
+
+        write(unit=iunit, fmt=*) "linspace((-5.0_dp, 100.0_dp), (20.0_dp, 13.0_dp), 20): "
+        write(unit=iunit,fmt='(70("="))')
+        do i = 1, n
+            write(unit=iunit,fmt=*) z(i)
+        end do
+        write(iunit,*)
+        write(iunit,*)
+
+    end subroutine
+
+    subroutine test_linspace_cmplx_sp
+
+        complex(sp) :: start = (0.5_sp, 5.0_sp)
+        complex(sp) :: end = (1.0_sp, -30.0_sp)
+        complex(sp) :: expected_interval
+
+        integer, parameter :: n = 10
+
+        complex(sp), allocatable :: z(:)
+
+        integer :: i
+
+        z = linspace(start, end, n)
+
+        expected_interval =( end - start ) / ( n - 1 )
+
+        call check(z(1) == start, msg="Initial value of array is not equal to the passed start parameter", warn=warn)
+        call check(z(n) == end, msg="Final array value is not equal to end parameter", warn=warn)
+        call check(size(z) == n, msg="Array not allocated to correct size", warn=warn)
+
+        ! Due to roundoff error, it is possible that the jump from x(n-1) to x(n) is slightly different than the expected interval
+        do i = 1, n-1
+
+            call check(abs( ( z(i + 1) - z(i) ) - expected_interval) < abs(expected_interval) * TOLERANCE_SP)
+
+        end do
+
+        write(unit=iunit, fmt=*) "linspace((0.5_sp, 5.0_sp), (1.0_sp, -30.0_sp), 10): "
+        write(unit=iunit,fmt='(70("="))')
+        do i = 1, n
+            write(unit=iunit,fmt=*) z(i)
+        end do
+        write(iunit,*)
+        write(iunit,*)
+
+    end subroutine
+
+    subroutine test_linspace_cmplx_sp_2
+
+        complex(sp) :: start = (50.0_sp, 500.0_sp)
+        complex(sp) :: end = (-100.0_sp, 2000.0_sp)
+        complex(sp) :: expected_interval
+        complex(sp) :: true_interval
+        real(sp) :: offset
+
+        integer, parameter :: n = DEFAULT_LINSPACE_LENGTH
+
+        complex(sp), allocatable :: z(:)
+
+        integer :: i
+
+        z = linspace(start, end)
+
+        expected_interval =( end - start ) / ( n - 1 )
+
+        call check(z(1) == start, msg="Initial value of array is not equal to the passed start parameter", warn=warn)
+        call check(z(n) == end, msg="Final array value is not equal to end parameter", warn=warn)
+        call check(size(z) == n, msg="Array not allocated to default size", warn=warn)
+
+        ! Due to roundoff error, it is possible that the jump from x(n-1) to x(n) is slightly different than the expected interval
+        do i = 1, n-1
+
+            true_interval = (z(i + 1) - z(i))
+            offset = abs(true_interval - expected_interval)
+            call check(abs( ( z(i + 1) - z(i) ) - expected_interval) < abs(expected_interval) * TOLERANCE_SP)
+            ! print *, i
+
+        end do
+
+        write(unit=iunit, fmt=*) "linspace((50.0_sp, 500.0_sp), (-100.0_sp, 2000.0_sp)): "
+        write(unit=iunit,fmt='(70("="))')
+        do i = 1, n
+            write(unit=iunit,fmt=*) z(i)
+        end do
+        write(iunit,*)
+        write(iunit,*)
+
+    end subroutine
+
+    subroutine test_linspace_int16
+
+        integer(int16) :: start = 5
+        integer(int16) :: end = 10
+        real(dp) :: expected_interval
+
+        integer, parameter :: n = 6
+
+        integer(int16), allocatable :: z(:)
+
+        integer :: i
+
+        z = linspace(start, end, n)
+
+        expected_interval =( end - start ) / ( n - 1 )
+
+        call check(z(1) == start, msg="Initial value of array is not equal to the passed start parameter", warn=warn)
+        call check(z(n) == end, msg="Final array value is not equal to end parameter", warn=warn)
+        call check(size(z) == n, msg="Array not allocated to correct size", warn=warn)
+
+        ! Due to roundoff error, it is possible that the jump from x(n-1) to x(n) is slightly different than the expected interval
+        do i = 1, n-1
+
+            call check(abs( ( z(i + 1) - z(i) ) - expected_interval) < abs(expected_interval) * TOLERANCE_DP)
+
+        end do
+
+        write(unit=iunit, fmt=*) "linspace(5_int16, 10_int16, 10): "
+        write(unit=iunit,fmt='(70("="))')
+        do i = 1, n
+            write(unit=iunit,fmt=*) z(i)
+        end do
+        write(iunit,*)
+        write(iunit,*)
+
+    end subroutine
+
+    subroutine test_linspace_int8
+
+        integer(int8) :: start = 20
+        integer(int8) :: end = 50
+
+        real(dp) :: expected_interval
+
+        integer, parameter :: n = 10
+
+        real(dp), allocatable :: z(:)
+        integer(int8) :: z_int(n)
+
+        integer :: i
+
+        z = linspace(start, end, n)
+        z_int = linspace(start, end, n)
+
+        expected_interval =real( end - start, dp ) / ( n - 1 )
+
+        call check(z(1) == start, msg="Initial value of array is not equal to the passed start parameter", warn=warn)
+        call check(z(n) == end, msg="Final array value is not equal to end parameter", warn=warn)
+        call check(size(z) == n, msg="Array not allocated to correct size", warn=warn)
+
+        ! Due to roundoff error, it is possible that the jump from x(n-1) to x(n) is slightly different than the expected interval
+        do i = 1, n-1
+
+            call check(abs( ( z(i + 1) - z(i) ) - expected_interval) < abs(expected_interval) * TOLERANCE_DP)
+
+        end do
+
+        write(unit=iunit, fmt=*) "linspace(5_int16, 10_int16, 10): "
+        write(unit=iunit,fmt='(70("="))')
+        do i = 1, n
+            write(unit=iunit,fmt=*) z(i)
+        end do
+        write(iunit,*)
+        write(iunit,*)
+
+    end subroutine
+
+
+
+end program

src/tests/math/test_logspace.f90

@@ -0,0 +1,243 @@
+program test_logspace
+
+    use stdlib_error, only: check
+    use stdlib_kinds, only: sp, dp, int8, int16, int32, int64
+    use stdlib_math, only: logspace, DEFAULT_LOGSPACE_BASE, DEFAULT_LOGSPACE_LENGTH
+
+    implicit none
+
+        logical :: warn = .false.
+
+        integer :: iunit
+
+        ! Testing logspace
+        !
+        ! logspace should return a rank 1 array of values equally logarithmically spaced
+        ! from the base**start to base**end, using 10 as the base. If no length
+        ! is specified, return a rank 1 array with 50 elements.
+        !
+        ! Also test to verify that the proportion between adjacent elements is constant within
+        ! a certain tolerance
+
+        real(sp), parameter :: TOLERANCE_SP = 1000 * epsilon(1.0_sp)
+        real(dp), parameter :: TOLERANCE_DP = 1000 * epsilon(1.0_dp) ! Percentage of the range for which the actual gap must not exceed
+
+        open(newunit=iunit, file="test_logspace_log.txt", status="unknown") ! Log the results of the function
+
+        write(iunit,*) "Writing to unit #: ", iunit
+
+        call test_logspace_sp
+        call test_logspace_dp
+        call test_logspace_default
+        call test_logspace_base_2
+        call test_logspace_base_2_cmplx_start
+        call test_logspace_base_i_int_start
+
+        close(unit=iunit)
+
+    contains
+
+        subroutine test_logspace_sp
+
+            integer :: n = 20
+            real(sp) :: start = 0.0_sp
+            real(sp) :: end = 2.0_sp
+
+            real(sp) :: expected_proportion
+            integer :: i = 1
+
+            real(sp), allocatable :: x(:)
+
+            x = logspace(start, end, n)
+
+            expected_proportion = 10 ** ( ( end - start ) / ( n - 1 ) )
+
+            call check(x(1) == DEFAULT_LOGSPACE_BASE ** start, msg="Initial value of array is not equal to 10^start", warn=warn)
+            call check(x(n) == DEFAULT_LOGSPACE_BASE ** end, msg="Final value of array is not equal to 10^end", warn=warn)
+            call check(size(x) == n, msg="Array not allocated to appropriate size", warn=warn)
+
+            do i = 1, n-1
+
+                call check(abs(x(i + 1) / x(i) - expected_proportion) < abs(expected_proportion) * TOLERANCE_SP)
+
+            end do
+
+
+            write(unit=iunit, fmt=*) "logspace(0.0_sp, 2.0_sp, 20): "
+            write(unit=iunit,fmt='(70("="))')
+            write(unit=iunit,fmt="(20(F7.3, 2X))") x
+            write(iunit,*)
+            write(iunit,*)
+
+        end subroutine
+
+        subroutine test_logspace_dp
+
+            integer :: n = 10
+            real(dp) :: start = 1.0_dp
+            real(dp) :: end = 0.0_dp
+            real(dp) :: expected_proportion
+            integer :: i = 1
+
+            real(dp), allocatable :: x(:)
+
+            x = logspace(start, end, n)
+
+            expected_proportion = 10 ** ( ( end - start ) / ( n - 1 ) )
+
+
+            call check(x(1) == DEFAULT_LOGSPACE_BASE ** start, msg="Initial value of array is not equal to 10^start", warn=warn)
+            call check(x(n) == DEFAULT_LOGSPACE_BASE ** end, msg="Final value of array is not equal to 10^end", warn=warn)
+            call check(size(x) == n, msg="Array not allocated to appropriate size", warn=warn)
+
+            do i = 1, n-1
+
+                call check(abs(x(i + 1) / x(i) - expected_proportion) < abs(expected_proportion) * TOLERANCE_DP)
+
+            end do
+
+            write(unit=iunit, fmt=*) "logspace(1.0_dp, 0.0_dp, 10): "
+            write(unit=iunit,fmt=99)
+            write(unit=iunit,fmt="(10(F7.3, 2X))") x
+            write(iunit,*)
+            write(iunit,*)
+
+            99 format(70("="))
+
+        end subroutine
+
+        subroutine test_logspace_default
+
+            real(dp) :: start = 0.0_dp
+            real(dp) :: end = 1.0_dp
+            integer :: n = DEFAULT_LOGSPACE_LENGTH
+            real(dp) :: expected_proportion
+            integer :: i
+
+            real(dp), allocatable :: x(:)
+
+            x = logspace(start, end)
+
+            expected_proportion = 10 ** ( ( end - start ) / ( n - 1 ) )
+
+
+            call check(x(1) == DEFAULT_LOGSPACE_BASE ** start, msg="Initial value of array is not equal to 10^start", warn=warn)
+            call check(x(n) == DEFAULT_LOGSPACE_BASE ** end, msg="Final value of array is not equal to 10^end", warn=warn)
+            call check(size(x) == n, msg="Array not allocated to appropriate size", warn=warn)
+
+            do i = 1, n-1
+
+                call check(abs(x(i + 1) / x(i) - expected_proportion) < abs(expected_proportion) * TOLERANCE_DP)
+
+            end do
+
+            write(unit=iunit, fmt=*) "logspace(0.0_dp, 1.0_dp): "
+            write(unit=iunit,fmt='(70("="))')
+            write(unit=iunit,fmt="(50(F7.3, 2X))") x
+            write(iunit,*)
+            write(iunit,*)
+
+        end subroutine
+
+        subroutine test_logspace_base_2
+
+            integer :: n = 10
+            real(dp) :: start = 1.0_dp
+            real(dp) :: end = 10.0_dp
+            integer :: base = 2
+            integer :: i
+            real(dp) :: expected_proportion
+
+            real(dp), allocatable :: x(:)
+
+            x = logspace(start, end, n, base)
+
+            expected_proportion = 2 ** ( ( end - start ) / ( n - 1 ) )
+
+            call check(x(1) == base ** start, msg="Initial value of array is not equal to 2^start", warn=warn)
+            call check(x(n) == base ** end, msg="Final value of array is not equal to 2^end", warn=warn)
+            call check(size(x) == n, msg="Array not allocated to appropriate size", warn=warn)
+
+            do i = 1, n-1
+
+                call check(abs(x(i + 1) / x(i) - expected_proportion) < abs(expected_proportion) * TOLERANCE_DP)
+
+            end do
+
+            write(unit=iunit, fmt=*) "logspace(1.0_dp, 10.0_dp, 10, 2): "
+            write(unit=iunit,fmt='(70("="))')
+            write(unit=iunit,fmt="(10(F9.3, 2X))") x
+            write(iunit,*)
+            write(iunit,*)
+
+        end subroutine
+
+        subroutine test_logspace_base_2_cmplx_start
+
+            integer :: n = 10
+            complex(dp) :: start = (1, 0)
+            complex(dp) :: end = (0, 1)
+            integer :: base = 2
+            complex(dp) :: expected_proportion
+            integer :: i
+
+            complex(dp), allocatable :: x(:)
+
+            x = logspace(start, end, n, base)
+
+            expected_proportion = 2 ** ( ( end - start ) / ( n - 1 ) )
+
+
+            call check(x(1) == base ** start, msg="Initial value of array is not equal to 2^start", warn=warn)
+            call check(x(n) == base ** end, msg="Final value of array is not equal to 2^end", warn=warn)
+            call check(size(x) == n, msg="Array not allocated to appropriate size", warn=warn)
+
+            do i = 1, n-1
+
+                call check(abs(x(i + 1) / x(i) - expected_proportion) < abs(expected_proportion) * TOLERANCE_DP)
+
+            end do
+
+            write(unit=iunit, fmt=*) "logspace(1, i, 10, 2): "
+            write(unit=iunit,fmt='(70("="))')
+            write(unit=iunit,fmt="(10('(', F6.3, ',', 1X, F6.3, ')', 2X))") x
+            write(iunit,*)
+            write(iunit,*)
+
+        end subroutine
+
+        subroutine test_logspace_base_i_int_start
+
+            integer :: n = 5
+            integer :: start = 1
+            integer :: end = 5
+            complex(dp) :: base = (0, 1) ! i
+            complex(dp) :: expected_proportion
+            integer :: i = 1
+
+            complex(dp), allocatable :: x(:)
+
+            x = logspace(start, end, n, base)
+
+            expected_proportion = base ** ( ( end - start ) / ( n - 1 ) )
+
+            call check(x(1) == base ** start, msg="Initial value of array is not equal to 2^start", warn=warn)
+            call check(x(n) == base ** end, msg="Final value of array is not equal to 2^end", warn=warn)
+            call check(size(x) == n, msg="Array not allocated to appropriate size", warn=warn)
+
+            do i = 1, n-1
+
+                call check(abs(x(i + 1) / x(i) - expected_proportion) < abs(expected_proportion) * TOLERANCE_DP)
+
+            end do
+
+            write(unit=iunit, fmt=*) "logspace(1, 5, 5, i): "
+            write(unit=iunit,fmt='(70("="))')
+            write(unit=iunit,fmt="(10('(', F6.3, ',', 1X, F6.3, ')', 2X))") x
+            write(iunit,*)
+            write(iunit,*)
+
+        end subroutine
+
+
+    end program