sync task with NumPy: no order-2, default=Frobenius

perazz · perazz · commit 753e31b3385d · 2024-11-06T04:33:16.000-06:00
diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
@@ -1585,12 +1585,12 @@ matrix norms are computed over specified dimensions.
 
 `a`: Shall be a rank-n `real` or `complex` array containing the data, where n � 2. It is an `intent(in)` argument.
 
-`order`: Shall be an `integer` value or a `character` flag that specifies the norm type, as follows. It is an `intent(in)` argument. 
+`order` (optional): Shall be an `integer` value or a `character` flag that specifies the norm type, as follows. It is an `intent(in)` argument. 
 
 | Integer input    | Character Input            | Norm type                                               |
 |------------------|---------------------------|---------------------------------------------------------|
 | `1`              | `'1'`                     | 1-norm (maximum column sum) \( \max_j \sum_i{ \left|a_{i,j}\right| } \) |
-| `2`              | `'2','Euclidean'`         | 2-norm/Euclidean/Frobenius norm \( \sqrt{\sum_{i,j}{ \left|a_{i,j}\right|^2 }} \) |
+| (not prov.)      | `'Euclidean','Frobenius','Fro'` | Frobenius norm \( \sqrt{\sum_{i,j}{ \left|a_{i,j}\right|^2 }} \) |
 | `huge(0)`        | `'inf', 'Inf', 'INF'`     | Infinity norm (maximum row sum) \( \max_i \sum_j{ \left|a_{i,j}\right| } \) |
 
 `dim` (optional): For arrays of rank > 2, shall be an integer array of size 2 specifying the dimensions over which to compute the matrix norm. Default value is `[1,2]`. It is an `intent(in)` argument.
diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
@@ -1333,7 +1333,8 @@ module stdlib_linalg
      !! that can have rank 2 or higher. For rank-2 arrays, the matrix norm is returned.
      !! If rank>2 and the optional input dimensions `dim` are specified, 
      !! a rank `n-2` array is returned with dimensions `dim(1),dim(2)` collapsed, containing all 
-     !! matrix norms evaluated over the specified dimensions only.
+     !! matrix norms evaluated over the specified dimensions only. `dim==[1,2]` are assumed as default
+     !! dimensions if not specified.
      !! 
      !!### Description
      !! 
@@ -1345,12 +1346,10 @@ module stdlib_linalg
      !! This can be provided as either an `integer` value or a `character` string. 
      !! Allowed metrics are: 
      !! - 1-norm: `order` = 1 or '1'    
-     !! - 2-norm/Euclidean: `order` = 2 or '2' or 'Euclidean'
-     !! - p-norm: `order` >= 3 
+     !! - Euclidean/Frobenius: `order` = 'Euclidean','Frobenius', or argument not specified
      !! - Infinity norm: `order` = huge(0) or 'Inf'
-     !! - Minus-infinity norm: `order` = '-Inf'
      !! 
-     !! If an invalid norm type is provided, the routine defaults to 1-norm and returns an error state.
+     !! If an invalid norm type is provided, the routine returns an error state.
      !!
      !!### Example
      !!
@@ -1359,14 +1358,15 @@ module stdlib_linalg
      !!    real(sp) :: b(3,3,4), nb(4)  ! Array of 4 3x3 matrices
      !!    a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
      !!    
-     !!    ! 2-norm/Euclidean norm of single matrix
+     !!    ! Euclidean/Frobenius norm of single matrix
+     !!    na = mnorm(a)
      !!    na = mnorm(a, 'Euclidean')
      !!   
      !!    ! 1-norm of each 3x3 matrix in b
      !!    nb = mnorm(b, 1, dim=[1,2])
      !!     
