benchmark/bench/spectralnorm.groovy

/*
 * The Computer Language Shootout
 * http://shootout.alioth.debian.org/
 *
 * contributed by Jochen Hinrichsen
 */

def approximate(n) {
    // create unit vector
    def u = [1.0D] * n

    // 20 steps of the power method
    def v = [0.0D] * n

    for (i in 1..10) {
        MultiplyAtAv(n, u, v)
        MultiplyAtAv(n, v, u)
    }

    // B=AtA         A multiplied by A transposed
    // v.Bv / (v.v)  eigenvalue of v
    double vBv = vv = 0.0D
    for (i in 0..<n) {
        vBv += u[i]*v[i]
        vv  += v[i]*v[i]
    }

    return Math.sqrt(vBv / vv)
}

/* return element i,j of infinite matrix A */
def A(i, j) {
    return 1.0D / ((i+j) * (i + j + 1.0D) / 2.0D  + i + 1.0D)
}

/* multiply vector v by matrix A */
def MultiplyAv(n, v, Av) {
    for (i in 0..<n) {
        Av[i] = 0.0D
        for (j in 0..<n) Av[i] += A(i,j) * v[j]
    }
}

/* multiply vector v by matrix A transposed */
def MultiplyAtv(n, v, Atv) {
    for (i in 0..<n) {
        Atv[i] = 0.0D
        for (j in 0..<n) Atv[i] += A(j,i) * v[j]
    }
}

/* multiply vector v by matrix A and then by matrix A transposed */
def MultiplyAtAv(n, v, AtAv) {
    double[] u = new double[n]
    MultiplyAv(n, v, u)
    MultiplyAtv(n, u, AtAv)
}

def n = (args.length == 0 ? 100 : args[0].toInteger())
def nf = java.text.NumberFormat.getInstance()
nf.setMaximumFractionDigits(9)
nf.setMinimumFractionDigits(9)
nf.setGroupingUsed(false)
println(nf.format(approximate(n)))