linalg: vector norms

[perazz]

Address #820 in joint effort with @jalvesz.
Compute several vector norms. Array $A$, real or complex, has general rank n>=1.

Proposed implementation

• x = norm(a, order=2 [, dim=dim] [, err=err]): (pure) function interface
• call get_norm(a, nrm=x, order=2 [, dim=dim] [, err=err]): pure subroutine interface

Key facts

• The following implementations are provided:

  • 1, '1': 1-norm, $\sum_i{ \left|a_i\right| }$
  • 2, '2', 'Euclidean': 2-norm, $\sqrt{\sum_i{ a_i^2 }}$
  • >=3: p-norm, $\left( \sum_i{ \left|a_i\right| ^p }\right) ^{1/p}$
  • 'inf', huge(0): $\infty$-norm, $\max_i{ \left|a_i\right| }$
  • '-inf', -huge(0): minimum norm, $\min_i{ \left|a_i\right| }$

• order can either be an integer (1, 2, ..., n, -huge(0), huge(0)) or a character input ('1', '2', '10', 'inf', '-inf', 'euclidean', ...)

• The implementation currently only uses Fortran intrinsics that handle all dim cases by default: would be hard to unroll all them to enable calls to BLAS backends for all ranks 1:15.

Progress

• base implementation
• tests
• documentation
• submodule
• examples
• all pure subroutine interfaces, pure function interfaces if no error flag is requested.

Prior art

• Numpy: linalg.norm(x, ord=None, axis=None, keepdims=False)
• Scipy: norm(a, ord=None, axis=None, keepdims=False, check_finite=True)

cc: @jvdp1 @jalvesz @loiseaujc

[jvdp1]

Overall LGTM. Thank you @perazz and @jalvesz .
It could be worthwhile to add a note in the specs to highlight the differences (or not) with the intrinsic norm2.

[perazz]

@jvdp1 @jalvesz in 437b96e I've introduced the intrinsic norm2 where possible (for real inputs). This is going to be one of the most used norms. So I'm wondering if we should replace it with a call to the BLAS backend. this:

• would allow to use a fast implementation, if available
• for the N-D case (with dim), some data sorting cannot be avoided.

So, should we introduce the BLAS backend at least for the 1D case evaluation?

[awvwgk]

I believe the BLAS API should allow also strided access to array elements, which would support dim != 1, the stride length would need to be computed based on the selected axis and the dimensions.

[perazz]

The compiler will create a temporary array passing from (:) (may be strided) to the (*) BLAS API, because the latter implies contiguous data. When that will happen is nontrivial. So, it is possible to infer the stride of an (:) array like I did, but there is no guarantee that the compiler will not create a temporary array in between the calls. See this example.

So, in the interest of safety, I suggest to avoid inferring strides and rather create the temporary ourselves, i.e. preprocessing it with a reshape(a, order=sorted_dims)

[perazz]

• unify calls to norm evaluation function
• $+\infty$ norm: replace maxval(abs(.)) with calls to BLAS, add test
• ND and dim-med versions: always call LAPACK/BLAS rather than Fortran intrinsics
• address contiguity via reshape(a, order=sorted_dims)

[jalvesz]

LGTM @perazz! On my end I have no further comments.