spectral norm benchmark

2011-12-07 15:36:07 +02:00
parent 6e7082ca5e
commit 41111c3bee
2 changed files with 297 additions and 0 deletions
@@ -0,0 +1,169 @@
+/*
+The Computer Language Benchmarks Game
+http://shootout.alioth.debian.org/
+
+Based on C# entry by Isaac Gouy
+contributed by Jarkko Miettinen
+Parallel by The Anh Tran
+*/
+
+import java.text.DecimalFormat;
+import java.text.NumberFormat;
+import java.util.concurrent.CyclicBarrier;
+
+public class SpectralNorm
+{
+    private static final NumberFormat formatter = new DecimalFormat ("#.000000000");
+
+    public static void main (String[] args)
+    {
+        int n = 5500;
+        if (args.length > 0) n = Integer.parseInt (args[0]);
+
+        final long millis = System.currentTimeMillis();
+        System.out.println (formatter.format (spectralnormGame (n)) );
+        long total = System.currentTimeMillis() - millis;
+        System.out.println("[SpectralNorm-" + System.getProperty("project.name")+ " Benchmark Result: " + total + "]");
+    }
+
+
+    private static final double spectralnormGame (int n)
+    {
+        // create unit vector
+        double[] u = new double[n];
+        double[] v = new double[n];
+        double[] tmp = new double[n];
+
+        for (int i = 0; i < n; i++)
+            u[i] = 1.0;
+
+        // get available processor, then set up syn object
+        int nthread = Runtime.getRuntime ().availableProcessors ();
+        Approximate.barrier = new CyclicBarrier (nthread);
+
+        int chunk = n / nthread;
+        Approximate[] ap = new Approximate[nthread];
+
+        for (int i = 0; i < nthread; i++)
+        {
+            int r1 = i * chunk;
+            int r2 = (i < (nthread -1)) ? r1 + chunk : n;
+
+            ap[i] = new Approximate (u, v, tmp, r1, r2);
+        }
+
+
+        double vBv = 0, vv = 0;
+        for (int i = 0; i < nthread; i++)
+        {
+            try
+            {
+                ap[i].join ();
+
+                vBv += ap[i].m_vBv;
+                vv += ap[i].m_vv;
+            }
+            catch (Exception e)
+            {
+                e.printStackTrace ();
+            }
+        }
+
+        return Math.sqrt (vBv/vv);
+    }
+
+
+    private static class Approximate extends Thread
+    {
+        private static CyclicBarrier barrier;
+
+        private double[] _u;
+        private double[] _v;
+        private double[] _tmp;
+
+        private int range_begin, range_end;
+        private double m_vBv = 0, m_vv = 0;
+
+
+        public Approximate (double[] u, double[] v, double[] tmp, int rbegin, int rend)
+        {
+            super ();
+
+            _u = u;
+            _v = v;
+            _tmp = tmp;
+
+            range_begin = rbegin;
+            range_end = rend;
+
+            start ();
+        }
+
+        public void run ()
+        {
+            // 20 steps of the power method
+            for (int i = 0; i < 10; i++)
+            {
+                MultiplyAtAv (_u, _tmp, _v);
+                MultiplyAtAv (_v, _tmp, _u);
+            }
+
+            for (int i = range_begin; i < range_end; i++)
+            {
+                m_vBv += _u[i] * _v[i];
+                m_vv  += _v[i] * _v[i];
+            }
+        }
+
+        /* return element i,j of infinite matrix A */
+        private final static double eval_A (int i, int j)
+        {
+            int div = ( ((i+j) * (i+j+1) >>> 1) +i+1 );
+            return 1.0 / div;
+        }
+
+        /* multiply vector v by matrix A, each thread evaluate its range only */
+        private final void MultiplyAv (final double[] v, double[] Av)
+        {
+            for (int i = range_begin; i < range_end; i++)
+            {
+                double sum = 0;
+                for (int j = 0; j < v.length; j++)
+                    sum += eval_A (i, j) * v[j];
+
+                Av[i] = sum;
+            }
+        }
+
+        /* multiply vector v by matrix A transposed */
+        private final void MultiplyAtv (final double[] v, double[] Atv)
+        {
+            for (int i = range_begin; i < range_end; i++)
+            {
+                double sum = 0;
+                for (int j = 0; j < v.length; j++)
+                    sum += eval_A (j, i) * v[j];
+
+                Atv[i] = sum;
+            }
+        }
+
+        /* multiply vector v by matrix A and then by matrix A transposed */
+        private final void MultiplyAtAv (final double[] v, double[] tmp, double[] AtAv)
+        {
+            try
+            {
+                MultiplyAv (v, tmp);
+                // all thread must syn at completion
+                barrier.await ();
+                MultiplyAtv (tmp, AtAv);
+                // all thread must syn at completion
+                barrier.await ();
+            }
+            catch (Exception e)
+            {
+                e.printStackTrace ();
+            }
+        }
+    }
+}
@@ -0,0 +1,128 @@
+namespace spectralnorm
+
+import java.text.DecimalFormat;
+import java.text.NumberFormat;
+import java.util.concurrent.CyclicBarrier;
+
+val formatter = DecimalFormat ("#.000000000");
+
+fun main (args : Array<String>) {
+    var n = 5500
+    if (args.size > 0)
+        n = Integer.parseInt (args[0]);
+
+    val millis = System.currentTimeMillis()
+    System.out?.println (formatter.format (spectralnormGame (n)) )
+    val total = System.currentTimeMillis() - millis;
+    System.out?.println("[SpectralNorm-" + System.getProperty("project.name")+ " Benchmark Result: " + total + "]");
+}
+
+fun spectralnormGame(n: Int) : Double {
+    val u = DoubleArray(n)
+    val v = DoubleArray(n)
+    val tmp = DoubleArray(n)
+
+    for(i in u.indices) {
+        u[i] = 1.0
+    }
+
+    val nthread = Runtime.getRuntime ().sure().availableProcessors ();
+    barrier = CyclicBarrier (nthread);
+
+    val chunk = n / nthread
+    val ap = Array<Approximate>(nthread,{
+        val r1 = it * chunk;
+        val r2 = if(it < (nthread -1)) r1 + chunk else nthread;
+
+        Approximate (u, v, tmp, r1, r2)
+    })
+
+    var vBv = 0.dbl
+    var vv = 0.dbl;
+    for (i in 0..nthread-1) {
+        try {
+            ap[i].join ();
+
+            vBv += ap[i].m_vBv;
+            vv += ap[i].m_vv;
+        }
+        catch (e: Exception )
+        {
+            e.printStackTrace ();
+        }
+    }
+
+    return Math.sqrt (vBv/vv);
+}
+
+fun eval_A (i: Int, j: Int) = 1.0 / ( ((i+j) * (i+j+1) shr 1) +i+1 )
+
+var barrier : CyclicBarrier? = null
+
+class Approximate(val u: DoubleArray, val v: DoubleArray, val _tmp: DoubleArray, val rbegin: Int, val rend: Int) : Thread() {
+    class object {
+    }
+
+    var m_vBv = 0.0
+    var m_vv  = 0.0
+
+    {
+        start()
+    }
+
+    override fun run () {
+        for(i in 0..10) {
+            MultiplyAtAv (u, _tmp, v);
+            MultiplyAtAv (v, _tmp, u);
+        }
+
+        for(i in rbegin..rend) {
+            m_vBv += u[i] * v[i];
+            m_vv  += v[i] * v[i];
+        }
+    }
+
+/* multiply vector v by matrix A, each thread evaluate its range only */
+    fun MultiplyAv (v: DoubleArray, Av: DoubleArray)
+    {
+        for (i in rbegin..rend)
+        {
+            var sum = 0.0;
+            for (j in v.indices)
+                sum += eval_A (i, j) * v[j];
+
+            Av[i] = sum;
+        }
+    }
+
+/* multiply vector v by matrix A transposed */
+    fun MultiplyAtv (v: DoubleArray, Atv: DoubleArray)
+    {
+        for (i in rbegin..rend)
+        {
+            var sum = 0.dbl
+            for (j in v.indices)
+                sum += eval_A (j, i) * v[j];
+
+            Atv[i] = sum;
+        }
+    }
+
+/* multiply vector v by matrix A and then by matrix A transposed */
+    fun MultiplyAtAv (v: DoubleArray, tmp: DoubleArray, AtAv: DoubleArray)
+    {
+        try
+        {
+            MultiplyAv (v, tmp);
+            // all thread must syn at completion
+            barrier?.await ();
+            MultiplyAtv (tmp, AtAv);
+            // all thread must syn at completion
+            barrier?.await ();
+        }
+        catch (e: Exception)
+        {
+            e.printStackTrace ();
+        }
+    }
+}