From 1dbf88e3fb108051ac404a1905ff8c6e98787346 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 13 Sep 2024 10:17:40 +0200 Subject: [PATCH 01/32] add `norms` module --- src/CMakeLists.txt | 1 + src/stdlib_linalg_norms.fypp | 335 +++++++++++++++++++++++++++++++++++ 2 files changed, 336 insertions(+) create mode 100644 src/stdlib_linalg_norms.fypp 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_norms.fypp b/src/stdlib_linalg_norms.fypp new file mode 100644 index 000000000..78db93d4c --- /dev/null +++ b/src/stdlib_linalg_norms.fypp @@ -0,0 +1,335 @@ +#:include "common.fypp" +#:set INPUT_TYPE = ["character(len=*)","integer(ilp)"] +#:set INPUT_SHORT = ["char","int"] +#:set INPUT_OPTIONS = list(zip(INPUT_TYPE,INPUT_SHORT)) +! Vector norms +module stdlib_linalg_norms + use stdlib_linalg_constants + 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 + implicit none(type,external) + private + + public :: norm, get_norm + + 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) + + !> Vector norm: function interface + interface norm + #:for rk,rt,ri in ALL_KINDS_TYPES + #:for it,ii in INPUT_OPTIONS + !> Scalar norms: ${rt}$ + #:for rank in range(1, MAXRANK + 1) + module procedure stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$ + module procedure stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$ + #:endfor + !> Array norms: ${rt}$ + #:for rank in range(2, MAXRANK + 1) + module procedure stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ + module procedure 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 + #:for rk,rt,ri in ALL_KINDS_TYPES + #:for it,ii in INPUT_OPTIONS + !> Scalar norms: ${rt}$ + #:for rank in range(1, MAXRANK + 1) + module procedure norm_${rank}$D_${ii}$_${ri}$ + #:endfor + !> Array norms: ${rt}$ + #:for rank in range(2, MAXRANK + 1) + module procedure norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ + #:endfor + #:endfor + #:endfor + end interface get_norm + + interface parse_norm_type + module procedure parse_norm_type_integer + module procedure parse_norm_type_character + end interface parse_norm_type + + 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:huge(0_ilp)-1_ilp) + norm_type = order + case (huge(0_ilp):) + norm_type = NORM_INF + case (:-huge(0_ilp)) + 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 + #: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 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 + 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 + pure subroutine norm_${rank}$D_${ii}$_${ri}$(a, nrm, order, err) + !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ + ${rt}$, intent(in) :: 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 + + 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 + + select case(norm_request) + case(NORM_ONE) + nrm = sum( abs(a) ) + case(NORM_TWO) + #:if rt.startswith('complex') + nrm = sqrt( real( sum( a * conjg(a) ), ${rk}$) ) + #:else + nrm = sqrt( sum( a ** 2 ) ) + #:endif + case(NORM_INF) + nrm = maxval( abs(a) ) + case(-NORM_INF) + nrm = minval( abs(a) ) + case (NORM_POW_FIRST:NORM_POW_LAST) + rorder = 1.0_${rk}$ / norm_request + nrm = sum( abs(a) ** norm_request ) ** rorder + case default + err_ = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid norm type after checking') + call linalg_error_handling(err_,err) + end select + + 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 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. + 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 subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err) + !> Input matrix a[..] + ${rt}$, intent(in), target :: 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,norm_request + real(${rk}$) :: rorder + + sze = size(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 + + select case(norm_request) + case(NORM_ONE) + nrm = sum( abs(a) , dim = dim ) + case(NORM_TWO) + #:if rt.startswith('complex') + nrm = sqrt( real( sum( a * conjg(a) , dim = dim ), ${rk}$) ) + #:else + nrm = sqrt( sum( a ** 2 , dim = dim ) ) + #:endif + case(NORM_INF) + nrm = maxval( abs(a) , dim = dim ) + case(-NORM_INF) + nrm = minval( abs(a) , dim = dim ) + case (NORM_POW_FIRST:NORM_POW_LAST) + rorder = 1.0_${rk}$ / norm_request + nrm = sum( abs(a) ** norm_request , dim = dim ) ** rorder + case default + err_ = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid norm type after checking') + call linalg_error_handling(err_,err) + end select + + end subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ + + #:endfor + #:endfor + #:endfor + +end module stdlib_linalg_norms From a9f4c7d1c58790b2024e2ae2ea06b2e2a3ffb907 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 13 Sep 2024 10:22:12 +0200 Subject: [PATCH 02/32] fix norms types --- src/stdlib_linalg_norms.fypp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 78db93d4c..4369441f0 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -1,4 +1,7 @@ #: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)) From 98b9b550f7216f44a617d3f0ac820622dcdc8409 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 13 Sep 2024 10:39:46 +0200 Subject: [PATCH 03/32] submodule --- src/stdlib_linalg.fypp | 98 ++++++++++++++++++++++++++++++++++++ src/stdlib_linalg_norms.fypp | 41 +-------------- 2 files changed, 100 insertions(+), 39 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 4574568e5..803a7c7bf 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,102 @@ 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 + #: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 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')}$ + 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 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')}$ + 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 + #: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) :: 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), target :: 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 + 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_norms.fypp b/src/stdlib_linalg_norms.fypp index 4369441f0..8df94ed42 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -6,17 +6,14 @@ #:set INPUT_SHORT = ["char","int"] #:set INPUT_OPTIONS = list(zip(INPUT_TYPE,INPUT_SHORT)) ! Vector norms -module stdlib_linalg_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_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR implicit none(type,external) - private - public :: norm, get_norm - character(*), parameter :: this = 'norm' !> List of internal norm flags @@ -27,40 +24,6 @@ module stdlib_linalg_norms integer(ilp), parameter :: NORM_POW_LAST = NORM_INF - 1_ilp integer(ilp), parameter :: NORM_MINUSINF = -huge(0_ilp) - !> Vector norm: function interface - interface norm - #:for rk,rt,ri in ALL_KINDS_TYPES - #:for it,ii in INPUT_OPTIONS - !> Scalar norms: ${rt}$ - #:for rank in range(1, MAXRANK + 1) - module procedure stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$ - module procedure stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$ - #:endfor - !> Array norms: ${rt}$ - #:for rank in range(2, MAXRANK + 1) - module procedure stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ - module procedure 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 - #:for rk,rt,ri in ALL_KINDS_TYPES - #:for it,ii in INPUT_OPTIONS - !> Scalar norms: ${rt}$ - #:for rank in range(1, MAXRANK + 1) - module procedure norm_${rank}$D_${ii}$_${ri}$ - #:endfor - !> Array norms: ${rt}$ - #:for rank in range(2, MAXRANK + 1) - module procedure norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ - #:endfor - #:endfor - #:endfor - end interface get_norm - interface parse_norm_type module procedure parse_norm_type_integer module procedure parse_norm_type_character @@ -335,4 +298,4 @@ module stdlib_linalg_norms #:endfor #:endfor -end module stdlib_linalg_norms +end submodule stdlib_linalg_norms From c4433cf7cb77578f7c410afd8bfc7574dd0e3fb0 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 13 Sep 2024 10:52:37 +0200 Subject: [PATCH 04/32] add examples --- example/linalg/CMakeLists.txt | 2 ++ example/linalg/example_get_norm.f90 | 51 +++++++++++++++++++++++++++++ example/linalg/example_norm.f90 | 40 ++++++++++++++++++++++ 3 files changed, 93 insertions(+) create mode 100644 example/linalg/example_get_norm.f90 create mode 100644 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..d58f655d6 --- /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..17a4d7625 --- /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 From 7d2fb852cf4667de59954ff7acb2350605567e02 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 13 Sep 2024 11:36:55 +0200 Subject: [PATCH 05/32] add tests --- test/linalg/CMakeLists.txt | 2 + test/linalg/test_linalg_norm.fypp | 202 ++++++++++++++++++++++++++++++ 2 files changed, 204 insertions(+) create mode 100644 test/linalg/test_linalg_norm.fypp 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..d1deffa0a --- /dev/null +++ b/test/linalg/test_linalg_norm.fypp @@ -0,0 +1,202 @@ +#: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 + #: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("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 + #:for rank in range(1, MAXRANK) + + !> 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 + #:for rank in range(2, MAXRANK) + #: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)) 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_norm + + From cd53be4c8b1be467281c95321ce894fa33a99f46 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 13 Sep 2024 12:14:00 +0200 Subject: [PATCH 06/32] document interfaces --- src/stdlib_linalg.fypp | 87 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 87 insertions(+) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 803a7c7bf..2dcc0b10d 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1074,6 +1074,49 @@ module stdlib_linalg #: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}$ @@ -1128,6 +1171,50 @@ module stdlib_linalg !> 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}$ From edf4757b965aa794d0ddc7cb45769d61496a84c8 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 13 Sep 2024 12:46:06 +0200 Subject: [PATCH 07/32] documentation --- doc/specs/stdlib_linalg.md | 108 +++++++++++++++++++++++++++++++++++++ 1 file changed, 108 insertions(+) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index b03fdd192..7c61650be 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1382,3 +1382,111 @@ If `err` is not present, exceptions trigger an `error stop`. {!example/linalg/example_inverse_function.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. + +Result `x` returns a `real`, scalar norm value for the whole array; if `dim` is specified, `x` is a rank-1 +array with the same shape as \(A \) and dimension `dim` collapsed, containing all norms evaluated along `dim`. + +### 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 collapsed. 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!} +``` + +## `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-1 +array with the same shape as \(A \) and dimension `dim` collapsed, 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!} +``` From 35fadc391b9a1a25b88f905ab198bac21d3eb85f Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 13 Sep 2024 12:59:51 +0200 Subject: [PATCH 08/32] attempt artifact v4 --- .github/workflows/ci_windows.yml | 2 +- .github/workflows/doc-deployment.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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/ From 5ba7ef0883db9d613c4109d6d1b28d6e3ee64d05 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 13 Sep 2024 13:09:54 +0200 Subject: [PATCH 09/32] The name of the module procedure conflicts with a name in the encompassing scoping unit --- src/stdlib_linalg_norms.fypp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 8df94ed42..06a0e37cf 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -102,7 +102,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms #:for rank in range(1, MAXRANK + 1) ! Pure function interface, with order specification. On error, the code will stop - pure function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$(a, order) result(nrm) + 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. @@ -115,7 +115,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms end function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$ ! Function interface with output error - function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$(a, order, err) result(nrm) + 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. @@ -130,7 +130,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms end function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$ ! Internal implementation - pure subroutine norm_${rank}$D_${ii}$_${ri}$(a, nrm, order, err) + pure module subroutine norm_${rank}$D_${ii}$_${ri}$(a, nrm, order, err) !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ ${rt}$, intent(in) :: a${ranksuffix(rank)}$ !> Norm of the matrix. @@ -196,7 +196,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms #:for rank in range(2, MAXRANK + 1) ! Pure function interface with DIM specifier. On error, the code will stop - pure function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, order, dim) result(nrm) + 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. @@ -211,7 +211,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ ! Function interface with DIM specifier and output error state. - function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$(a, order, dim, err) result(nrm) + 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. @@ -228,7 +228,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$ ! Internal implementation - pure subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err) + pure module subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err) !> Input matrix a[..] ${rt}$, intent(in), target :: a${ranksuffix(rank)}$ !> Dimension to collapse by computing the norm w.r.t other dimensions From fe23d1e2e2405084f67dcc3bd6bf2adb778e2df2 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 15 Sep 2024 14:42:59 +0200 Subject: [PATCH 10/32] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 7c61650be..b44cea784 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1393,8 +1393,9 @@ Experimental 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. -Result `x` returns a `real`, scalar norm value for the whole array; if `dim` is specified, `x` is a rank-1 -array with the same shape as \(A \) and dimension `dim` collapsed, containing all norms evaluated along `dim`. +Result `x` returns a `real` scalar norm value for the whole array; if `dim` is specified, `x` is an array of rank n-1, +where n equals the rank of ARRAY, and a shape similar to that of \( A \) with dimension `dim` dropped, +containing all norms evaluated along `dim`. ### Syntax From 7ad3ea84c1622330197779803cc3510420eb8386 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 15 Sep 2024 14:43:07 +0200 Subject: [PATCH 11/32] 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 b44cea784..0d0fcce49 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1423,7 +1423,7 @@ containing all norms evaluated along `dim`. 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 collapsed. Each element of `x` contains the 1D norm of the elements of \( A \), +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. From 4e772e134c61499996086e1edcdeb052065d8182 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 15 Sep 2024 14:43:17 +0200 Subject: [PATCH 12/32] 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 0d0fcce49..3d7a29906 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1447,7 +1447,7 @@ Experimental 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-1 +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` collapsed, containing all norms evaluated along `dim`. ### Syntax From b405671d4443916d6673df1e44a37157fef92aa3 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 15 Sep 2024 14:43:24 +0200 Subject: [PATCH 13/32] 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 3d7a29906..d7e4c6641 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1448,7 +1448,7 @@ This `pure subroutine` interface computes one of several vector norms of `real` 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` collapsed, containing all norms evaluated along `dim`. +array with the same shape as \(A \) and dimension `dim` dropped, containing all norms evaluated along `dim`. ### Syntax From aa734de2707b808ccfcc2fc8a390c63543d37d5f Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 15 Sep 2024 14:46:57 +0200 Subject: [PATCH 14/32] improve argument descriptions --- src/stdlib_linalg.fypp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 2dcc0b10d..8f3ef30a9 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1147,9 +1147,9 @@ module stdlib_linalg ${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 + !> Dimension the norm is computed along integer(ilp), intent(in) :: dim - !> Norm of the matrix. + !> 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) @@ -1157,11 +1157,11 @@ module stdlib_linalg ${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 + !> 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. + !> 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 @@ -1235,10 +1235,9 @@ module stdlib_linalg pure module subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err) !> Input matrix a[..] ${rt}$, intent(in), target :: 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`) + !> Dimension the norm is computed along integer(ilp), intent(in) :: dim - !> Norm of the matrix. + !> 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 From 437b96eb020516d25048b810833884f69689ac8d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 15 Sep 2024 14:52:38 +0200 Subject: [PATCH 15/32] use intrinsic `norm2` where possible --- src/stdlib_linalg_norms.fypp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 06a0e37cf..4400579dd 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -144,6 +144,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms integer(ilp) :: sze,norm_request real(${rk}$) :: rorder + intrinsic :: abs, sum, sqrt, norm2, maxval, minval, conjg sze = size(a,kind=ilp) @@ -171,7 +172,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms #:if rt.startswith('complex') nrm = sqrt( real( sum( a * conjg(a) ), ${rk}$) ) #:else - nrm = sqrt( sum( a ** 2 ) ) + nrm = norm2( a ) #:endif case(NORM_INF) nrm = maxval( abs(a) ) @@ -244,6 +245,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms type(linalg_state_type) :: err_ integer(ilp) :: sze,norm_request real(${rk}$) :: rorder + intrinsic :: abs, sum, sqrt, norm2, maxval, minval, conjg sze = size(a,kind=ilp) @@ -278,7 +280,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms #:if rt.startswith('complex') nrm = sqrt( real( sum( a * conjg(a) , dim = dim ), ${rk}$) ) #:else - nrm = sqrt( sum( a ** 2 , dim = dim ) ) + nrm = norm2( a , dim = dim ) #:endif case(NORM_INF) nrm = maxval( abs(a) , dim = dim ) From de4f8d386eef0c45535993b8f2c05b293ef58e06 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 17 Sep 2024 13:36:36 +0200 Subject: [PATCH 16/32] 2-norm: use BLAS on contiguous or strided arrays if possible - add nonstandard-named `complex` norms to `nrm2` interface - test sliced and reshaped 2-norm --- src/stdlib_linalg.fypp | 2 +- src/stdlib_linalg_blas.fypp | 26 +++++++++++++++++- src/stdlib_linalg_norms.fypp | 44 ++++++++++++++++++++++++++----- test/linalg/test_linalg_norm.fypp | 36 +++++++++++++++++++++++++ 4 files changed, 99 insertions(+), 9 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 8f3ef30a9..9ec5213e1 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1221,7 +1221,7 @@ module stdlib_linalg #: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) :: 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. diff --git a/src/stdlib_linalg_blas.fypp b/src/stdlib_linalg_blas.fypp index 504a15108..38e1ca87c 100644 --- a/src/stdlib_linalg_blas.fypp +++ b/src/stdlib_linalg_blas.fypp @@ -974,12 +974,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 +1003,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 index 4400579dd..858832cbd 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -8,10 +8,11 @@ ! Vector norms submodule(stdlib_linalg) stdlib_linalg_norms use stdlib_linalg_constants - use stdlib_linalg_blas, only: nrm2 + 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 implicit none(type,external) character(*), parameter :: this = 'norm' @@ -29,6 +30,13 @@ submodule(stdlib_linalg) stdlib_linalg_norms 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 @@ -93,6 +101,25 @@ submodule(stdlib_linalg) stdlib_linalg_norms 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}$ + + #:for it,ii in INPUT_OPTIONS !============================================== @@ -132,7 +159,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms ! Internal implementation pure module subroutine norm_${rank}$D_${ii}$_${ri}$(a, nrm, order, err) !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ - ${rt}$, intent(in) :: 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. @@ -142,9 +169,10 @@ submodule(stdlib_linalg) stdlib_linalg_norms type(linalg_state_type) :: err_ - integer(ilp) :: sze,norm_request + integer(ilp) :: sze,norm_request,str real(${rk}$) :: rorder - intrinsic :: abs, sum, sqrt, norm2, maxval, minval, conjg + ${rt}$, pointer :: a1d(:) + intrinsic :: abs, sum, sqrt, maxval, minval, conjg sze = size(a,kind=ilp) @@ -169,10 +197,12 @@ submodule(stdlib_linalg) stdlib_linalg_norms case(NORM_ONE) nrm = sum( abs(a) ) case(NORM_TWO) - #:if rt.startswith('complex') - nrm = sqrt( real( sum( a * conjg(a) ), ${rk}$) ) + #:if rank==1 + nrm = nrm2(sze,a,incx=stride_1d(a)) + #:elif rt.startswith('complex') + nrm = sqrt( real( sum( a * conjg(a) ), ${rk}$) ) #:else - nrm = norm2( a ) + nrm = norm2(a) #:endif case(NORM_INF) nrm = maxval( abs(a) ) diff --git a/test/linalg/test_linalg_norm.fypp b/test/linalg/test_linalg_norm.fypp index d1deffa0a..e24c324e1 100644 --- a/test/linalg/test_linalg_norm.fypp +++ b/test/linalg/test_linalg_norm.fypp @@ -32,6 +32,7 @@ module test_linalg_norm 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 @@ -46,6 +47,41 @@ module test_linalg_norm 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 From 37616ad68c7402fae24afac94ac9737f69f06755 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 18 Sep 2024 12:27:24 +0200 Subject: [PATCH 17/32] do not use strides --- src/stdlib_linalg_norms.fypp | 6 +++--- test/linalg/test_linalg_norm.fypp | 15 +++++++++++++++ 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 858832cbd..1793ab919 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -8,7 +8,7 @@ ! Vector norms submodule(stdlib_linalg) stdlib_linalg_norms use stdlib_linalg_constants - use stdlib_linalg_blas!, only: nrm2 + 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 @@ -169,7 +169,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms type(linalg_state_type) :: err_ - integer(ilp) :: sze,norm_request,str + integer(ilp) :: sze,norm_request real(${rk}$) :: rorder ${rt}$, pointer :: a1d(:) intrinsic :: abs, sum, sqrt, maxval, minval, conjg @@ -198,7 +198,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms nrm = sum( abs(a) ) case(NORM_TWO) #:if rank==1 - nrm = nrm2(sze,a,incx=stride_1d(a)) + nrm = nrm2(sze,a,incx=1) #:elif rt.startswith('complex') nrm = sqrt( real( sum( a * conjg(a) ), ${rk}$) ) #:else diff --git a/test/linalg/test_linalg_norm.fypp b/test/linalg/test_linalg_norm.fypp index e24c324e1..3ae08822b 100644 --- a/test/linalg/test_linalg_norm.fypp +++ b/test/linalg/test_linalg_norm.fypp @@ -80,6 +80,21 @@ module test_linalg_norm '2d-reshaped ${rt}$ norm(a,2)') if (allocated(error)) return + ! Test row norm (strided access) + slice = twod(3,:) + call check(error,abs(norm(twod(3,:),2)-norm(slice,2)) Date: Thu, 19 Sep 2024 13:00:49 +0200 Subject: [PATCH 18/32] interface to, use and test 1-norm `asum` from BLAS --- src/stdlib_linalg_blas.fypp | 54 +++++++++++++++++++++++++++++++ src/stdlib_linalg_norms.fypp | 9 ++++-- test/linalg/test_linalg_norm.fypp | 7 ++++ 3 files changed, 67 insertions(+), 3 deletions(-) diff --git a/src/stdlib_linalg_blas.fypp b/src/stdlib_linalg_blas.fypp index 38e1ca87c..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 diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 1793ab919..79a3536fb 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -8,7 +8,7 @@ ! Vector norms submodule(stdlib_linalg) stdlib_linalg_norms use stdlib_linalg_constants - use stdlib_linalg_blas, only: nrm2 + use stdlib_linalg_blas, only: nrm2,asum 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 @@ -171,7 +171,6 @@ submodule(stdlib_linalg) stdlib_linalg_norms integer(ilp) :: sze,norm_request real(${rk}$) :: rorder - ${rt}$, pointer :: a1d(:) intrinsic :: abs, sum, sqrt, maxval, minval, conjg sze = size(a,kind=ilp) @@ -195,10 +194,14 @@ submodule(stdlib_linalg) stdlib_linalg_norms select case(norm_request) case(NORM_ONE) + #:if rank==1 + nrm = asum(sze,a,incx=1_ilp) + #:else nrm = sum( abs(a) ) + #:endif case(NORM_TWO) #:if rank==1 - nrm = nrm2(sze,a,incx=1) + nrm = nrm2(sze,a,incx=1_ilp) #:elif rt.startswith('complex') nrm = sqrt( real( sum( a * conjg(a) ), ${rk}$) ) #:else diff --git a/test/linalg/test_linalg_norm.fypp b/test/linalg/test_linalg_norm.fypp index 3ae08822b..17495962c 100644 --- a/test/linalg/test_linalg_norm.fypp +++ b/test/linalg/test_linalg_norm.fypp @@ -212,6 +212,13 @@ module test_linalg_norm if (allocated(error)) return end do + + ! Compare ND whole vector norm with unrolled vector norm + write(msg,"('Unrolled (1d) vs ${rank}$d order=',i0,' norm')") order + call check(error,abs(norm(a,order)-norm(b,order)) Date: Mon, 23 Sep 2024 10:48:48 +0200 Subject: [PATCH 19/32] drop duplicate descr --- doc/specs/stdlib_linalg.md | 4 ---- 1 file changed, 4 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index d7e4c6641..22569246d 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1393,10 +1393,6 @@ Experimental 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. -Result `x` returns a `real` scalar norm value for the whole array; if `dim` is specified, `x` is an array of rank n-1, -where n equals the rank of ARRAY, and a shape similar to that of \( A \) with dimension `dim` dropped, -containing all norms evaluated along `dim`. - ### Syntax `x = ` [[stdlib_linalg(module):norm(interface)]] `(a, order, [, dim, err])` From 9e020a84ae0eee442d3a2c3d9e125f1c3d956e0b Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 23 Sep 2024 10:49:54 +0200 Subject: [PATCH 20/32] move `get_norm` before `norm` --- doc/specs/stdlib_linalg.md | 42 ++++++++++++++++++++------------------ 1 file changed, 22 insertions(+), 20 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 22569246d..5f956b234 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1382,7 +1382,7 @@ If `err` is not present, exceptions trigger an `error stop`. {!example/linalg/example_inverse_function.f90!} ``` -## `norm` - Computes the vector norm of a generic-rank array. +## `get_norm` - Computes the vector norm of a generic-rank array. ### Status @@ -1390,17 +1390,23 @@ Experimental ### Description -This function computes one of several vector norms of `real` or `complex` array \( A \), depending on +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 -`x = ` [[stdlib_linalg(module):norm(interface)]] `(a, order, [, dim, err])` +`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 | @@ -1417,9 +1423,9 @@ the value of the `order` input argument. \( A \) may be an array of any rank. ### 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 \), +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. @@ -1429,10 +1435,10 @@ If `err` is not present, exceptions trigger an `error stop`. ### Example ```fortran -{!example/linalg/example_norm.f90!} +{!example/linalg/example_get_norm.f90!} ``` -## `get_norm` - Computes the vector norm of a generic-rank array. +## `norm` - Computes the vector norm of a generic-rank array. ### Status @@ -1440,23 +1446,17 @@ Experimental ### Description -This `pure subroutine` interface computes one of several vector norms of `real` or `complex` array \( A \), depending on +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. -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])` +`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. -`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 | @@ -1473,9 +1473,9 @@ to that of `a` with dimension `dim` dropped. ### 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 \), +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. @@ -1485,5 +1485,7 @@ If `err` is not present, exceptions trigger an `error stop`. ### Example ```fortran -{!example/linalg/example_get_norm.f90!} +{!example/linalg/example_norm.f90!} ``` + + From 0dabeac40e61180790ca157b1f8a8c37ff0589f7 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 23 Sep 2024 10:50:03 +0200 Subject: [PATCH 21/32] Update src/stdlib_linalg_norms.fypp Co-authored-by: jalvesz <102541118+jalvesz@users.noreply.github.com> --- src/stdlib_linalg_norms.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 79a3536fb..12b01bcf8 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -53,7 +53,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms norm_type = NORM_ONE case (2_ilp) norm_type = NORM_TWO - case (3_ilp:huge(0_ilp)-1_ilp) + case (3_ilp:NORM_POW_LAST) norm_type = order case (huge(0_ilp):) norm_type = NORM_INF From 54caa7f561e204d4e58e07f8d8cfde5266cdd427 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 23 Sep 2024 10:50:18 +0200 Subject: [PATCH 22/32] Update src/stdlib_linalg_norms.fypp Co-authored-by: jalvesz <102541118+jalvesz@users.noreply.github.com> --- src/stdlib_linalg_norms.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 12b01bcf8..fed64f536 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -55,7 +55,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms norm_type = NORM_TWO case (3_ilp:NORM_POW_LAST) norm_type = order - case (huge(0_ilp):) + case (NORM_INF:) norm_type = NORM_INF case (:-huge(0_ilp)) norm_type = NORM_MINUSINF From 07f05254f01c83cd30f8004f57592deb052150f7 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 23 Sep 2024 10:50:25 +0200 Subject: [PATCH 23/32] Update src/stdlib_linalg_norms.fypp Co-authored-by: jalvesz <102541118+jalvesz@users.noreply.github.com> --- src/stdlib_linalg_norms.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index fed64f536..c3678999f 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -57,7 +57,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms norm_type = order case (NORM_INF:) norm_type = NORM_INF - case (:-huge(0_ilp)) + case (:NORM_MINUSINF) norm_type = NORM_MINUSINF case default From a94647c297fc15ec4af871cd41fc24104b8ce4cb Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 23 Sep 2024 10:50:35 +0200 Subject: [PATCH 24/32] Update src/stdlib_linalg_norms.fypp Co-authored-by: jalvesz <102541118+jalvesz@users.noreply.github.com> --- src/stdlib_linalg_norms.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index c3678999f..1a8fe6b38 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -209,7 +209,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms #:endif case(NORM_INF) nrm = maxval( abs(a) ) - case(-NORM_INF) + case(NORM_MINUSINF) nrm = minval( abs(a) ) case (NORM_POW_FIRST:NORM_POW_LAST) rorder = 1.0_${rk}$ / norm_request From 13636a90aa3828457b0fc56eb2e8115e95f870a9 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 23 Sep 2024 10:50:43 +0200 Subject: [PATCH 25/32] Update example/linalg/example_get_norm.f90 Co-authored-by: jalvesz <102541118+jalvesz@users.noreply.github.com> --- example/linalg/example_get_norm.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/linalg/example_get_norm.f90 b/example/linalg/example_get_norm.f90 index d58f655d6..1d1e779bb 100644 --- a/example/linalg/example_get_norm.f90 +++ b/example/linalg/example_get_norm.f90 @@ -44,7 +44,7 @@ program example_get_norm 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] + ! [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() From 92ab520c3a2388b56e420c0fcac066db017c25b4 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 23 Sep 2024 10:50:50 +0200 Subject: [PATCH 26/32] Update src/stdlib_linalg_norms.fypp Co-authored-by: jalvesz <102541118+jalvesz@users.noreply.github.com> --- src/stdlib_linalg_norms.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 1a8fe6b38..398912d6b 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -317,7 +317,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms #:endif case(NORM_INF) nrm = maxval( abs(a) , dim = dim ) - case(-NORM_INF) + case(NORM_MINUSINF) nrm = minval( abs(a) , dim = dim ) case (NORM_POW_FIRST:NORM_POW_LAST) rorder = 1.0_${rk}$ / norm_request From 91882fa9bf68871f2259033f01a26fff1389b8e4 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 23 Sep 2024 10:50:56 +0200 Subject: [PATCH 27/32] Update example/linalg/example_get_norm.f90 Co-authored-by: jalvesz <102541118+jalvesz@users.noreply.github.com> --- example/linalg/example_get_norm.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/linalg/example_get_norm.f90 b/example/linalg/example_get_norm.f90 index 1d1e779bb..f150a9033 100644 --- a/example/linalg/example_get_norm.f90 +++ b/example/linalg/example_get_norm.f90 @@ -3,7 +3,7 @@ program example_get_norm use stdlib_linalg, only: get_norm, linalg_state_type implicit none - real :: a(3,3),nrm,nrmd(3) + real :: a(3,3), nrm, nrmd(3) integer :: j type(linalg_state_type) :: err From 1a79128d26fd55c103023d5a901fe59bc775b6d9 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 23 Sep 2024 10:51:03 +0200 Subject: [PATCH 28/32] Update example/linalg/example_norm.f90 Co-authored-by: jalvesz <102541118+jalvesz@users.noreply.github.com> --- example/linalg/example_norm.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/linalg/example_norm.f90 b/example/linalg/example_norm.f90 index 17a4d7625..43dc7e3eb 100644 --- a/example/linalg/example_norm.f90 +++ b/example/linalg/example_norm.f90 @@ -33,7 +33,7 @@ program example_norm 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] + ! [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() From 5f6b5e894745fc23b2e06e05ba414fe75eadec87 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 27 Sep 2024 10:21:08 +0200 Subject: [PATCH 29/32] ND arrays: use LAPACK internals via loops --- include/common.fypp | 101 ++++++++++++++++++++++++++++++ src/stdlib_linalg.fypp | 2 +- src/stdlib_linalg_norms.fypp | 86 ++++++++++++++++++------- test/linalg/test_linalg_norm.fypp | 6 +- 4 files changed, 167 insertions(+), 28 deletions(-) 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/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 9ec5213e1..a26a10630 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1234,7 +1234,7 @@ module stdlib_linalg #: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), target :: a${ranksuffix(rank)}$ + ${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). diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 398912d6b..75c607ec8 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -264,7 +264,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms ! 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), target :: a${ranksuffix(rank)}$ + ${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 @@ -276,11 +276,17 @@ submodule(stdlib_linalg) stdlib_linalg_norms type(linalg_state_type), intent(out), optional :: err type(linalg_state_type) :: err_ - integer(ilp) :: sze,norm_request + 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 - sze = size(a,kind=ilp) + ! Input matrix properties + sze = size (a,kind=ilp) + spe = shape(a,kind=ilp) ! Initialize norm to zero nrm = 0.0_${rk}$ @@ -304,28 +310,60 @@ submodule(stdlib_linalg) stdlib_linalg_norms if (err_%error()) then call linalg_error_handling(err_,err) return - endif + endif - select case(norm_request) - case(NORM_ONE) - nrm = sum( abs(a) , dim = dim ) - case(NORM_TWO) - #:if rt.startswith('complex') - nrm = sqrt( real( sum( a * conjg(a) , dim = dim ), ${rk}$) ) - #:else - nrm = norm2( a , dim = dim ) - #:endif - case(NORM_INF) - nrm = maxval( abs(a) , dim = dim ) - case(NORM_MINUSINF) - nrm = minval( abs(a) , dim = dim ) - case (NORM_POW_FIRST:NORM_POW_LAST) - rorder = 1.0_${rk}$ / norm_request - nrm = sum( abs(a) ** norm_request , dim = dim ) ** rorder - case default - err_ = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid norm type after checking') - call linalg_error_handling(err_,err) - end select + ! The norm's leading dimension + lda = spe(dim) + + ! Check if input column data is contiguous + contiguous_data = dim==1 .or. all(norm_request/=[NORM_ONE,NORM_TWO]) + + ! 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)}$ + select case(norm_request) + case(NORM_ONE) + nrm(${loop_variables('j',rank-1,1)}$) = & + asum(lda,apack(:, ${loop_variables('j',rank-1,1)}$),incx=1_ilp) + case(NORM_TWO) + nrm(${loop_variables('j',rank-1,1)}$) = & + nrm2(lda,apack(:, ${loop_variables('j',rank-1,1)}$),incx=1_ilp) + end select +${loop_variables_end(rank-1," "*12)}$ + + else + + select case(norm_request) + case(NORM_ONE) +${loop_variables_start('j', 'a', rank-1, 1," "*20)}$ + nrm(${loop_variables('j',rank-1,1)}$) = & + asum(lda,a(:, ${loop_variables('j',rank-1,1)}$),incx=1_ilp) +${loop_variables_end(rank-1," "*20)}$ + case(NORM_TWO) +${loop_variables_start('j', 'a', rank-1, 1," "*20)}$ + nrm(${loop_variables('j',rank-1,1)}$) = & + nrm2(lda,a(:, ${loop_variables('j',rank-1,1)}$),incx=1_ilp) +${loop_variables_end(rank-1," "*20)}$ + case(NORM_INF) + nrm = maxval( abs(a) , dim = dim ) + case(NORM_MINUSINF) + nrm = minval( abs(a) , dim = dim ) + case (NORM_POW_FIRST:NORM_POW_LAST) + rorder = 1.0_${rk}$ / norm_request + nrm = sum( abs(a) ** norm_request , dim = dim ) ** rorder + case default + err_ = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid norm type after checking') + call linalg_error_handling(err_,err) + end select + + endif end subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ diff --git a/test/linalg/test_linalg_norm.fypp b/test/linalg/test_linalg_norm.fypp index 17495962c..ef660a1ff 100644 --- a/test/linalg/test_linalg_norm.fypp +++ b/test/linalg/test_linalg_norm.fypp @@ -159,17 +159,17 @@ module test_linalg_norm ! Test some norms call check(error,abs(norm(a,2) - norm2(a)) Date: Fri, 27 Sep 2024 11:23:42 +0200 Subject: [PATCH 30/32] streamline ND->1D norm functions --- src/stdlib_linalg_norms.fypp | 81 ++++++++++++++++++++++-------------- 1 file changed, 49 insertions(+), 32 deletions(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index 75c607ec8..c387ba3b7 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,asum - 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 @@ -119,6 +119,47 @@ submodule(stdlib_linalg) stdlib_linalg_norms end function stride_1d_${ri}$ + ! Private internal implementation: 1D + pure subroutine internal_norm_1D_${ri}$(sze, a, nrm, norm_request, err) + !> 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 + !> State return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(inout) :: err + + 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 @@ -155,8 +196,8 @@ submodule(stdlib_linalg) stdlib_linalg_norms 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 + + ! 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)}$ @@ -168,7 +209,6 @@ submodule(stdlib_linalg) stdlib_linalg_norms 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 @@ -190,34 +230,11 @@ submodule(stdlib_linalg) stdlib_linalg_norms if (err_%error()) then call linalg_error_handling(err_,err) return - endif + endif - select case(norm_request) - case(NORM_ONE) - #:if rank==1 - nrm = asum(sze,a,incx=1_ilp) - #:else - nrm = sum( abs(a) ) - #:endif - case(NORM_TWO) - #:if rank==1 - nrm = nrm2(sze,a,incx=1_ilp) - #:elif rt.startswith('complex') - nrm = sqrt( real( sum( a * conjg(a) ), ${rk}$) ) - #:else - nrm = norm2(a) - #:endif - case(NORM_INF) - nrm = maxval( abs(a) ) - 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 - case default - err_ = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid norm type after checking') - call linalg_error_handling(err_,err) - end select + ! Get norm + call internal_norm_1D_${ri}$(sze, a, nrm, norm_request, err_) + call linalg_error_handling(err_,err) end subroutine norm_${rank}$D_${ii}$_${ri}$ From a65e7713921fe34c1943850587e616110d203d8e Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 27 Sep 2024 11:28:20 +0200 Subject: [PATCH 31/32] streamline `dim`-med norm functions --- src/stdlib_linalg_norms.fypp | 47 +++++++++--------------------------- 1 file changed, 11 insertions(+), 36 deletions(-) diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp index c387ba3b7..8e3fba4a4 100644 --- a/src/stdlib_linalg_norms.fypp +++ b/src/stdlib_linalg_norms.fypp @@ -119,8 +119,9 @@ submodule(stdlib_linalg) stdlib_linalg_norms end function stride_1d_${ri}$ - ! Private internal implementation: 1D - pure subroutine internal_norm_1D_${ri}$(sze, a, nrm, norm_request, err) + ! 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(*) @@ -129,8 +130,6 @@ submodule(stdlib_linalg) stdlib_linalg_norms real(${rk}$), intent(out) :: nrm !> Internal matrix request integer(ilp), intent(in) :: norm_request - !> State return flag. On error if not requested, the code will stop - type(linalg_state_type), intent(inout) :: err integer(ilp) :: i real(${rk}$) :: rorder @@ -233,7 +232,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms endif ! Get norm - call internal_norm_1D_${ri}$(sze, a, nrm, norm_request, err_) + call internal_norm_1D_${ri}$(sze, a, nrm, norm_request) call linalg_error_handling(err_,err) end subroutine norm_${rank}$D_${ii}$_${ri}$ @@ -333,7 +332,7 @@ submodule(stdlib_linalg) stdlib_linalg_norms lda = spe(dim) ! Check if input column data is contiguous - contiguous_data = dim==1 .or. all(norm_request/=[NORM_ONE,NORM_TWO]) + contiguous_data = dim==1 ! Get packed data with the norm dimension as the first dimension if (.not.contiguous_data) then @@ -345,40 +344,16 @@ submodule(stdlib_linalg) stdlib_linalg_norms apack = reshape(a, shape=spack, order=iperm) ${loop_variables_start('j', 'apack', rank-1, 1," "*12)}$ - select case(norm_request) - case(NORM_ONE) - nrm(${loop_variables('j',rank-1,1)}$) = & - asum(lda,apack(:, ${loop_variables('j',rank-1,1)}$),incx=1_ilp) - case(NORM_TWO) - nrm(${loop_variables('j',rank-1,1)}$) = & - nrm2(lda,apack(:, ${loop_variables('j',rank-1,1)}$),incx=1_ilp) - end select + 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 - select case(norm_request) - case(NORM_ONE) -${loop_variables_start('j', 'a', rank-1, 1," "*20)}$ - nrm(${loop_variables('j',rank-1,1)}$) = & - asum(lda,a(:, ${loop_variables('j',rank-1,1)}$),incx=1_ilp) -${loop_variables_end(rank-1," "*20)}$ - case(NORM_TWO) -${loop_variables_start('j', 'a', rank-1, 1," "*20)}$ - nrm(${loop_variables('j',rank-1,1)}$) = & - nrm2(lda,a(:, ${loop_variables('j',rank-1,1)}$),incx=1_ilp) -${loop_variables_end(rank-1," "*20)}$ - case(NORM_INF) - nrm = maxval( abs(a) , dim = dim ) - case(NORM_MINUSINF) - nrm = minval( abs(a) , dim = dim ) - case (NORM_POW_FIRST:NORM_POW_LAST) - rorder = 1.0_${rk}$ / norm_request - nrm = sum( abs(a) ** norm_request , dim = dim ) ** rorder - case default - err_ = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid norm type after checking') - call linalg_error_handling(err_,err) - end select +${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 From 39d6daaebea66b2665845899402ab279e1b85c0d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 27 Sep 2024 11:32:39 +0200 Subject: [PATCH 32/32] test infinity norm vs. `maxval(abs(.))` --- test/linalg/test_linalg_norm.fypp | 44 +++++++++++++++++++++++++++++-- 1 file changed, 42 insertions(+), 2 deletions(-) diff --git a/test/linalg/test_linalg_norm.fypp b/test/linalg/test_linalg_norm.fypp index ef660a1ff..0ec88c2ec 100644 --- a/test/linalg/test_linalg_norm.fypp +++ b/test/linalg/test_linalg_norm.fypp @@ -40,6 +40,7 @@ module test_linalg_norm #: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 @@ -135,9 +136,9 @@ module test_linalg_norm end subroutine test_norm_${ri}$_${rank}$d #:endfor - - !> Test Euclidean norm; compare with Fortran intrinsic norm2 for reals + #:for rank in range(2, MAXRANK) + !> 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 @@ -178,6 +179,45 @@ module test_linalg_norm end subroutine test_norm2_${ri}$_${rank}$d #:endif + !> 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)))