From mboxrd@z Thu Jan  1 00:00:00 1970
Received: (from weis@localhost) by pauillac.inria.fr (8.7.6/8.7.3) id WAA03412 for caml-red; Fri, 26 Jan 2001 22:03:16 +0100 (MET)
Received: from concorde.inria.fr (concorde.inria.fr [192.93.2.39]) by pauillac.inria.fr (8.7.6/8.7.3) with ESMTP id LAA29805 for <caml-list@pauillac.inria.fr>; Thu, 25 Jan 2001 11:34:36 +0100 (MET)
Received: from mailgw1.fhg.de (mailgw1.fhg.de [153.96.1.2])
	by concorde.inria.fr (8.11.1/8.10.0) with ESMTP id f0PAYZL08533
	for <caml-list@inria.fr>; Thu, 25 Jan 2001 11:34:35 +0100 (MET)
Received: from ilt.fhg.de (ilt.ilt.fhg.de [153.96.180.2])
	by mailgw1.fhg.de (8.11.1/8.11.1) with ESMTP id f0PAYWH15498
	for <caml-list@inria.fr>; Thu, 25 Jan 2001 11:34:32 +0100 (MET)
Received: from eent230 (eelta136.llt.RWTH-Aachen.DE [137.226.47.136])
	by ilt.fhg.de (8.9.3/8.9.3) with SMTP id LAA06543
	for <caml-list@inria.fr>; Thu, 25 Jan 2001 11:34:31 +0100 (MET)
Message-Id: <200101251034.LAA06543@ilt.fhg.de>
From: wester@ilt.fhg.de
To: caml-list@inria.fr
Date: Thu, 25 Jan 2001 11:34:32 +0100
MIME-Version: 1.0
Content-type: text/plain; charset=US-ASCII
Content-transfer-encoding: 7BIT
Subject: Floating point array computations
Priority: normal
X-mailer: Pegasus Mail for Win32 (v3.01d)
Sender: weis@pauillac.inria.fr

Hi,

I'm a OCmal beginner and would like to use it for some numerical
computations. Speed is not the only criterion but for number crunching it
is an important point. So I would like to make my OCaml code as fast 
as possible. I have tried some small examples (source attached below) 
and have some questions. 

In sample code form the ICFP 200 contest I found:

type v = float array

is array a build in type?

The functions main1 .. main4 below do the same but take different time for it
(compiled with ocamlopt). I expected main2 to be the fastest but main3 is 
fastest. main3 300 300 1000 takes on my PC (450 Mhz PI, 520 MB RAM) 
26 sec, this is about twice the time of C++ and about the same as Java.

Can somone tell me why main3 is faster than main2 (36 sec) although
main3 makes function calls to step?

I also tried Array.unsafe_get/set in main4. I expected this to be even faster
because I have seen it in the ICFP contest sample code. But actually 
it is much slower (about 100 sec). Did I make anything wrong?

Would Bigarray perform better? 

I would be very appriciative for your help to clarify this points.

Thank you in advance.

Rolf


-----------------------------------------------------------------------------------------------------------

let main1 nx ny nt =
    let matrix1 = Array.make_matrix nx ny 0.0 in
    let matrix2 = Array.make_matrix nx ny 0.0 in
    let sum   = ref 0  in
    let start = ref 0. in
    let upper_idx = Array.length matrix1 - 1 in
    for r = 1 to nt do
        for i = 1 to upper_idx - 1  do
		  let mat2  = matrix2.(i)   and
		      mat1m = matrix1.(i-1) and
			  mat10 = matrix1.(i)   and
			  mat1p = matrix1.(i+1) in
            for j = 1 to upper_idx - 1 do
                mat2.(j) <- mat10.(j) +. 1.0  +.
				                0.1*.(mat1p.(j)   +.
                                      mat1m.(j)   +.
                                      mat10.(j+1) +.
                                      mat10.(j-1) -.
                               4.0 *. mat10.(j));
            done;
        done;
        for i = 1 to upper_idx - 1  do
            for j = 1 to upper_idx - 1 do
                matrix1.(i).(j) <- matrix2.(i).(j) +. 1.0 +.
				0.1*.(matrix2.(i+1).(j) +.
                                      matrix2.(i-1).(j) +.
                                      matrix2.(i).(j+1) +.
                                      matrix2.(i).(j-1) -.
                                   4.0 *. matrix2.(i).(j));
            done;
        done;
	Printf.printf "%f\n" (matrix1.(150).(150));
	flush stdout;
    done;
;;


