diff --git a/config/fypp_deployment.py b/config/fypp_deployment.py index 9f3549e4f..ee7cb02cd 100644 --- a/config/fypp_deployment.py +++ b/config/fypp_deployment.py @@ -122,11 +122,11 @@ def fpm_build(args,unknown): flags= flags + unknown[idx+1] #========================================== # build with fpm - subprocess.run(["fpm build"]+ - [" --compiler "]+[FPM_FC]+ - [" --c-compiler "]+[FPM_CC]+ - [" --cxx-compiler "]+[FPM_CXX]+ - [" --flag "]+[flags], shell=True, check=True) + subprocess.run("fpm build"+ + " --compiler "+FPM_FC+ + " --c-compiler "+FPM_CC+ + " --cxx-compiler "+FPM_CXX+ + " --flag \"{}\"".format(flags), shell=True, check=True) return if __name__ == "__main__": @@ -137,7 +137,7 @@ def fpm_build(args,unknown): parser.add_argument("--vpatch", type=int, default=0, help="Project Version Patch") parser.add_argument("--njob", type=int, default=4, help="Number of parallel jobs for preprocessing") - parser.add_argument("--maxrank",type=int, default=7, help="Set the maximum allowed rank for arrays") + parser.add_argument("--maxrank",type=int, default=4, help="Set the maximum allowed rank for arrays") parser.add_argument("--with_qp",action='store_true', help="Include WITH_QP in the command") parser.add_argument("--with_xdp",action='store_true', help="Include WITH_XDP in the command") parser.add_argument("--lnumbering",action='store_true', help="Add line numbering in preprocessed files") diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 01a19e624..2f03f0fad 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1459,3 +1459,110 @@ If `err` is not present, exceptions trigger an `error stop`. {!example/linalg/example_inverse_function.f90!} ``` +## `get_norm` - Computes the vector norm of a generic-rank array. + +### Status + +Experimental + +### Description + +This `pure subroutine` interface computes one of several vector norms of `real` or `complex` array \( A \), depending on +the value of the `order` input argument. \( A \) may be an array of any rank. + +Result `nrm` returns a `real`, scalar norm value for the whole array; if `dim` is specified, `nrm` is a rank n-1 +array with the same shape as \(A \) and dimension `dim` dropped, containing all norms evaluated along `dim`. + +### Syntax + +`call ` [[stdlib_linalg(module):get_norm(interface)]] `(a, nrm, order, [, dim, err])` + +### Arguments + +`a`: Shall be a rank-n `real` or `complex` array containing the data. It is an `intent(in)` argument. + +`nrm`: if `dim` is absent, shall be a scalar with the norm evaluated over all the elements of the array. Otherwise, an array of rank `n-1`, and a shape similar +to that of `a` with dimension `dim` dropped. + +`order`: Shall be an `integer` value or a `character` flag that specifies the norm type, as follows. It is an `intent(in)` argument. + +| Integer input | Character Input | Norm type | +|------------------|------------------|---------------------------------------------------------| +| `-huge(0)` | `'-inf', '-Inf'` | Minimum absolute value \( \min_i{ \left|a_i\right| } \) | +| `1` | `'1'` | 1-norm \( \sum_i{ \left|a_i\right| } \) | +| `2` | `'2'` | Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \) | +| `>=3` | `'3','4',...` | p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \) | +| `huge(0)` | `'inf', 'Inf'` | Maximum absolute value \( \max_i{ \left|a_i\right| } \) | + +`dim` (optional): Shall be a scalar `integer` value with a value in the range from `1` to `n`, where `n` is the rank of the array. It is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. If `err` is not present, the function is `pure`. + +### Return value + +By default, the return value `nrm` is a scalar, and contains the norm as evaluated over all elements of the generic-rank array \( A \). +If the optional `dim` argument is present, `nrm` is a rank `n-1` array with the same shape as \( A \) except +for dimension `dim`, that is collapsed. Each element of `nrm` contains the 1D norm of the elements of \( A \), +evaluated along dimension `dim` only. + +Raises `LINALG_ERROR` if the requested norm type is invalid. +Raises `LINALG_VALUE_ERROR` if any of the arguments has an invalid size. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_get_norm.f90!} +``` + +## `norm` - Computes the vector norm of a generic-rank array. + +### Status + +Experimental + +### Description + +This function computes one of several vector norms of `real` or `complex` array \( A \), depending on +the value of the `order` input argument. \( A \) may be an array of any rank. + +### Syntax + +`x = ` [[stdlib_linalg(module):norm(interface)]] `(a, order, [, dim, err])` + +### Arguments + +`a`: Shall be a rank-n `real` or `complex` array containing the data. It is an `intent(in)` argument. + +`order`: Shall be an `integer` value or a `character` flag that specifies the norm type, as follows. It is an `intent(in)` argument. + +| Integer input | Character Input | Norm type | +|------------------|------------------|---------------------------------------------------------| +| `-huge(0)` | `'-inf', '-Inf'` | Minimum absolute value \( \min_i{ \left|a_i\right| } \) | +| `1` | `'1'` | 1-norm \( \sum_i{ \left|a_i\right| } \) | +| `2` | `'2'` | Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \) | +| `>=3` | `'3','4',...` | p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \) | +| `huge(0)` | `'inf', 'Inf'` | Maximum absolute value \( \max_i{ \left|a_i\right| } \) | + +`dim` (optional): Shall be a scalar `integer` value with a value in the range from `1` to `n`, where `n` is the rank of the array. It is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. If `err` is not present, the function is `pure`. + +### Return value + +By default, the return value `x` is a scalar, and contains the norm as evaluated over all elements of the generic-rank array \( A \). +If the optional `dim` argument is present, `x` is a rank `n-1` array with the same shape as \( A \) except +for dimension `dim`, that is dropped. Each element of `x` contains the 1D norm of the elements of \( A \), +evaluated along dimension `dim` only. + +Raises `LINALG_ERROR` if the requested norm type is invalid. +Raises `LINALG_VALUE_ERROR` if any of the arguments has an invalid size. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_norm.f90!} +``` + + diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 487e2197b..e12236375 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -28,6 +28,8 @@ ADD_EXAMPLE(blas_gemv) ADD_EXAMPLE(lapack_getrf) ADD_EXAMPLE(lstsq1) ADD_EXAMPLE(lstsq2) +ADD_EXAMPLE(norm) +ADD_EXAMPLE(get_norm) ADD_EXAMPLE(solve1) ADD_EXAMPLE(solve2) ADD_EXAMPLE(solve3) diff --git a/example/linalg/example_get_norm.f90 b/example/linalg/example_get_norm.f90 new file mode 100644 index 000000000..f150a9033 --- /dev/null +++ b/example/linalg/example_get_norm.f90 @@ -0,0 +1,51 @@ +! Vector norm: demonstrate usage of the function interface +program example_get_norm + use stdlib_linalg, only: get_norm, linalg_state_type + implicit none + + real :: a(3,3), nrm, nrmd(3) + integer :: j + type(linalg_state_type) :: err + + ! a = [ -3.00000000 0.00000000 3.00000000 + ! -2.00000000 1.00000000 4.00000000 + ! -1.00000000 2.00000000 5.00000000 ] + a = reshape([(j-4,j=1,9)], [3,3]) + + print "(' a = [ ',3(g0,3x),2(/9x,3(g0,3x)),']')", transpose(a) + + ! Norm with integer input + call get_norm(a, nrm, 2) + print *, 'Euclidean norm = ',nrm ! 8.30662346 + + ! Norm with character input + call get_norm(a, nrm, '2') + print *, 'Euclidean norm = ',nrm ! 8.30662346 + + ! Euclidean norm of row arrays, a(i,:) + call get_norm(a, nrmd, 2, dim=2) + print *, 'Rows norms = ',nrmd ! 4.24264050 4.58257580 5.47722578 + + ! Euclidean norm of columns arrays, a(:,i) + call get_norm(a, nrmd, 2, dim=1) + print *, 'Columns norms = ',nrmd ! 3.74165750 2.23606801 7.07106781 + + ! Infinity norms + call get_norm(a, nrm, 'inf') + print *, 'maxval(||a||) = ',nrm ! 5.00000000 + + call get_norm(a, nrmd, 'inf', dim=2) + print *, 'maxval(||a(i,:)||) = ',nrmd ! 3.00000000 4.00000000 5.00000000 + + call get_norm(a, nrm, '-inf') + print *, 'minval(||a||) = ',nrm ! 0.00000000 + + call get_norm(a, nrmd, '-inf', dim=1) + print *, 'minval(||a(:,i)||) = ',nrmd ! 1.00000000 0.00000000 3.00000000 + + ! Catch Error: + ! [norm] returned Value Error: dimension 4 is out of rank for shape(a)= [3, 3] + call get_norm(a, nrmd, 'inf', dim=4, err=err) + print *, 'invalid: ',err%print() + +end program example_get_norm diff --git a/example/linalg/example_norm.f90 b/example/linalg/example_norm.f90 new file mode 100644 index 000000000..43dc7e3eb --- /dev/null +++ b/example/linalg/example_norm.f90 @@ -0,0 +1,40 @@ +! Vector norm: demonstrate usage of the function interface +program example_norm + use stdlib_linalg, only: norm, linalg_state_type + implicit none + + real :: a(3,3),na + integer :: j + type(linalg_state_type) :: err + + ! a = [ -3.00000000 0.00000000 3.00000000 + ! -2.00000000 1.00000000 4.00000000 + ! -1.00000000 2.00000000 5.00000000 ] + a = reshape([(j-4,j=1,9)], [3,3]) + + print "(' a = [ ',3(g0,3x),2(/9x,3(g0,3x)),']')", transpose(a) + + ! Norm with integer input + print *, 'Euclidean norm = ',norm(a, 2) ! 8.30662346 + + ! Norm with character input + print *, 'Euclidean norm = ',norm(a, '2') ! 8.30662346 + + ! Euclidean norm of row arrays, a(i,:) + print *, 'Rows norms = ',norm(a, 2, dim=2) ! 4.24264050 4.58257580 5.47722578 + + ! Euclidean norm of columns arrays, a(:,i) + print *, 'Columns norms = ',norm(a, 2, dim=1) ! 3.74165750 2.23606801 7.07106781 + + ! Infinity norms + print *, 'maxval(||a||) = ',norm(a, 'inf') ! 5.00000000 + print *, 'maxval(||a(i,:)||) = ',norm(a, 'inf', dim=2) ! 3.00000000 4.00000000 5.00000000 + print *, 'minval(||a||) = ',norm(a, '-inf') ! 0.00000000 + print *, 'minval(||a(:,i)||) = ',norm(a, '-inf', dim=1) ! 1.00000000 0.00000000 3.00000000 + + ! Catch Error: + ! [norm] returned Value Error: dimension 4 is out of rank for shape(a)= [3, 3] + print *, 'invalid: ',norm(a,'inf', dim=4, err=err) + print *, 'error = ',err%print() + +end program example_norm diff --git a/include/common.fypp b/include/common.fypp index bd71f7394..d43708cb1 100644 --- a/include/common.fypp +++ b/include/common.fypp @@ -321,4 +321,105 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #:endif #:enddef +#! +#! Generates a list of loop variables +#! +#! Args: +#! varname(str): Name of the variable to be used as prefix +#! n (int): Number of loop variables to be created +#! offset (int): Optional index offset +#! +#! Returns: +#! Variable definition string +#! +#! E.g., +#! loop_variables('j', 5) +#! -> "j1, j2, j3, j4, j5 +#! +#:def loop_variables(varname, n, offset=0) + #:assert n > 0 + #:call join_lines(joinstr=", ") + #:for i in range(1, n + 1) + ${varname}$${i+offset}$ + #:endfor + #:endcall +#:enddef + +#! Generates an array shape specifier from an N-D array size +#! +#! Args: +#! name (str): Name of the original variable +#! rank (int): Rank of the original variable +#! offset(int): optional offset of the dimension loop (default = 0) +#! +#! Returns: +#! Array rank suffix string enclosed in braces +#! +#! E.g., +#! shape_from_array_size('mat', 5)}$ +#! -> (size(mat,1),size(mat,2),size(mat,3),size(mat,4),size(mat,5)) +#! shape_from_array_size('mat', 5, 2)}$ +#! -> (size(mat,3),size(mat,4),size(mat,5),size(mat,6),size(mat,7)) +#! +#:def shape_from_array_size(name, rank, offset=0) + #:assert rank > 0 + #:call join_lines(joinstr=", ", prefix="(", suffix=")") + #:for i in range(1, rank + 1) + size(${name}$,${i+offset}$) + #:endfor + #:endcall +#:enddef + +#! Generates an array shape specifier from an N-D array of sizes +#! +#! Args: +#! name (str): Name of the original variable +#! rank (int): Rank of the original variable +#! offset(int): optional offset of the dimension loop (default = 0) +#! +#! Returns: +#! Array rank suffix string enclosed in braces +#! +#! E.g., +#! shape_from_array_data('mat', 5)}$ +#! -> (1:mat(1),1:mat(2),1:mat(3),1:mat(4),1:mat(5)) +#! shape_from_array_data('mat', 5, 2)}$ +#! -> (1:mat(3),1:mat(4),1:mat(5),1:mat(6),1:mat(7)) +#! +#:def shape_from_array_data(name, rank, offset=0) + #:assert rank > 0 + #:call join_lines(joinstr=", ", prefix="(", suffix=")") + #:for i in range(1, rank + 1) + 1:${name}$(${i+offset}$) + #:endfor + #:endcall +#:enddef + +#! +#! Start a sequence of loop with indexed variables over an N-D array +#! +#! Args: +#! varname (str): Name of the variable to be used as prefix +#! matname (str): Name of the variable to be used as array +#! n (int): Number of nested loops to be created (1=innermost; n=outermost) +#! dim_offset (int): Optional dimension offset (1st loop is over dimension 1+dim_offset) +#! intent (str): Optional indentation. Default: 8 spaces +#! +#! +#:def loop_variables_start(varname, matname, n, dim_offset=0, indent=" "*8) + #:assert n > 0 + #:for i in range(1, n + 1) +${indent}$do ${varname}$${n+1+dim_offset-i}$ = lbound(${matname}$, ${n+1+dim_offset-i}$), ubound(${matname}$, ${n+1+dim_offset-i}$) + #:endfor +#:enddef + +#:def loop_variables_end(n, indent=" "*8) +#:assert n > 0 + #:call join_lines(joinstr="; ",prefix=indent) + #:for i in range(1, n + 1) + enddo + #:endfor + #:endcall +#:enddef + #:endmute diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8090b5775..8b8efaccb 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -32,6 +32,7 @@ set(fppFiles stdlib_linalg_determinant.fypp stdlib_linalg_qr.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_io.fypp b/src/stdlib_io.fypp index 7aceae2e2..ac207a944 100644 --- a/src/stdlib_io.fypp +++ b/src/stdlib_io.fypp @@ -9,7 +9,6 @@ module stdlib_io use, intrinsic :: iso_fortran_env, only : input_unit use stdlib_kinds, only: sp, dp, xdp, qp, & int8, int16, int32, int64 - use stdlib_error, only: error_stop use stdlib_optval, only: optval use stdlib_ascii, only: is_blank use stdlib_string_type, only : string_type @@ -120,7 +119,8 @@ contains !! ... !! integer :: s - integer :: nrow, ncol, i, skiprows_, max_rows_ + integer :: nrow, ncol, i, ios, skiprows_, max_rows_ + character(len=128) :: iomsg, msgout skiprows_ = max(optval(skiprows, 0), 0) max_rows_ = optval(max_rows, -1) @@ -142,56 +142,51 @@ contains allocate(d(max_rows_, ncol)) do i = 1, skiprows_ - read(s, *) + read(s, *, iostat=ios, iomsg=iomsg) + + if (ios/=0) then + write(msgout,1) trim(iomsg),i,trim(filename) + error stop trim(msgout) + end if + end do - - #:if 'real' in t1 + ! Default to format used for savetxt if fmt not specified. - fmt_ = optval(fmt, "(*"//FMT_REAL_${k1}$(1:len(FMT_REAL_${k1}$)-1)//",1x))") - - if ( fmt_ == '*' ) then - ! Use list directed read if user has specified fmt='*' - do i = 1, max_rows_ - read (s,*) d(i, :) - enddo - else - ! Otherwise pass default or user specified fmt string. - do i = 1, max_rows_ - read (s,fmt_) d(i, :) - enddo - endif + #:if 'real' in t1 + fmt_ = optval(fmt, "(*"//FMT_REAL_${k1}$(1:len(FMT_REAL_${k1}$)-1)//",:,1x))") #:elif 'complex' in t1 - ! Default to format used for savetxt if fmt not specified. - fmt_ = optval(fmt, "(*"//FMT_COMPLEX_${k1}$(1:len(FMT_COMPLEX_${k1}$)-1)//",1x))") - if ( fmt_ == '*' ) then - ! Use list directed read if user has specified fmt='*' - do i = 1, max_rows_ - read (s,*) d(i, :) - enddo - else - ! Otherwise pass default or user specified fmt string. - do i = 1, max_rows_ - read (s,fmt_) d(i, :) - enddo - endif + fmt_ = optval(fmt, "(*"//FMT_COMPLEX_${k1}$(1:len(FMT_COMPLEX_${k1}$)-1)//",:,1x))") #:else - ! Default to list directed for integer fmt_ = optval(fmt, "*") - ! Use list directed read if user has specified fmt='*' + #:endif + if ( fmt_ == '*' ) then + ! Use list directed read if user has specified fmt='*' do i = 1, max_rows_ - read (s,*) d(i, :) + read (s,*,iostat=ios,iomsg=iomsg) d(i, :) + + if (ios/=0) then + write(msgout,1) trim(iomsg),i,trim(filename) + error stop trim(msgout) + end if + enddo else - ! Otherwise pass default user specified fmt string. + ! Otherwise pass default or user specified fmt string. do i = 1, max_rows_ - read (s,fmt_) d(i, :) + read (s,fmt_,iostat=ios,iomsg=iomsg) d(i, :) + + if (ios/=0) then + write(msgout,1) trim(iomsg),i,trim(filename) + error stop trim(msgout) + end if + enddo endif - #:endif - close(s) + + 1 format('loadtxt: error <',a,'> reading line ',i0,' of ',a,'.') end subroutine loadtxt_${t1[0]}$${k1}$ #:endfor @@ -218,20 +213,31 @@ contains !!``` !! - integer :: s, i + integer :: s, i, ios + character(len=128) :: iomsg, msgout s = open(filename, "w") do i = 1, size(d, 1) #:if 'real' in t1 - write(s, "(*"//FMT_REAL_${k1}$(1:len(FMT_REAL_${k1}$)-1)//",1x))") d(i, :) + write(s, "(*"//FMT_REAL_${k1}$(1:len(FMT_REAL_${k1}$)-1)//",:,1x))", & #:elif 'complex' in t1 - write(s, "(*"//FMT_COMPLEX_${k1}$(1:len(FMT_COMPLEX_${k1}$)-1)//",1x))") d(i, :) + write(s, "(*"//FMT_COMPLEX_${k1}$(1:len(FMT_COMPLEX_${k1}$)-1)//",:,1x))", & #:elif 'integer' in t1 - write(s, "(*"//FMT_INT(1:len(FMT_INT)-1)//",1x))") d(i, :) + write(s, "(*"//FMT_INT(1:len(FMT_INT)-1)//",:,1x))", & #:else - write(s, *) d(i, :) + write(s, *, & #:endif + iostat=ios,iomsg=iomsg) d(i, :) + + if (ios/=0) then + write(msgout,1) trim(iomsg),i,trim(filename) + error stop trim(msgout) + end if + end do close(s) + + 1 format('savetxt: error <',a,'> writing line ',i0,' of ',a,'.') + end subroutine savetxt_${t1[0]}$${k1}$ #:endfor @@ -360,7 +366,7 @@ contains position_='asis' status_='new' case default - call error_stop("Unsupported mode: "//mode_(1:2)) + error stop "Unsupported mode: "//mode_(1:2) end select select case (mode_(3:3)) @@ -369,7 +375,7 @@ contains case('b') form_='unformatted' case default - call error_stop("Unsupported mode: "//mode_(3:3)) + error stop "Unsupported mode: "//mode_(3:3) end select access_ = 'stream' @@ -415,9 +421,9 @@ contains else if (a(i:i) == ' ') then cycle else if(any(.not.lfirst)) then - call error_stop("Wrong mode: "//trim(a)) + error stop "Wrong mode: "//trim(a) else - call error_stop("Wrong character: "//a(i:i)) + error stop "Wrong character: "//a(i:i) endif end do @@ -466,7 +472,7 @@ contains if (present(iostat)) then iostat = stat else if (stat /= 0) then - call error_stop(trim(msg)) + error stop trim(msg) end if end subroutine getline_char diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index cd4de131c..ad36ad677 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 @@ -1136,6 +1138,188 @@ module stdlib_linalg #:endfor end interface svdvals + + #! Allow for integer or character norm input: i.e., norm(a,2) or norm(a, '2') + #:set NORM_INPUT_TYPE = ["character(len=*)","integer(ilp)"] + #:set NORM_INPUT_SHORT = ["char","int"] + #:set NORM_INPUT_OPTIONS = list(zip(NORM_INPUT_TYPE,NORM_INPUT_SHORT)) + ! Vector norms: function interface + interface norm + !! version: experimental + !! + !! Computes the vector norm of a generic-rank array \( A \). + !! ([Specification](../page/specs/stdlib_linalg.html#norm-computes-the-vector-norm-of-a-generic-rank-array)) + !! + !!### Summary + !! Return one of several scalar norm metrics of a `real` or `complex` input array \( A \), + !! that can have any rank. For generic rank-n arrays, the scalar norm over the whole + !! array is returned by default. If `n>=2` and the optional input dimension `dim` is specified, + !! a rank `n-1` array is returned with dimension `dim` collapsed, containing all 1D array norms + !! evaluated along dimension `dim` only. + !! + !! + !!### Description + !! + !! This interface provides methods for computing the vector norm(s) of an array. + !! Supported data types include `real` and `complex`. + !! Input arrays may have generic rank from 1 to ${MAXRANK}$. + !! + !! Norm type input is mandatory, and it is provided via the `order` argument. + !! This can be provided as either an `integer` value or a `character` string. + !! Allowed metrics are: + !! - 1-norm \( \sum_i{ \left|a_i\right| } \): `order` = 1 or '1' + !! - Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \): `order` = 2 or '2' + !! - p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \): `integer` `order`, order>=3 + !! - Infinity norm \( \max_i{ \left|a_i\right| } \): order = huge(0) or 'inf' + !! - Minus-infinity norm \( \min_i{ \left|a_i\right| } \): order = -huge(0) or '-inf' + !! + !!### Example + !! + !!```fortran + !! + !! real(sp) :: a(3,3), na, rown(3) + !! a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3]) + !! + !! ! L2 norm: whole matrix + !! na = norm(a, 2) + !! + !! ! Infinity norm of each row + !! rown = norm(a, 'inf', dim=2) + !! + !!``` + !! + #:for rk,rt,ri in RC_KINDS_TYPES + #:for it,ii in NORM_INPUT_OPTIONS + !> Scalar norms: ${rt}$ + #:for rank in range(1, MAXRANK + 1) + pure module function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$(a, order) result(nrm) + !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ + ${rt}$, intent(in) :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Norm of the matrix. + real(${rk}$) :: nrm + end function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$ + module function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$(a, order, err) result(nrm) + !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ + ${rt}$, intent(in) :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Output state return flag. + type(linalg_state_type), intent(out) :: err + !> Norm of the matrix. + real(${rk}$) :: nrm + end function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$ + #:endfor + !> Array norms: ${rt}$ + #:for rank in range(2, MAXRANK + 1) + pure module function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, order, dim) result(nrm) + !> Input matrix a[..] + ${rt}$, intent(in), target :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Dimension the norm is computed along + integer(ilp), intent(in) :: dim + !> Norm of the matrix. (Same shape as `a`, with `dim` dropped). + real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$ + end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ + module function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$(a, order, dim, err) result(nrm) + !> Input matrix a[..] + ${rt}$, intent(in), target :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Dimension the norm is computed along + integer(ilp), intent(in) :: dim + !> Output state return flag. + type(linalg_state_type), intent(out) :: err + !> Norm of the matrix. (Same shape as `a`, with `dim` dropped). + real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$ + end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$ + #:endfor + #:endfor + #:endfor + end interface norm + + !> Vector norm: subroutine interface + interface get_norm + !! version: experimental + !! + !! Computes the vector norm of a generic-rank array \( A \). + !! ([Specification](../page/specs/stdlib_linalg.html#get-norm-computes-the-vector-norm-of-a-generic-rank-array)) + !! + !!### Summary + !! Subroutine interface that returns one of several scalar norm metrics of a `real` or `complex` + !! input array \( A \), that can have any rank. For generic rank-n arrays, the scalar norm over + !! the whole array is returned by default. If `n>=2` and the optional input dimension `dim` is + !! specified, a rank `n-1` array is returned with dimension `dim` collapsed, containing all 1D + !! array norms evaluated along dimension `dim` only. + !! + !! + !!### Description + !! + !! This `pure subroutine `interface provides methods for computing the vector norm(s) of an array. + !! Supported data types include `real` and `complex`. + !! Input arrays may have generic rank from 1 to ${MAXRANK}$. + !! + !! Norm type input is mandatory, and it is provided via the `order` argument. + !! This can be provided as either an `integer` value or a `character` string. + !! Allowed metrics are: + !! - 1-norm \( \sum_i{ \left|a_i\right| } \): `order` = 1 or '1' + !! - Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \): `order` = 2 or '2' + !! - p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \): `integer` `order`, order>=3 + !! - Infinity norm \( \max_i{ \left|a_i\right| } \): order = huge(0) or 'inf' + !! - Minus-infinity norm \( \min_i{ \left|a_i\right| } \): order = -huge(0) or '-inf' + !! + !!### Example + !! + !!```fortran + !! + !! real(sp) :: a(3,3), na, rown(3) + !! type(linalg_state_type) :: err + !! a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3]) + !! + !! ! L2 norm: whole matrix + !! call get_norm(a, na, 2) + !! + !! ! Infinity norms of each row, with error control + !! call get_norm(a, rown, 'inf', dim=2, err=err) + !! + !!``` + !! + #:for rk,rt,ri in RC_KINDS_TYPES + #:for it,ii in NORM_INPUT_OPTIONS + !> Scalar norms: ${rt}$ + #:for rank in range(1, MAXRANK + 1) + pure module subroutine norm_${rank}$D_${ii}$_${ri}$(a, nrm, order, err) + !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ + ${rt}$, intent(in), target :: a${ranksuffix(rank)}$ + !> Norm of the matrix. + real(${rk}$), intent(out) :: nrm + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out), optional :: err + end subroutine norm_${rank}$D_${ii}$_${ri}$ + #:endfor + !> Array norms: ${rt}$ + #:for rank in range(2, MAXRANK + 1) + pure module subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err) + !> Input matrix a[..] + ${rt}$, intent(in) :: a${ranksuffix(rank)}$ + !> Dimension the norm is computed along + integer(ilp), intent(in) :: dim + !> Norm of the matrix. (Same shape as `a`, with `dim` dropped). + real(${rk}$), intent(out) :: nrm${reduced_shape('a', rank, 'dim')}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out), optional :: err + end subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ + #:endfor + #:endfor + #:endfor + end interface get_norm + contains diff --git a/src/stdlib_linalg_blas.fypp b/src/stdlib_linalg_blas.fypp index 504a15108..2c0618e8d 100644 --- a/src/stdlib_linalg_blas.fypp +++ b/src/stdlib_linalg_blas.fypp @@ -9,6 +9,60 @@ module stdlib_linalg_blas implicit none(type,external) public + interface asum + !! ASUM takes the sum of the absolute values. +#ifdef STDLIB_EXTERNAL_BLAS + pure real(dp) function dasum( n, x, incx ) + import sp,dp,qp,ilp,lk + implicit none(type,external) + integer(ilp), intent(in) :: incx,n + real(dp), intent(in) :: x(*) + end function dasum +#else + module procedure stdlib_dasum +#endif +#ifdef STDLIB_EXTERNAL_BLAS + pure real(dp) function dzasum( n, x, incx ) + import sp,dp,qp,ilp,lk + implicit none(type,external) + integer(ilp), intent(in) :: incx,n + complex(dp), intent(in) :: x(*) + end function dzasum +#else + module procedure stdlib_dzasum +#endif +#:for rk,rt,ri in REAL_KINDS_TYPES +#:if not rk in ["sp","dp"] + module procedure stdlib_${ri}$asum +#:endif +#:endfor +#:for rk,rt,ri in CMPLX_KINDS_TYPES +#:if not rk in ["sp","dp"] + module procedure stdlib_${c2ri(ri)}$zasum +#:endif +#:endfor +#ifdef STDLIB_EXTERNAL_BLAS + pure real(sp) function sasum( n, x, incx ) + import sp,dp,qp,ilp,lk + implicit none(type,external) + integer(ilp), intent(in) :: incx,n + real(sp), intent(in) :: x(*) + end function sasum +#else + module procedure stdlib_sasum +#endif +#ifdef STDLIB_EXTERNAL_BLAS + pure real(sp) function scasum( n, x, incx ) + import sp,dp,qp,ilp,lk + implicit none(type,external) + integer(ilp), intent(in) :: incx,n + complex(sp), intent(in) :: x(*) + end function scasum +#else + module procedure stdlib_scasum +#endif + end interface asum + interface axpy !! AXPY constant times a vector plus a vector. #ifdef STDLIB_EXTERNAL_BLAS @@ -974,12 +1028,26 @@ module stdlib_linalg_blas #else module procedure stdlib_dnrm2 #endif +#ifdef STDLIB_EXTERNAL_BLAS + pure real(dp) function dznrm2( n, x, incx ) + import sp,dp,qp,ilp,lk + implicit none(type,external) + integer(ilp), intent(in) :: incx,n + complex(dp), intent(in) :: x(*) + end function dznrm2 +#else + module procedure stdlib_dznrm2 +#endif #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] module procedure stdlib_${ri}$nrm2 - #:endif #:endfor +#:for rk,rt,ri in CMPLX_KINDS_TYPES +#:if not rk in ["sp","dp"] + module procedure stdlib_${c2ri(ri)}$znrm2 +#:endif +#:endfor #ifdef STDLIB_EXTERNAL_BLAS pure real(sp) function snrm2( n, x, incx ) import sp,dp,qp,ilp,lk @@ -989,6 +1057,16 @@ module stdlib_linalg_blas end function snrm2 #else module procedure stdlib_snrm2 +#endif +#ifdef STDLIB_EXTERNAL_BLAS + pure real(sp) function scnrm2( n, x, incx ) + import sp,dp,qp,ilp,lk + implicit none(type,external) + integer(ilp), intent(in) :: incx,n + complex(sp), intent(in) :: x(*) + end function scnrm2 +#else + module procedure stdlib_scnrm2 #endif end interface nrm2 diff --git a/src/stdlib_linalg_inverse.fypp b/src/stdlib_linalg_inverse.fypp index 38ad3dda7..20567b2fa 100644 --- a/src/stdlib_linalg_inverse.fypp +++ b/src/stdlib_linalg_inverse.fypp @@ -121,7 +121,7 @@ submodule (stdlib_linalg) stdlib_linalg_inverse inva = a !> Compute matrix inverse - call stdlib_linalg_invert_inplace_${ri}$(inva,err=err0) + call stdlib_linalg_invert_inplace_${ri}$(inva,pivot=pivot,err=err0) end if @@ -164,7 +164,7 @@ submodule (stdlib_linalg) stdlib_linalg_inverse if (allocated(inva)) deallocate(inva) allocate(inva(size(a,1,kind=ilp),size(a,2,kind=ilp))) - #:if rt.startswith('complex') + #:if rt.startswith('real') inva = ieee_value(1.0_${rk}$,ieee_quiet_nan) #:else inva = cmplx(ieee_value(1.0_${rk}$,ieee_quiet_nan), & diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp new file mode 100644 index 000000000..8e3fba4a4 --- /dev/null +++ b/src/stdlib_linalg_norms.fypp @@ -0,0 +1,366 @@ +#:include "common.fypp" +#:set ALL_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + +#! Allow for integer or character norm input: i.e., norm(a,2) or norm(a, '2') +#:set INPUT_TYPE = ["character(len=*)","integer(ilp)"] +#:set INPUT_SHORT = ["char","int"] +#:set INPUT_OPTIONS = list(zip(INPUT_TYPE,INPUT_SHORT)) +! Vector norms +submodule(stdlib_linalg) stdlib_linalg_norms + use stdlib_linalg_constants + use stdlib_linalg_blas + use stdlib_linalg_lapack + use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & + LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR + use iso_c_binding, only: c_intptr_t,c_char,c_loc + implicit none(type,external) + + character(*), parameter :: this = 'norm' + + !> List of internal norm flags + integer(ilp), parameter :: NORM_ONE = 1_ilp + integer(ilp), parameter :: NORM_TWO = 2_ilp + integer(ilp), parameter :: NORM_POW_FIRST = 3_ilp + integer(ilp), parameter :: NORM_INF = +huge(0_ilp) ! infinity norm + integer(ilp), parameter :: NORM_POW_LAST = NORM_INF - 1_ilp + integer(ilp), parameter :: NORM_MINUSINF = -huge(0_ilp) + + interface parse_norm_type + module procedure parse_norm_type_integer + module procedure parse_norm_type_character + end interface parse_norm_type + + + interface stride_1d + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stride_1d_${ri}$ + #:endfor + end interface stride_1d + + contains + + !> Parse norm type from an integer user input + pure subroutine parse_norm_type_integer(order,norm_type,err) + !> User input value + integer(ilp), intent(in) :: order + !> Return value: norm type + integer(ilp), intent(out) :: norm_type + !> State return flag + type(linalg_state_type), intent(out) :: err + + select case (order) + case (1_ilp) + norm_type = NORM_ONE + case (2_ilp) + norm_type = NORM_TWO + case (3_ilp:NORM_POW_LAST) + norm_type = order + case (NORM_INF:) + norm_type = NORM_INF + case (:NORM_MINUSINF) + norm_type = NORM_MINUSINF + + case default + norm_type = NORM_ONE + err = linalg_state_type(this,LINALG_ERROR,'Input norm type ',order,' is not recognized.') + end select + + end subroutine parse_norm_type_integer + + pure subroutine parse_norm_type_character(order,norm_type,err) + !> User input value + character(len=*), intent(in) :: order + !> Return value: norm type + integer(ilp), intent(out) :: norm_type + !> State return flag + type(linalg_state_type), intent(out) :: err + + integer(ilp) :: int_order,read_err + + select case (order) + case ('inf','Inf','INF') + norm_type = NORM_INF + case ('-inf','-Inf','-INF') + norm_type = NORM_MINUSINF + case ('Euclidean','euclidean','EUCLIDEAN') + norm_type = NORM_TWO + case default + + ! Check if this input can be read as an integer + read(order,*,iostat=read_err) int_order + if (read_err/=0) then + ! Cannot read as an integer + norm_type = NORM_ONE + err = linalg_state_type(this,LINALG_ERROR,'Input norm type ',order,' is not recognized.') + else + call parse_norm_type_integer(int_order,norm_type,err) + endif + + end select + + end subroutine parse_norm_type_character + + #:for rk,rt,ri in ALL_KINDS_TYPES + + ! Compute stride of a 1d array + pure integer(ilp) function stride_1d_${ri}$(a) result(stride) + !> Input 1-d array + ${rt}$, intent(in), target :: a(:) + + integer(c_intptr_t) :: a1,a2 + + if (size(a,kind=ilp)<=1_ilp) then + stride = 1_ilp + else + a1 = transfer(c_loc(a(1)),a1) + a2 = transfer(c_loc(a(2)),a2) + stride = bit_size(0_c_char)*int(a2-a1, ilp)/storage_size(a, kind=ilp) + endif + + end function stride_1d_${ri}$ + + ! Private internal 1D implementation. This has to be used only internally, + ! when all inputs are already checked for correctness. + pure subroutine internal_norm_1D_${ri}$(sze, a, nrm, norm_request) + !> Input matrix length + integer(ilp), intent(in) :: sze + !> Input contiguous 1-d matrix a(*) + ${rt}$, intent(in) :: a(sze) + !> Norm of the matrix. + real(${rk}$), intent(out) :: nrm + !> Internal matrix request + integer(ilp), intent(in) :: norm_request + + integer(ilp) :: i + real(${rk}$) :: rorder + intrinsic :: abs, sum, sqrt, maxval, minval, conjg + + ! Initialize norm to zero + nrm = 0.0_${rk}$ + + select case(norm_request) + case(NORM_ONE) + nrm = asum(sze,a,incx=1_ilp) + case(NORM_TWO) + nrm = nrm2(sze,a,incx=1_ilp) + case(NORM_INF) + #:if rt.startswith('complex') + ! The default BLAS stdlib_i${ri}$amax uses |Re(.)|+|Im(.)|. Do not use it. + i = stdlib_i${ri}$max1(sze,a,incx=1_ilp) + #:else + i = stdlib_i${ri}$amax(sze,a,incx=1_ilp) + #:endif + nrm = abs(a(i)) + case(NORM_MINUSINF) + nrm = minval( abs(a) ) + case (NORM_POW_FIRST:NORM_POW_LAST) + rorder = 1.0_${rk}$ / norm_request + nrm = sum( abs(a) ** norm_request ) ** rorder + end select + + end subroutine internal_norm_1D_${ri}$ + + #:for it,ii in INPUT_OPTIONS + + !============================================== + ! Norms : any rank to scalar + !============================================== + + #:for rank in range(1, MAXRANK + 1) + + ! Pure function interface, with order specification. On error, the code will stop + pure module function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$(a, order) result(nrm) + !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ + ${rt}$, intent(in) :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Norm of the matrix. + real(${rk}$) :: nrm + + call norm_${rank}$D_${ii}$_${ri}$(a, nrm=nrm, order=order) + + end function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$ + + ! Function interface with output error + module function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$(a, order, err) result(nrm) + !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ + ${rt}$, intent(in) :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Output state return flag. + type(linalg_state_type), intent(out) :: err + !> Norm of the matrix. + real(${rk}$) :: nrm + + call norm_${rank}$D_${ii}$_${ri}$(a, nrm=nrm, order=order, err=err) + + end function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$ + + ! Internal implementation: ${rank}$-d + pure module subroutine norm_${rank}$D_${ii}$_${ri}$(a, nrm, order, err) + !> Input ${rank}$-d matrix a${ranksuffix(rank)}$ + ${rt}$, intent(in), target :: a${ranksuffix(rank)}$ + !> Norm of the matrix. + real(${rk}$), intent(out) :: nrm + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out), optional :: err + + type(linalg_state_type) :: err_ + integer(ilp) :: sze,norm_request + real(${rk}$) :: rorder + intrinsic :: abs, sum, sqrt, maxval, minval, conjg + + sze = size(a,kind=ilp) + + ! Initialize norm to zero + nrm = 0.0_${rk}$ + + ! Check matrix size + if (sze<=0) then + err_ = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix shape: a=',shape(a,kind=ilp)) + call linalg_error_handling(err_,err) + return + end if + + ! Check norm request + call parse_norm_type(order,norm_request,err_) + if (err_%error()) then + call linalg_error_handling(err_,err) + return + endif + + ! Get norm + call internal_norm_1D_${ri}$(sze, a, nrm, norm_request) + call linalg_error_handling(err_,err) + + end subroutine norm_${rank}$D_${ii}$_${ri}$ + + #:endfor + + !==================================================================== + ! Norms : any rank to rank-1, with DIM specifier and ${ii}$ input + !==================================================================== + + #:for rank in range(2, MAXRANK + 1) + + ! Pure function interface with DIM specifier. On error, the code will stop + pure module function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, order, dim) result(nrm) + !> Input matrix a[..] + ${rt}$, intent(in), target :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Dimension to collapse by computing the norm w.r.t other dimensions + integer(ilp), intent(in) :: dim + !> Norm of the matrix. + real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$ + + call norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim) + + end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ + + ! Function interface with DIM specifier and output error state. + module function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$(a, order, dim, err) result(nrm) + !> Input matrix a[..] + ${rt}$, intent(in), target :: a${ranksuffix(rank)}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> Dimension to collapse by computing the norm w.r.t other dimensions + integer(ilp), intent(in) :: dim + !> Output state return flag. + type(linalg_state_type), intent(out) :: err + !> Norm of the matrix. + real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$ + + call norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err) + + end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$ + + ! Internal implementation + pure module subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err) + !> Input matrix a[..] + ${rt}$, intent(in) :: a${ranksuffix(rank)}$ + !> Dimension to collapse by computing the norm w.r.t other dimensions + ! (dim must be defined before it is used for `nrm`) + integer(ilp), intent(in) :: dim + !> Norm of the matrix. + real(${rk}$), intent(out) :: nrm${reduced_shape('a', rank, 'dim')}$ + !> Order of the matrix norm being computed. + ${it}$, intent(in) :: order + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out), optional :: err + + type(linalg_state_type) :: err_ + integer(ilp) :: sze,lda,norm_request,${loop_variables('j',rank-1,1)}$ + logical :: contiguous_data + real(${rk}$) :: rorder + integer(ilp), dimension(${rank}$) :: spe,spack,perm,iperm + integer(ilp), dimension(${rank}$), parameter :: dim_range = [(lda,lda=1_ilp,${rank}$_ilp)] + ${rt}$, allocatable :: apack${ranksuffix(rank)}$ + intrinsic :: abs, sum, sqrt, norm2, maxval, minval, conjg + + ! Input matrix properties + sze = size (a,kind=ilp) + spe = shape(a,kind=ilp) + + ! Initialize norm to zero + nrm = 0.0_${rk}$ + + if (sze<=0) then + err_ = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix shape: a=',shape(a,kind=ilp)) + call linalg_error_handling(err_,err) + return + end if + + ! Check dimension choice + if (dim<1 .or. dim>${rank}$) then + err_ = linalg_state_type(this,LINALG_VALUE_ERROR,'dimension ',dim, & + 'is out of rank for shape(a)=',shape(a,kind=ilp)) + call linalg_error_handling(err_,err) + return + end if + + ! Check norm request + call parse_norm_type(order,norm_request,err_) + if (err_%error()) then + call linalg_error_handling(err_,err) + return + endif + + ! The norm's leading dimension + lda = spe(dim) + + ! Check if input column data is contiguous + contiguous_data = dim==1 + + ! Get packed data with the norm dimension as the first dimension + if (.not.contiguous_data) then + + ! Permute array to map dim to 1 + perm = [dim,pack(dim_range,dim_range/=dim)] + iperm(perm) = dim_range + spack = spe(perm) + apack = reshape(a, shape=spack, order=iperm) + +${loop_variables_start('j', 'apack', rank-1, 1," "*12)}$ + call internal_norm_1D_${ri}$(lda, apack(:, ${loop_variables('j',rank-1,1)}$), & + nrm(${loop_variables('j',rank-1,1)}$), norm_request) +${loop_variables_end(rank-1," "*12)}$ + + else + +${loop_variables_start('j', 'a', rank-1, 1," "*12)}$ + call internal_norm_1D_${ri}$(lda, a(:, ${loop_variables('j',rank-1,1)}$), & + nrm(${loop_variables('j',rank-1,1)}$), norm_request) +${loop_variables_end(rank-1," "*12)}$ + + endif + + end subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$ + + #:endfor + #:endfor + #:endfor + +end submodule stdlib_linalg_norms diff --git a/src/stdlib_linalg_svd.fypp b/src/stdlib_linalg_svd.fypp index 8910bd3ce..6fa991d49 100644 --- a/src/stdlib_linalg_svd.fypp +++ b/src/stdlib_linalg_svd.fypp @@ -269,7 +269,11 @@ submodule(stdlib_linalg) stdlib_linalg_svd if (info==0) then !> Prepare working storage - lwork = nint(real(work_dummy(1),kind=${rk}$), kind=ilp) + ! Check if the returned working storage space is smaller than the largest value + ! allowed by lwork + lwork = merge(nint(real(work_dummy(1),kind=${rk}$), kind=ilp) & + , huge(lwork) & + , real(work_dummy(1),kind=${rk}$) < real(huge(lwork),kind=${rk}$) ) allocate(work(lwork)) !> Compute SVD diff --git a/src/stdlib_specialfunctions_gamma.fypp b/src/stdlib_specialfunctions_gamma.fypp index e91f14954..8fff966c2 100644 --- a/src/stdlib_specialfunctions_gamma.fypp +++ b/src/stdlib_specialfunctions_gamma.fypp @@ -528,9 +528,9 @@ contains ! Log(n!) ! ${t1}$, intent(in) :: n - real :: res + real(dp) :: res ${t1}$, parameter :: zero = 0_${k1}$, one = 1_${k1}$, two = 2_${k1}$ - real, parameter :: zero_k2 = 0.0 + real(dp), parameter :: zero_dp = 0.0_dp if(n < zero) call error_stop("Error(l_factorial): Logarithm of" & //" factorial function argument must be non-negative") @@ -539,15 +539,15 @@ contains case (zero) - res = zero_k2 + res = zero_dp case (one) - res = zero_k2 + res = zero_dp case (two : ) - res = l_gamma(n + 1, 1.0D0) + res = l_gamma(n + 1, 1.0_dp) end select end function l_factorial_${t1[0]}$${k1}$ diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index adf8b5db2..f835872df 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_qr.fypp" "test_linalg_svd.fypp" @@ -20,6 +21,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_qr) diff --git a/test/linalg/test_linalg_norm.fypp b/test/linalg/test_linalg_norm.fypp new file mode 100644 index 000000000..0ec88c2ec --- /dev/null +++ b/test/linalg/test_linalg_norm.fypp @@ -0,0 +1,300 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + +#! Generate an array rank suffix with the same fixed size for all dimensions. +#! +#! Args: +#! rank (int): Rank of the variable +#! size (int): Size along each dimension +#! +#! Returns: +#! Array rank suffix string (e.g. (4,4,4) if rank = 3 and size = 4) +#! +#:def fixedranksuffix(rank,size) +#{if rank > 0}#(${str(size) + (","+str(size)) * (rank - 1)}$)#{endif}# +#:enddef +! Test vector norms +module test_linalg_norm + use testdrive, only: error_type, check, new_unittest, unittest_type + use stdlib_linalg_constants + use stdlib_linalg, only: norm, linalg_state_type + use stdlib_linalg_state, only: linalg_state_type + + implicit none (type,external) + + contains + + !> Vector norm tests + subroutine test_vector_norms(tests) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + allocate(tests(0)) + + #:for rk,rt,ri in RC_KINDS_TYPES + tests = [tests,new_unittest("strided_1d_norm_${ri}$",test_strided_1d_${ri}$)] + #:for rank in range(1, MAXRANK) + tests = [tests,new_unittest("norm_${ri}$_${rank}$d",test_norm_${ri}$_${rank}$d)] + #:endfor + #:for rank in range(2, MAXRANK) + #:if rt.startswith('real') + tests = [tests,new_unittest("norm2_${ri}$_${rank}$d",test_norm2_${ri}$_${rank}$d)] + #:endif + tests = [tests,new_unittest("maxabs_${ri}$_${rank}$d",test_maxabs_${ri}$_${rank}$d)] + tests = [tests,new_unittest("norm_dimmed_${ri}$_${rank}$d",test_norm_dimmed_${ri}$_${rank}$d)] + #:endfor + #:endfor + + end subroutine test_vector_norms + + #:for rk,rt,ri in RC_KINDS_TYPES + + !> Test strided norm + subroutine test_strided_1d_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp), parameter :: m = 8_ilp + integer(ilp), parameter :: n = m**2 + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + ${rt}$, target :: a(n) + ${rt}$, allocatable :: slice(:) + ${rt}$, pointer :: twod(:,:) + real(${rk}$) :: rea(n),ima(n) + + call random_number(rea) + #:if rt.startswith('real') + a = rea + #:else + call random_number(ima) + a = cmplx(rea,ima,kind=${rk}$) + #:endif + + ! Test sliced array results + slice = a(4:7:59) + call check(error,abs(norm(a(4:7:59),2)-norm(slice,2)) a + call check(error,abs(norm(twod,2)-norm(a,2)) Test several norms with different dimensions + subroutine test_norm_${ri}$_${rank}$d(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp) :: j,order + integer(ilp), parameter :: n = 2_ilp**${rank}$ + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$ + character(64) :: msg + + allocate(a(n), b${fixedranksuffix(rank,2)}$) + + ! Init as a range,but with small elements such that all power norms will + ! never overflow, even in single precision + a = [(0.01_${rk}$*(j-n/2_ilp), j=1_ilp,n)] + b = reshape(a, shape(b)) + + ! Test some norms + do order = 1, 10 + write(msg,"('reshaped order-',i0,' p-norm is the same')") order + call check(error,abs(norm(a,order)-norm(b,order)) Test Euclidean norm; compare with Fortran intrinsic norm2 for reals + #:if rt.startswith('real') + subroutine test_norm2_${ri}$_${rank}$d(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp) :: j,dim + integer(ilp), parameter :: ndim = ${rank}$ + integer(ilp), parameter :: n = 2_ilp**ndim + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$ + intrinsic :: norm2 + character(64) :: msg + + allocate(a(n), b${fixedranksuffix(rank,2)}$) + + ! Init as a range,but with small elements such that all power norms will + ! never overflow, even in single precision + a = [(0.01_${rk}$*(j-n/2_ilp), j=1_ilp,n)] + b = reshape(a, shape(b)) + + ! Test some norms + call check(error,abs(norm(a,2) - norm2(a)) Test Infinity norm; compare with Fortran intrinsic max(abs(a)) + subroutine test_maxabs_${ri}$_${rank}$d(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp) :: j,dim + integer(ilp), parameter :: ndim = ${rank}$ + integer(ilp), parameter :: n = 2_ilp**ndim + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$ + intrinsic :: maxval, abs + character(128) :: msg + + allocate(a(n), b${fixedranksuffix(rank,2)}$) + + ! Init as a range,but with small elements such that all power norms will + ! never overflow, even in single precision + a = [(0.01_${rk}$*(j-n/2_ilp), j=1_ilp,n)] + b = reshape(a, shape(b)) + + ! Test some norms + call check(error,abs(norm(a,'inf') - maxval(abs(a))) 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_norm + + diff --git a/test/specialfunctions/test_specialfunctions_gamma.fypp b/test/specialfunctions/test_specialfunctions_gamma.fypp index 62ee4c1f9..c29cf7a60 100644 --- a/test/specialfunctions/test_specialfunctions_gamma.fypp +++ b/test/specialfunctions/test_specialfunctions_gamma.fypp @@ -89,36 +89,33 @@ contains ${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, & 5_${k1}$, 100_${k1}$] - real(sp), parameter :: ans(n) = [0.0, 0.0, 0.693147180, 3.17805383, & - 4.78749174, 3.63739376e2] + real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, & + 4.78749174278204_dp, 3.637393755555e2_dp] #:elif k1 == "int16" ${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, & 7_${k1}$, 500_${k1}$] - real(sp), parameter :: ans(n) = [0.0, 0.0, 0.693147180, 3.17805383, & - 8.52516136, 2.61133046e3] - + real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, & + 8.52516136106541_dp, 2.611330458460e3_dp] #:elif k1 == "int32" ${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, & 12_${k1}$, 7000_${k1}$] - real(sp), parameter :: ans(n) = [0.0, 0.0, 0.693147180, 3.17805383, & - 1.99872145e1, 5.49810038e4] - + real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, & + 1.998721449566e1_dp, 5.498100377941e4_dp] #:elif k1 == "int64" ${t1}$, parameter :: x(n) = [0_${k1}$, 1_${k1}$, 2_${k1}$, 4_${k1}$, & 20_${k1}$, 90000_${k1}$] - real(sp), parameter :: ans(n) = [0.0, 0.0, 0.693147180, 3.17805383, & - 4.23356165e1, 9.36687468e5] - + real(dp), parameter :: ans(n) = [0.0_dp, 0.0_dp, 0.69314718055994_dp, 3.17805383034794_dp, & + 4.233561646075e1_dp, 9.366874681600e5_dp] #:endif do i = 1, n call check(error, log_factorial(x(i)), ans(i), "Integer kind " & - //"${k1}$ failed", thr = tol_sp, rel = .true.) + //"${k1}$ failed", thr = tol_dp, rel = .true.) end do end subroutine test_logfact_${t1[0]}$${k1}$