Skip to content

Commit 1dbf88e

Browse files
perazzperazz
committed
add norms module
1 parent 1d27d27 commit 1dbf88e

File tree

2 files changed

+336
-0
lines changed

2 files changed

+336
-0
lines changed

‎src/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@ set(fppFiles
3131
stdlib_linalg_solve.fypp
3232
stdlib_linalg_determinant.fypp
3333
stdlib_linalg_inverse.fypp
34+
stdlib_linalg_norms.fypp
3435
stdlib_linalg_state.fypp
3536
stdlib_linalg_svd.fypp
3637
stdlib_linalg_cholesky.fypp

‎src/stdlib_linalg_norms.fypp

Lines changed: 335 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,335 @@
1+
#:include "common.fypp"
2+
#:set INPUT_TYPE = ["character(len=*)","integer(ilp)"]
3+
#:set INPUT_SHORT = ["char","int"]
4+
#:set INPUT_OPTIONS = list(zip(INPUT_TYPE,INPUT_SHORT))
5+
! Vector norms
6+
module stdlib_linalg_norms
7+
use stdlib_linalg_constants
8+
use stdlib_linalg_blas, only: nrm2
9+
use stdlib_linalg_lapack, only: lange
10+
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
11+
LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
12+
implicit none(type,external)
13+
private
14+
15+
public :: norm, get_norm
16+
17+
character(*), parameter :: this = 'norm'
18+
19+
!> List of internal norm flags
20+
integer(ilp), parameter :: NORM_ONE = 1_ilp
21+
integer(ilp), parameter :: NORM_TWO = 2_ilp
22+
integer(ilp), parameter :: NORM_POW_FIRST = 3_ilp
23+
integer(ilp), parameter :: NORM_INF = +huge(0_ilp) ! infinity norm
24+
integer(ilp), parameter :: NORM_POW_LAST = NORM_INF - 1_ilp
25+
integer(ilp), parameter :: NORM_MINUSINF = -huge(0_ilp)
26+
27+
!> Vector norm: function interface
28+
interface norm
29+
#:for rk,rt,ri in ALL_KINDS_TYPES
30+
#:for it,ii in INPUT_OPTIONS
31+
!> Scalar norms: ${rt}$
32+
#:for rank in range(1, MAXRANK + 1)
33+
module procedure stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$
34+
module procedure stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$
35+
#:endfor
36+
!> Array norms: ${rt}$
37+
#:for rank in range(2, MAXRANK + 1)
38+
module procedure stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$
39+
module procedure stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$
40+
#:endfor
41+
#:endfor
42+
#:endfor
43+
end interface norm
44+
45+
!> Vector norm: subroutine interface
46+
interface get_norm
47+
#:for rk,rt,ri in ALL_KINDS_TYPES
48+
#:for it,ii in INPUT_OPTIONS
49+
!> Scalar norms: ${rt}$
50+
#:for rank in range(1, MAXRANK + 1)
51+
module procedure norm_${rank}$D_${ii}$_${ri}$
52+
#:endfor
53+
!> Array norms: ${rt}$
54+
#:for rank in range(2, MAXRANK + 1)
55+
module procedure norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$
56+
#:endfor
57+
#:endfor
58+
#:endfor
59+
end interface get_norm
60+
61+
interface parse_norm_type
62+
module procedure parse_norm_type_integer
63+
module procedure parse_norm_type_character
64+
end interface parse_norm_type
65+
66+
contains
67+
68+
!> Parse norm type from an integer user input
69+
pure subroutine parse_norm_type_integer(order,norm_type,err)
70+
!> User input value
71+
integer(ilp), intent(in) :: order
72+
!> Return value: norm type
73+
integer(ilp), intent(out) :: norm_type
74+
!> State return flag
75+
type(linalg_state_type), intent(out) :: err
76+
77+
select case (order)
78+
case (1_ilp)
79+
norm_type = NORM_ONE
80+
case (2_ilp)
81+
norm_type = NORM_TWO
82+
case (3_ilp:huge(0_ilp)-1_ilp)
83+
norm_type = order
84+
case (huge(0_ilp):)
85+
norm_type = NORM_INF
86+
case (:-huge(0_ilp))
87+
norm_type = NORM_MINUSINF
88+
89+
case default
90+
norm_type = NORM_ONE
91+
err = linalg_state_type(this,LINALG_ERROR,'Input norm type ',order,' is not recognized.')
92+
end select
93+
94+
end subroutine parse_norm_type_integer
95+
96+
pure subroutine parse_norm_type_character(order,norm_type,err)
97+
!> User input value
98+
character(len=*), intent(in) :: order
99+
!> Return value: norm type
100+
integer(ilp), intent(out) :: norm_type
101+
!> State return flag
102+
type(linalg_state_type), intent(out) :: err
103+
104+
integer(ilp) :: int_order,read_err
105+
106+
select case (order)
107+
case ('inf','Inf','INF')
108+
norm_type = NORM_INF
109+
case ('-inf','-Inf','-INF')
110+
norm_type = NORM_MINUSINF
111+
case ('Euclidean','euclidean','EUCLIDEAN')
112+
norm_type = NORM_TWO
113+
case default
114+
115+
! Check if this input can be read as an integer
116+
read(order,*,iostat=read_err) int_order
117+
if (read_err/=0) then
118+
! Cannot read as an integer
119+
norm_type = NORM_ONE
120+
err = linalg_state_type(this,LINALG_ERROR,'Input norm type ',order,' is not recognized.')
121+
else
122+
call parse_norm_type_integer(int_order,norm_type,err)
123+
endif
124+
125+
end select
126+
127+
end subroutine parse_norm_type_character
128+
129+
#:for rk,rt,ri in ALL_KINDS_TYPES
130+
#:for it,ii in INPUT_OPTIONS
131+
132+
!==============================================
133+
! Norms : any rank to scalar
134+
!==============================================
135+
136+
#:for rank in range(1, MAXRANK + 1)
137+
138+
! Pure function interface, with order specification. On error, the code will stop
139+
pure function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$(a, order) result(nrm)
140+
!> Input ${rank}$-d matrix a${ranksuffix(rank)}$
141+
${rt}$, intent(in) :: a${ranksuffix(rank)}$
142+
!> Order of the matrix norm being computed.
143+
${it}$, intent(in) :: order
144+
!> Norm of the matrix.
145+
real(${rk}$) :: nrm
146+
147+
call norm_${rank}$D_${ii}$_${ri}$(a, nrm=nrm, order=order)
148+
149+
end function stdlib_linalg_norm_${rank}$D_order_${ii}$_${ri}$
150+
151+
! Function interface with output error
152+
function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$(a, order, err) result(nrm)
153+
!> Input ${rank}$-d matrix a${ranksuffix(rank)}$
154+
${rt}$, intent(in) :: a${ranksuffix(rank)}$
155+
!> Order of the matrix norm being computed.
156+
${it}$, intent(in) :: order
157+
!> Output state return flag.
158+
type(linalg_state_type), intent(out) :: err
159+
!> Norm of the matrix.
160+
real(${rk}$) :: nrm
161+
162+
call norm_${rank}$D_${ii}$_${ri}$(a, nrm=nrm, order=order, err=err)
163+
164+
end function stdlib_linalg_norm_${rank}$D_order_err_${ii}$_${ri}$
165+
166+
! Internal implementation
167+
pure subroutine norm_${rank}$D_${ii}$_${ri}$(a, nrm, order, err)
168+
!> Input ${rank}$-d matrix a${ranksuffix(rank)}$
169+
${rt}$, intent(in) :: a${ranksuffix(rank)}$
170+
!> Norm of the matrix.
171+
real(${rk}$), intent(out) :: nrm
172+
!> Order of the matrix norm being computed.
173+
${it}$, intent(in) :: order
174+
!> [optional] state return flag. On error if not requested, the code will stop
175+
type(linalg_state_type), intent(out), optional :: err
176+
177+
type(linalg_state_type) :: err_
178+
179+
integer(ilp) :: sze,norm_request
180+
real(${rk}$) :: rorder
181+
182+
sze = size(a,kind=ilp)
183+
184+
! Initialize norm to zero
185+
nrm = 0.0_${rk}$
186+
187+
! Check matrix size
188+
if (sze<=0) then
189+
err_ = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix shape: a=',shape(a,kind=ilp))
190+
call linalg_error_handling(err_,err)
191+
return
192+
end if
193+
194+
! Check norm request
195+
call parse_norm_type(order,norm_request,err_)
196+
if (err_%error()) then
197+
call linalg_error_handling(err_,err)
198+
return
199+
endif
200+
201+
select case(norm_request)
202+
case(NORM_ONE)
203+
nrm = sum( abs(a) )
204+
case(NORM_TWO)
205+
#:if rt.startswith('complex')
206+
nrm = sqrt( real( sum( a * conjg(a) ), ${rk}$) )
207+
#:else
208+
nrm = sqrt( sum( a ** 2 ) )
209+
#:endif
210+
case(NORM_INF)
211+
nrm = maxval( abs(a) )
212+
case(-NORM_INF)
213+
nrm = minval( abs(a) )
214+
case (NORM_POW_FIRST:NORM_POW_LAST)
215+
rorder = 1.0_${rk}$ / norm_request
216+
nrm = sum( abs(a) ** norm_request ) ** rorder
217+
case default
218+
err_ = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid norm type after checking')
219+
call linalg_error_handling(err_,err)
220+
end select
221+
222+
end subroutine norm_${rank}$D_${ii}$_${ri}$
223+
224+
#:endfor
225+
226+
!====================================================================
227+
! Norms : any rank to rank-1, with DIM specifier and ${ii}$ input
228+
!====================================================================
229+
230+
#:for rank in range(2, MAXRANK + 1)
231+
232+
! Pure function interface with DIM specifier. On error, the code will stop
233+
pure function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, order, dim) result(nrm)
234+
!> Input matrix a[..]
235+
${rt}$, intent(in), target :: a${ranksuffix(rank)}$
236+
!> Order of the matrix norm being computed.
237+
${it}$, intent(in) :: order
238+
!> Dimension to collapse by computing the norm w.r.t other dimensions
239+
integer(ilp), intent(in) :: dim
240+
!> Norm of the matrix.
241+
real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$
242+
243+
call norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim)
244+
245+
end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$
246+
247+
! Function interface with DIM specifier and output error state.
248+
function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$(a, order, dim, err) result(nrm)
249+
!> Input matrix a[..]
250+
${rt}$, intent(in), target :: a${ranksuffix(rank)}$
251+
!> Order of the matrix norm being computed.
252+
${it}$, intent(in) :: order
253+
!> Dimension to collapse by computing the norm w.r.t other dimensions
254+
integer(ilp), intent(in) :: dim
255+
!> Output state return flag.
256+
type(linalg_state_type), intent(out) :: err
257+
!> Norm of the matrix.
258+
real(${rk}$) :: nrm${reduced_shape('a', rank, 'dim')}$
259+
260+
call norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err)
261+
262+
end function stdlib_linalg_norm_${rank}$D_to_${rank-1}$D_err_${ii}$_${ri}$
263+
264+
! Internal implementation
265+
pure subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$(a, nrm, order, dim, err)
266+
!> Input matrix a[..]
267+
${rt}$, intent(in), target :: a${ranksuffix(rank)}$
268+
!> Dimension to collapse by computing the norm w.r.t other dimensions
269+
! (dim must be defined before it is used for `nrm`)
270+
integer(ilp), intent(in) :: dim
271+
!> Norm of the matrix.
272+
real(${rk}$), intent(out) :: nrm${reduced_shape('a', rank, 'dim')}$
273+
!> Order of the matrix norm being computed.
274+
${it}$, intent(in) :: order
275+
!> [optional] state return flag. On error if not requested, the code will stop
276+
type(linalg_state_type), intent(out), optional :: err
277+
278+
type(linalg_state_type) :: err_
279+
integer(ilp) :: sze,norm_request
280+
real(${rk}$) :: rorder
281+
282+
sze = size(a,kind=ilp)
283+
284+
! Initialize norm to zero
285+
nrm = 0.0_${rk}$
286+
287+
if (sze<=0) then
288+
err_ = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix shape: a=',shape(a,kind=ilp))
289+
call linalg_error_handling(err_,err)
290+
return
291+
end if
292+
293+
! Check dimension choice
294+
if (dim<1 .or. dim>${rank}$) then
295+
err_ = linalg_state_type(this,LINALG_VALUE_ERROR,'dimension ',dim, &
296+
'is out of rank for shape(a)=',shape(a,kind=ilp))
297+
call linalg_error_handling(err_,err)
298+
return
299+
end if
300+
301+
! Check norm request
302+
call parse_norm_type(order,norm_request,err_)
303+
if (err_%error()) then
304+
call linalg_error_handling(err_,err)
305+
return
306+
endif
307+
308+
select case(norm_request)
309+
case(NORM_ONE)
310+
nrm = sum( abs(a) , dim = dim )
311+
case(NORM_TWO)
312+
#:if rt.startswith('complex')
313+
nrm = sqrt( real( sum( a * conjg(a) , dim = dim ), ${rk}$) )
314+
#:else
315+
nrm = sqrt( sum( a ** 2 , dim = dim ) )
316+
#:endif
317+
case(NORM_INF)
318+
nrm = maxval( abs(a) , dim = dim )
319+
case(-NORM_INF)
320+
nrm = minval( abs(a) , dim = dim )
321+
case (NORM_POW_FIRST:NORM_POW_LAST)
322+
rorder = 1.0_${rk}$ / norm_request
323+
nrm = sum( abs(a) ** norm_request , dim = dim ) ** rorder
324+
case default
325+
err_ = linalg_state_type(this,LINALG_INTERNAL_ERROR,'invalid norm type after checking')
326+
call linalg_error_handling(err_,err)
327+
end select
328+
329+
end subroutine norm_${rank}$D_to_${rank-1}$D_${ii}$_${ri}$
330+
331+
#:endfor
332+
#:endfor
333+
#:endfor
334+
335+
end module stdlib_linalg_norms

0 commit comments

Comments
 (0)