linalg: vector norms

.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

.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/

doc/specs/stdlib_linalg.md

@@ -1382,3 +1382,110 @@ If `err` is not present, exceptions trigger an `error stop`.
 {!example/linalg/example_inverse_function.f90!}
 ```
 
+## `get_norm` - Computes the vector norm of a generic-rank array.
+
+### Status
+
+Experimental
+
+### Description
+
+This `pure subroutine` interface computes one of several vector norms of `real` or `complex` array \( A \), depending on 
+the value of the `order` input argument. \( A \) may be an array of any rank. 
+
+Result `nrm` returns a `real`, scalar norm value for the whole array; if `dim` is specified, `nrm` is a rank n-1 
+array with the same shape as \(A \) and dimension `dim` dropped, containing all norms evaluated along `dim`.
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):get_norm(interface)]] `(a, nrm, order, [, dim, err])`
+
+### Arguments
+
+`a`: Shall be a rank-n `real` or `complex` array containing the data. It is an `intent(in)` argument.
+
+`nrm`: if `dim` is absent, shall be a scalar with the norm evaluated over all the elements of the array. Otherwise, an array of rank `n-1`, and a shape similar
+to that of `a` with dimension `dim` dropped.
+
+`order`: Shall be an `integer` value or a `character` flag that specifies the norm type, as follows. It is an `intent(in)` argument. 
+
+| Integer input    | Character Input  | Norm type                                               |
+|------------------|------------------|---------------------------------------------------------|
+| `-huge(0)`       | `'-inf', '-Inf'` | Minimum absolute value \( \min_i{ \left|a_i\right| } \) |
+| `1`              | `'1'`            | 1-norm \( \sum_i{ \left|a_i\right| } \)                 |
+| `2`              | `'2'`            | Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \)             |
+| `>=3`            | `'3','4',...`    | p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \) |
+| `huge(0)`        | `'inf', 'Inf'`   | Maximum absolute value \( \max_i{ \left|a_i\right| } \) |
+
+`dim` (optional): Shall be a scalar `integer` value with a value in the range from `1` to `n`, where `n` is the rank of the array. It is an `intent(in)` argument.
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. If `err` is not present, the function is `pure`. 
+
+### Return value
+
+By default, the return value `nrm` is a scalar, and contains the norm as evaluated over all elements of the generic-rank array \( A \). 
+If the optional `dim` argument is present, `nrm` is a rank `n-1` array with the same shape as \( A \) except 
+for dimension `dim`, that is collapsed. Each element of `nrm` contains the 1D norm of the elements of \( A \), 
+evaluated along dimension `dim` only.
+
+Raises `LINALG_ERROR` if the requested norm type is invalid.
+Raises `LINALG_VALUE_ERROR` if any of the arguments has an invalid size.
+If `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_get_norm.f90!}
+```
+
+## `norm` - Computes the vector norm of a generic-rank array.
+
+### Status
+
+Experimental
+
+### Description
+
+This function computes one of several vector norms of `real` or `complex` array \( A \), depending on 
+the value of the `order` input argument. \( A \) may be an array of any rank. 
+
+### Syntax
+
+`x = ` [[stdlib_linalg(module):norm(interface)]] `(a, order, [, dim, err])`
+
+### Arguments
+
+`a`: Shall be a rank-n `real` or `complex` array containing the data. It is an `intent(in)` argument.
+
+`order`: Shall be an `integer` value or a `character` flag that specifies the norm type, as follows. It is an `intent(in)` argument. 
+
+| Integer input    | Character Input  | Norm type                                               |
+|------------------|------------------|---------------------------------------------------------|
+| `-huge(0)`       | `'-inf', '-Inf'` | Minimum absolute value \( \min_i{ \left|a_i\right| } \) |
+| `1`              | `'1'`            | 1-norm \( \sum_i{ \left|a_i\right| } \)                 |
+| `2`              | `'2'`            | Euclidean norm \( \sqrt{\sum_i{ a_i^2 }} \)             |
+| `>=3`            | `'3','4',...`    | p-norm \( \left( \sum_i{ \left|a_i\right|^p }\right) ^{1/p} \) |
+| `huge(0)`        | `'inf', 'Inf'`   | Maximum absolute value \( \max_i{ \left|a_i\right| } \) |
+
+`dim` (optional): Shall be a scalar `integer` value with a value in the range from `1` to `n`, where `n` is the rank of the array. It is an `intent(in)` argument.
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. If `err` is not present, the function is `pure`. 
+
+### Return value
+
+By default, the return value `x` is a scalar, and contains the norm as evaluated over all elements of the generic-rank array \( A \). 
+If the optional `dim` argument is present, `x` is a rank `n-1` array with the same shape as \( A \) except 
+for dimension `dim`, that is dropped. Each element of `x` contains the 1D norm of the elements of \( A \), 
+evaluated along dimension `dim` only.
+
+Raises `LINALG_ERROR` if the requested norm type is invalid.
+Raises `LINALG_VALUE_ERROR` if any of the arguments has an invalid size.
+If `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_norm.f90!}
+```
+
+

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)

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

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

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

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