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