From 5532ded5f55017397fbf4c99071165ecfd7ff955 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 5 Nov 2024 05:41:50 -0600 Subject: [PATCH 01/23] interpret `*LANGE` inputs --- src/stdlib_linalg_norms.fypp | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 8e3fba4a4..82afeaec3 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -24,6 +24,12 @@ submodule(stdlib_linalg) stdlib_linalg_norms 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) + + !> Wrappers to LAPACK *LANGE matrix norm flags + character, parameter :: LANGE_NORM_MAT = 'M' ! maxval(sum(abs(a))) ! over whole matrix: unused + character, parameter :: LANGE_NORM_ONE = '1' ! maxval(sum(abs(a),1)) ! over columns + character, parameter :: LANGE_NORM_INF = 'I' ! maxval(sum(abs(a),2)) ! over rows + character, parameter :: LANGE_NORM_TWO = 'E' ! "Euclidean" or "Frobenius" interface parse_norm_type module procedure parse_norm_type_integer @@ -100,6 +106,27 @@ submodule(stdlib_linalg) stdlib_linalg_norms end subroutine parse_norm_type_character + !> From a user norm request, generate a *LANGE task command + pure subroutine lange_task_request(norm_type,lange_task,err) + !> Parsed matrix norm type + integer(ilp), intent(in) :: norm_type + !> LANGE task + character, intent(out) :: lange_task + !> Error flag + type(linalg_state_type), intent(inout) :: err + + select case (norm_type) + case (NORM_INF) + lange_task = LANGE_NORM_INF + case (NORM_ONE) + lange_task = LANGE_NORM_ONE + case (NORM_TWO) + lange_task = LANGE_NORM_TWO + case default + err = linalg_state_type(this,LINALG_VALUE_ERROR,'Order ',norm_type,' is not a valid matrix norm input.') + end select + end subroutine lange_task_request + #:for rk,rt,ri in ALL_KINDS_TYPES ! Compute stride of a 1d array From 7b3b14059a99682d614ae8f7e8fa3060ac300a75 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 5 Nov 2024 05:44:28 -0600 Subject: [PATCH 02/23] introduce internal implementation --- src/stdlib_linalg_norms.fypp | 156 +++++++++++++++++++++++++++++++++++ 1 file changed, 156 insertions(+) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 82afeaec3..08fcc8aa9 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -387,6 +387,162 @@ ${loop_variables_end(rank-1," "*12)}$ end subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ #:endfor + + !==================================================================== + ! Matrix norms + !==================================================================== + + ! Internal implementation + module function matrix_norm_${ii}$_${ri}$(a, order, err) result(nrm) + !> Input matrix a(m,n) + ${rt}$, intent(in), target :: a(:,:) + !> Norm of the matrix. + real(${rk}$) :: 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) :: m,n,norm_request + character :: lange_task + real(${rk}$), target :: work1(1) + real(${rk}$), pointer :: work(:) + + m = size(a,dim=1,kind=ilp) + n = size(a,dim=2,kind=ilp) + + ! Initialize norm to zero + nrm = 0.0_${rk}$ + + if (m<=0 .or. n<=0) then + err_ = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix shape: a=',[m,n]) + call linalg_error_handling(err_,err) + return + end if + + ! Check norm request: user + *LANGE support + call parse_norm_type(order,norm_request,err_) + call lange_task_request(norm_request,lange_task,err_) + if (err_%error()) then + call linalg_error_handling(err_,err) + return + endif + + if (lange_task==LANGE_NORM_INF) then + allocate(work(m)) + else + work => work1 + endif + + ! LAPACK interface + nrm = lange(lange_task,m,n,a,m,work) + + if (lange_task==LANGE_NORM_INF) deallocate(work) + + end function matrix_norm_${ii}$_${ri}$ + + #:for rank in range(3, MAXRANK + 1) + + ! Pure function interface with DIM specifier. On error, the code will stop + module function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$(a, order, dim, err) result(nrm) + !> Input matrix a(m,n) + ${rt}$, intent(in), contiguous, target :: a${ranksuffix(rank)}$ + !> Norm of the matrix. + real(${rk}$), allocatable :: nrm${ranksuffix(rank-2)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> [optional] dimensions of the sub-matrices the norms should be evaluated at (default = [1,2]) + integer(ilp), optional, intent(in) :: dim(2) + !> [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) :: j,m,n,lda,dims(2),norm_request + integer(ilp), dimension(${rank}$) :: s,spack,perm + integer(ilp), dimension(${rank}$), parameter :: dim_range = [(m,m=1_ilp,${rank}$_ilp)] + integer(ilp) :: ${loop_variables('j',rank-2,2)}$ + logical :: contiguous_data + character :: lange_task + real(${rk}$), target :: work1(1) + real(${rk}$), pointer :: work(:) + ${rt}$, pointer :: apack${ranksuffix(rank)}$ + + ! Get dimensions + if (present(dim)) then + dims = dim + else + dims = [1,2] + endif + + nullify(apack) + + if (dims(1)==dims(2) .or. .not.all(dims>0 .and. dims<=${rank}$)) then + err_ = linalg_state_type(this,LINALG_VALUE_ERROR,'Rank-',${rank}$,' matrix norm has invalid dim=',dims) + allocate(nrm${emptyranksuffix(rank-2)}$) + call linalg_error_handling(err_,err) + return + endif + + ! Check norm request: user + *LANGE support + call parse_norm_type(order,norm_request,err_) + call lange_task_request(norm_request,lange_task,err_) + if (err_%error()) then + call linalg_error_handling(err_,err) + return + endif + + ! Input matrix properties + s = shape(a,kind=ilp) + + ! Check if input column data is contiguous + contiguous_data = dims(1)==1 + + ! Matrix norm size + m = s(dims(1)) + n = s(dims(2)) + + ! Get packed data with norm dimensions as 1:2 + if (contiguous_data) then + + ! Collapse everything before the 1st dimension as apack's dim #1 + ! Set size==1 for all unused trailing dimensions + spack = [product(s(1:dims(2)-1)),s(dims(2):),(1_ilp,j=1,dims(2)-2)] + + ! Reshape without moving data + apack${shape_from_array_data('spack',rank)}$ => a + + else + + ! Dimension permutations to map dims(1),dims(2) => 1:2 + perm = [dims,pack(dim_range, dim_range/=dims(1) .and. dim_range/=dims(2))] + spack = s(perm) + apack = reshape(a, shape=spack, order=perm) + + endif + + if (lange_task==LANGE_NORM_INF) then + allocate(work(m)) + else + work => work1 + endif + + ! Allocate norm + allocate(nrm${shape_from_array_size('apack',rank-2, 2)}$) + + lda = size(apack,dim=1,kind=ilp) + + ! LAPACK interface +${loop_variables_start('j', 'apack', rank-2, 2)}$ + nrm(${loop_variables('j',rank-2,2)}$) = & + lange(lange_task,m,n,apack(:,:,${loop_variables('j',rank-2,2)}$),lda,work) +${loop_variables_end(rank-2)}$ + + if (lange_task==LANGE_NORM_INF) deallocate(work) + if (.not.contiguous_data) deallocate(apack) + + end function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$ + #:endfor #:endfor From c2beec82ace1ce57831cbc1e34ffbf2db5ef78c1 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 5 Nov 2024 05:57:08 -0600 Subject: [PATCH 03/23] introduce module interface --- src/stdlib_linalg.fypp | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index ad36ad677..0aee323e7 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -32,6 +32,7 @@ module stdlib_linalg public :: lstsq public :: lstsq_space public :: norm + public :: mnorm public :: get_norm public :: solve public :: solve_lu @@ -1320,6 +1321,40 @@ module stdlib_linalg #:endfor end interface get_norm + !> Matrix norm: function interface + interface mnorm + #:for rk,rt,ri in RC_KINDS_TYPES + #:for it,ii in NORM_INPUT_OPTIONS + module function matrix_norm_${ii}$_${ri}$(a, order, err) result(nrm) + !> Input matrix a(m,n) + ${rt}$, intent(in), target :: a(:,:) + !> Norm of the matrix. + real(${rk}$) :: 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 function matrix_norm_${ii}$_${ri}$ + + #:for rank in range(3, MAXRANK + 1) + !> N-D array with optional `dim` specifier + module function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$(a, order, dim, err) result(nrm) + !> Input matrix a(m,n) + ${rt}$, intent(in), contiguous, target :: a${ranksuffix(rank)}$ + !> Norm of the matrix. + real(${rk}$), allocatable :: nrm${ranksuffix(rank-2)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> [optional] dimensions of the sub-matrices the norms should be evaluated at (default = [1,2]) + integer(ilp), optional, intent(in) :: dim(2) + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out), optional :: err + end function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$ + #:endfor + #:endfor + #:endfor + end interface mnorm + contains From e3c85bfe87e62407f97dba7debdd38748f0ae677 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 5 Nov 2024 05:57:17 -0600 Subject: [PATCH 04/23] fix `reshape` permutation index --- src/stdlib_linalg_norms.fypp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 08fcc8aa9..e07be6a77 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -8,8 +8,8 @@ ! Vector norms submodule(stdlib_linalg) stdlib_linalg_norms use stdlib_linalg_constants - use stdlib_linalg_blas - use stdlib_linalg_lapack + use stdlib_linalg_blas, only: nrm2 + use stdlib_linalg_lapack, only: lange 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 @@ -459,7 +459,7 @@ ${loop_variables_end(rank-1," "*12)}$ type(linalg_state_type) :: err_ integer(ilp) :: j,m,n,lda,dims(2),norm_request - integer(ilp), dimension(${rank}$) :: s,spack,perm + integer(ilp), dimension(${rank}$) :: s,spack,perm,iperm integer(ilp), dimension(${rank}$), parameter :: dim_range = [(m,m=1_ilp,${rank}$_ilp)] integer(ilp) :: ${loop_variables('j',rank-2,2)}$ logical :: contiguous_data @@ -515,9 +515,10 @@ ${loop_variables_end(rank-1," "*12)}$ else ! Dimension permutations to map dims(1),dims(2) => 1:2 - perm = [dims,pack(dim_range, dim_range/=dims(1) .and. dim_range/=dims(2))] + perm = [dims,pack(dim_range, dim_range/=dims(1) .and. dim_range/=dims(2))] + iperm(perm) = dim_range spack = s(perm) - apack = reshape(a, shape=spack, order=perm) + apack = reshape(a, shape=spack, order=iperm) endif @@ -543,6 +544,8 @@ ${loop_variables_end(rank-2)}$ end function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$ + #:endfor + #:endfor #:endfor From d4d1e858dbd7dfc0852d73280b848b88234dd764 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 5 Nov 2024 14:26:34 +0100 Subject: [PATCH 05/23] docs: add interface explanation --- src/stdlib_linalg.fypp | 52 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 50 insertions(+), 2 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 0aee323e7..a413cc5f5 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1321,10 +1321,58 @@ module stdlib_linalg #:endfor end interface get_norm - !> Matrix norm: function interface + !> Matrix norms: function interface interface mnorm + !! version: experimental + !! + !! Computes the matrix norm of a generic-rank array \( A \). + !! ([Specification](../page/specs/stdlib_linalg.html#mnorm-computes-the-matrix-norm-of-a-generic-rank-array)) + !! + !!### Summary + !! Return one of several matrix norm metrics of a `real` or `complex` input array \( A \), + !! that can have rank 2 or higher. For rank-2 arrays, the matrix norm is returned. + !! If rank>2 and the optional input dimensions `dim` are specified, + !! a rank `n-2` array is returned with dimensions `dim(1),dim(2)` collapsed, containing all + !! matrix norms evaluated over the specified dimensions only. + !! + !!### Description + !! + !! This interface provides methods for computing the matrix norm(s) of an array. + !! Supported data types include `real` and `complex`. + !! Input arrays must have rank >= 2. + !! + !! 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: `order` = 1 or '1' + !! - 2-norm/Euclidean: `order` = 2 or '2' or 'Euclidean' + !! - p-norm: `order` >= 3 + !! - Infinity norm: `order` = huge(0) or 'Inf' + !! - Minus-infinity norm: `order` = '-Inf' + !! + !! If an invalid norm type is provided, the routine defaults to 1-norm and returns an error state. + !! + !!### Example + !! + !!```fortran + !! real(sp) :: a(3,3), na + !! real(sp) :: b(3,3,4), nb(4) ! Array of 4 3x3 matrices + !! a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3]) + !! + !! ! 2-norm/Euclidean norm of single matrix + !! na = mnorm(a, 'Euclidean') + !! + !! ! 1-norm of each 3x3 matrix in b + !! nb = mnorm(b, 1, dim=[1,2]) + !! + !! ! p-norm with p=3 + !! na = mnorm(a, 3) + !!``` + !! #:for rk,rt,ri in RC_KINDS_TYPES #:for it,ii in NORM_INPUT_OPTIONS + + !> Matrix norms: ${rt}$ rank-2 arrays module function matrix_norm_${ii}$_${ri}$(a, order, err) result(nrm) !> Input matrix a(m,n) ${rt}$, intent(in), target :: a(:,:) @@ -1336,8 +1384,8 @@ module stdlib_linalg type(linalg_state_type), intent(out), optional :: err end function matrix_norm_${ii}$_${ri}$ + !> Matrix norms: ${rt}$ higher rank arrays #:for rank in range(3, MAXRANK + 1) - !> N-D array with optional `dim` specifier module function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$(a, order, dim, err) result(nrm) !> Input matrix a(m,n) ${rt}$, intent(in), contiguous, target :: a${ranksuffix(rank)}$ From 48de116b7840d32ec7f0d7ecf1a10e1601fcfc08 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 5 Nov 2024 14:40:43 +0100 Subject: [PATCH 06/23] add specs --- doc/specs/stdlib_linalg.md | 48 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 2f03f0fad..6df881ad8 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1565,4 +1565,52 @@ If `err` is not present, exceptions trigger an `error stop`. {!example/linalg/example_norm.f90!} ``` +## `mnorm` - Computes the matrix norm of a generic-rank array. + +### Status + +Experimental + +### Description + +This function computes one of several matrix norms of `real` or `complex` array \( A \), depending on +the value of the `order` input argument. \( A \) must be an array of rank 2 or higher. For arrays of rank > 2, +matrix norms are computed over specified dimensions. + +### Syntax + +`x = ` [[stdlib_linalg(module):mnorm(interface)]] `(a, order [, dim, err])` + +### Arguments + +`a`: Shall be a rank-n `real` or `complex` array containing the data, where n ³ 2. 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 | +|------------------|---------------------------|---------------------------------------------------------| +| `1` | `'1'` | 1-norm (maximum column sum) \( \max_j \sum_i{ \left|a_{i,j}\right| } \) | +| `2` | `'2','Euclidean'` | 2-norm/Euclidean/Frobenius norm \( \sqrt{\sum_{i,j}{ \left|a_{i,j}\right|^2 }} \) | +| `huge(0)` | `'inf', 'Inf', 'INF'` | Infinity norm (maximum row sum) \( \max_i \sum_j{ \left|a_{i,j}\right| } \) | + +`dim` (optional): For arrays of rank > 2, shall be an integer array of size 2 specifying the dimensions over which to compute the matrix norm. Default value is `[1,2]`. It is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +For rank-2 input arrays, the return value `x` is a scalar containing the matrix norm. +For arrays of rank > 2, if the optional `dim` argument is present, `x` is a rank `n-2` array with the same shape as \( A \) except +for dimensions `dim(1)` and `dim(2)`, which are dropped. Each element of `x` contains the matrix norm of the corresponding submatrix of \( A \), +evaluated over the specified dimensions only, with the given order. + +If an invalid norm type is provided, defaults to 1-norm and raises `LINALG_ERROR`. +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_mnorm.f90!} +``` From 04a61b7d80c837f3be276266eb44b42fde3ce359 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 5 Nov 2024 14:40:50 +0100 Subject: [PATCH 07/23] add example --- example/linalg/CMakeLists.txt | 1 + example/linalg/example_mnorm.f90 | 24 ++++++++++++++++++++++++ 2 files changed, 25 insertions(+) create mode 100644 example/linalg/example_mnorm.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index e12236375..8c34b300b 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -29,6 +29,7 @@ ADD_EXAMPLE(lapack_getrf) ADD_EXAMPLE(lstsq1) ADD_EXAMPLE(lstsq2) ADD_EXAMPLE(norm) +ADD_EXAMPLE(mnorm) ADD_EXAMPLE(get_norm) ADD_EXAMPLE(solve1) ADD_EXAMPLE(solve2) diff --git a/example/linalg/example_mnorm.f90 b/example/linalg/example_mnorm.f90 new file mode 100644 index 000000000..24b2036e2 --- /dev/null +++ b/example/linalg/example_mnorm.f90 @@ -0,0 +1,24 @@ +program example_mnorm + use stdlib_linalg, only: mnorm + use stdlib_kinds, only: sp + implicit none + real(sp) :: a(3,3), na + real(sp) :: b(3,3,4), nb(4) ! Array of 4 3x3 matrices + + ! Initialize example matrix + a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3]) + + ! Compute Euclidean norm of single matrix + na = mnorm(a, 'Euclidean') + print *, "Euclidean norm of matrix a:", na + + ! Initialize array of matrices + b(:,:,1) = a + b(:,:,2) = 2*a + b(:,:,3) = 3*a + b(:,:,4) = 4*a + + ! Compute infinity norm of each 3x3 matrix in b + nb = mnorm(b, 'inf', dim=[1,2]) + print *, "Infinity norms of matrices in b:", nb +end program example_mnorm From 16657dfdf184cb7826064b4cd228e3916c9b701e Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 5 Nov 2024 14:56:37 +0100 Subject: [PATCH 08/23] fix `use` --- src/stdlib_linalg_norms.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index e07be6a77..ac417c63d 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -8,8 +8,8 @@ ! Vector norms submodule(stdlib_linalg) stdlib_linalg_norms use stdlib_linalg_constants - use stdlib_linalg_blas, only: nrm2 - use stdlib_linalg_lapack, only: lange + 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 From 53de5e6166187d2aed0a2695306aa776aa3b2c07 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 5 Nov 2024 15:10:02 +0100 Subject: [PATCH 09/23] first test program --- test/linalg/CMakeLists.txt | 2 + test/linalg/test_linalg_mnorm.fypp | 122 +++++++++++++++++++++++++++++ 2 files changed, 124 insertions(+) create mode 100644 test/linalg/test_linalg_mnorm.fypp diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index f835872df..a723b6ae1 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -8,6 +8,7 @@ set( "test_linalg_inverse.fypp" "test_linalg_lstsq.fypp" "test_linalg_norm.fypp" + "test_linalg_mnorm.fypp" "test_linalg_determinant.fypp" "test_linalg_qr.fypp" "test_linalg_svd.fypp" @@ -22,6 +23,7 @@ ADDTEST(linalg_eigenvalues) ADDTEST(linalg_matrix_property_checks) ADDTEST(linalg_inverse) ADDTEST(linalg_norm) +ADDTEST(linalg_mnorm) ADDTEST(linalg_solve) ADDTEST(linalg_lstsq) ADDTEST(linalg_qr) diff --git a/test/linalg/test_linalg_mnorm.fypp b/test/linalg/test_linalg_mnorm.fypp new file mode 100644 index 000000000..c79324373 --- /dev/null +++ b/test/linalg/test_linalg_mnorm.fypp @@ -0,0 +1,122 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +module test_linalg_mnorm + use testdrive, only: error_type, check, new_unittest, unittest_type + use stdlib_linalg_constants + use stdlib_linalg, only: mnorm, linalg_state_type + use stdlib_linalg_state, only: linalg_state_type + + implicit none (type,external) + + contains + + !> Matrix norm tests + subroutine test_matrix_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("test_matrix_norms_${ri}$",test_matrix_norms_${ri}$)] + #:endfor + + end subroutine test_matrix_norms + + #:for rk,rt,ri in RC_KINDS_TYPES + !> Test 1-norm, 2-norm (Euclidean), and infinity norm for ${rt}$ matrices + subroutine test_matrix_norms_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp) :: i,j,m,n + integer(ilp), parameter :: mtx_dim = 5 + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + ${rt}$, allocatable :: A(:,:) + type(linalg_state_type) :: err + character(64) :: msg + + allocate(A(mtx_dim,mtx_dim)) + + ! Initialize matrix with small values to avoid overflow + A = reshape([(0.01_${rk}$*(i-mtx_dim/2_ilp), i=1_ilp,mtx_dim*mtx_dim)], [mtx_dim,mtx_dim]) + + ! 1-norm (Maximum absolute column sum) + call check(error, abs(mnorm(A, '1', err) - maxval(sum(abs(A), dim=1),1)) < tol*mnorm(A, '1', err), & + 'Matrix 1-norm does not match expected value') + if (allocated(error)) return + + ! 2-norm (Frobenius norm) + call check(error, abs(mnorm(A, '2', err) - sqrt(sum(A**2)) < tol*mnorm(A, '2', err), & + 'Matrix Frobenius norm does not match expected value') + if (allocated(error)) return + + ! Inf-norm (Maximum absolute row sum) + call check(error, abs(mnorm(A, 'Inf', err) - maxval(sum(abs(A), dim=2),1)) < tol*mnorm(A, 'Inf', err), & + 'Matrix Infinity norm does not match expected value') + if (allocated(error)) return + + end subroutine test_matrix_norms_${ri}$ + + #:for rank in range(3, MAXRANK) + !> Test N-D norms + subroutine test_mnorm_${ri}$_${rank}$d(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp) :: j,dim(2) + 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)}$ + + 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 norm as collapsed around dimension + do dim = 1, ndim + write(msg,"('Not all dim=',i0,' Euclidean norms match ${rt}$ `norm2` intrinsic')") dim + call check(error,all(abs(mnorm(b,2,dim)-norm2(b,dim)) 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_mnorm From 8323e14343cf201347a321d57f235e40f66f1ab1 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 5 Nov 2024 15:10:02 +0100 Subject: [PATCH 10/23] first test program --- include/common.fypp | 18 +++- test/linalg/CMakeLists.txt | 2 + test/linalg/test_linalg_mnorm.fypp | 138 +++++++++++++++++++++++++++++ 3 files changed, 157 insertions(+), 1 deletion(-) create mode 100644 test/linalg/test_linalg_mnorm.fypp diff --git a/include/common.fypp b/include/common.fypp index b1169fcdb..256efa447 100644 --- a/include/common.fypp +++ b/include/common.fypp @@ -152,6 +152,22 @@ $:"s" if cmplx=="c" else "d" if cmplx=="z" else "x" if cmplx=="y" else "q" if cm #{if rank > 0}#(${"0" + ",0" * (rank - 1)}$)#{endif}# #:enddef +#! Generates an array rank suffix with a fixed integer size for all dimensions. +#! +#! Args: +#! rank (int): Rank of the variable +#! size (int): Size along each dimension +#! +#! Returns: +#! Array rank suffix string +#! E.g., +#! fixedranksuffix(3,4) +#! -> (4,4,4) +#! +#:def fixedranksuffix(rank,size) +#{if rank > 0}#(${str(size) + (","+str(size)) * (rank - 1)}$)#{endif}# +#:enddef + #! Joins stripped lines with given character string #! #! Args: @@ -222,7 +238,7 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #! Array rank suffix string enclosed in braces #! #! E.g., -#! select_subarray(5 , [(4, 'i'), (5, 'j')])}$ +#! select_subarray(5 , [(4, 'i'), (5, 'j')]) #! -> (:, :, :, i, j) #! #:def select_subarray(rank, selectors) diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index f835872df..a723b6ae1 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -8,6 +8,7 @@ set( "test_linalg_inverse.fypp" "test_linalg_lstsq.fypp" "test_linalg_norm.fypp" + "test_linalg_mnorm.fypp" "test_linalg_determinant.fypp" "test_linalg_qr.fypp" "test_linalg_svd.fypp" @@ -22,6 +23,7 @@ ADDTEST(linalg_eigenvalues) ADDTEST(linalg_matrix_property_checks) ADDTEST(linalg_inverse) ADDTEST(linalg_norm) +ADDTEST(linalg_mnorm) ADDTEST(linalg_solve) ADDTEST(linalg_lstsq) ADDTEST(linalg_qr) diff --git a/test/linalg/test_linalg_mnorm.fypp b/test/linalg/test_linalg_mnorm.fypp new file mode 100644 index 000000000..371f2eba3 --- /dev/null +++ b/test/linalg/test_linalg_mnorm.fypp @@ -0,0 +1,138 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +module test_linalg_mnorm + use testdrive, only: error_type, check, new_unittest, unittest_type + use stdlib_linalg_constants + use stdlib_linalg, only: mnorm, linalg_state_type + use stdlib_linalg_state, only: linalg_state_type + + implicit none (type,external) + + contains + + !> Matrix norm tests + subroutine test_matrix_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("test_matrix_norms_${ri}$",test_matrix_norms_${ri}$)] + #:for rank in range(3, MAXRANK) + tests = [tests,new_unittest("test_mnorm_${ri}$_${rank}$d",test_mnorm_${ri}$_${rank}$d)] + #:endfor + #:endfor + + end subroutine test_matrix_norms + + #:for rk,rt,ri in RC_KINDS_TYPES + !> Test 1-norm, 2-norm (Euclidean), and infinity norm for ${rt}$ matrices + subroutine test_matrix_norms_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp) :: i,j,m,n + integer(ilp), parameter :: mtx_dim = 5 + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + ${rt}$, allocatable :: A(:,:) + type(linalg_state_type) :: err + character(64) :: msg + + allocate(A(mtx_dim,mtx_dim)) + + ! Initialize matrix with small values to avoid overflow + A = reshape([(0.01_${rk}$*(i-mtx_dim/2_ilp), i=1_ilp,mtx_dim*mtx_dim)], [mtx_dim,mtx_dim]) + + ! 1-norm (Maximum absolute column sum) + call check(error, abs(mnorm(A, '1', err) - maxval(sum(abs(A), dim=1),1)) < tol*mnorm(A, '1', err), & + 'Matrix 1-norm does not match expected value') + if (allocated(error)) return + + ! 2-norm (Frobenius norm) + call check(error, abs(mnorm(A, '2', err) - sqrt(sum(A**2))) < tol*mnorm(A, '2', err), & + 'Matrix Frobenius norm does not match expected value') + if (allocated(error)) return + + ! Inf-norm (Maximum absolute row sum) + call check(error, abs(mnorm(A, 'Inf', err) - maxval(sum(abs(A), dim=2),1)) < tol*mnorm(A, 'Inf', err), & + 'Matrix Infinity norm does not match expected value') + if (allocated(error)) return + + end subroutine test_matrix_norms_${ri}$ + + #:for rank in range(3, MAXRANK) + !> Test N-D norms + subroutine test_mnorm_${ri}$_${rank}$d(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp) :: j,dim1,dim2,dim(2),order + integer(ilp), parameter :: orders(*) = [1_ilp,2_ilp,huge(0_ilp)] + integer(ilp), parameter :: ndim = ${rank}$ + integer(ilp), parameter :: n = 2_ilp**ndim + integer(ilp), parameter :: dims(*) = [(dim1, dim1=1,ndim)] + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + real(${rk}$), allocatable :: bnrm${ranksuffix(rank-2)}$ + ${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 norm as collapsed around dimensions + do j = 1, size(orders) + order = orders(j) + do dim1 = 1, ndim + do dim2 = dim1+1, ndim + + dim = [dim1,dim2] + + ! Get norms collapsed on these dims + bnrm = mnorm(b,order,dim) + + ! Assert size + write(msg,"('dim=[',i0,',',i0,'] order=',i0,' ${rk}$ norm returned wrong shape')") dim, order + call check(error,all(shape(bnrm)==pack(shape(b),dims/=dim1 .and. dims/=dim2) ), trim(msg)) + if (allocated(error)) return + + end do + end do + end do + + end subroutine test_mnorm_${ri}$_${rank}$d + #:endfor + + + #:endfor + +end module test_linalg_mnorm + +program test_mnorm + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_linalg_mnorm, only : test_matrix_norms + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("matrix_norms", test_matrix_norms) & + ] + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_mnorm From 3f30ed8404114d6503e8de3192b61fe3e4c1a014 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 6 Nov 2024 03:12:36 -0600 Subject: [PATCH 11/23] fix return size --- src/stdlib_linalg_norms.fypp | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index ac417c63d..dc643b0b0 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -488,6 +488,7 @@ ${loop_variables_end(rank-1," "*12)}$ call parse_norm_type(order,norm_request,err_) call lange_task_request(norm_request,lange_task,err_) if (err_%error()) then + allocate(nrm${emptyranksuffix(rank-2)}$) call linalg_error_handling(err_,err) return endif @@ -496,7 +497,7 @@ ${loop_variables_end(rank-1," "*12)}$ s = shape(a,kind=ilp) ! Check if input column data is contiguous - contiguous_data = dims(1)==1 + contiguous_data = all(dims==[1,2]) ! Matrix norm size m = s(dims(1)) @@ -505,12 +506,9 @@ ${loop_variables_end(rank-1," "*12)}$ ! Get packed data with norm dimensions as 1:2 if (contiguous_data) then - ! Collapse everything before the 1st dimension as apack's dim #1 - ! Set size==1 for all unused trailing dimensions - spack = [product(s(1:dims(2)-1)),s(dims(2):),(1_ilp,j=1,dims(2)-2)] - ! Reshape without moving data - apack${shape_from_array_data('spack',rank)}$ => a + spack = s + apack => a else @@ -518,7 +516,7 @@ ${loop_variables_end(rank-1," "*12)}$ perm = [dims,pack(dim_range, dim_range/=dims(1) .and. dim_range/=dims(2))] iperm(perm) = dim_range spack = s(perm) - apack = reshape(a, shape=spack, order=iperm) + allocate(apack,source=reshape(a, shape=spack, order=iperm)) endif From 5d4964566683335b0bac9b01d3371c037fd2163f Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 6 Nov 2024 03:13:59 -0600 Subject: [PATCH 12/23] Update test_linalg_mnorm.fypp --- test/linalg/test_linalg_mnorm.fypp | 1 + 1 file changed, 1 insertion(+) diff --git a/test/linalg/test_linalg_mnorm.fypp b/test/linalg/test_linalg_mnorm.fypp index 5d5f3c3f9..43a7e7e94 100644 --- a/test/linalg/test_linalg_mnorm.fypp +++ b/test/linalg/test_linalg_mnorm.fypp @@ -22,6 +22,7 @@ module test_linalg_mnorm #:for rank in range(3, MAXRANK) tests = [tests,new_unittest("test_mnorm_${ri}$_${rank}$d",test_mnorm_${ri}$_${rank}$d)] #:endfor + #:endfor end subroutine test_matrix_norms From ba596940a445d34d684eb89e82928bfa3923d3de Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 6 Nov 2024 03:40:09 -0600 Subject: [PATCH 13/23] add generic-rank tests --- include/common.fypp | 28 +++++++++++++++++++++ test/linalg/test_linalg_mnorm.fypp | 39 ++++++++++++++++++++++++------ 2 files changed, 60 insertions(+), 7 deletions(-) diff --git a/include/common.fypp b/include/common.fypp index 256efa447..87b861bf2 100644 --- a/include/common.fypp +++ b/include/common.fypp @@ -338,6 +338,34 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #:endcall #:enddef +#! +#! Generates a list of loop variables from an array +#! +#! Args: +#! varname(str): Name of the array 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_array_variables('j', 5) +#! -> "j(1), j(2), j(3), j(4), j(5) +#! +#! loop_array_variables('j', 5, 2) +#! -> "j(3), j(4), j(5), j(6), j(7) +#! +#:def loop_array_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: diff --git a/test/linalg/test_linalg_mnorm.fypp b/test/linalg/test_linalg_mnorm.fypp index 43a7e7e94..3c4031e34 100644 --- a/test/linalg/test_linalg_mnorm.fypp +++ b/test/linalg/test_linalg_mnorm.fypp @@ -65,14 +65,15 @@ module test_linalg_mnorm subroutine test_mnorm_${ri}$_${rank}$d(error) type(error_type), allocatable, intent(out) :: error - integer(ilp) :: j,dim1,dim2,dim(2),order + integer(ilp) :: i,j,k,l,dim1,dim2,dim(2),dim_sizes(2),order,ptr(${rank}$) integer(ilp), parameter :: orders(*) = [1_ilp,2_ilp,huge(0_ilp)] integer(ilp), parameter :: ndim = ${rank}$ - integer(ilp), parameter :: n = 2_ilp**ndim + integer(ilp), parameter :: n = 2_ilp**ndim integer(ilp), parameter :: dims(*) = [(dim1, dim1=1,ndim)] real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + real(${rk}$) :: one_nrm real(${rk}$), allocatable :: bnrm${ranksuffix(rank-2)}$ - ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$ + ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$, one_mat(:,:) character(64) :: msg @@ -84,12 +85,13 @@ module test_linalg_mnorm b = reshape(a, shape(b)) ! Test norm as collapsed around dimensions - do j = 1, size(orders) - order = orders(j) + do k = 1, size(orders) + order = orders(k) do dim1 = 1, ndim do dim2 = dim1+1, ndim - dim = [dim1,dim2] + dim = [dim1,dim2] + dim_sizes = [size(b,dim1,kind=ilp),size(b,dim2,kind=ilp)] ! Get norms collapsed on these dims bnrm = mnorm(b,order,dim) @@ -97,7 +99,30 @@ module test_linalg_mnorm ! Assert size write(msg,"('dim=[',i0,',',i0,'] order=',i0,' ${rk}$ norm returned wrong shape')") dim, order call check(error,all(shape(bnrm)==pack(shape(b),dims/=dim1 .and. dims/=dim2) ), trim(msg)) - if (allocated(error)) return + if (allocated(error)) return + + ! Assert some matrix results: check that those on same index i.e. (l,l,l,:,l,l,:) etc. + ! are equal to the corresponding 2d-array result + do l = 1, minval(shape(b)) + + ptr = l + + allocate(one_mat(dim_sizes(1),dim_sizes(2))) + do j = 1, dim_sizes(2) + ptr(dim(2)) = j + do i = 1, dim_sizes(1) + ptr(dim(1)) = i + one_mat(i,j) = b(${loop_array_variables('ptr',rank)}$) + end do + end do + one_nrm = mnorm(one_mat,order) + + write(msg,"('dim=[',i0,',',i0,'] order=',i0,' ${rk}$ ',i0,'-th norm is wrong')") dim, order, l + call check(error, abs(one_nrm-bnrm(${fixedranksuffix(rank-2,'l')}$)) Date: Wed, 6 Nov 2024 03:41:36 -0600 Subject: [PATCH 14/23] better example --- example/linalg/example_mnorm.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/example/linalg/example_mnorm.f90 b/example/linalg/example_mnorm.f90 index 24b2036e2..4af32b857 100644 --- a/example/linalg/example_mnorm.f90 +++ b/example/linalg/example_mnorm.f90 @@ -20,5 +20,7 @@ program example_mnorm ! Compute infinity norm of each 3x3 matrix in b nb = mnorm(b, 'inf', dim=[1,2]) + + ! 18.0000000 36.0000000 54.0000000 72.0000000 print *, "Infinity norms of matrices in b:", nb end program example_mnorm From 753e31b3385dda38d4bb004030d81d06f14ae471 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 6 Nov 2024 04:09:56 -0600 Subject: [PATCH 15/23] sync task with NumPy: no order-2, default=Frobenius --- doc/specs/stdlib_linalg.md | 4 +- src/stdlib_linalg.fypp | 20 ++++---- src/stdlib_linalg_norms.fypp | 79 ++++++++++++++++++++++-------- test/linalg/test_linalg_mnorm.fypp | 13 ++--- 4 files changed, 78 insertions(+), 38 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 6df881ad8..e3a79de02 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1585,12 +1585,12 @@ matrix norms are computed over specified dimensions. `a`: Shall be a rank-n `real` or `complex` array containing the data, where n ³ 2. 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. +`order` (optional): 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 | |------------------|---------------------------|---------------------------------------------------------| | `1` | `'1'` | 1-norm (maximum column sum) \( \max_j \sum_i{ \left|a_{i,j}\right| } \) | -| `2` | `'2','Euclidean'` | 2-norm/Euclidean/Frobenius norm \( \sqrt{\sum_{i,j}{ \left|a_{i,j}\right|^2 }} \) | +| (not prov.) | `'Euclidean','Frobenius','Fro'` | Frobenius norm \( \sqrt{\sum_{i,j}{ \left|a_{i,j}\right|^2 }} \) | | `huge(0)` | `'inf', 'Inf', 'INF'` | Infinity norm (maximum row sum) \( \max_i \sum_j{ \left|a_{i,j}\right| } \) | `dim` (optional): For arrays of rank > 2, shall be an integer array of size 2 specifying the dimensions over which to compute the matrix norm. Default value is `[1,2]`. It is an `intent(in)` argument. diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index a413cc5f5..ea120f005 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1333,7 +1333,8 @@ module stdlib_linalg !! that can have rank 2 or higher. For rank-2 arrays, the matrix norm is returned. !! If rank>2 and the optional input dimensions `dim` are specified, !! a rank `n-2` array is returned with dimensions `dim(1),dim(2)` collapsed, containing all - !! matrix norms evaluated over the specified dimensions only. + !! matrix norms evaluated over the specified dimensions only. `dim==[1,2]` are assumed as default + !! dimensions if not specified. !! !!### Description !! @@ -1345,12 +1346,10 @@ module stdlib_linalg !! This can be provided as either an `integer` value or a `character` string. !! Allowed metrics are: !! - 1-norm: `order` = 1 or '1' - !! - 2-norm/Euclidean: `order` = 2 or '2' or 'Euclidean' - !! - p-norm: `order` >= 3 + !! - Euclidean/Frobenius: `order` = 'Euclidean','Frobenius', or argument not specified !! - Infinity norm: `order` = huge(0) or 'Inf' - !! - Minus-infinity norm: `order` = '-Inf' !! - !! If an invalid norm type is provided, the routine defaults to 1-norm and returns an error state. + !! If an invalid norm type is provided, the routine returns an error state. !! !!### Example !! @@ -1359,14 +1358,15 @@ module stdlib_linalg !! real(sp) :: b(3,3,4), nb(4) ! Array of 4 3x3 matrices !! a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3]) !! - !! ! 2-norm/Euclidean norm of single matrix + !! ! Euclidean/Frobenius norm of single matrix + !! na = mnorm(a) !! na = mnorm(a, 'Euclidean') !! !! ! 1-norm of each 3x3 matrix in b !! nb = mnorm(b, 1, dim=[1,2]) !! - !! ! p-norm with p=3 - !! na = mnorm(a, 3) + !! ! Infinity-norm + !! na = mnorm(b, 'inf', dim=[3,2]) !!``` !! #:for rk,rt,ri in RC_KINDS_TYPES @@ -1379,7 +1379,7 @@ module stdlib_linalg !> Norm of the matrix. real(${rk}$) :: nrm !> Order of the matrix norm being computed. - ${it}$, intent(in) :: order + ${it}$, #{if 'integer' in it}#optional, #{endif}#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 function matrix_norm_${ii}$_${ri}$ @@ -1392,7 +1392,7 @@ module stdlib_linalg !> Norm of the matrix. real(${rk}$), allocatable :: nrm${ranksuffix(rank-2)}$ !> Order of the matrix norm being computed. - ${it}$, intent(in) :: order + ${it}$, #{if 'integer' in it}#optional, #{endif}#intent(in) :: order !> [optional] dimensions of the sub-matrices the norms should be evaluated at (default = [1,2]) integer(ilp), optional, intent(in) :: dim(2) !> [optional] state return flag. On error if not requested, the code will stop diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index dc643b0b0..7be06ffff 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -29,13 +29,18 @@ submodule(stdlib_linalg) stdlib_linalg_norms character, parameter :: LANGE_NORM_MAT = 'M' ! maxval(sum(abs(a))) ! over whole matrix: unused character, parameter :: LANGE_NORM_ONE = '1' ! maxval(sum(abs(a),1)) ! over columns character, parameter :: LANGE_NORM_INF = 'I' ! maxval(sum(abs(a),2)) ! over rows - character, parameter :: LANGE_NORM_TWO = 'E' ! "Euclidean" or "Frobenius" + character, parameter :: LANGE_NORM_FRO = 'E' ! sqrt(sum(a**2)) ! "Euclidean" or "Frobenius" interface parse_norm_type module procedure parse_norm_type_integer module procedure parse_norm_type_character end interface parse_norm_type + interface lange_task_request + module procedure lange_task_request_integer + module procedure lange_task_request_character + end interface lange_task_request + interface stride_1d #:for rk,rt,ri in ALL_KINDS_TYPES @@ -107,25 +112,61 @@ submodule(stdlib_linalg) stdlib_linalg_norms end subroutine parse_norm_type_character !> From a user norm request, generate a *LANGE task command - pure subroutine lange_task_request(norm_type,lange_task,err) + pure subroutine lange_task_request_integer(order,lange_task,err) !> Parsed matrix norm type - integer(ilp), intent(in) :: norm_type + integer(ilp), optional, intent(in) :: order !> LANGE task character, intent(out) :: lange_task !> Error flag type(linalg_state_type), intent(inout) :: err - select case (norm_type) - case (NORM_INF) - lange_task = LANGE_NORM_INF - case (NORM_ONE) + if (present(order)) then + + select case (order) + case (NORM_INF) + lange_task = LANGE_NORM_INF + case (NORM_ONE) + lange_task = LANGE_NORM_ONE + case default + err = linalg_state_type(this,LINALG_VALUE_ERROR,'Integer order ',order,' is not a valid matrix norm input.') + end select + + else + + ! No user input: Frobenius norm + lange_task = LANGE_NORM_FRO + + endif + end subroutine lange_task_request_integer + + pure subroutine lange_task_request_character(order,lange_task,err) + !> User input value + character(len=*), intent(in) :: order + !> Return value: norm type + character, intent(out) :: lange_task + !> State return flag + type(linalg_state_type), intent(out) :: err + + integer(ilp) :: int_order,read_err + + select case (order) + case ('inf','Inf','INF') + lange_task = LANGE_NORM_INF + case ('Euclidean','euclidean','EUCLIDEAN','Frobenius','frobenius','FROBENIUS','Fro','fro','frob') + lange_task = LANGE_NORM_FRO + case default + + ! Check if this input can be read as an integer + read(order,*,iostat=read_err) int_order + if (read_err/=0 .or. int_order/=1) then + ! Cannot read as an integer + err = linalg_state_type(this,LINALG_ERROR,'Matrix norm input',order,' is not recognized.') + endif lange_task = LANGE_NORM_ONE - case (NORM_TWO) - lange_task = LANGE_NORM_TWO - case default - err = linalg_state_type(this,LINALG_VALUE_ERROR,'Order ',norm_type,' is not a valid matrix norm input.') - end select - end subroutine lange_task_request + + end select + + end subroutine lange_task_request_character #:for rk,rt,ri in ALL_KINDS_TYPES @@ -399,12 +440,12 @@ ${loop_variables_end(rank-1," "*12)}$ !> Norm of the matrix. real(${rk}$) :: nrm !> Order of the matrix norm being computed. - ${it}$, intent(in) :: order + ${it}$, #{if 'integer' in it}#optional, #{endif}#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) :: m,n,norm_request + integer(ilp) :: m,n character :: lange_task real(${rk}$), target :: work1(1) real(${rk}$), pointer :: work(:) @@ -422,8 +463,7 @@ ${loop_variables_end(rank-1," "*12)}$ end if ! Check norm request: user + *LANGE support - call parse_norm_type(order,norm_request,err_) - call lange_task_request(norm_request,lange_task,err_) + call lange_task_request(order,lange_task,err_) if (err_%error()) then call linalg_error_handling(err_,err) return @@ -451,7 +491,7 @@ ${loop_variables_end(rank-1," "*12)}$ !> Norm of the matrix. real(${rk}$), allocatable :: nrm${ranksuffix(rank-2)}$ !> Order of the matrix norm being computed. - ${it}$, intent(in) :: order + ${it}$, #{if 'integer' in it}#optional, #{endif}#intent(in) :: order !> [optional] dimensions of the sub-matrices the norms should be evaluated at (default = [1,2]) integer(ilp), optional, intent(in) :: dim(2) !> [optional] state return flag. On error if not requested, the code will stop @@ -485,8 +525,7 @@ ${loop_variables_end(rank-1," "*12)}$ endif ! Check norm request: user + *LANGE support - call parse_norm_type(order,norm_request,err_) - call lange_task_request(norm_request,lange_task,err_) + call lange_task_request(order,lange_task,err_) if (err_%error()) then allocate(nrm${emptyranksuffix(rank-2)}$) call linalg_error_handling(err_,err) diff --git a/test/linalg/test_linalg_mnorm.fypp b/test/linalg/test_linalg_mnorm.fypp index 3c4031e34..cd59f63a7 100644 --- a/test/linalg/test_linalg_mnorm.fypp +++ b/test/linalg/test_linalg_mnorm.fypp @@ -49,7 +49,7 @@ module test_linalg_mnorm if (allocated(error)) return ! 2-norm (Frobenius norm) - call check(error, abs(mnorm(A, '2', err) - sqrt(sum(A**2))) < tol*mnorm(A, '2', err), & + call check(error, abs(mnorm(A, err=err) - sqrt(sum(A**2))) < tol*mnorm(A, err=err), & 'Matrix Frobenius norm does not match expected value') if (allocated(error)) return @@ -65,8 +65,8 @@ module test_linalg_mnorm subroutine test_mnorm_${ri}$_${rank}$d(error) type(error_type), allocatable, intent(out) :: error - integer(ilp) :: i,j,k,l,dim1,dim2,dim(2),dim_sizes(2),order,ptr(${rank}$) - integer(ilp), parameter :: orders(*) = [1_ilp,2_ilp,huge(0_ilp)] + integer(ilp) :: i,j,k,l,dim1,dim2,dim(2),dim_sizes(2),ptr(${rank}$) + character(3), parameter :: orders(*) = ['1 ','fro','inf'] integer(ilp), parameter :: ndim = ${rank}$ integer(ilp), parameter :: n = 2_ilp**ndim integer(ilp), parameter :: dims(*) = [(dim1, dim1=1,ndim)] @@ -74,6 +74,7 @@ module test_linalg_mnorm real(${rk}$) :: one_nrm real(${rk}$), allocatable :: bnrm${ranksuffix(rank-2)}$ ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$, one_mat(:,:) + character(:), allocatable :: order character(64) :: msg @@ -86,7 +87,7 @@ module test_linalg_mnorm ! Test norm as collapsed around dimensions do k = 1, size(orders) - order = orders(k) + order = trim(orders(k)) do dim1 = 1, ndim do dim2 = dim1+1, ndim @@ -97,7 +98,7 @@ module test_linalg_mnorm bnrm = mnorm(b,order,dim) ! Assert size - write(msg,"('dim=[',i0,',',i0,'] order=',i0,' ${rk}$ norm returned wrong shape')") dim, order + write(msg,"('dim=[',i0,',',i0,'] order=',a,' ${rk}$ norm returned wrong shape')") dim, order call check(error,all(shape(bnrm)==pack(shape(b),dims/=dim1 .and. dims/=dim2) ), trim(msg)) if (allocated(error)) return @@ -117,7 +118,7 @@ module test_linalg_mnorm end do one_nrm = mnorm(one_mat,order) - write(msg,"('dim=[',i0,',',i0,'] order=',i0,' ${rk}$ ',i0,'-th norm is wrong')") dim, order, l + write(msg,"('dim=[',i0,',',i0,'] order=',a,' ${rk}$ ',i0,'-th norm is wrong')") dim, order, l call check(error, abs(one_nrm-bnrm(${fixedranksuffix(rank-2,'l')}$)) Date: Thu, 7 Nov 2024 08:26:39 +0100 Subject: [PATCH 16/23] replace `LANGE_*` with `MAT_*`, add `max(svdvals(a))` option --- src/stdlib_linalg_norms.fypp | 73 ++++++++++++++++++++---------------- 1 file changed, 41 insertions(+), 32 deletions(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 7be06ffff..615969480 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -26,20 +26,22 @@ submodule(stdlib_linalg) stdlib_linalg_norms integer(ilp), parameter :: NORM_MINUSINF = -huge(0_ilp) !> Wrappers to LAPACK *LANGE matrix norm flags - character, parameter :: LANGE_NORM_MAT = 'M' ! maxval(sum(abs(a))) ! over whole matrix: unused - character, parameter :: LANGE_NORM_ONE = '1' ! maxval(sum(abs(a),1)) ! over columns - character, parameter :: LANGE_NORM_INF = 'I' ! maxval(sum(abs(a),2)) ! over rows - character, parameter :: LANGE_NORM_FRO = 'E' ! sqrt(sum(a**2)) ! "Euclidean" or "Frobenius" + character, parameter :: MAT_NORM_MAT = 'M' ! maxval(sum(abs(a))) ! over whole matrix: unused + character, parameter :: MAT_NORM_ONE = '1' ! maxval(sum(abs(a),1)) ! over columns + character, parameter :: MAT_NORM_INF = 'I' ! maxval(sum(abs(a),2)) ! over rows + character, parameter :: MAT_NORM_FRO = 'E' ! sqrt(sum(a**2)) ! "Euclidean" or "Frobenius" + !> Other wrappers to matrix norms + character, parameter :: MAT_NORM_SVD = '2' ! maxval(svdvals(a)) ! Maximum singular value interface parse_norm_type module procedure parse_norm_type_integer module procedure parse_norm_type_character end interface parse_norm_type - interface lange_task_request - module procedure lange_task_request_integer - module procedure lange_task_request_character - end interface lange_task_request + interface mat_task_request + module procedure mat_task_request_integer + module procedure mat_task_request_character + end interface mat_task_request interface stride_1d @@ -112,11 +114,11 @@ submodule(stdlib_linalg) stdlib_linalg_norms end subroutine parse_norm_type_character !> From a user norm request, generate a *LANGE task command - pure subroutine lange_task_request_integer(order,lange_task,err) + pure subroutine mat_task_request_integer(order,mat_task,err) !> Parsed matrix norm type integer(ilp), optional, intent(in) :: order !> LANGE task - character, intent(out) :: lange_task + character, intent(out) :: mat_task !> Error flag type(linalg_state_type), intent(inout) :: err @@ -124,9 +126,11 @@ submodule(stdlib_linalg) stdlib_linalg_norms select case (order) case (NORM_INF) - lange_task = LANGE_NORM_INF + mat_task = MAT_NORM_INF + case (NORM_TWO) + mat_task = MAT_NORM_SVD case (NORM_ONE) - lange_task = LANGE_NORM_ONE + mat_task = MAT_NORM_ONE case default err = linalg_state_type(this,LINALG_VALUE_ERROR,'Integer order ',order,' is not a valid matrix norm input.') end select @@ -134,16 +138,16 @@ submodule(stdlib_linalg) stdlib_linalg_norms else ! No user input: Frobenius norm - lange_task = LANGE_NORM_FRO + mat_task = MAT_NORM_FRO endif - end subroutine lange_task_request_integer + end subroutine mat_task_request_integer - pure subroutine lange_task_request_character(order,lange_task,err) + pure subroutine mat_task_request_character(order,mat_task,err) !> User input value character(len=*), intent(in) :: order !> Return value: norm type - character, intent(out) :: lange_task + character, intent(out) :: mat_task !> State return flag type(linalg_state_type), intent(out) :: err @@ -151,22 +155,27 @@ submodule(stdlib_linalg) stdlib_linalg_norms select case (order) case ('inf','Inf','INF') - lange_task = LANGE_NORM_INF + mat_task = MAT_NORM_INF case ('Euclidean','euclidean','EUCLIDEAN','Frobenius','frobenius','FROBENIUS','Fro','fro','frob') - lange_task = LANGE_NORM_FRO + mat_task = MAT_NORM_FRO case default ! Check if this input can be read as an integer read(order,*,iostat=read_err) int_order - if (read_err/=0 .or. int_order/=1) then + if (read_err/=0 .or. all(int_order/=[1,2]) then ! Cannot read as an integer err = linalg_state_type(this,LINALG_ERROR,'Matrix norm input',order,' is not recognized.') endif - lange_task = LANGE_NORM_ONE - + + select case (int_order) + case (1); mat_task = MAT_NORM_ONE + case (2); mat_task = MAT_NORM_SVD + case default; mat_task = MAT_NORM_ONE + end select + end select - end subroutine lange_task_request_character + end subroutine mat_task_request_character #:for rk,rt,ri in ALL_KINDS_TYPES @@ -446,7 +455,7 @@ ${loop_variables_end(rank-1," "*12)}$ type(linalg_state_type) :: err_ integer(ilp) :: m,n - character :: lange_task + character :: mat_task real(${rk}$), target :: work1(1) real(${rk}$), pointer :: work(:) @@ -463,22 +472,22 @@ ${loop_variables_end(rank-1," "*12)}$ end if ! Check norm request: user + *LANGE support - call lange_task_request(order,lange_task,err_) + call mat_task_request(order,mat_task,err_) if (err_%error()) then call linalg_error_handling(err_,err) return endif - if (lange_task==LANGE_NORM_INF) then + if (mat_task==MAT_NORM_INF) then allocate(work(m)) else work => work1 endif ! LAPACK interface - nrm = lange(lange_task,m,n,a,m,work) + nrm = lange(mat_task,m,n,a,m,work) - if (lange_task==LANGE_NORM_INF) deallocate(work) + if (mat_task==MAT_NORM_INF) deallocate(work) end function matrix_norm_${ii}$_${ri}$ @@ -503,7 +512,7 @@ ${loop_variables_end(rank-1," "*12)}$ integer(ilp), dimension(${rank}$), parameter :: dim_range = [(m,m=1_ilp,${rank}$_ilp)] integer(ilp) :: ${loop_variables('j',rank-2,2)}$ logical :: contiguous_data - character :: lange_task + character :: mat_task real(${rk}$), target :: work1(1) real(${rk}$), pointer :: work(:) ${rt}$, pointer :: apack${ranksuffix(rank)}$ @@ -525,7 +534,7 @@ ${loop_variables_end(rank-1," "*12)}$ endif ! Check norm request: user + *LANGE support - call lange_task_request(order,lange_task,err_) + call mat_task_request(order,mat_task,err_) if (err_%error()) then allocate(nrm${emptyranksuffix(rank-2)}$) call linalg_error_handling(err_,err) @@ -559,7 +568,7 @@ ${loop_variables_end(rank-1," "*12)}$ endif - if (lange_task==LANGE_NORM_INF) then + if (mat_task==MAT_NORM_INF) then allocate(work(m)) else work => work1 @@ -573,10 +582,10 @@ ${loop_variables_end(rank-1," "*12)}$ ! LAPACK interface ${loop_variables_start('j', 'apack', rank-2, 2)}$ nrm(${loop_variables('j',rank-2,2)}$) = & - lange(lange_task,m,n,apack(:,:,${loop_variables('j',rank-2,2)}$),lda,work) + lange(mat_task,m,n,apack(:,:,${loop_variables('j',rank-2,2)}$),lda,work) ${loop_variables_end(rank-2)}$ - if (lange_task==LANGE_NORM_INF) deallocate(work) + if (mat_task==MAT_NORM_INF) deallocate(work) if (.not.contiguous_data) deallocate(apack) end function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$ From 106bf232bb10459d521af0522818ddc1686cb8bd Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 7 Nov 2024 08:43:55 +0100 Subject: [PATCH 17/23] implement 2-norm `maxval(svdvals(a))` --- doc/specs/stdlib_linalg.md | 1 + src/stdlib_linalg.fypp | 1 + src/stdlib_linalg_norms.fypp | 44 ++++++++++++++++++++++++++++-------- 3 files changed, 37 insertions(+), 9 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index e3a79de02..eaf9d6110 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1590,6 +1590,7 @@ matrix norms are computed over specified dimensions. | Integer input | Character Input | Norm type | |------------------|---------------------------|---------------------------------------------------------| | `1` | `'1'` | 1-norm (maximum column sum) \( \max_j \sum_i{ \left|a_{i,j}\right| } \) | +| `2` | `'2'` | 2-norm (largest singular value) | | (not prov.) | `'Euclidean','Frobenius','Fro'` | Frobenius norm \( \sqrt{\sum_{i,j}{ \left|a_{i,j}\right|^2 }} \) | | `huge(0)` | `'inf', 'Inf', 'INF'` | Infinity norm (maximum row sum) \( \max_i \sum_j{ \left|a_{i,j}\right| } \) | diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index ea120f005..a85a9022e 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1346,6 +1346,7 @@ module stdlib_linalg !! This can be provided as either an `integer` value or a `character` string. !! Allowed metrics are: !! - 1-norm: `order` = 1 or '1' + !! - 2-norm: `order` = 2 or '2' !! - Euclidean/Frobenius: `order` = 'Euclidean','Frobenius', or argument not specified !! - Infinity norm: `order` = huge(0) or 'Inf' !! diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 615969480..85bc8a248 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -162,7 +162,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms ! Check if this input can be read as an integer read(order,*,iostat=read_err) int_order - if (read_err/=0 .or. all(int_order/=[1,2]) then + if (read_err/=0 .or. all(int_order/=[1,2])) then ! Cannot read as an integer err = linalg_state_type(this,LINALG_ERROR,'Matrix norm input',order,' is not recognized.') endif @@ -482,10 +482,15 @@ ${loop_variables_end(rank-1," "*12)}$ allocate(work(m)) else work => work1 - endif - - ! LAPACK interface - nrm = lange(mat_task,m,n,a,m,work) + end if + + if (mat_task==MAT_NORM_SVD) then + nrm = maxval(svdvals(a,err_),1) + call linalg_error_handling(err_,err) + else + ! LAPACK interface + nrm = lange(mat_task,m,n,a,m,work) + end if if (mat_task==MAT_NORM_INF) deallocate(work) @@ -507,7 +512,7 @@ ${loop_variables_end(rank-1," "*12)}$ type(linalg_state_type), intent(out), optional :: err type(linalg_state_type) :: err_ - integer(ilp) :: j,m,n,lda,dims(2),norm_request + integer(ilp) :: j,m,n,lda,dims(2),norm_request,svd_errors integer(ilp), dimension(${rank}$) :: s,spack,perm,iperm integer(ilp), dimension(${rank}$), parameter :: dim_range = [(m,m=1_ilp,${rank}$_ilp)] integer(ilp) :: ${loop_variables('j',rank-2,2)}$ @@ -525,6 +530,7 @@ ${loop_variables_end(rank-1," "*12)}$ endif nullify(apack) + svd_errors = 0 if (dims(1)==dims(2) .or. .not.all(dims>0 .and. dims<=${rank}$)) then err_ = linalg_state_type(this,LINALG_VALUE_ERROR,'Rank-',${rank}$,' matrix norm has invalid dim=',dims) @@ -551,6 +557,13 @@ ${loop_variables_end(rank-1," "*12)}$ m = s(dims(1)) n = s(dims(2)) + if (m<=0 .or. n<=0) then + err_ = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix shape: a=',[m,n]) + allocate(nrm${emptyranksuffix(rank-2)}$) + call linalg_error_handling(err_,err) + return + end if + ! Get packed data with norm dimensions as 1:2 if (contiguous_data) then @@ -570,6 +583,8 @@ ${loop_variables_end(rank-1," "*12)}$ if (mat_task==MAT_NORM_INF) then allocate(work(m)) + elseif (mat_task==MAT_NORM_SVD) then + allocate(work(min(m,n))) else work => work1 endif @@ -581,12 +596,23 @@ ${loop_variables_end(rank-1," "*12)}$ ! LAPACK interface ${loop_variables_start('j', 'apack', rank-2, 2)}$ - nrm(${loop_variables('j',rank-2,2)}$) = & - lange(mat_task,m,n,apack(:,:,${loop_variables('j',rank-2,2)}$),lda,work) + if (mat_task==MAT_NORM_SVD) then + work(:) = svdvals(apack(:,:,${loop_variables('j',rank-2,2)}$),err_) + nrm(${loop_variables('j',rank-2,2)}$) = maxval(work,1) + if (err_%error()) svd_errors = svd_errors+1 + else + nrm(${loop_variables('j',rank-2,2)}$) = & + lange(mat_task,m,n,apack(:,:,${loop_variables('j',rank-2,2)}$),lda,work) + end if ${loop_variables_end(rank-2)}$ - if (mat_task==MAT_NORM_INF) deallocate(work) + if (any(mat_task==[MAT_NORM_INF,MAT_NORM_SVD])) deallocate(work) if (.not.contiguous_data) deallocate(apack) + + if (svd_errors>0) then + err_ = linalg_state_type(this,LINALG_VALUE_ERROR,svd_errors,'failed SVDs') + call linalg_error_handling(err_,err) + endif end function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$ From a2dd310b0a33360175784c72b3caf95e56918e28 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 7 Nov 2024 08:44:30 +0100 Subject: [PATCH 18/23] add 2-norm to the tests --- test/linalg/test_linalg_mnorm.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/linalg/test_linalg_mnorm.fypp b/test/linalg/test_linalg_mnorm.fypp index cd59f63a7..03a5a627e 100644 --- a/test/linalg/test_linalg_mnorm.fypp +++ b/test/linalg/test_linalg_mnorm.fypp @@ -66,7 +66,7 @@ module test_linalg_mnorm type(error_type), allocatable, intent(out) :: error integer(ilp) :: i,j,k,l,dim1,dim2,dim(2),dim_sizes(2),ptr(${rank}$) - character(3), parameter :: orders(*) = ['1 ','fro','inf'] + character(3), parameter :: orders(*) = ['1 ','2 ','fro','inf'] integer(ilp), parameter :: ndim = ${rank}$ integer(ilp), parameter :: n = 2_ilp**ndim integer(ilp), parameter :: dims(*) = [(dim1, dim1=1,ndim)] From 89f2fa4d01ee1744d4be3728491a6f3dd4342e80 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 7 Nov 2024 08:47:12 +0100 Subject: [PATCH 19/23] test flipped dimensions e.g. [2,1] --- test/linalg/test_linalg_mnorm.fypp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/linalg/test_linalg_mnorm.fypp b/test/linalg/test_linalg_mnorm.fypp index 03a5a627e..890ca1885 100644 --- a/test/linalg/test_linalg_mnorm.fypp +++ b/test/linalg/test_linalg_mnorm.fypp @@ -89,7 +89,9 @@ module test_linalg_mnorm do k = 1, size(orders) order = trim(orders(k)) do dim1 = 1, ndim - do dim2 = dim1+1, ndim + do dim2 = 1, ndim + + if (dim1==dim2) cycle dim = [dim1,dim2] dim_sizes = [size(b,dim1,kind=ilp),size(b,dim2,kind=ilp)] From 65d06f7b6f656bd207b0db353999382adc11d613 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 7 Nov 2024 08:49:13 +0100 Subject: [PATCH 20/23] align specs table --- doc/specs/stdlib_linalg.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index eaf9d6110..f9284f587 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1587,12 +1587,12 @@ matrix norms are computed over specified dimensions. `order` (optional): 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 | -|------------------|---------------------------|---------------------------------------------------------| -| `1` | `'1'` | 1-norm (maximum column sum) \( \max_j \sum_i{ \left|a_{i,j}\right| } \) | -| `2` | `'2'` | 2-norm (largest singular value) | -| (not prov.) | `'Euclidean','Frobenius','Fro'` | Frobenius norm \( \sqrt{\sum_{i,j}{ \left|a_{i,j}\right|^2 }} \) | -| `huge(0)` | `'inf', 'Inf', 'INF'` | Infinity norm (maximum row sum) \( \max_i \sum_j{ \left|a_{i,j}\right| } \) | +| Integer input | Character Input | Norm type | +|------------------|---------------------------------|-----------------------------------------------------------------------------| +| `1` | `'1'` | 1-norm (maximum column sum) \( \max_j \sum_i{ \left|a_{i,j}\right| } \) | +| `2` | `'2'` | 2-norm (largest singular value) | +| (not prov.) | `'Euclidean','Frobenius','Fro'` | Frobenius norm \( \sqrt{\sum_{i,j}{ \left|a_{i,j}\right|^2 }} \) | +| `huge(0)` | `'inf', 'Inf', 'INF'` | Infinity norm (maximum row sum) \( \max_i \sum_j{ \left|a_{i,j}\right| } \) | `dim` (optional): For arrays of rank > 2, shall be an integer array of size 2 specifying the dimensions over which to compute the matrix norm. Default value is `[1,2]`. It is an `intent(in)` argument. From 706ea1c1b9719a34d83c1c2f7d3acf07083a210f Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 18 Nov 2024 09:20:53 +0100 Subject: [PATCH 21/23] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index f9284f587..b7d00e828 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1583,7 +1583,7 @@ matrix norms are computed over specified dimensions. ### Arguments -`a`: Shall be a rank-n `real` or `complex` array containing the data, where n ³ 2. It is an `intent(in)` argument. +`a`: Shall be a rank-n `real` or `complex` array containing the data, where n >= 2. It is an `intent(in)` argument. `order` (optional): Shall be an `integer` value or a `character` flag that specifies the norm type, as follows. It is an `intent(in)` argument. From f7ec7da09ccb5dce206dc62a2f56cbce7f21ea6b Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 18 Nov 2024 09:24:22 +0100 Subject: [PATCH 22/23] Update doc/specs/stdlib_linalg.md --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index b7d00e828..cb34c0cc0 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1579,7 +1579,7 @@ matrix norms are computed over specified dimensions. ### Syntax -`x = ` [[stdlib_linalg(module):mnorm(interface)]] `(a, order [, dim, err])` +`x = ` [[stdlib_linalg(module):mnorm(interface)]] `(a [, order, dim, err])` ### Arguments From f52d57a78a92d2398f1932053eba325bdf6d0578 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 18 Nov 2024 10:46:19 +0100 Subject: [PATCH 23/23] Update src/stdlib_linalg.fypp --- src/stdlib_linalg.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index a85a9022e..ad8ca3fb8 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1342,7 +1342,7 @@ module stdlib_linalg !! Supported data types include `real` and `complex`. !! Input arrays must have rank >= 2. !! - !! Norm type input is mandatory, and it is provided via the `order` argument. + !! Norm type input is optional, 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: `order` = 1 or '1'