/*
 * Copyright (C) 2015 Linaro Limited. Ported from:
 *  https://github.com/WebKit/webkit/blob/master/PerformanceTests/SunSpider/tests/sunspider-1.0.2/math-spectral-norm.js
 *
 * Description:     Martrix spectral norm calculation.
 * Main Focus:      Floating-Point operations.
 *
 */

// The Great Computer Language Shootout
// http://shootout.alioth.debian.org/
//
// contributed by Ian Osgood

// http://benchmarksgame.alioth.debian.org/license.html  (BSD 3-clause license)
// See NOTICE file for license.

package org.linaro.benchmarks.math;

import java.lang.Math;
import java.lang.System;
import org.openjdk.jmh.annotations.*;
import java.util.concurrent.TimeUnit;

@BenchmarkMode(Mode.AverageTime)
@OutputTimeUnit(TimeUnit.MICROSECONDS)
@State(Scope.Benchmark)

public class MathSpectralNorm {
  private static final double SPECTRAL_NORM_EXPECTED = 5.086694231303284d;

  private double aa(int i, int j) {
    return 1.0 / ((i + j) * (i + j + 1) / 2 + i + 1);
  }

  private void au(double[] u, double[] v) {
    for (int i = 0; i < u.length; ++i) {
      double t = 0.0;
      for (int j = 0; j < u.length; ++j) {
        t += aa(i, j) * u[j];
      }
      v[i] = t;
    }
  }

  private void atu(double[] u, double[] v) {
    for (int i = 0; i < u.length; ++i) {
      double t = 0.0;
      for (int j = 0; j < u.length; ++j) {
        t += aa(j, i) * u[j];
      }
      v[i] = t;
    }
  }

  private void atAu(double[] u, double[] v, double[] w) {
    au(u, w);
    atu(w, v);
  }

  double spectralNorm(int n) {
    double[] u = new double[n];
    double[] v = new double[n];
    double[] w = new double[n];
    double vv;
    double vBv;

    vv = vBv = 0.0;
    for (int i = 0; i < n; ++i) {
      u[i] = 1.0;
    }

    for (int i = 0; i < 10; ++i) {
      atAu(u, v, w);
      atAu(v, u, w);
    }

    for (int i = 0; i < n; ++i) {
      vBv += u[i] * v[i];
      vv += v[i] * v[i];
    }
    return Math.sqrt(vBv / vv);
  }
  @Benchmark
  public void jmhTimeMathSpectralNorm() {
    for (int j = 6; j <= 48; j *= 2) {
      spectralNorm(j);
    }
  }
}