1;+
2; NAME: NORM
3;
4; PURPOSE:
5;       For vectors, returns a norm
6;       For matrix, returns the L0, L1 or L2 norm (only)
7;
8; CATEGORY:
9;       Algebra
10;
11; CALLING SEQUENCE:
12;       result=norm(a [, /double] [, /lnorm] )
13;
14; INPUTS:
15;       a       Vector or Matrix (can be complex)
16;
17; OUTPUTS:
18;       result  Norm (see Purpose)
19;
20; IDL DIFFERENCES:
21;       scalars are treated as vectors
22;       double=0 does not convert the result to float (like total)
23;
24; MODIFICATION HISTORY:
25;  12-Jan-2006 : written by Pierre Chanial
26;  12-Jul-2017 : AlainC., except the case (Lnorm=2 & matrix & Complex
27;  input), try to be complete.
28;
29; LICENCE:
30; Copyright (C) 2006, P. Chanial, 2017, A. Coulais
31; This program is free software; you can redistribute it and/or modify
32; it under the terms of the GNU General Public License as published by
33; the Free Software Foundation; either version 2 of the License, or
34; (at your option) any later version.
35;
36;-
37;
38function NORM, array, double=double, lnorm=lnormOut, test=test, help=help
39;
40ON_ERROR, 2
41;
42if KEYWORD_SET(help) then begin
43   print, "function NORM, array, double=double, lnorm=lnormOut, $"
44   print, "               test=test, help=help"
45   return, -1
46endif
47;
48type=SIZE(array, /type)
49ndim=SIZE(array, /n_dimensions)
50;
51if ((ndim LT 1) OR (ndim GT 2)) then begin
52   MESSAGE, 'Input must be an N-element vector or an M by N array.'
53endif
54;
55; "lnorm" value must be carefuly processed ...
56; single element of numerical type but not complex
57;
58if (N_ELEMENTS(lnormOut) GT 0) then begin
59   if (N_ELEMENTS(lnormOut) GT 1) then $
60      MESSAGE, 'LNORM must be a scalar or a 1-element array.'
61   lnorm=lnormOut[0]
62endif else begin
63   if (ndim EQ 1) then lnorm=2
64   if (ndim EQ 2) then lnorm=0
65endelse
66; "lnorm" must be numeric : integer or float but no complex
67if ~(ISA(lnorm,/integer) OR ISA(lnorm,/float)) then $
68   MESSAGE,/continue, 'bad type for LNORM, converted into Double'
69;
70; "lnorm" must be a double
71lnorm=DOUBLE(lnorm)
72;
73if KEYWORD_SET(test) then $
74   print, "input norm value, used norm value : ", lnormout, lnorm
75;
76; Do we have NaN in input array ?? O if no NaN, 1 otherwise
77nan=1-ARRAY_EQUAL(FINITE(array), 1)
78;
79; First we process the Vector case (ndim == 1),
80; no peculiar problem :)
81;
82if (ndim EQ 1) then begin
83   case lnorm of
84      ;; infinity norm
85      0: result=MAX(ABS(array), NAN=nan)
86      ;; L(1) norm
87      1: result=TOTAL(ABS(array), /DOUBLE, NAN=nan)
88      ;; L(2) norm (Euclidean)
89      2: result=SQRT(TOTAL(ABS(array)^2d, /DOUBLE, NAN=nan))
90      ;; General L(n) norm (tested with negatives values)
91      ;; print, norm(findgen(10)+1, lnorm=-10,/doubl)
92      ;; 0.99990060
93      else: result=(TOTAL(ABS(array)^lnorm, /DOUBLE, NAN=nan))^(1/lnorm)
94   endcase
95endif
96;
97if (ndim EQ 2) then begin
98   case lnorm of
99      ;; infinity norm (maximum absolute *row* sum norm)
100      0: result=MAX(TOTAL(ABS(array), 1, /DOUBLE, NAN=nan))
101      ;; L(1) norm (maximum absolute *column* sum norm)
102      1: result=MAX(TOTAL(ABS(array), 2, /DOUBLE, NAN=nan))
103      ;; L(2) norm (spectral norm)
104      2: begin
105         ;; we are not ready for Complex inputs in GDL
106         if ISA(array, /complex) then begin
107            MESSAGE, /continue, 'Complex Input not ready in GDL !'
108            MESSAGE, 'please contribute : LA_SVD missing :('
109         endif else begin
110            SVDC, array, W, U, V, /DOUBLE
111            result = MAX(W)
112         endelse
113      end
114      ;; no other case for Matrix
115      else: MESSAGE, 'For MATRIX input, LNORM must be equal to 0, 1, or 2'
116      endcase
117endif
118;
119; Which type for the result ? Float or Double !
120; Default is "float"
121result_type=4
122if KEYWORD_SET(double) then result_type=5
123if (type EQ 5) OR (type EQ 9) then result_type=5
124;
125if KEYWORD_SET(test) then STOP
126;
127return, FIX(result, TYPE=result_type)
128;
129end
130