-     !!    ! p-norm with p=3
-     !!    na = mnorm(a, 3)
+     !!    ! Infinity-norm 
+     !!    na = mnorm(b, 'inf', dim=[3,2])
      !!```     
      !!
       #:for rk,rt,ri in RC_KINDS_TYPES
@@ -1379,7 +1379,7 @@ module stdlib_linalg
         !> Norm of the matrix.        
         real(${rk}$) :: nrm
         !> Order of the matrix norm being computed.
-        ${it}$, intent(in) :: order
+        ${it}$, #{if 'integer' in it}#optional, #{endif}#intent(in) :: order
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), intent(out), optional :: err      
       end function matrix_norm_${ii}$_${ri}$
@@ -1392,7 +1392,7 @@ module stdlib_linalg
           !> Norm of the matrix.        
           real(${rk}$), allocatable :: nrm${ranksuffix(rank-2)}$
           !> Order of the matrix norm being computed.
-          ${it}$, intent(in) :: order
+          ${it}$, #{if 'integer' in it}#optional, #{endif}#intent(in) :: order
           !> [optional] dimensions of the sub-matrices the norms should be evaluated at (default = [1,2])
           integer(ilp), optional, intent(in) :: dim(2)
           !> [optional] state return flag. On error if not requested, the code will stop
diff --git a/src/stdlib_linalg_norms.fypp b/src/stdlib_linalg_norms.fypp
@@ -29,13 +29,18 @@ submodule(stdlib_linalg) stdlib_linalg_norms
      character, parameter :: LANGE_NORM_MAT = 'M' ! maxval(sum(abs(a)))   ! over whole matrix: unused
      character, parameter :: LANGE_NORM_ONE = '1' ! maxval(sum(abs(a),1)) ! over columns
      character, parameter :: LANGE_NORM_INF = 'I' ! maxval(sum(abs(a),2)) ! over rows
-     character, parameter :: LANGE_NORM_TWO = 'E' ! "Euclidean" or "Frobenius"     
+     character, parameter :: LANGE_NORM_FRO = 'E' ! sqrt(sum(a**2))       ! "Euclidean" or "Frobenius"     
      
      interface parse_norm_type
         module procedure parse_norm_type_integer
         module procedure parse_norm_type_character
      end interface parse_norm_type
      
+     interface lange_task_request
+        module procedure lange_task_request_integer
+        module procedure lange_task_request_character
+     end interface lange_task_request
+     
      
      interface stride_1d
         #:for rk,rt,ri in ALL_KINDS_TYPES
@@ -107,25 +112,61 @@ submodule(stdlib_linalg) stdlib_linalg_norms
      end subroutine parse_norm_type_character
 
      !> From a user norm request, generate a *LANGE task command
-     pure subroutine lange_task_request(norm_type,lange_task,err)
+     pure subroutine lange_task_request_integer(order,lange_task,err)
         !> Parsed matrix norm type
-        integer(ilp), intent(in) :: norm_type
+        integer(ilp), optional, intent(in) :: order
         !> LANGE task
         character, intent(out) :: lange_task
         !> Error flag
         type(linalg_state_type), intent(inout) :: err
         
-        select case (norm_type)
-           case (NORM_INF)
-              lange_task = LANGE_NORM_INF 
-           case (NORM_ONE)
+        if (present(order)) then 
+        
+            select case (order)
+               case (NORM_INF)
+                  lange_task = LANGE_NORM_INF 
+               case (NORM_ONE)
+                  lange_task = LANGE_NORM_ONE
+               case default 
+                  err = linalg_state_type(this,LINALG_VALUE_ERROR,'Integer order ',order,' is not a valid matrix norm input.')
+            end select
+            
+        else
+            
+            ! No user input: Frobenius norm
+            lange_task = LANGE_NORM_FRO
+            
+        endif
+     end subroutine lange_task_request_integer
+
+     pure subroutine lange_task_request_character(order,lange_task,err)
+        !> User input value
+        character(len=*), intent(in) :: order
+        !> Return value: norm type
+        character, intent(out) :: lange_task
+        !> State return flag
+        type(linalg_state_type), intent(out) :: err
+        
+        integer(ilp) :: int_order,read_err
+        
+        select case (order)
+           case ('inf','Inf','INF')
+              lange_task = LANGE_NORM_INF
+           case ('Euclidean','euclidean','EUCLIDEAN','Frobenius','frobenius','FROBENIUS','Fro','fro','frob')
+              lange_task = LANGE_NORM_FRO
+           case default
+            
+              ! Check if this input can be read as an integer
+              read(order,*,iostat=read_err) int_order
+              if (read_err/=0 .or. int_order/=1) then 
+                 ! Cannot read as an integer                 
+                 err = linalg_state_type(this,LINALG_ERROR,'Matrix norm input',order,' is not recognized.')                 
+              endif
               lange_task = LANGE_NORM_ONE
-           case (NORM_TWO)
-              lange_task = LANGE_NORM_TWO
-           case default 
-              err = linalg_state_type(this,LINALG_VALUE_ERROR,'Order ',norm_type,' is not a valid matrix norm input.')
-        end select
-     end subroutine lange_task_request
+
+        end select    
+        
+     end subroutine lange_task_request_character
 
     #:for rk,rt,ri in ALL_KINDS_TYPES
     
@@ -399,12 +440,12 @@ ${loop_variables_end(rank-1," "*12)}$
         !> Norm of the matrix.        
         real(${rk}$) :: nrm
         !> Order of the matrix norm being computed.
-        ${it}$, intent(in) :: order
+        ${it}$, #{if 'integer' in it}#optional, #{endif}#intent(in) :: order
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), intent(out), optional :: err           
         
         type(linalg_state_type) :: err_
-        integer(ilp) :: m,n,norm_request
+        integer(ilp) :: m,n
         character :: lange_task
         real(${rk}$), target :: work1(1)
         real(${rk}$), pointer :: work(:)        
@@ -422,8 +463,7 @@ ${loop_variables_end(rank-1," "*12)}$
         end if
 
         ! Check norm request: user + *LANGE support
-        call parse_norm_type(order,norm_request,err_)
-        call lange_task_request(norm_request,lange_task,err_)
+        call lange_task_request(order,lange_task,err_)
         if (err_%error()) then 
             call linalg_error_handling(err_,err)
             return
@@ -451,7 +491,7 @@ ${loop_variables_end(rank-1," "*12)}$
         !> Norm of the matrix.        
         real(${rk}$), allocatable :: nrm${ranksuffix(rank-2)}$
         !> Order of the matrix norm being computed.
-        ${it}$, intent(in) :: order
+        ${it}$, #{if 'integer' in it}#optional, #{endif}#intent(in) :: order
         !> [optional] dimensions of the sub-matrices the norms should be evaluated at (default = [1,2])
         integer(ilp), optional, intent(in) :: dim(2)
         !> [optional] state return flag. On error if not requested, the code will stop
@@ -485,8 +525,7 @@ ${loop_variables_end(rank-1," "*12)}$
         endif
         
         ! Check norm request: user + *LANGE support
-        call parse_norm_type(order,norm_request,err_)
-        call lange_task_request(norm_request,lange_task,err_)
+        call lange_task_request(order,lange_task,err_)
         if (err_%error()) then 
             allocate(nrm${emptyranksuffix(rank-2)}$)
             call linalg_error_handling(err_,err)
diff --git a/test/linalg/test_linalg_mnorm.fypp b/test/linalg/test_linalg_mnorm.fypp
@@ -49,7 +49,7 @@ module test_linalg_mnorm
         if (allocated(error)) return
 
         ! 2-norm (Frobenius norm)
-        call check(error, abs(mnorm(A, '2', err) - sqrt(sum(A**2))) < tol*mnorm(A, '2', err), &
+        call check(error, abs(mnorm(A, err=err) - sqrt(sum(A**2))) < tol*mnorm(A, err=err), &
                    'Matrix Frobenius norm does not match expected value')
         if (allocated(error)) return
 
@@ -65,15 +65,16 @@ module test_linalg_mnorm
     subroutine test_mnorm_${ri}$_${rank}$d(error)
         type(error_type), allocatable, intent(out) :: error
 
-        integer(ilp) :: i,j,k,l,dim1,dim2,dim(2),dim_sizes(2),order,ptr(${rank}$)
-        integer(ilp), parameter :: orders(*) = [1_ilp,2_ilp,huge(0_ilp)]
+        integer(ilp) :: i,j,k,l,dim1,dim2,dim(2),dim_sizes(2),ptr(${rank}$)
+        character(3), parameter :: orders(*) = ['1  ','fro','inf']
         integer(ilp), parameter :: ndim = ${rank}$
         integer(ilp), parameter :: n    = 2_ilp**ndim
         integer(ilp), parameter :: dims(*) = [(dim1, dim1=1,ndim)]
         real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$))
         real(${rk}$) :: one_nrm
         real(${rk}$), allocatable :: bnrm${ranksuffix(rank-2)}$
         ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$, one_mat(:,:)
+        character(:), allocatable :: order
         
         character(64) :: msg
         
@@ -86,7 +87,7 @@ module test_linalg_mnorm
         
         ! Test norm as collapsed around dimensions
         do k = 1, size(orders)
-            order = orders(k)
+            order = trim(orders(k))
             do dim1 = 1, ndim
                 do dim2 = dim1+1, ndim
                                     
@@ -97,7 +98,7 @@ module test_linalg_mnorm
                     bnrm = mnorm(b,order,dim)
                                    
                     ! Assert size
-                    write(msg,"('dim=[',i0,',',i0,'] order=',i0,' ${rk}$ norm returned wrong shape')") dim, order
+                    write(msg,"('dim=[',i0,',',i0,'] order=',a,' ${rk}$ norm returned wrong shape')") dim, order
                     call check(error,all(shape(bnrm)==pack(shape(b),dims/=dim1 .and. dims/=dim2) ), trim(msg))
                     if (allocated(error)) return    
                     
@@ -117,7 +118,7 @@ module test_linalg_mnorm
                        end do
                        one_nrm = mnorm(one_mat,order)
                     
-                       write(msg,"('dim=[',i0,',',i0,'] order=',i0,' ${rk}$ ',i0,'-th norm is wrong')") dim, order, l
+                       write(msg,"('dim=[',i0,',',i0,'] order=',a,' ${rk}$ ',i0,'-th norm is wrong')") dim, order, l
                        call check(error, abs(one_nrm-bnrm(${fixedranksuffix(rank-2,'l')}$))<tol*one_nrm, trim(msg))
                        if (allocated(error)) return    
                        deallocate(one_mat) 
