[Shootout-list] Numerical medium-sized benchmarks

[Sebastien Loisel]

------=_Part_454_31391994.1111745285716
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

The included ASCII source got munged by either gmail or the web
archive of the mailing list, so I'm resending with attachments
instead.

------=_Part_454_31391994.1111745285716
Content-Type: text/plain; name=numbench-matrix-norm.c; charset=us-ascii
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment; filename="numbench-matrix-norm.c"

#include <stdio.h>
#include <math.h>

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

void eval_A_times_u(int N, const double u[], double Au[])
{
  int i,j;
  for(i=0;i<N;i++)
    {
      Au[i]=0;
      for(j=0;j<N;j++) Au[i]+=eval_A(i,j)*u[j];
    }
}

void eval_At_times_u(int N, const double u[], double Au[])
{
  int i,j;
  for(i=0;i<N;i++)
    {
      Au[i]=0;
      for(j=0;j<N;j++) Au[i]+=eval_A(j,i)*u[j];
    }
}

void eval_AtA_times_u(int N, const double u[], double AtAu[])
{ double v[N]; eval_A_times_u(N,u,v); eval_At_times_u(N,v,AtAu); }

int main()
{
  int N=2000,i;
  double u[N],v[N],vBv,vv;
  for(i=0;i<N;i++) u[i]=1;
  for(i=0;i<10;i++)
    {
      eval_AtA_times_u(N,u,v);
      eval_AtA_times_u(N,v,u);
    }
  vBv=vv=0;
  for(i=0;i<N;i++) { vBv+=u[i]*v[i]; vv+=v[i]*v[i]; }
  printf("%0.20f\n",sqrt(vBv/vv));
  return 0;
}


------=_Part_454_31391994.1111745285716
Content-Type: application/octet-stream; name=numbench-matrix-norm.ml
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment; filename="numbench-matrix-norm.ml"

let eval_A i j = 1.0 /. (float_of_int ((i+j)*(i+j+1)/2+i+1))
let eval_A_times_u u v =
  for i = 0 to (Array.length v) - 1 do
    v.(i) <- 0.0 ;
    for j = 0 to (Array.length u) - 1 do
      v.(i) <- v.(i) +. (eval_A i j) *. u.(j)
    done
  done
let eval_At_times_u u v =
  for i = 0 to (Array.length v)-1 do
    v.(i)<-0.0 ;
    for j = 0 to (Array.length u)-1 do
      v.(i)<-v.(i) +. (eval_A j i) *. u.(j)
    done
  done
let print_array u =
  for i = 0 to (Array.length u)-1 do
    Printf.printf " %f" u.(i)
  done;
  Printf.printf "\n"

let eval_AtA_times_u u v =
  let w = Array.create (Array.length u) 0.0 in
  eval_A_times_u u w; eval_At_times_u w v

let n=2000
let u = Array.create n 1.0
let v = Array.create n 0.0 ;;

for i = 0 to 9 do
  eval_AtA_times_u u v; eval_AtA_times_u v u
done;;

let vv = ref 0.0
let vBv = ref 0.0 ;;

for i=0 to n-1 do
  vv := !vv +. v.(i) *. v.(i);
  vBv := !vBv +. u.(i) *. v.(i)
done;;

Printf.printf "%.20f\n" (sqrt (!vBv /. !vv));;

------=_Part_454_31391994.1111745285716
Content-Type: text/plain; name=numbench-matrix-norm.py; charset=us-ascii
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment; filename="numbench-matrix-norm.py"

#!/usr/bin/python
#numarray version left as an exercise

from math import sqrt
def eval_A(i,j):
	return 1.0/((i+j)*(i+j+1)/2+i+1)

def eval_A_times_u(u):
	v=[0]*len(u)
	for i in range(len(u)):
		for j in range(len(u)):
			v[i]=v[i]+eval_A(i,j)*u[j]
	return v

def eval_At_times_u(u):
	v=[0]*len(u)
	for i in range(len(u)):
		for j in range(len(u)):
			v[i]=v[i]+eval_A(j,i)*u[j]
	return v

def eval_AtA_times_u(u):
	return eval_At_times_u(eval_A_times_u(u))

n=2000
u=[1]*n
for i in range(10):
	v=eval_AtA_times_u(u)
	u=eval_AtA_times_u(v)
vBv=0
vv=0
for i in range(n):
	vBv=vBv+u[i]*v[i]
	vv=vv+v[i]*v[i]
print sqrt(vBv/vv)


------=_Part_454_31391994.1111745285716--