Bug #9204: The Frobenius norm of a complex vector returns a real.
[scilab.git] / scilab / modules / linear_algebra / macros / norm.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) ????-2008 - INRIA
3 // Copyright (C) 2009 - INRIA Michael Baudin
4 // Copyright (C) 2012 - Scilab Enterprises - Adeline CARNIS
5 //
6 // This file must be used under the terms of the CeCILL.
7 // This source file is licensed as described in the file COPYING, which
8 // you should have received as part of this distribution.  The terms
9 // are also available at
10 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11
12 //
13 // norm --
14 //   Returns the norm of the given vector/matrix A.
15 //   Uses scaling to improve accuracy for Pythagorean sums.
16 // References
17 //   Moler C, Morrison D. 
18 //   Replacing square roots by pythagorean sums. 
19 //   IBM Journal of Research and Development 1983; 27(6):577-581. 
20 //
21 function y=norm(A,flag)
22     //compute various matrix norms
23     if argn(2)==1 then flag=2,end
24
25     if type(A)==1 then
26         if A==[] then y=0,return,end
27         if or(size(A)==1) then // vector norm
28             if type(flag)==10 then //'inf' or 'fro'
29                 select convstr(part(flag,1))
30                 case 'i' then //'inf'
31                     y=max(abs(A))
32                 case 'f' then //'fro'
33                     A=A(:)
34                     //
35                     // Scaling for better floating point accuracy.
36                     //
37                     //
38                     s = max(abs(A));
39                     if s==0.0 then
40                         y=sqrt(A'*A);
41                     else
42                         sA = A/s;
43                         // return real result
44                         y = s * sqrt(abs(sA'*sA));
45                     end
46                 else
47                     error("invalid value for flag")
48                 end
49             elseif type(flag)==1 then //p_norm
50                 p=flag;
51                 if ~isreal(p) then
52                     error('flag must be real')
53                 end
54                 if p==%inf then
55                     y=max(abs(A))
56                 elseif p==1 then
57                     y=sum(abs(A))
58                 elseif p==-%inf then
59                     y=min(abs(A))
60                 elseif isnan(p) then
61                     y=%nan
62                 elseif p==0 then
63                     y=%inf
64                 else
65                     //
66                     // Scaling for better floating point accuracy.
67                     //
68                     s = max(abs(A));
69                     if s==0.0 then
70                         y=sum(abs(A).^p)^(1/p);
71                     else
72                         sA = A/s;
73                         y = s * sum(abs(sA).^p)^(1/p);
74                     end
75                 end
76             else
77                 error("invalid value for flag")
78             end
79         else //matrix norm
80             if type(flag)==10 then //'inf' or 'fro'
81                 select convstr(part(flag,1))
82                 case 'i' then //'inf'
83                     y=max(sum(abs(A),2))  
84                 case 'f' then //'fro'
85                     //
86                     // Scaling for better floating point accuracy.
87                     //
88                     s = max(abs(A));
89                     if s==0.0 then
90                         if size(A,1)>size(A,2) then
91                             y=sqrt(sum(diag(A'*A))) 
92                         else
93                             y=sqrt(sum(diag(A*A'))) 
94                         end
95                     else
96                         sA = A/s;
97                         if size(A,1)>size(A,2) then
98                             // return real result
99                             y = s * sqrt(sum(abs(diag(sA'*sA))))
100                         else
101                             y = s * sqrt(sum(abs(diag(sA*sA'))))
102                         end
103                     end
104                 else
105                     error("invalid value for flag")
106                 end
107             elseif type(flag)==1 then //p_norm
108                 p=flag;
109                 select p
110                 case 1 then 
111                     y=max(sum(abs(A),1))
112                 case 2 then
113                     y=max(svd(A))
114                 case %inf then
115                     y=max(sum(abs(A),2))  
116                 else
117                     error('flag must be 1 2 or inf')
118                 end
119             else
120                 error("invalid value for flag")
121             end    
122         end
123     else
124         if type(A)==16|type(A)==17 then
125             n=getfield(1,A);n=n(1)
126         else
127             [t,n]=typename()
128             n=stripblanks(n(find(t==type(A))))
129         end
130         fun='%'+n+'_norm'
131         if exists(fun)==1 then
132             execstr('y='+fun+'(A,flag)')
133         else
134             error('norm not defined for type ""'+n+'"" .'+..
135             'Check argument or define function '+fun)
136         end
137     end
138 endfunction
139