linalg: Matrix norms (#885)

Author: perazz

doc/specs/stdlib_linalg.md

@@ -1565,4 +1565,53 @@ If `err` is not present, exceptions trigger an `error stop`.
 {!example/linalg/example_norm.f90!}
 ```
 
+## `mnorm` - Computes the matrix norm of a generic-rank array.
+
+### Status
+
+Experimental
+
+### Description
+
+This function computes one of several matrix norms of `real` or `complex` array \( A \), depending on 
+the value of the `order` input argument. \( A \) must be an array of rank 2 or higher. For arrays of rank > 2,
+matrix norms are computed over specified dimensions.
+
+### Syntax
+
+`x = ` [[stdlib_linalg(module):mnorm(interface)]] `(a [, order, dim, err])`
+
+### Arguments
+
+`a`: Shall be a rank-n `real` or `complex` array containing the data, where n >= 2. It is an `intent(in)` argument.
+
+`order` (optional): Shall be an `integer` value or a `character` flag that specifies the norm type, as follows. It is an `intent(in)` argument. 
+
+| Integer input    | Character Input                 | Norm type                                                                   |
+|------------------|---------------------------------|-----------------------------------------------------------------------------|
+| `1`              | `'1'`                           | 1-norm (maximum column sum) \( \max_j \sum_i{ \left|a_{i,j}\right| } \)     |
+| `2`              | `'2'`                           | 2-norm (largest singular value)                                             |
+| (not prov.)      | `'Euclidean','Frobenius','Fro'` | Frobenius norm \( \sqrt{\sum_{i,j}{ \left|a_{i,j}\right|^2 }} \)            |
+| `huge(0)`        | `'inf', 'Inf', 'INF'`           | Infinity norm (maximum row sum) \( \max_i \sum_j{ \left|a_{i,j}\right| } \) |
+
+`dim` (optional): For arrays of rank > 2, shall be an integer array of size 2 specifying the dimensions over which to compute the matrix norm. Default value is `[1,2]`. It is an `intent(in)` argument.
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. 
+
+### Return value
+
+For rank-2 input arrays, the return value `x` is a scalar containing the matrix norm.
+For arrays of rank > 2, if the optional `dim` argument is present, `x` is a rank `n-2` array with the same shape as \( A \) except 
+for dimensions `dim(1)` and `dim(2)`, which are dropped. Each element of `x` contains the matrix norm of the corresponding submatrix of \( A \),
+evaluated over the specified dimensions only, with the given order.
+
+If an invalid norm type is provided, defaults to 1-norm and raises `LINALG_ERROR`.
+Raises `LINALG_VALUE_ERROR` if any of the arguments has an invalid size.
+If `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_mnorm.f90!}
+```
 

example/linalg/CMakeLists.txt

@@ -29,6 +29,7 @@ ADD_EXAMPLE(lapack_getrf)
 ADD_EXAMPLE(lstsq1)
 ADD_EXAMPLE(lstsq2)
 ADD_EXAMPLE(norm)
+ADD_EXAMPLE(mnorm)
 ADD_EXAMPLE(get_norm)
 ADD_EXAMPLE(solve1)
 ADD_EXAMPLE(solve2)

example/linalg/example_mnorm.f90

@@ -0,0 +1,26 @@
+program example_mnorm
+    use stdlib_linalg, only: mnorm
+    use stdlib_kinds, only: sp
+    implicit none
+    real(sp) :: a(3,3), na
+    real(sp) :: b(3,3,4), nb(4)  ! Array of 4 3x3 matrices
+    
+    ! Initialize example matrix
+    a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
+    
+    ! Compute Euclidean norm of single matrix
+    na = mnorm(a, 'Euclidean')
+    print *, "Euclidean norm of matrix a:", na
+    
+    ! Initialize array of matrices
+    b(:,:,1) = a
+    b(:,:,2) = 2*a
+    b(:,:,3) = 3*a
+    b(:,:,4) = 4*a
+    
+    ! Compute infinity norm of each 3x3 matrix in b
+    nb = mnorm(b, 'inf', dim=[1,2])
+
+    ! 18.0000000       36.0000000       54.0000000       72.0000000
+    print *, "Infinity norms of matrices in b:", nb
+end program example_mnorm

include/common.fypp

@@ -157,6 +157,22 @@ $:"s" if cmplx=="c" else "d" if cmplx=="z" else "x" if cmplx=="y" else "q" if cm
 #{if rank > 0}#(${"0" + ",0" * (rank - 1)}$)#{endif}#
 #:enddef
 
+#! Generates an array rank suffix with a fixed integer size for all dimensions.
+#!
+#! Args:
+#!     rank (int): Rank of the variable
+#!     size (int): Size along each dimension
+#!
+#! Returns:
+#!     Array rank suffix string
+#! E.g.,
+#!   fixedranksuffix(3,4)
+#!   -> (4,4,4) 
+#!
+#:def fixedranksuffix(rank,size)
+#{if rank > 0}#(${str(size) + (","+str(size)) * (rank - 1)}$)#{endif}#
+#:enddef
+
 #! Joins stripped lines with given character string
 #!
 #! Args:
@@ -227,7 +243,7 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$
 #!   Array rank suffix string enclosed in braces
 #!
 #! E.g.,
-#!   select_subarray(5 , [(4, 'i'), (5, 'j')])}$
+#!   select_subarray(5 , [(4, 'i'), (5, 'j')])
 #!   -> (:, :, :, i, j)
 #!
 #:def select_subarray(rank, selectors)
@@ -327,6 +343,34 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$
     #:endcall
 #:enddef
 
+#!
+#! Generates a list of loop variables from an array
+#!
+#! Args:
+#!   varname(str): Name of the array variable to be used as prefix
+#!   n      (int): Number of loop variables to be created
+#!   offset (int): Optional index offset
+#!
+#! Returns:
+#!   Variable definition string 
+#!
+#! E.g.,
+#!  loop_array_variables('j', 5)
+#!   -> "j(1), j(2), j(3), j(4), j(5)
+#!
+#!  loop_array_variables('j', 5, 2)
+#!   -> "j(3), j(4), j(5), j(6), j(7)
+#!
+#:def loop_array_variables(varname, n, offset=0)
+  #:assert n > 0
+    #:call join_lines(joinstr=", ")
+      #:for i in range(1, n + 1)
+        ${varname}$(${i+offset}$)
+      #:endfor
+    #:endcall
+#:enddef
+
+
 #! Generates an array shape specifier from an N-D array size
 #!
 #! Args:

src/stdlib_linalg.fypp

@@ -32,6 +32,7 @@ module stdlib_linalg
   public :: lstsq
   public :: lstsq_space
   public :: norm
+  public :: mnorm
   public :: get_norm
   public :: solve
   public :: solve_lu  
@@ -1320,6 +1321,89 @@ module stdlib_linalg
      #:endfor
   end interface get_norm
 
+  !> Matrix norms: function interface
+  interface mnorm
+     !! version: experimental 
+     !!
+     !! Computes the matrix norm of a generic-rank array \( A \). 
+     !! ([Specification](../page/specs/stdlib_linalg.html#mnorm-computes-the-matrix-norm-of-a-generic-rank-array))
+     !! 
+     !!### Summary 
+     !! Return one of several matrix norm metrics of a `real` or `complex` input array \( A \), 
+     !! that can have rank 2 or higher. For rank-2 arrays, the matrix norm is returned.
+     !! If rank>2 and the optional input dimensions `dim` are specified, 
+     !! a rank `n-2` array is returned with dimensions `dim(1),dim(2)` collapsed, containing all 
+     !! matrix norms evaluated over the specified dimensions only. `dim==[1,2]` are assumed as default
+     !! dimensions if not specified.
+     !! 
+     !!### Description
+     !! 
+     !! This interface provides methods for computing the matrix norm(s) of an array.  
+     !! Supported data types include `real` and `complex`. 
+     !! Input arrays must have rank >= 2.
+     !!
+     !! Norm type input is optional, and it is provided via the `order` argument. 
+     !! This can be provided as either an `integer` value or a `character` string. 
+     !! Allowed metrics are: 
+     !! - 1-norm: `order` = 1 or '1'    
+     !! - 2-norm: `order` = 2 or '2'
+     !! - Euclidean/Frobenius: `order` = 'Euclidean','Frobenius', or argument not specified
+     !! - Infinity norm: `order` = huge(0) or 'Inf'
+     !! 
+     !! If an invalid norm type is provided, the routine returns an error state.
+     !!
+     !!### Example
+     !!
+     !!```fortran
+     !!    real(sp) :: a(3,3), na
+     !!    real(sp) :: b(3,3,4), nb(4)  ! Array of 4 3x3 matrices
+     !!    a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
+     !!    
+     !!    ! Euclidean/Frobenius norm of single matrix
+     !!    na = mnorm(a)
+     !!    na = mnorm(a, 'Euclidean')
+     !!   
+     !!    ! 1-norm of each 3x3 matrix in b
+     !!    nb = mnorm(b, 1, dim=[1,2])
+     !!     
+     !!    ! Infinity-norm 
+     !!    na = mnorm(b, 'inf', dim=[3,2])
+     !!```     
+     !!
+      #:for rk,rt,ri in RC_KINDS_TYPES
+      #:for it,ii in NORM_INPUT_OPTIONS   
+      
+      !> Matrix norms: ${rt}$ rank-2 arrays
+      module function matrix_norm_${ii}$_${ri}$(a, order, err) result(nrm)
+        !> Input matrix a(m,n)
+        ${rt}$, intent(in), target :: a(:,:)
+        !> Norm of the matrix.        
+        real(${rk}$) :: nrm
+        !> Order of the matrix norm being computed.
+        ${it}$, #{if 'integer' in it}#optional, #{endif}#intent(in) :: order
+        !> [optional] state return flag. On error if not requested, the code will stop
+        type(linalg_state_type), intent(out), optional :: err      
+      end function matrix_norm_${ii}$_${ri}$
+      
+      !> Matrix norms: ${rt}$ higher rank arrays
+      #:for rank in range(3, MAXRANK + 1)
+      module function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$(a, order, dim, err) result(nrm)
+          !> Input matrix a(m,n)
+          ${rt}$, intent(in), contiguous, target :: a${ranksuffix(rank)}$
+          !> Norm of the matrix.        
+          real(${rk}$), allocatable :: nrm${ranksuffix(rank-2)}$
+          !> Order of the matrix norm being computed.
+          ${it}$, #{if 'integer' in it}#optional, #{endif}#intent(in) :: order
+          !> [optional] dimensions of the sub-matrices the norms should be evaluated at (default = [1,2])
+          integer(ilp), optional, intent(in) :: dim(2)
+          !> [optional] state return flag. On error if not requested, the code will stop
+          type(linalg_state_type), intent(out), optional :: err        
+      end function matrix_norm_${rank}$D_to_${rank-2}$D_${ii}$_${ri}$
+      #:endfor            
+      #:endfor
+      #:endfor            
+  end interface mnorm
+
 contains