Bug #10178: norm returned error message: 'division by zero' for some sparse matrices.
[scilab.git] / scilab / modules / overloading / macros / %sp_norm.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA
3 // Copyright (C) Scilab Enterprises - Adeline CARNIS
4 // 
5 // This file must be used under the terms of the CeCILL.
6 // This source file is licensed as described in the file COPYING, which
7 // you should have received as part of this distribution.  The terms
8 // are also available at    
9 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10
11 function res=%sp_norm(S,flag)
12
13     [lhs,rhs]=argn(0)
14     if rhs==1 then flag=2;end //norm(S)=norm(S,2)
15     [m,n]=size(S)
16
17     if m==1|n==1 then //vector norm
18         [ij,v]=spget(S);
19         res=norm(v,flag);
20         return
21     end
22
23     select flag
24     case 1 then
25         res=max(ones(1,m)*abs(S))
26     case 2 then
27         if m<n then
28             S1=S*S'
29         elseif m>n then
30             S1=S'*S
31         else
32             S1 = S;
33         end
34
35         tol=%eps;
36         x = sum(abs(S1),1)';
37         res = norm(x);
38         if res==0 then return,end
39         x = x/res;res0 = 0;
40         while abs(res-res0) > tol*res
41             res0 = res;   Sx = S1*x;
42             
43             // Bug #10178: norm failed for some sparse matrix
44             // If Sx = 0, we had "division by zero" with x/norm(x)
45             // So, use to sum(abs(S).^2).^(1/2)
46             if Sx == 0 then
47                 res = sum(abs(S).^2).^(1/2);
48                 return
49             end
50             // End Bug #10178
51             
52             res = norm(Sx);
53             x = S1'*Sx;
54             x = x/norm(x);
55         end
56         if m<>n then res=sqrt(res),end
57     case %inf then
58         res=max(abs(S)*ones(n,1))
59     case 'inf' then
60         res=max(abs(S)*ones(n,1))
61     case 'fro' then
62         [ij,v]=spget(S);
63         res=sqrt(sum(abs(v.*v)))
64     else
65         [ij,v]=spget(S);
66         res=norm(v,flag);
67     end
68 endfunction