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
125
126
127
128
129
|
// Licensed to the .NET Foundation under one or more agreements.
// The .NET Foundation licenses this file to you under the MIT license.
// See the LICENSE file in the project root for more information.
/* The Computer Language Benchmarks Game
http://benchmarksgame.alioth.debian.org/
contributed by Isaac Gouy
modified for use with xunit-performance
*/
using Microsoft.Xunit.Performance;
using System;
[assembly: OptimizeForBenchmarks]
[assembly: MeasureInstructionsRetired]
namespace BenchmarksGame
{
public class SpectralNorm
{
#if DEBUG
public const int Iterations = 1;
#else
public const int Iterations = 300;
#endif
public static int Main(String[] args)
{
int n = 100;
if (args.Length > 0) n = Int32.Parse(args[0]);
double norm = new SpectralNorm().Approximate(n);
Console.WriteLine("Norm={0:f9}", norm);
double expected = 1.274219991;
bool result = Math.Abs(norm - expected) < 1e-4;
return (result ? 100 : -1);
}
[Benchmark]
public static void Bench()
{
int n = 100;
foreach (var iteration in Benchmark.Iterations)
{
double a = 0;
using (iteration.StartMeasurement())
{
for (int i = 0; i < Iterations; i++)
{
SpectralNorm s = new SpectralNorm();
a += s.Approximate(n);
}
}
double norm = a / Iterations;
double expected = 1.274219991;
bool valid = Math.Abs(norm - expected) < 1e-4;
if (!valid)
{
throw new Exception("Benchmark failed to validate");
}
}
}
private double Approximate(int n)
{
// create unit vector
double[] u = new double[n];
for (int i = 0; i < n; i++) u[i] = 1;
// 20 steps of the power method
double[] v = new double[n];
for (int i = 0; i < n; i++) v[i] = 0;
for (int i = 0; i < 10; i++)
{
MultiplyAtAv(n, u, v);
MultiplyAtAv(n, v, u);
}
// B=AtA A multiplied by A transposed
// v.Bv /(v.v) eigenvalue of v
double vBv = 0, vv = 0;
for (int i = 0; i < n; i++)
{
vBv += u[i] * v[i];
vv += v[i] * v[i];
}
return Math.Sqrt(vBv / vv);
}
/* return element i,j of infinite matrix A */
private double A(int i, int j)
{
return 1.0 / ((i + j) * (i + j + 1) / 2 + i + 1);
}
/* multiply vector v by matrix A */
private void MultiplyAv(int n, double[] v, double[] Av)
{
for (int i = 0; i < n; i++)
{
Av[i] = 0;
for (int j = 0; j < n; j++) Av[i] += A(i, j) * v[j];
}
}
/* multiply vector v by matrix A transposed */
private void MultiplyAtv(int n, double[] v, double[] Atv)
{
for (int i = 0; i < n; i++)
{
Atv[i] = 0;
for (int j = 0; j < n; j++) Atv[i] += A(j, i) * v[j];
}
}
/* multiply vector v by matrix A and then by matrix A transposed */
private void MultiplyAtAv(int n, double[] v, double[] AtAv)
{
double[] u = new double[n];
MultiplyAv(n, v, u);
MultiplyAtv(n, u, AtAv);
}
}
}
|
