File tree Expand file tree Collapse file tree 1 file changed +70
-0
lines changed Expand file tree Collapse file tree 1 file changed +70
-0
lines changed Original file line number Diff line number Diff line change
1
+ // Based on spectalnorm.gcc by Sebastien Loisel
2
+
3
+ use std;
4
+
5
+ fn eval_A ( i : uint , j : uint ) -> float {
6
+ 1.0 /( ( ( i+j) * ( i+j+1 u) /2 u+i+1 u) as float )
7
+ }
8
+
9
+ fn eval_A_times_u ( u : [ const float ] , Au : [ mutable float] ) {
10
+ let N = vec:: len ( u) ;
11
+ let i = 0 u;
12
+ while i < N {
13
+ Au [ i] = 0.0 ;
14
+ let j = 0 u;
15
+ while j < N {
16
+ Au [ i] += eval_A ( i, j) * u[ j] ;
17
+ j += 1 u;
18
+ }
19
+ i += 1 u;
20
+ }
21
+ }
22
+
23
+ fn eval_At_times_u ( u : [ const float ] , Au : [ mutable float] ) {
24
+ let N = vec:: len ( u) ;
25
+ let i = 0 u;
26
+ while i < N {
27
+ Au [ i] = 0.0 ;
28
+ let j = 0 u;
29
+ while j < N {
30
+ Au [ i] += eval_A ( j, i) * u[ j] ;
31
+ j += 1 u;
32
+ }
33
+ i += 1 u;
34
+ }
35
+ }
36
+
37
+ fn eval_AtA_times_u ( u : [ const float ] , AtAu : [ mutable float] ) {
38
+ let v = vec:: init_elt_mut ( vec:: len ( u) , 0.0 ) ;
39
+ eval_A_times_u ( u, v) ;
40
+ eval_At_times_u ( v, AtAu ) ;
41
+ }
42
+
43
+ fn main ( args : [ str ] ) {
44
+
45
+ let N = if vec:: len ( args) == 2 u {
46
+ uint:: from_str ( args[ 1 ] )
47
+ } else {
48
+ 1000 u
49
+ } ;
50
+
51
+ let u = vec:: init_elt_mut ( N , 1.0 ) ;
52
+ let v = vec:: init_elt_mut ( N , 0.0 ) ;
53
+ let i = 0 u;
54
+ while i < 10 u {
55
+ eval_AtA_times_u ( u, v) ;
56
+ eval_AtA_times_u ( v, u) ;
57
+ i += 1 u;
58
+ }
59
+
60
+ let vBv = 0.0 ;
61
+ let vv = 0.0 ;
62
+ let i = 0 u;
63
+ while i < N {
64
+ vBv += u[ i] * v[ i] ;
65
+ vv += v[ i] * v[ i] ;
66
+ i += 1 u;
67
+ }
68
+
69
+ std:: io:: println ( #fmt ( "%0.9f\n " , float:: sqrt ( vBv / vv) ) ) ;
70
+ }
You can’t perform that action at this time.
0 commit comments