1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
// =============================================================================
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) ????-2008 - INRIA Michael Baudin
//
// This file is distributed under the same license as the Scilab package.
// =============================================================================
// TODO : use relative error criteria instead of absolute error
eps=100*%eps;
//
// assert_close --
// Returns 1 if the two real matrices computed and expected are close,
// i.e. if the relative distance between computed and expected is lesser than epsilon.
// Arguments
// computed, expected : the two matrices to compare
// epsilon : a small number
//
function flag = assert_close ( computed , expected , epsilon )
if expected==0.0 then
shift = norm(computed-expected);
else
shift = norm(computed-expected)/norm(expected);
end
if shift < epsilon then
flag = 1;
else
flag = 0;
end
if flag <> 1 then bugmes();quit;end
endfunction
// inf
if norm([1,2,3,-1,-2,-3],0)<>%inf then bugmes();quit;end
if ~isnan(norm([1,2,3,-1,-2,-3],%nan)) then bugmes();quit;end
if norm([])<>0 then bugmes();quit;end
// vector
x=[1,2,3,-4];
if abs(norm(x,1) - sum(abs(x))) > eps then bugmes();quit;end
if abs(norm(x,2) - sqrt(sum(abs(x).*abs(x)))) > eps then bugmes();quit;end
if abs(norm(x,2) - norm(x)) > eps then bugmes();quit;end
p=0.5;
if abs(norm(x,p) - sum(abs(x)^p)^(1/p)) > eps then bugmes();quit;end
p=2.5;
if abs(norm(x,p) - sum(abs(x)^p)^(1/p)) > eps then bugmes();quit;end
if abs(norm(x,'inf') -maxi(abs(x))) > eps then bugmes();quit;end
if abs(norm(x,'inf') -norm(x,%inf)) > eps then bugmes();quit;end
if abs(norm(x,'fro') -norm(x,2)) > eps then bugmes();quit;end
// complex
x=x+%i*x;
if abs(norm(x,1) - sum(abs(x))) > eps then bugmes();quit;end
if abs(norm(x,2) - sqrt(sum(abs(x).*abs(x)))) > eps then bugmes();quit;end
if abs(norm(x,2) - norm(x)) > eps then bugmes();quit;end
p=0.5;
// 100*%eps is needed for linux
if abs(norm(x,p) - maxi(abs(x))*sum((abs(x)/maxi(abs(x)))^p)^(1/p))> 100*%eps then bugmes();quit;end
p=2.5;
if abs(norm(x,p) - maxi(abs(x))*sum((abs(x)/maxi(abs(x)))^p)^(1/p))> 100*%eps then bugmes();quit;end
if abs(norm(x,'inf') -maxi(abs(x)))> eps then bugmes();quit;end
if abs(norm(x,'inf') -norm(x,%inf)) > eps then bugmes();quit;end
if abs(norm(x,'fro') -norm(x,2))> eps then bugmes();quit;end
// scalar
x=[1.23];
if abs(norm(x,1) - sum(abs(x))) > eps then bugmes();quit;end
if abs(norm(x,2) - sqrt(sum(abs(x).*abs(x)))) > eps then bugmes();quit;end
if abs(norm(x,2) - norm(x)) > eps then bugmes();quit;end
p=0.5;
if abs(norm(x,p) - sum(abs(x)^p)^(1/p)) > eps then bugmes();quit;end
p=2.5;
if abs(norm(x,p) - sum(abs(x)^p)^(1/p)) > eps then bugmes();quit;end
if abs(norm(x,'inf') -maxi(abs(x))) > eps then bugmes();quit;end
if abs(norm(x,'inf') -norm(x,%inf)) > eps then bugmes();quit;end
if abs(norm(x,'fro') -norm(x,2)) > eps then bugmes();quit;end
// Matrices
a=rand(10,10,'g');
if abs(norm(a,1) - maxi(sum(abs(a),'r'))) > eps then bugmes();quit;end
if abs(norm(a,'inf') - maxi(sum(abs(a),'c'))) > eps then bugmes();quit;end
if abs(norm(a,%inf) - maxi(sum(abs(a),'c'))) > eps then bugmes();quit;end
if abs(norm(a,2) - maxi(svd(a))) > eps then bugmes();quit;end
if abs(norm(a,'fro') - norm(matrix(a,1,size(a,'*')),2)) > eps then bugmes();quit;end
a=a+%i*a;
if abs(norm(a,1) - maxi(sum(abs(a),'r'))) > eps then bugmes();quit;end
if abs(norm(a,'inf') - maxi(sum(abs(a),'c'))) > eps then bugmes();quit;end
if abs(norm(a,%inf) - maxi(sum(abs(a),'c'))) > eps then bugmes();quit;end
if abs(norm(a,2) - maxi(svd(a))) > eps then bugmes();quit;end
if abs(norm(a,'fro') - norm(matrix(a,1,size(a,'*')),2)) > eps then bugmes();quit;end
//
// Difficult cases for large/small vectors
//
//norm 2
x = 1.e307 * [1 1];
assert_close ( norm(x) , sqrt(2) * 1.e307 , %eps );
x = 1.e-307 * [1 1];
assert_close ( norm(x) , sqrt(2) * 1.e-307 , %eps );
// norm f
x = 1.e307 * [1 1];
assert_close ( norm(x,"f") , sqrt(2) * 1.e307 , %eps );
x = 1.e-307 * [1 1];
assert_close ( norm(x,"f") , sqrt(2) * 1.e-307 , %eps );
//
// Difficult cases for large/small matrices
//
// norm f - case 1 : n < m
x = 1.e307 * ones(10,20);
assert_close ( norm(x,"f") , sqrt(200) * 1.e307 , %eps );
x = 1.e-307 * ones(10,20);
assert_close ( norm(x,"f") , sqrt(200) * 1.e-307 , %eps );
// norm f - case 2 : n > m
x = 1.e307 * ones(20,10);
assert_close ( norm(x,"f") , sqrt(200) * 1.e307 , %eps );
x = 1.e-307 * ones(20,10);
assert_close ( norm(x,"f") , sqrt(200) * 1.e-307 , %eps );
//
// Special cases for zero vectors
//
// 2-norm of a zero vector
x=[0 0 0];
assert_close ( norm(x,2) , 0.0 , %eps );
// f-norm of a zero vector
x=zeros(4,1);
assert_close ( norm(x,"f") , 0.0 , %eps );
// f-norm of a zero matrix, case 1 n>m
x=zeros(4,2);
assert_close ( norm(x,"f") , 0.0 , %eps );
// f-norm of a zero matrix, case 2 m>n
x=zeros(2,4);
assert_close ( norm(x,"f") , 0.0 , %eps );
|
