diff --git a/.github/workflows/ci_windows.yml b/.github/workflows/ci_windows.yml index 66551fdd2..e55aa6e3b 100644 --- a/.github/workflows/ci_windows.yml +++ b/.github/workflows/ci_windows.yml @@ -59,7 +59,7 @@ jobs: - name: CTest run: PATH=$PATH:/mingw64/bin/ ctest --test-dir build --output-on-failure --parallel -V -LE quadruple_precision - - uses: actions/upload-artifact@v1 + - uses: actions/upload-artifact@v4 if: failure() with: name: WindowsCMakeTestlog diff --git a/.github/workflows/doc-deployment.yml b/.github/workflows/doc-deployment.yml index d48149407..f41817eeb 100644 --- a/.github/workflows/doc-deployment.yml +++ b/.github/workflows/doc-deployment.yml @@ -40,7 +40,7 @@ jobs: ford -r $(git describe --always) --debug ${MAYBE_SKIP_SEARCH} "${FORD_FILE}" - name: Upload Documentation - uses: actions/upload-artifact@v2 + uses: actions/upload-artifact@v4 with: name: FORD-API-docs path: ./API-doc/ diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index b03fdd192..5f956b234 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1382,3 +1382,110 @@ If `err` is not present, exceptions trigger an `error stop`. {!example/linalg/example_inverse_function.f90!} ``` +## `get_norm` - Computes the vector norm of a generic-rank array. + +### Status + +Experimental + +### Description + +This `pure subroutine` interface computes one of several vector norms of `real` or `complex` array \( A \), depending on +the value of the `order` input argument. \( A \) may be an array of any rank. + +Result `nrm` returns a `real`, scalar norm value for the whole array; if `dim` is specified, `nrm` is a rank n-1 +array with the same shape as \(A \) and dimension `dim` dropped, containing all norms evaluated along `dim`. + +### Syntax + +`call ` [[stdlib_linalg(module):get_norm(interface)]] `(a, nrm, order, [, dim, err])` + +### Arguments + +`a`: Shall be a rank-n `real` or `complex` array containing the data. It is an `intent(in)` argument. + +`nrm`: if `dim` is absent, shall be a scalar with the norm evaluated over all the elements of the array. Otherwise, an array of rank `n-1`, and a shape similar +to that of `a` with dimension `dim` dropped. + +`order`: Shall be an `integer` value or a `character` flag that specifies the norm type, as follows. It is an `intent(in)` argument. + +| Integer input | Character Input | Norm type | +|------------------|------------------|---------------------------------------------------------| +| `-huge(0)` | `'-inf', '-Inf'` | Minimum absolute value \( \min_i{ \left|a_i\right| } \) | +| `1` | `'1'` | 1-norm \( \sum_i{ \left|a_i\right| } \) | +| `2` | `'2'` | Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \) | +| `>=3` | `'3','4',...` | p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \) | +| `huge(0)` | `'inf', 'Inf'` | Maximum absolute value \( \max_i{ \left|a_i\right| } \) | + +`dim` (optional): Shall be a scalar `integer` value with a value in the range from `1` to `n`, where `n` is the rank of the array. It is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. If `err` is not present, the function is `pure`. + +### Return value + +By default, the return value `nrm` is a scalar, and contains the norm as evaluated over all elements of the generic-rank array \( A \). +If the optional `dim` argument is present, `nrm` is a rank `n-1` array with the same shape as \( A \) except +for dimension `dim`, that is collapsed. Each element of `nrm` contains the 1D norm of the elements of \( A \), +evaluated along dimension `dim` only. + +Raises `LINALG_ERROR` if the requested norm type is invalid. +Raises `LINALG_VALUE_ERROR` if any of the arguments has an invalid size. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_get_norm.f90!} +``` + +## `norm` - Computes the vector norm of a generic-rank array. + +### Status + +Experimental + +### Description + +This function computes one of several vector norms of `real` or `complex` array \( A \), depending on +the value of the `order` input argument. \( A \) may be an array of any rank. + +### Syntax + +`x = ` [[stdlib_linalg(module):norm(interface)]] `(a, order, [, dim, err])` + +### Arguments + +`a`: Shall be a rank-n `real` or `complex` array containing the data. It is an `intent(in)` argument. + +`order`: Shall be an `integer` value or a `character` flag that specifies the norm type, as follows. It is an `intent(in)` argument. + +| Integer input | Character Input | Norm type | +|------------------|------------------|---------------------------------------------------------| +| `-huge(0)` | `'-inf', '-Inf'` | Minimum absolute value \( \min_i{ \left|a_i\right| } \) | +| `1` | `'1'` | 1-norm \( \sum_i{ \left|a_i\right| } \) | +| `2` | `'2'` | Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \) | +| `>=3` | `'3','4',...` | p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \) | +| `huge(0)` | `'inf', 'Inf'` | Maximum absolute value \( \max_i{ \left|a_i\right| } \) | + +`dim` (optional): Shall be a scalar `integer` value with a value in the range from `1` to `n`, where `n` is the rank of the array. It is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. If `err` is not present, the function is `pure`. + +### Return value + +By default, the return value `x` is a scalar, and contains the norm as evaluated over all elements of the generic-rank array \( A \). +If the optional `dim` argument is present, `x` is a rank `n-1` array with the same shape as \( A \) except +for dimension `dim`, that is dropped. Each element of `x` contains the 1D norm of the elements of \( A \), +evaluated along dimension `dim` only. + +Raises `LINALG_ERROR` if the requested norm type is invalid. +Raises `LINALG_VALUE_ERROR` if any of the arguments has an invalid size. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_norm.f90!} +``` + + diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index c8c676043..20ca284f8 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -28,6 +28,8 @@ ADD_EXAMPLE(blas_gemv) ADD_EXAMPLE(lapack_getrf) ADD_EXAMPLE(lstsq1) ADD_EXAMPLE(lstsq2) +ADD_EXAMPLE(norm) +ADD_EXAMPLE(get_norm) ADD_EXAMPLE(solve1) ADD_EXAMPLE(solve2) ADD_EXAMPLE(solve3) diff --git a/example/linalg/example_get_norm.f90 b/example/linalg/example_get_norm.f90 new file mode 100644 index 000000000..f150a9033 --- /dev/null +++ b/example/linalg/example_get_norm.f90 @@ -0,0 +1,51 @@ +! Vector norm: demonstrate usage of the function interface +program example_get_norm + use stdlib_linalg, only: get_norm, linalg_state_type + implicit none + + real :: a(3,3), nrm, nrmd(3) + integer :: j + type(linalg_state_type) :: err + + ! a = [ -3.00000000 0.00000000 3.00000000 + ! -2.00000000 1.00000000 4.00000000 + ! -1.00000000 2.00000000 5.00000000 ] + a = reshape([(j-4,j=1,9)], [3,3]) + + print "(' a = [ ',3(g0,3x),2(/9x,3(g0,3x)),']')", transpose(a) + + ! Norm with integer input + call get_norm(a, nrm, 2) + print *, 'Euclidean norm = ',nrm ! 8.30662346 + + ! Norm with character input + call get_norm(a, nrm, '2') + print *, 'Euclidean norm = ',nrm ! 8.30662346 + + ! Euclidean norm of row arrays, a(i,:) + call get_norm(a, nrmd, 2, dim=2) + print *, 'Rows norms = ',nrmd ! 4.24264050 4.58257580 5.47722578 + + ! Euclidean norm of columns arrays, a(:,i) + call get_norm(a, nrmd, 2, dim=1) + print *, 'Columns norms = ',nrmd ! 3.74165750 2.23606801 7.07106781 + + ! Infinity norms + call get_norm(a, nrm, 'inf') + print *, 'maxval(||a||) = ',nrm ! 5.00000000 + + call get_norm(a, nrmd, 'inf', dim=2) + print *, 'maxval(||a(i,:)||) = ',nrmd ! 3.00000000 4.00000000 5.00000000 + + call get_norm(a, nrm, '-inf') + print *, 'minval(||a||) = ',nrm ! 0.00000000 + + call get_norm(a, nrmd, '-inf', dim=1) + print *, 'minval(||a(:,i)||) = ',nrmd ! 1.00000000 0.00000000 3.00000000 + + ! Catch Error: + ! [norm] returned Value Error: dimension 4 is out of rank for shape(a)= [3, 3] + call get_norm(a, nrmd, 'inf', dim=4, err=err) + print *, 'invalid: ',err%print() + +end program example_get_norm diff --git a/example/linalg/example_norm.f90 b/example/linalg/example_norm.f90 new file mode 100644 index 000000000..43dc7e3eb --- /dev/null +++ b/example/linalg/example_norm.f90 @@ -0,0 +1,40 @@ +! Vector norm: demonstrate usage of the function interface +program example_norm + use stdlib_linalg, only: norm, linalg_state_type + implicit none + + real :: a(3,3),na + integer :: j + type(linalg_state_type) :: err + + ! a = [ -3.00000000 0.00000000 3.00000000 + ! -2.00000000 1.00000000 4.00000000 + ! -1.00000000 2.00000000 5.00000000 ] + a = reshape([(j-4,j=1,9)], [3,3]) + + print "(' a = [ ',3(g0,3x),2(/9x,3(g0,3x)),']')", transpose(a) + + ! Norm with integer input + print *, 'Euclidean norm = ',norm(a, 2) ! 8.30662346 + + ! Norm with character input + print *, 'Euclidean norm = ',norm(a, '2') ! 8.30662346 + + ! Euclidean norm of row arrays, a(i,:) + print *, 'Rows norms = ',norm(a, 2, dim=2) ! 4.24264050 4.58257580 5.47722578 + + ! Euclidean norm of columns arrays, a(:,i) + print *, 'Columns norms = ',norm(a, 2, dim=1) ! 3.74165750 2.23606801 7.07106781 + + ! Infinity norms + print *, 'maxval(||a||) = ',norm(a, 'inf') ! 5.00000000 + print *, 'maxval(||a(i,:)||) = ',norm(a, 'inf', dim=2) ! 3.00000000 4.00000000 5.00000000 + print *, 'minval(||a||) = ',norm(a, '-inf') ! 0.00000000 + print *, 'minval(||a(:,i)||) = ',norm(a, '-inf', dim=1) ! 1.00000000 0.00000000 3.00000000 + + ! Catch Error: + ! [norm] returned Value Error: dimension 4 is out of rank for shape(a)= [3, 3] + print *, 'invalid: ',norm(a,'inf', dim=4, err=err) + print *, 'error = ',err%print() + +end program example_norm diff --git a/include/common.fypp b/include/common.fypp index 0d861aead..b1169fcdb 100644 --- a/include/common.fypp +++ b/include/common.fypp @@ -298,4 +298,105 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #:endif #:enddef +#! +#! Generates a list of loop variables +#! +#! Args: +#! varname(str): Name of the variable to be used as prefix +#! n (int): Number of loop variables to be created +#! offset (int): Optional index offset +#! +#! Returns: +#! Variable definition string +#! +#! E.g., +#! loop_variables('j', 5) +#! -> "j1, j2, j3, j4, j5 +#! +#:def loop_variables(varname, n, offset=0) + #:assert n > 0 + #:call join_lines(joinstr=", ") + #:for i in range(1, n + 1) + ${varname}$${i+offset}$ + #:endfor + #:endcall +#:enddef + +#! Generates an array shape specifier from an N-D array size +#! +#! Args: +#! name (str): Name of the original variable +#! rank (int): Rank of the original variable +#! offset(int): optional offset of the dimension loop (default = 0) +#! +#! Returns: +#! Array rank suffix string enclosed in braces +#! +#! E.g., +#! shape_from_array_size('mat', 5)}$ +#! -> (size(mat,1),size(mat,2),size(mat,3),size(mat,4),size(mat,5)) +#! shape_from_array_size('mat', 5, 2)}$ +#! -> (size(mat,3),size(mat,4),size(mat,5),size(mat,6),size(mat,7)) +#! +#:def shape_from_array_size(name, rank, offset=0) + #:assert rank > 0 + #:call join_lines(joinstr=", ", prefix="(", suffix=")") + #:for i in range(1, rank + 1) + size(${name}$,${i+offset}$) + #:endfor + #:endcall +#:enddef + +#! Generates an array shape specifier from an N-D array of sizes +#! +#! Args: +#! name (str): Name of the original variable +#! rank (int): Rank of the original variable +#! offset(int): optional offset of the dimension loop (default = 0) +#! +#! Returns: +#! Array rank suffix string enclosed in braces +#! +#! E.g., +#! shape_from_array_data('mat', 5)}$ +#! -> (1:mat(1),1:mat(2),1:mat(3),1:mat(4),1:mat(5)) +#! shape_from_array_data('mat', 5, 2)}$ +#! -> (1:mat(3),1:mat(4),1:mat(5),1:mat(6),1:mat(7)) +#! +#:def shape_from_array_data(name, rank, offset=0) + #:assert rank > 0 + #:call join_lines(joinstr=", ", prefix="(", suffix=")") + #:for i in range(1, rank + 1) + 1:${name}$(${i+offset}$) + #:endfor + #:endcall +#:enddef + +#! +#! Start a sequence of loop with indexed variables over an N-D array +#! +#! Args: +#! varname (str): Name of the variable to be used as prefix +#! matname (str): Name of the variable to be used as array +#! n (int): Number of nested loops to be created (1=innermost; n=outermost) +#! dim_offset (int): Optional dimension offset (1st loop is over dimension 1+dim_offset) +#! intent (str): Optional indentation. Default: 8 spaces +#! +#! +#:def loop_variables_start(varname, matname, n, dim_offset=0, indent=" "*8) + #:assert n > 0 + #:for i in range(1, n + 1) +${indent}$do ${varname}$${n+1+dim_offset-i}$ = lbound(${matname}$, ${n+1+dim_offset-i}$), ubound(${matname}$, ${n+1+dim_offset-i}$) + #:endfor +#:enddef + +#:def loop_variables_end(n, indent=" "*8) +#:assert n > 0 + #:call join_lines(joinstr="; ",prefix=indent) + #:for i in range(1, n + 1) + enddo + #:endfor + #:endcall +#:enddef + #:endmute diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ef11b642e..ec491eeb8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -31,6 +31,7 @@ set(fppFiles stdlib_linalg_solve.fypp stdlib_linalg_determinant.fypp stdlib_linalg_inverse.fypp + stdlib_linalg_norms.fypp stdlib_linalg_state.fypp stdlib_linalg_svd.fypp stdlib_linalg_cholesky.fypp diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 4574568e5..a26a10630 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -31,6 +31,8 @@ module stdlib_linalg public :: operator(.inv.) public :: lstsq public :: lstsq_space + public :: norm + public :: get_norm public :: solve public :: solve_lu public :: solve_lstsq @@ -1065,6 +1067,188 @@ module stdlib_linalg #:endfor end interface svdvals + + #! Allow for integer or character norm input: i.e., norm(a,2) or norm(a, '2') + #:set NORM_INPUT_TYPE = ["character(len=*)","integer(ilp)"] + #:set NORM_INPUT_SHORT = ["char","int"] + #:set NORM_INPUT_OPTIONS = list(zip(NORM_INPUT_TYPE,NORM_INPUT_SHORT)) + ! Vector norms: function interface + interface norm + !! version: experimental + !! + !! Computes the vector norm of a generic-rank array \( A \). + !! ([Specification](../page/specs/stdlib_linalg.html#norm-computes-the-vector-norm-of-a-generic-rank-array)) + !! + !!### Summary + !! Return one of several scalar norm metrics of a `real` or `complex` input array \( A \), + !! that can have any rank. For generic rank-n arrays, the scalar norm over the whole + !! array is returned by default. If `n>=2` and the optional input dimension `dim` is specified, + !! a rank `n-1` array is returned with dimension `dim` collapsed, containing all 1D array norms + !! evaluated along dimension `dim` only. + !! + !! + !!### Description + !! + !! This interface provides methods for computing the vector norm(s) of an array. + !! Supported data types include `real` and `complex`. + !! Input arrays may have generic rank from 1 to ${MAXRANK}$. + !! + !! Norm type input is mandatory, and it is provided via the `order` argument. + !! This can be provided as either an `integer` value or a `character` string. + !! Allowed metrics are: + !! - 1-norm \( \sum_i{ \left|a_i\right| } \): `order` = 1 or '1' + !! - Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \): `order` = 2 or '2' + !! - p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \): `integer` `order`, order>=3 + !! - Infinity norm \( \max_i{ \left|a_i\right| } \): order = huge(0) or 'inf' + !! - Minus-infinity norm \( \min_i{ \left|a_i\right| } \): order = -huge(0) or '-inf' + !! + !!### Example + !! + !!```fortran + !! + !! real(sp) :: a(3,3), na, rown(3) + !! a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3]) + !! + !! ! L2 norm: whole matrix + !! na = norm(a, 2) + !! + !! ! Infinity norm of each row + !! rown = norm(a, 'inf', dim=2) + !! + !!``` + !! + #:for rk,rt,ri in RC_KINDS_TYPES + #:for it,ii in NORM_INPUT_OPTIONS + !> Scalar norms: ${rt}$ + #:for rank in range(1, MAXRANK + 1) + pure module function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$(a, order) result(nrm) + !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ + ${rt}$, intent(in) :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Norm of the matrix. + real(${rk}$) :: nrm + end function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$ + module function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$(a, order, err) result(nrm) + !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ + ${rt}$, intent(in) :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Output state return flag. + type(linalg_state_type), intent(out) :: err + !> Norm of the matrix. + real(${rk}$) :: nrm + end function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$ + #:endfor + !> Array norms: ${rt}$ + #:for rank in range(2, MAXRANK + 1) + pure module function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, order, dim) result(nrm) + !> Input matrix a[..] + ${rt}$, intent(in), target :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Dimension the norm is computed along + integer(ilp), intent(in) :: dim + !> Norm of the matrix. (Same shape as `a`, with `dim` dropped). + real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$ + end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ + module function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$(a, order, dim, err) result(nrm) + !> Input matrix a[..] + ${rt}$, intent(in), target :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Dimension the norm is computed along + integer(ilp), intent(in) :: dim + !> Output state return flag. + type(linalg_state_type), intent(out) :: err + !> Norm of the matrix. (Same shape as `a`, with `dim` dropped). + real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$ + end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$ + #:endfor + #:endfor + #:endfor + end interface norm + + !> Vector norm: subroutine interface + interface get_norm + !! version: experimental + !! + !! Computes the vector norm of a generic-rank array \( A \). + !! ([Specification](../page/specs/stdlib_linalg.html#get-norm-computes-the-vector-norm-of-a-generic-rank-array)) + !! + !!### Summary + !! Subroutine interface that returns one of several scalar norm metrics of a `real` or `complex` + !! input array \( A \), that can have any rank. For generic rank-n arrays, the scalar norm over + !! the whole array is returned by default. If `n>=2` and the optional input dimension `dim` is + !! specified, a rank `n-1` array is returned with dimension `dim` collapsed, containing all 1D + !! array norms evaluated along dimension `dim` only. + !! + !! + !!### Description + !! + !! This `pure subroutine `interface provides methods for computing the vector norm(s) of an array. + !! Supported data types include `real` and `complex`. + !! Input arrays may have generic rank from 1 to ${MAXRANK}$. + !! + !! Norm type input is mandatory, and it is provided via the `order` argument. + !! This can be provided as either an `integer` value or a `character` string. + !! Allowed metrics are: + !! - 1-norm \( \sum_i{ \left|a_i\right| } \): `order` = 1 or '1' + !! - Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \): `order` = 2 or '2' + !! - p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \): `integer` `order`, order>=3 + !! - Infinity norm \( \max_i{ \left|a_i\right| } \): order = huge(0) or 'inf' + !! - Minus-infinity norm \( \min_i{ \left|a_i\right| } \): order = -huge(0) or '-inf' + !! + !!### Example + !! + !!```fortran + !! + !! real(sp) :: a(3,3), na, rown(3) + !! type(linalg_state_type) :: err + !! a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3]) + !! + !! ! L2 norm: whole matrix + !! call get_norm(a, na, 2) + !! + !! ! Infinity norms of each row, with error control + !! call get_norm(a, rown, 'inf', dim=2, err=err) + !! + !!``` + !! + #:for rk,rt,ri in RC_KINDS_TYPES + #:for it,ii in NORM_INPUT_OPTIONS + !> Scalar norms: ${rt}$ + #:for rank in range(1, MAXRANK + 1) + pure module subroutine norm_${rank}$D_${ii}$_${ri}$(a, nrm, order, err) + !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ + ${rt}$, intent(in), target :: a${ranksuffix(rank)}$ + !> Norm of the matrix. + real(${rk}$), intent(out) :: nrm + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out), optional :: err + end subroutine norm_${rank}$D_${ii}$_${ri}$ + #:endfor + !> Array norms: ${rt}$ + #:for rank in range(2, MAXRANK + 1) + pure module subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err) + !> Input matrix a[..] + ${rt}$, intent(in) :: a${ranksuffix(rank)}$ + !> Dimension the norm is computed along + integer(ilp), intent(in) :: dim + !> Norm of the matrix. (Same shape as `a`, with `dim` dropped). + real(${rk}$), intent(out) :: nrm${reduced_shape('a', rank, 'dim')}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out), optional :: err + end subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ + #:endfor + #:endfor + #:endfor + end interface get_norm + contains diff --git a/src/stdlib_linalg_blas.fypp b/src/stdlib_linalg_blas.fypp index 504a15108..2c0618e8d 100644 --- a/src/stdlib_linalg_blas.fypp +++ b/src/stdlib_linalg_blas.fypp @@ -9,6 +9,60 @@ module stdlib_linalg_blas implicit none(type,external) public + interface asum + !! ASUM takes the sum of the absolute values. +#ifdef STDLIB_EXTERNAL_BLAS + pure real(dp) function dasum( n, x, incx ) + import sp,dp,qp,ilp,lk + implicit none(type,external) + integer(ilp), intent(in) :: incx,n + real(dp), intent(in) :: x(*) + end function dasum +#else + module procedure stdlib_dasum +#endif +#ifdef STDLIB_EXTERNAL_BLAS + pure real(dp) function dzasum( n, x, incx ) + import sp,dp,qp,ilp,lk + implicit none(type,external) + integer(ilp), intent(in) :: incx,n + complex(dp), intent(in) :: x(*) + end function dzasum +#else + module procedure stdlib_dzasum +#endif +#:for rk,rt,ri in REAL_KINDS_TYPES +#:if not rk in ["sp","dp"] + module procedure stdlib_${ri}$asum +#:endif +#:endfor +#:for rk,rt,ri in CMPLX_KINDS_TYPES +#:if not rk in ["sp","dp"] + module procedure stdlib_${c2ri(ri)}$zasum +#:endif +#:endfor +#ifdef STDLIB_EXTERNAL_BLAS + pure real(sp) function sasum( n, x, incx ) + import sp,dp,qp,ilp,lk + implicit none(type,external) + integer(ilp), intent(in) :: incx,n + real(sp), intent(in) :: x(*) + end function sasum +#else + module procedure stdlib_sasum +#endif +#ifdef STDLIB_EXTERNAL_BLAS + pure real(sp) function scasum( n, x, incx ) + import sp,dp,qp,ilp,lk + implicit none(type,external) + integer(ilp), intent(in) :: incx,n + complex(sp), intent(in) :: x(*) + end function scasum +#else + module procedure stdlib_scasum +#endif + end interface asum + interface axpy !! AXPY constant times a vector plus a vector. #ifdef STDLIB_EXTERNAL_BLAS @@ -974,12 +1028,26 @@ module stdlib_linalg_blas #else module procedure stdlib_dnrm2 #endif +#ifdef STDLIB_EXTERNAL_BLAS + pure real(dp) function dznrm2( n, x, incx ) + import sp,dp,qp,ilp,lk + implicit none(type,external) + integer(ilp), intent(in) :: incx,n + complex(dp), intent(in) :: x(*) + end function dznrm2 +#else + module procedure stdlib_dznrm2 +#endif #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] module procedure stdlib_${ri}$nrm2 - #:endif #:endfor +#:for rk,rt,ri in CMPLX_KINDS_TYPES +#:if not rk in ["sp","dp"] + module procedure stdlib_${c2ri(ri)}$znrm2 +#:endif +#:endfor #ifdef STDLIB_EXTERNAL_BLAS pure real(sp) function snrm2( n, x, incx ) import sp,dp,qp,ilp,lk @@ -989,6 +1057,16 @@ module stdlib_linalg_blas end function snrm2 #else module procedure stdlib_snrm2 +#endif +#ifdef STDLIB_EXTERNAL_BLAS + pure real(sp) function scnrm2( n, x, incx ) + import sp,dp,qp,ilp,lk + implicit none(type,external) + integer(ilp), intent(in) :: incx,n + complex(sp), intent(in) :: x(*) + end function scnrm2 +#else + module procedure stdlib_scnrm2 #endif end interface nrm2 diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp new file mode 100644 index 000000000..8e3fba4a4 --- /dev/null +++ b/src/stdlib_linalg_norms.fypp @@ -0,0 +1,366 @@ +#:include "common.fypp" +#:set ALL_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + +#! Allow for integer or character norm input: i.e., norm(a,2) or norm(a, '2') +#:set INPUT_TYPE = ["character(len=*)","integer(ilp)"] +#:set INPUT_SHORT = ["char","int"] +#:set INPUT_OPTIONS = list(zip(INPUT_TYPE,INPUT_SHORT)) +! Vector norms +submodule(stdlib_linalg) stdlib_linalg_norms + use stdlib_linalg_constants + use stdlib_linalg_blas + use stdlib_linalg_lapack + use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & + LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR + use iso_c_binding, only: c_intptr_t,c_char,c_loc + implicit none(type,external) + + character(*), parameter :: this = 'norm' + + !> List of internal norm flags + integer(ilp), parameter :: NORM_ONE = 1_ilp + integer(ilp), parameter :: NORM_TWO = 2_ilp + integer(ilp), parameter :: NORM_POW_FIRST = 3_ilp + integer(ilp), parameter :: NORM_INF = +huge(0_ilp) ! infinity norm + integer(ilp), parameter :: NORM_POW_LAST = NORM_INF - 1_ilp + integer(ilp), parameter :: NORM_MINUSINF = -huge(0_ilp) + + interface parse_norm_type + module procedure parse_norm_type_integer + module procedure parse_norm_type_character + end interface parse_norm_type + + + interface stride_1d + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stride_1d_${ri}$ + #:endfor + end interface stride_1d + + contains + + !> Parse norm type from an integer user input + pure subroutine parse_norm_type_integer(order,norm_type,err) + !> User input value + integer(ilp), intent(in) :: order + !> Return value: norm type + integer(ilp), intent(out) :: norm_type + !> State return flag + type(linalg_state_type), intent(out) :: err + + select case (order) + case (1_ilp) + norm_type = NORM_ONE + case (2_ilp) + norm_type = NORM_TWO + case (3_ilp:NORM_POW_LAST) + norm_type = order + case (NORM_INF:) + norm_type = NORM_INF + case (:NORM_MINUSINF) + norm_type = NORM_MINUSINF + + case default + norm_type = NORM_ONE + err = linalg_state_type(this,LINALG_ERROR,'Input norm type ',order,' is not recognized.') + end select + + end subroutine parse_norm_type_integer + + pure subroutine parse_norm_type_character(order,norm_type,err) + !> User input value + character(len=*), intent(in) :: order + !> Return value: norm type + integer(ilp), intent(out) :: norm_type + !> State return flag + type(linalg_state_type), intent(out) :: err + + integer(ilp) :: int_order,read_err + + select case (order) + case ('inf','Inf','INF') + norm_type = NORM_INF + case ('-inf','-Inf','-INF') + norm_type = NORM_MINUSINF + case ('Euclidean','euclidean','EUCLIDEAN') + norm_type = NORM_TWO + case default + + ! Check if this input can be read as an integer + read(order,*,iostat=read_err) int_order + if (read_err/=0) then + ! Cannot read as an integer + norm_type = NORM_ONE + err = linalg_state_type(this,LINALG_ERROR,'Input norm type ',order,' is not recognized.') + else + call parse_norm_type_integer(int_order,norm_type,err) + endif + + end select + + end subroutine parse_norm_type_character + + #:for rk,rt,ri in ALL_KINDS_TYPES + + ! Compute stride of a 1d array + pure integer(ilp) function stride_1d_${ri}$(a) result(stride) + !> Input 1-d array + ${rt}$, intent(in), target :: a(:) + + integer(c_intptr_t) :: a1,a2 + + if (size(a,kind=ilp)<=1_ilp) then + stride = 1_ilp + else + a1 = transfer(c_loc(a(1)),a1) + a2 = transfer(c_loc(a(2)),a2) + stride = bit_size(0_c_char)*int(a2-a1, ilp)/storage_size(a, kind=ilp) + endif + + end function stride_1d_${ri}$ + + ! Private internal 1D implementation. This has to be used only internally, + ! when all inputs are already checked for correctness. + pure subroutine internal_norm_1D_${ri}$(sze, a, nrm, norm_request) + !> Input matrix length + integer(ilp), intent(in) :: sze + !> Input contiguous 1-d matrix a(*) + ${rt}$, intent(in) :: a(sze) + !> Norm of the matrix. + real(${rk}$), intent(out) :: nrm + !> Internal matrix request + integer(ilp), intent(in) :: norm_request + + integer(ilp) :: i + real(${rk}$) :: rorder + intrinsic :: abs, sum, sqrt, maxval, minval, conjg + + ! Initialize norm to zero + nrm = 0.0_${rk}$ + + select case(norm_request) + case(NORM_ONE) + nrm = asum(sze,a,incx=1_ilp) + case(NORM_TWO) + nrm = nrm2(sze,a,incx=1_ilp) + case(NORM_INF) + #:if rt.startswith('complex') + ! The default BLAS stdlib_i${ri}$amax uses |Re(.)|+|Im(.)|. Do not use it. + i = stdlib_i${ri}$max1(sze,a,incx=1_ilp) + #:else + i = stdlib_i${ri}$amax(sze,a,incx=1_ilp) + #:endif + nrm = abs(a(i)) + case(NORM_MINUSINF) + nrm = minval( abs(a) ) + case (NORM_POW_FIRST:NORM_POW_LAST) + rorder = 1.0_${rk}$ / norm_request + nrm = sum( abs(a) ** norm_request ) ** rorder + end select + + end subroutine internal_norm_1D_${ri}$ + + #:for it,ii in INPUT_OPTIONS + + !============================================== + ! Norms : any rank to scalar + !============================================== + + #:for rank in range(1, MAXRANK + 1) + + ! Pure function interface, with order specification. On error, the code will stop + pure module function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$(a, order) result(nrm) + !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ + ${rt}$, intent(in) :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Norm of the matrix. + real(${rk}$) :: nrm + + call norm_${rank}$D_${ii}$_${ri}$(a, nrm=nrm, order=order) + + end function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$ + + ! Function interface with output error + module function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$(a, order, err) result(nrm) + !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ + ${rt}$, intent(in) :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Output state return flag. + type(linalg_state_type), intent(out) :: err + !> Norm of the matrix. + real(${rk}$) :: nrm + + call norm_${rank}$D_${ii}$_${ri}$(a, nrm=nrm, order=order, err=err) + + end function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$ + + ! Internal implementation: ${rank}$-d + pure module subroutine norm_${rank}$D_${ii}$_${ri}$(a, nrm, order, err) + !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ + ${rt}$, intent(in), target :: a${ranksuffix(rank)}$ + !> Norm of the matrix. + real(${rk}$), intent(out) :: nrm + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out), optional :: err + + type(linalg_state_type) :: err_ + integer(ilp) :: sze,norm_request + real(${rk}$) :: rorder + intrinsic :: abs, sum, sqrt, maxval, minval, conjg + + sze = size(a,kind=ilp) + + ! Initialize norm to zero + nrm = 0.0_${rk}$ + + ! Check matrix size + if (sze<=0) then + err_ = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix shape: a=',shape(a,kind=ilp)) + call linalg_error_handling(err_,err) + return + end if + + ! Check norm request + call parse_norm_type(order,norm_request,err_) + if (err_%error()) then + call linalg_error_handling(err_,err) + return + endif + + ! Get norm + call internal_norm_1D_${ri}$(sze, a, nrm, norm_request) + call linalg_error_handling(err_,err) + + end subroutine norm_${rank}$D_${ii}$_${ri}$ + + #:endfor + + !==================================================================== + ! Norms : any rank to rank-1, with DIM specifier and ${ii}$ input + !==================================================================== + + #:for rank in range(2, MAXRANK + 1) + + ! Pure function interface with DIM specifier. On error, the code will stop + pure module function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, order, dim) result(nrm) + !> Input matrix a[..] + ${rt}$, intent(in), target :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Dimension to collapse by computing the norm w.r.t other dimensions + integer(ilp), intent(in) :: dim + !> Norm of the matrix. + real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$ + + call norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim) + + end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ + + ! Function interface with DIM specifier and output error state. + module function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$(a, order, dim, err) result(nrm) + !> Input matrix a[..] + ${rt}$, intent(in), target :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Dimension to collapse by computing the norm w.r.t other dimensions + integer(ilp), intent(in) :: dim + !> Output state return flag. + type(linalg_state_type), intent(out) :: err + !> Norm of the matrix. + real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$ + + call norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err) + + end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$ + + ! Internal implementation + pure module subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err) + !> Input matrix a[..] + ${rt}$, intent(in) :: a${ranksuffix(rank)}$ + !> Dimension to collapse by computing the norm w.r.t other dimensions + ! (dim must be defined before it is used for `nrm`) + integer(ilp), intent(in) :: dim + !> Norm of the matrix. + real(${rk}$), intent(out) :: nrm${reduced_shape('a', rank, 'dim')}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out), optional :: err + + type(linalg_state_type) :: err_ + integer(ilp) :: sze,lda,norm_request,${loop_variables('j',rank-1,1)}$ + logical :: contiguous_data + real(${rk}$) :: rorder + integer(ilp), dimension(${rank}$) :: spe,spack,perm,iperm + integer(ilp), dimension(${rank}$), parameter :: dim_range = [(lda,lda=1_ilp,${rank}$_ilp)] + ${rt}$, allocatable :: apack${ranksuffix(rank)}$ + intrinsic :: abs, sum, sqrt, norm2, maxval, minval, conjg + + ! Input matrix properties + sze = size (a,kind=ilp) + spe = shape(a,kind=ilp) + + ! Initialize norm to zero + nrm = 0.0_${rk}$ + + if (sze<=0) then + err_ = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix shape: a=',shape(a,kind=ilp)) + call linalg_error_handling(err_,err) + return + end if + + ! Check dimension choice + if (dim<1 .or. dim>${rank}$) then + err_ = linalg_state_type(this,LINALG_VALUE_ERROR,'dimension ',dim, & + 'is out of rank for shape(a)=',shape(a,kind=ilp)) + call linalg_error_handling(err_,err) + return + end if + + ! Check norm request + call parse_norm_type(order,norm_request,err_) + if (err_%error()) then + call linalg_error_handling(err_,err) + return + endif + + ! The norm's leading dimension + lda = spe(dim) + + ! Check if input column data is contiguous + contiguous_data = dim==1 + + ! Get packed data with the norm dimension as the first dimension + if (.not.contiguous_data) then + + ! Permute array to map dim to 1 + perm = [dim,pack(dim_range,dim_range/=dim)] + iperm(perm) = dim_range + spack = spe(perm) + apack = reshape(a, shape=spack, order=iperm) + +${loop_variables_start('j', 'apack', rank-1, 1," "*12)}$ + call internal_norm_1D_${ri}$(lda, apack(:, ${loop_variables('j',rank-1,1)}$), & + nrm(${loop_variables('j',rank-1,1)}$), norm_request) +${loop_variables_end(rank-1," "*12)}$ + + else + +${loop_variables_start('j', 'a', rank-1, 1," "*12)}$ + call internal_norm_1D_${ri}$(lda, a(:, ${loop_variables('j',rank-1,1)}$), & + nrm(${loop_variables('j',rank-1,1)}$), norm_request) +${loop_variables_end(rank-1," "*12)}$ + + endif + + end subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ + + #:endfor + #:endfor + #:endfor + +end submodule stdlib_linalg_norms diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 0cdbbaa3c..0f0036f2c 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -7,6 +7,7 @@ set( "test_linalg_solve.fypp" "test_linalg_inverse.fypp" "test_linalg_lstsq.fypp" + "test_linalg_norm.fypp" "test_linalg_determinant.fypp" "test_linalg_svd.fypp" "test_linalg_matrix_property_checks.fypp" @@ -19,6 +20,7 @@ ADDTEST(linalg_determinant) ADDTEST(linalg_eigenvalues) ADDTEST(linalg_matrix_property_checks) ADDTEST(linalg_inverse) +ADDTEST(linalg_norm) ADDTEST(linalg_solve) ADDTEST(linalg_lstsq) ADDTEST(linalg_svd) diff --git a/test/linalg/test_linalg_norm.fypp b/test/linalg/test_linalg_norm.fypp new file mode 100644 index 000000000..0ec88c2ec --- /dev/null +++ b/test/linalg/test_linalg_norm.fypp @@ -0,0 +1,300 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + +#! Generate an array rank suffix with the same fixed size for all dimensions. +#! +#! Args: +#! rank (int): Rank of the variable +#! size (int): Size along each dimension +#! +#! Returns: +#! Array rank suffix string (e.g. (4,4,4) if rank = 3 and size = 4) +#! +#:def fixedranksuffix(rank,size) +#{if rank > 0}#(${str(size) + (","+str(size)) * (rank - 1)}$)#{endif}# +#:enddef +! Test vector norms +module test_linalg_norm + use testdrive, only: error_type, check, new_unittest, unittest_type + use stdlib_linalg_constants + use stdlib_linalg, only: norm, linalg_state_type + use stdlib_linalg_state, only: linalg_state_type + + implicit none (type,external) + + contains + + !> Vector norm tests + subroutine test_vector_norms(tests) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + allocate(tests(0)) + + #:for rk,rt,ri in RC_KINDS_TYPES + tests = [tests,new_unittest("strided_1d_norm_${ri}$",test_strided_1d_${ri}$)] + #:for rank in range(1, MAXRANK) + tests = [tests,new_unittest("norm_${ri}$_${rank}$d",test_norm_${ri}$_${rank}$d)] + #:endfor + #:for rank in range(2, MAXRANK) + #:if rt.startswith('real') + tests = [tests,new_unittest("norm2_${ri}$_${rank}$d",test_norm2_${ri}$_${rank}$d)] + #:endif + tests = [tests,new_unittest("maxabs_${ri}$_${rank}$d",test_maxabs_${ri}$_${rank}$d)] + tests = [tests,new_unittest("norm_dimmed_${ri}$_${rank}$d",test_norm_dimmed_${ri}$_${rank}$d)] + #:endfor + #:endfor + + end subroutine test_vector_norms + + #:for rk,rt,ri in RC_KINDS_TYPES + + !> Test strided norm + subroutine test_strided_1d_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp), parameter :: m = 8_ilp + integer(ilp), parameter :: n = m**2 + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + ${rt}$, target :: a(n) + ${rt}$, allocatable :: slice(:) + ${rt}$, pointer :: twod(:,:) + real(${rk}$) :: rea(n),ima(n) + + call random_number(rea) + #:if rt.startswith('real') + a = rea + #:else + call random_number(ima) + a = cmplx(rea,ima,kind=${rk}$) + #:endif + + ! Test sliced array results + slice = a(4:7:59) + call check(error,abs(norm(a(4:7:59),2)-norm(slice,2)) a + call check(error,abs(norm(twod,2)-norm(a,2)) Test several norms with different dimensions + subroutine test_norm_${ri}$_${rank}$d(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp) :: j,order + integer(ilp), parameter :: n = 2_ilp**${rank}$ + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$ + character(64) :: msg + + allocate(a(n), b${fixedranksuffix(rank,2)}$) + + ! Init as a range,but with small elements such that all power norms will + ! never overflow, even in single precision + a = [(0.01_${rk}$*(j-n/2_ilp), j=1_ilp,n)] + b = reshape(a, shape(b)) + + ! Test some norms + do order = 1, 10 + write(msg,"('reshaped order-',i0,' p-norm is the same')") order + call check(error,abs(norm(a,order)-norm(b,order)) Test Euclidean norm; compare with Fortran intrinsic norm2 for reals + #:if rt.startswith('real') + subroutine test_norm2_${ri}$_${rank}$d(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp) :: j,dim + integer(ilp), parameter :: ndim = ${rank}$ + integer(ilp), parameter :: n = 2_ilp**ndim + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$ + intrinsic :: norm2 + character(64) :: msg + + allocate(a(n), b${fixedranksuffix(rank,2)}$) + + ! Init as a range,but with small elements such that all power norms will + ! never overflow, even in single precision + a = [(0.01_${rk}$*(j-n/2_ilp), j=1_ilp,n)] + b = reshape(a, shape(b)) + + ! Test some norms + call check(error,abs(norm(a,2) - norm2(a)) Test Infinity norm; compare with Fortran intrinsic max(abs(a)) + subroutine test_maxabs_${ri}$_${rank}$d(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp) :: j,dim + integer(ilp), parameter :: ndim = ${rank}$ + integer(ilp), parameter :: n = 2_ilp**ndim + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$ + intrinsic :: maxval, abs + character(128) :: msg + + allocate(a(n), b${fixedranksuffix(rank,2)}$) + + ! Init as a range,but with small elements such that all power norms will + ! never overflow, even in single precision + a = [(0.01_${rk}$*(j-n/2_ilp), j=1_ilp,n)] + b = reshape(a, shape(b)) + + ! Test some norms + call check(error,abs(norm(a,'inf') - maxval(abs(a))) 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_norm + +