/*
* Licensed to the Apache Software Foundation (ASF) under one
* or more contributor license agreements. See the NOTICE file
* distributed with this work for additional information
* regarding copyright ownership. The ASF licenses this file
* to you under the Apache License, Version 2.0 (the
* "License"); you may not use this file except in compliance
* with the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing,
* software distributed under the License is distributed on an
* "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
* KIND, either express or implied. See the License for the
* specific language governing permissions and limitations
* under the License.
*/
package org.codehaus.groovy.benchmarks.vm5;
/*
The Great Computer Language Shootout
http://shootout.alioth.debian.org/
contributed by Java novice Jarkko Miettinen
modified ~3 lines of the original C#-version
by Isaac Gouy
*/
import java.text.DecimalFormat;
import java.text.NumberFormat;
public class spectralnorm
{
private static final NumberFormat formatter = new DecimalFormat("#.000000000");
public static void main(String[] args) {
long start = System.currentTimeMillis ();
int n = 100;
if (args.length > 0) n = Integer.parseInt(args[0]);
System.out.println(formatter.format(new spectralnorm().Approximate(n)));
System.out.println((System.currentTimeMillis () - start) + "ms");
}
private final 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 final double A(int i, int j){
return 1.0/((i+j)*(i+j+1)/2 +i+1);
}
/* multiply vector v by matrix A */
private final 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 final 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 final void MultiplyAtAv(int n, double[] v, double[] AtAv){
double[] u = new double[n];
MultiplyAv(n,v,u);
MultiplyAtv(n,u,AtAv);
}
}