← Back to branch summary

~ubuntu-branches/ubuntu/raring/scilab/raring-proposed

« back to all changes in this revision

Viewing changes to modules/overloading/macros/%sp_norm.sci

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2012-08-30 14:42:38 UTC
  • mfrom: (1.4.7)
  • Revision ID: package-import@ubuntu.com-20120830144238-c1y2og7dbm7m9nig
Tags: 5.4.0-beta-3-1~exp1
* New upstream release
* Update the scirenderer dep
* Get ride of libjhdf5-java dependency

Show diffs side-by-side

added added

removed removed

Lines of Context:
modules/overloading/macros/%sp_norm.sci
1
1
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2
2
// Copyright (C) INRIA
 
3
// Copyright (C) Scilab Enterprises - Adeline CARNIS
3
4
// 
4
5
// This file must be used under the terms of the CeCILL.
5
6
// This source file is licensed as described in the file COPYING, which
9
10
 
10
11
function res=%sp_norm(S,flag)
11
12
 
12
 
[lhs,rhs]=argn(0)
13
 
if rhs==1 then flag=2;end //norm(S)=norm(S,2)
14
 
[m,n]=size(S)
15
 
 
16
 
if m==1|n==1 then //vector norm
17
 
  [ij,v]=spget(S);
18
 
  res=norm(v,flag);
19
 
  return
20
 
end
21
 
 
22
 
select flag
23
 
case 1 then
24
 
  res=max(ones(1,m)*abs(S))
25
 
case 2 then
26
 
  if m<n then
27
 
    S=S*S'
28
 
  elseif m>n then
29
 
    S=S'*S
30
 
  end
31
 
  tol=1.d-5;
32
 
  x = sum(abs(S),1)';
33
 
  res = norm(x);
34
 
  if res==0 then return,end
35
 
  x = x/res;res0 = 0;
36
 
  while abs(res-res0) > tol*res
37
 
    res0 = res;   Sx = S*x; res = norm(Sx);
38
 
    x = S'*Sx;
39
 
    x = x/norm(x);
40
 
  end
41
 
  if m<>n then res=sqrt(res),end
42
 
case %inf then
43
 
  res=max(abs(S)*ones(n,1))
44
 
case 'inf' then
45
 
  res=max(abs(S)*ones(n,1))
46
 
case 'fro' then
47
 
  [ij,v]=spget(S);
48
 
  res=sqrt(sum(abs(v.*v)))
49
 
else
50
 
  [ij,v]=spget(S);
51
 
  res=norm(v,flag);
52
 
end
 
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
53
68
endfunction