let main2 nx ny nt =
    let matrix1 = Array.make_matrix nx ny 0.0 in
    let matrix2 = Array.make_matrix nx ny 0.0 in
    let sum   = ref 0  in
    let start = ref 0. in
    let upper_idx = Array.length matrix1 - 1 in
    for r = 1 to nt do
        for i = 1 to upper_idx - 1  do
		  let mat2  = matrix2.(i)   and
		      mat1m = matrix1.(i-1) and
			  mat10 = matrix1.(i)   and
			  mat1p = matrix1.(i+1) in
              for j = 1 to upper_idx - 1 do
                mat2.(j) <- mat10.(j) +. 1.0  +.
				                0.1*.(mat1p.(j)   +.
                                      mat1m.(j)   +.
                                      mat10.(j+1) +.
                                      mat10.(j-1) -.
                               4.0 *. mat10.(j));
            done;
        done;
        for i = 1 to upper_idx - 1  do
		  let mat2  = matrix1.(i)   and
		      mat1m = matrix2.(i-1) and
			  mat10 = matrix2.(i)   and
			  mat1p = matrix2.(i+1) in
              for j = 1 to upper_idx - 1 do
                mat2.(j) <- mat10.(j) +. 1.0  +.
				                0.1*.(mat1p.(j)   +.
                                      mat1m.(j)   +.
                                      mat10.(j+1) +.
                                      mat10.(j-1) -.
                               4.0 *. mat10.(j));
            done;
        done;
	Printf.printf "%f\n" (matrix1.(nx / 2).(ny / 2));
	flush stdout;
    done;
;;

let main3 nx ny nt =
  let step nmax m1p m1m m10 m2 =
     for j = 1 to nmax - 1 do
          m2.(j) <- m10.(j) +. 1.0  +.
				   0.1*.(m1p.(j)   +.
                         m1m.(j)   +.
                         m10.(j+1) +.
                         m10.(j-1) -.
                  4.0 *. m10.(j));
      done;
	in
    let matrix1 = Array.make_matrix nx ny 0.0 in
    let matrix2 = Array.make_matrix nx ny 0.0 in
    let sum   = ref 0  in
    let start = ref 0. in
    let upper_idx = Array.length matrix1 - 1 in
    for r = 1 to nt do
        for i = 1 to upper_idx - 1  do
		  let mat2  = matrix2.(i)   and
		      mat1m = matrix1.(i-1) and
			  mat10 = matrix1.(i)   and
			  mat1p = matrix1.(i+1) in
		  step upper_idx mat1p mat1m mat10 mat2;
        done;
        for i = 1 to upper_idx - 1  do
		  let mat2  = matrix1.(i)   and
		      mat1m = matrix2.(i-1) and
			  mat10 = matrix2.(i)   and
			  mat1p = matrix2.(i+1) in
		  step upper_idx mat1p mat1m mat10 mat2;
        done;
	Printf.printf "%f\n" (matrix1.(nx / 2).(ny / 2));
	flush stdout;
    done;
;;


let uset = Array.unsafe_set;;
let uget = Array.unsafe_get;;

let main4 nx ny nt =
  let step nmax m1p m1m m10 m2 =
     for j = 1 to nmax - 1 do
          uset m2 j  ((uget m10 j) +. 1.0  +.
				             0.1*.((uget m1p j    ) +.
                                   (uget m1m j    ) +.
                                   (uget m10 (j+1)) +.
                                   (uget m10 (j-1)) -.
                            4.0 *. (uget m10 j    )    ));
      done;
	in
    let matrix1 = Array.make_matrix nx ny 0.0 in
    let matrix2 = Array.make_matrix nx ny 0.0 in
    let sum   = ref 0  in
    let start = ref 0. in
    let upper_idx = Array.length matrix1 - 1 in
    for r = 1 to nt do
        for i = 1 to upper_idx - 1  do
		  let mat2  = matrix2.(i)   and
		      mat1m = matrix1.(i-1) and
			  mat10 = matrix1.(i)   and
			  mat1p = matrix1.(i+1) in
		  step upper_idx mat1p mat1m mat10 mat2;
        done;
        for i = 1 to upper_idx - 1  do
		  let mat2  = matrix1.(i)   and
		      mat1m = matrix2.(i-1) and
			  mat10 = matrix2.(i)   and
			  mat1p = matrix2.(i+1) in
		  step upper_idx mat1p mat1m mat10 mat2;
        done;
	Printf.printf "%f\n" (matrix1.(nx / 2).(ny / 2));
	flush stdout;
    done;
;;

let _ = main3 300 300 1000;;


-------------------------------------
Rolf Wester
wester@ilt.fhg.de