← Back to branch summary

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

« back to all changes in this revision

Viewing changes to modules/linear_algebra/macros/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/linear_algebra/macros/norm.sci
1
 
 
2
1
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3
2
// Copyright (C) ????-2008 - INRIA
4
3
// Copyright (C) 2009 - INRIA Michael Baudin
 
4
// Copyright (C) 2012 - Scilab Enterprises - Adeline CARNIS
5
5
//
6
6
// This file must be used under the terms of the CeCILL.
7
7
// This source file is licensed as described in the file COPYING, which
19
19
//   IBM Journal of Research and Development 1983; 27(6):577-581. 
20
20
//
21
21
function y=norm(A,flag)
22
 
//compute various matrix norms
23
 
if argn(2)==1 then flag=2,end
 
22
    //compute various matrix norms
 
23
    if argn(2)==1 then flag=2,end
24
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
 
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
75
123
    else
76
 
      error("invalid value for flag")
 
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
77
137
    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
138
endfunction
137
139