1 |
|
2 |
|
3 | import java.text.DecimalFormat;
|
4 | import java.text.NumberFormat;
|
5 | import java.util.concurrent.CyclicBarrier;
|
6 |
|
7 | public class spectralnorm
|
8 | {
|
9 | private static final NumberFormat formatter = new DecimalFormat ("#.000000000");
|
10 |
|
11 | public static void main (String[] args)
|
12 | {
|
13 | long start = System.currentTimeMillis();
|
14 | System.out.println (formatter.format (spectralnormGame (5500)) );
|
15 | System.out.println(String.format("%dms", System.currentTimeMillis() - start));
|
16 | }
|
17 |
|
18 |
|
19 | private static final double spectralnormGame (int n)
|
20 | {
|
21 |
|
22 | double[] u = new double[n];
|
23 | double[] v = new double[n];
|
24 | double[] tmp = new double[n];
|
25 |
|
26 | for (int i = 0; i < n; i++)
|
27 | u[i] = 1.0;
|
28 |
|
29 |
|
30 | int nthread = Runtime.getRuntime ().availableProcessors ();
|
31 | Approximate.barrier = new CyclicBarrier (nthread);
|
32 |
|
33 | int chunk = n / nthread;
|
34 | Approximate[] ap = new Approximate[nthread];
|
35 |
|
36 | for (int i = 0; i < nthread; i++)
|
37 | {
|
38 | int r1 = i * chunk;
|
39 | int r2 = (i < (nthread -1)) ? r1 + chunk : n;
|
40 |
|
41 | ap[i] = new Approximate (u, v, tmp, r1, r2);
|
42 | }
|
43 |
|
44 |
|
45 | double vBv = 0, vv = 0;
|
46 | for (int i = 0; i < nthread; i++)
|
47 | {
|
48 | try
|
49 | {
|
50 | ap[i].join ();
|
51 |
|
52 | vBv += ap[i].m_vBv;
|
53 | vv += ap[i].m_vv;
|
54 | }
|
55 | catch (Exception e)
|
56 | {
|
57 | e.printStackTrace ();
|
58 | }
|
59 | }
|
60 |
|
61 | return Math.sqrt (vBv/vv);
|
62 | }
|
63 |
|
64 |
|
65 | private static class Approximate extends Thread
|
66 | {
|
67 | private static CyclicBarrier barrier;
|
68 |
|
69 | private double[] _u;
|
70 | private double[] _v;
|
71 | private double[] _tmp;
|
72 |
|
73 | private int range_begin, range_end;
|
74 | private double m_vBv = 0, m_vv = 0;
|
75 |
|
76 |
|
77 | public Approximate (double[] u, double[] v, double[] tmp, int rbegin, int rend)
|
78 | {
|
79 | super ();
|
80 |
|
81 | _u = u;
|
82 | _v = v;
|
83 | _tmp = tmp;
|
84 |
|
85 | range_begin = rbegin;
|
86 | range_end = rend;
|
87 |
|
88 | start ();
|
89 | }
|
90 |
|
91 | public void run ()
|
92 | {
|
93 |
|
94 | for (int i = 0; i < 10; i++)
|
95 | {
|
96 | MultiplyAtAv (_u, _tmp, _v);
|
97 | MultiplyAtAv (_v, _tmp, _u);
|
98 | }
|
99 |
|
100 | for (int i = range_begin; i < range_end; i++)
|
101 | {
|
102 | m_vBv += _u[i] * _v[i];
|
103 | m_vv += _v[i] * _v[i];
|
104 | }
|
105 | }
|
106 |
|
107 |
|
108 | private final static double eval_A (int i, int j)
|
109 | {
|
110 | int div = ( ((i+j) * (i+j+1) >>> 1) +i+1 );
|
111 | return 1.0 / div;
|
112 | }
|
113 |
|
114 |
|
115 | private final void MultiplyAv (final double[] v, double[] Av)
|
116 | {
|
117 | for (int i = range_begin; i < range_end; i++)
|
118 | {
|
119 | double sum = 0;
|
120 | for (int j = 0; j < v.length; j++)
|
121 | sum += eval_A (i, j) * v[j];
|
122 |
|
123 | Av[i] = sum;
|
124 | }
|
125 | }
|
126 |
|
127 |
|
128 | private final void MultiplyAtv (final double[] v, double[] Atv)
|
129 | {
|
130 | for (int i = range_begin; i < range_end; i++)
|
131 | {
|
132 | double sum = 0;
|
133 | for (int j = 0; j < v.length; j++)
|
134 | sum += eval_A (j, i) * v[j];
|
135 |
|
136 | Atv[i] = sum;
|
137 | }
|
138 | }
|
139 |
|
140 |
|
141 | private final void MultiplyAtAv (final double[] v, double[] tmp, double[] AtAv)
|
142 | {
|
143 | try
|
144 | {
|
145 | MultiplyAv (v, tmp);
|
146 |
|
147 | barrier.await ();
|
148 | MultiplyAtv (tmp, AtAv);
|
149 |
|
150 | barrier.await ();
|
151 | }
|
152 | catch (Exception e)
|
153 | {
|
154 | e.printStackTrace ();
|
155 | }
|
156 | }
|
157 | }
|
158 | }
|