4829dc80f9ee0d1d6cb07b481e14f1cec834456e
[scilab.git] / scilab / modules / linear_algebra / macros / norm.sci
1
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) ????-2008 - INRIA
4 // Copyright (C) 2009 - INRIA Michael Baudin
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           y = s * sqrt(sA'*sA);
44         end
45       else
46         error("invalid value for flag")
47       end
48     elseif type(flag)==1 then //p_norm
49       p=flag;
50       if ~isreal(p) then
51         error('flag must be real')
52       end
53       if p==%inf then
54         y=max(abs(A))
55       elseif p==1 then
56         y=sum(abs(A))
57       elseif p==-%inf then
58         y=min(abs(A))
59       elseif isnan(p) then
60         y=%nan
61       elseif p==0 then
62         y=%inf
63       else
64         //
65         // Scaling for better floating point accuracy.
66         //
67         s = max(abs(A));
68         if s==0.0 then
69           y=sum(abs(A).^p)^(1/p);
70         else
71           sA = A/s;
72           y = s * sum(abs(sA).^p)^(1/p);
73         end
74       end
75     else
76       error("invalid value for flag")
77     end
78   else //matrix norm
79     if type(flag)==10 then //'inf' or 'fro'
80       select convstr(part(flag,1))
81       case 'i' then //'inf'
82         y=max(sum(abs(A),2))  
83         case 'f' then //'fro'
84         //
85         // Scaling for better floating point accuracy.
86         //
87         s = max(abs(A));
88         if s==0.0 then
89           if size(A,1)>size(A,2) then
90             y=sqrt(sum(diag(A'*A))) 
91           else
92             y=sqrt(sum(diag(A*A'))) 
93           end
94         else
95           sA = A/s;
96           if size(A,1)>size(A,2) then
97             y = s * sqrt(sum(diag(sA'*sA)))
98           else
99             y = s * sqrt(sum(diag(sA*sA')))
100           end
101         end
102       else
103         error("invalid value for flag")
104       end
105     elseif type(flag)==1 then //p_norm
106       p=flag;
107       select p
108         case 1 then 
109         y=max(sum(abs(A),1))
110         case 2 then
111         y=max(svd(A))
112         case %inf then
113         y=max(sum(abs(A),2))  
114       else
115         error('flag must be 1 2 or inf')
116       end
117     else
118       error("invalid value for flag")
119     end    
120   end
121 else
122   if type(A)==16|type(A)==17 then
123     n=getfield(1,A);n=n(1)
124   else
125     [t,n]=typename()
126     n=stripblanks(n(find(t==type(A))))
127   end
128   fun='%'+n+'_norm'
129   if exists(fun)==1 then
130     execstr('y='+fun+'(A,flag)')
131   else
132     error('norm not defined for type ""'+n+'"" .'+..
133           'Check argument or define function '+fun)
134   end
135 end
136 endfunction
137