diff --git a/doc/specs/stdlib_math.md b/doc/specs/stdlib_math.md
index 3e45b6696..2e09796de 100644
--- a/doc/specs/stdlib_math.md
+++ b/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
+```
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 9e8bee894..81f4a1cda 100644
--- a/src/CMakeLists.txt
+++ b/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
)
diff --git a/src/Makefile.manual b/src/Makefile.manual
index 09c7c27e2..2bb75c319 100644
--- a/src/Makefile.manual
+++ b/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
diff --git a/src/stdlib_math.fypp b/src/stdlib_math.fypp
index fc00d94f8..eae02ccc7 100644
--- a/src/stdlib_math.fypp
+++ b/src/stdlib_math.fypp
@@ -1,12 +1,24 @@
#: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
@@ -14,8 +26,243 @@ module stdlib_math
#: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
diff --git a/src/stdlib_math_linspace.fypp b/src/stdlib_math_linspace.fypp
new file mode 100644
index 000000000..655051b41
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/src/stdlib_math_logspace.fypp b/src/stdlib_math_logspace.fypp
new file mode 100644
index 000000000..87544ffcc
--- /dev/null
+++ b/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
diff --git a/src/tests/math/CMakeLists.txt b/src/tests/math/CMakeLists.txt
index ed5f32894..c4335e9ce 100644
--- a/src/tests/math/CMakeLists.txt
+++ b/src/tests/math/CMakeLists.txt
@@ -1 +1,3 @@
ADDTEST(stdlib_math)
+ADDTEST(linspace)
+ADDTEST(logspace)
\ No newline at end of file
diff --git a/src/tests/math/Makefile.manual b/src/tests/math/Makefile.manual
index de5f87d26..209990732 100644
--- a/src/tests/math/Makefile.manual
+++ b/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
diff --git a/src/tests/math/test_linspace.f90 b/src/tests/math/test_linspace.f90
new file mode 100644
index 000000000..25e5b9202
--- /dev/null
+++ b/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
diff --git a/src/tests/math/test_logspace.f90 b/src/tests/math/test_logspace.f90
new file mode 100644
index 000000000..57196ead6
--- /dev/null
+++ b/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
\ No newline at end of file