* Funny effect computation time
@ 2005-12-18 14:17 Andries Hekstra
2005-12-18 15:04 ` [Caml-list] " malc
0 siblings, 1 reply; 2+ messages in thread
From: Andries Hekstra @ 2005-12-18 14:17 UTC (permalink / raw)
To: caml-list
[-- Attachment #1.1: Type: text/plain, Size: 800 bytes --]
Dear OCaml-list,
I would like to thank several people of the ocaml list for constructive
comments to the LDPC decoder program I posted earlier. I attached the
revised Main.ml to this email.
I noticed that when compiling Main.ml with -inline 100 -noassert -unsafe
on a Pentium M, and I removing the unused argument messageCheckNodeOut of
the function processCheckNode (resulting in Main2.ml), the computation
time increases ca. 1.5-fold. I wonder whether this has to do with
(un)boxing?
=Andries
------------------------------------------------------------------------
Dr. Ir. Andries P. Hekstra
Philips Research
High Tech Campus 37
5656 AG Eindhoven
Tel./Fax/Secr. +31 40 27 42048/42566/44051
If you need anti-RSI break software : Look at http://www.workrave.org
(good and for free)
[-- Attachment #1.2: Type: text/html, Size: 1056 bytes --]
[-- Attachment #2: Main.ml --]
[-- Type: application/octet-stream, Size: 13567 bytes --]
open Printf (** Allows direct access to Printf functions *)
(* Parameters of the LDPC code, left and right degree*)
let ldpc_j = 3
let ldpc_k = 30
let codeRate = 1. -. float_of_int ldpc_j /. float_of_int ldpc_k
let targetCodewordLength = (ldpc_k * 300)
let blockSize = targetCodewordLength / ldpc_k
let codewordLength = blockSize * ldpc_k
let numberOfParityChecks = blockSize * ldpc_j
let fudgeFactor = 0.75
let maxNumberOfIterations = 30
(** With « let () = » you make sure your code returns a unit type --
thus you did not forgot an argument. Moreover, the variables
declared with a "let" binding are clearly local to the computation
instead of being global (see (+)). Printf is useful! *)
let () =
printf "LDPC decoder for rate 1 - %i/%i Gallager code with length %i\n%!"
ldpc_j ldpc_k codewordLength
(* Construction of a sparse parity check matrix*)
type vec = int array
type matrix = vec array (** An easy way to apply Gerd Stolpmann "trick" *)
(** Imperative style may be more readable and slightly faster (?) *)
let applyToMatrix f (a:matrix) =
let w = Array.length a.(0) - 1 in
for i = 0 to Array.length a - 1 do
let ai = a.(i) in
for j = 0 to w do
ai.(j) <- f i j ai.(j)
done
done
let order = Array.create_matrix ldpc_j codewordLength 0
let permutation = Array.create_matrix ldpc_j codewordLength 0
let () =
applyToMatrix (fun _ _ _ -> (Random.int (100 * codewordLength))) order;
applyToMatrix (fun _ i _ -> i) permutation
let swapElement i j (a:vec) = let x = a.(i) in a.(i) <- a.(j); a.(j) <- x
let qsort op (a:vec) =
let l = Array.length a in
let rec qsortAux i j op a =
if (i < j) then (
let x = a.((i + j) / 2) in
let i1 = ref i in
let j1 = ref j in
while (!i1 <= !j1) do
while (op a.(!i1) x) do incr i1 done;
while (op x a.(!j1)) do decr j1 done;
if (i1 <= j1) then (
swapElement !i1 !j1 a;
incr i1;
decr j1
)
done;
qsortAux !i1 j op a;
qsortAux i !j1 op a;
)
in qsortAux 0 (l - 1) op a;;
let qsort2 op (a:vec) (b:vec) =
let l = Array.length a in
let rec qsort2Aux i j op (a:vec) (b:vec) =
if (i < j) then (
let x = a.((i + j) / 2) in
let i1 = ref i in
let j1 = ref j in
while (!i1 <= !j1) do
while (op a.(!i1) x) do incr i1 done;
while (op x a.(!j1)) do decr j1 done;
if (i1 <= j1) then (
swapElement !i1 !j1 a;
swapElement !i1 !j1 b;
incr i1;
decr j1
)
done;
qsort2Aux !i1 j op a b;
qsort2Aux i !j1 op a b;
)
in qsort2Aux 0 (l - 1) op a b;;
let sortRows a =
for i = 0 to Array.length a - 1 do qsort (<) a.(i) done
let sortRows2 a b =
for i = 0 to Array.length a - 1 do qsort2 (<) a.(i) b.(i) done
(** This is a matter of taste but this is mostly imperative, thus a
for loop is appropriate (still it is good to have the choice).
For float arrays (not the case here), for loops are slightly
better as references to floats are "unboxed" (not as arguments of
a recursive function). *)
let oneColumnInRow = Array.create_matrix numberOfParityChecks ldpc_k 0
let jOneColumnInRow = Array.create_matrix numberOfParityChecks ldpc_k 0
let () =
sortRows2 order permutation;
(** (+) The following variable is local to the current computation.
This is clear (and devoid of error) now that it is inside a
« let () = ... ». *)
let inversePermutation = Array.create_matrix ldpc_j codewordLength 0 in
for j = 0 to ldpc_j - 1 do
for i = 0 to codewordLength - 1 do
inversePermutation.(j).(permutation.(j).(i)) <- i
done
done;
(* Here, i indexes parity checks *)
applyToMatrix (fun i k x ->
(** Local variable to help readability *)
let cwl = k + ldpc_k * (i mod blockSize) in
inversePermutation.(i / blockSize).(cwl)) oneColumnInRow;
applyToMatrix (fun i _ _ -> i / blockSize) jOneColumnInRow;
sortRows2 oneColumnInRow jOneColumnInRow
let () =
let oneRowInColumn = Array.create_matrix codewordLength ldpc_j 0
and kOneRowInColumn = Array.create_matrix codewordLength ldpc_j 0 in
for j = 0 to ldpc_j-1 do
for r = 0 to blockSize -1 do
let i = r + j * blockSize in
for k = 0 to ldpc_k-1 do
oneRowInColumn .(oneColumnInRow.(i).(k)).(j) <- i;
kOneRowInColumn.(oneColumnInRow.(i).(k)).(j) <- k;
done
done
done;
sortRows2 oneRowInColumn kOneRowInColumn;
let numberOfColumnPairsWithOverlap = ref 0 in
for i1 = 0 to codewordLength - 1 do
let oi1 = oneRowInColumn.(i1) in
for i2 = 0 to i1-1 do
let oi2 = oneRowInColumn.(i2) in
let overlap = ref 0 in
for j1 = 0 to ldpc_j-1 do
let oi1j1 = oi1.(j1) in
for j2 = 0 to ldpc_j-1 do
if oi1j1 = oi2.(j2) then
incr overlap;
done
done;
if !overlap > 1 then incr numberOfColumnPairsWithOverlap;
done
done;
printf "numberOfColumnPairsWithOverlap = %i\n%!"
!numberOfColumnPairsWithOverlap;;
(* Random bit vector *) (* First I used 'app)' of array : App)ing to an array is very inefficient *)
let bitArrayFromInt n iii =
let iii = ref iii in
Array.init n (fun _ -> let v = !iii mod 2 in iii := !iii / 2; v)
let makeRandomBit a =
let i = ref(Array.length a - 1) in
while !i > 0 do
let b = bitArrayFromInt 8 (Random.int 256) in
let count = ref 0 in
while !i >= 0 && !count < 8 do
a.(!i) <- b.(!count);
decr i;
incr count;
done;
done
(* Normal distribution *)
let pi = 2. *. asin 1.
let randomGaussianInStore = ref false
let storedRandomGaussian = ref 0.
let randomGaussian sigma =
if !randomGaussianInStore then (
randomGaussianInStore := false;
(sigma *. !storedRandomGaussian)
) else
let u1 = Random.float 1. in
let u2 = Random.float 1. in
(* log(.) is a slow computation *)
let r = sqrt (-2. *. log (u1)) in
let phi = 2. *. pi *. u2 in (
randomGaussianInStore := true;
storedRandomGaussian := r *. sin (phi);
sigma *. r *. cos (phi)
)
(* Message passing functions *)
let sum a =
let r = ref 0. in
for i = 0 to Array.length a - 1 do
r := !r +. a.(i)
done;
!r
let hardBit softBit = if (softBit >= 0.) then 0 else 1
let reliability softBit = if (softBit >= 0.) then softBit else (-. softBit)
(* LDPC decode *)
let softBitOut channelLLR1 messageSymbolNodeIn1 =
channelLLR1 +. sum messageSymbolNodeIn1
let hardBitOut channelLLR1 messageSymbolNodeIn1 =
hardBit (softBitOut channelLLR1 messageSymbolNodeIn1)
let processCheckNode channelLLR
oneColumnInRow
jOneColumnInRow
syndrome
messagesSymbolNodeIn
messageCheckNodeIn
messageCheckNodeOut =
let messageSymbolNodeOut channelLLR1 messageSymbolNodeIn1 j =
(softBitOut channelLLR1 messageSymbolNodeIn1) -. messageSymbolNodeIn1.(j) in
(* Get the input messages to the current check node *)
for k = 0 to ldpc_k-1 do
messageCheckNodeIn.(k) <- fudgeFactor *.
(messageSymbolNodeOut channelLLR.(oneColumnInRow.(k))
messagesSymbolNodeIn.(oneColumnInRow.(k))
jOneColumnInRow.(k))
done;
(* Compute the total EXOR, the min. and second min. reliability and
the index *)
let k1 = ref 0 in
let minReliability1 = ref 1e10 in
let minReliability2 = ref 1e10 in
let totalExor = ref syndrome in (** ( : unnecessary and cause
extra identation *)
for k = 0 to ldpc_k-1 do
totalExor := !totalExor + hardBit messageCheckNodeIn.(k);
let h = reliability messageCheckNodeIn.(k) in
if !minReliability1 >= h then (
minReliability2 := !minReliability1;
minReliability1 := h;
k1 := k;
) else (
if !minReliability2 >= h then minReliability2 := h
)
done;
totalExor := !totalExor mod 2;
(* Compute the output messages *)
for k = 0 to ldpc_k-1 do
let h = if k = !k1 then !minReliability2 else !minReliability1 in
let ok = oneColumnInRow.(k) in
let jk = jOneColumnInRow.(k) in
if ((!totalExor + (hardBit messageCheckNodeIn.(k))) mod 2) = 0 then
messagesSymbolNodeIn.(ok).(jk) <- h
else
messagesSymbolNodeIn.(ok).(jk) <- (-. h)
done;
!totalExor = 0
let ldpcDecode transmittedCodeword
syndrome
channelLLR
decodedCodeword
oneColumnInRow
jOneColumnInRow
messagesSymbolNodeIn
numberOfIterations =
let messageCheckNodeIn = Array.create ldpc_k 0. in
let messageCheckNodeOut = Array.create ldpc_k 0. in
let allZeroesSyndrome = ref false in
numberOfIterations := 0;
while (!numberOfIterations < maxNumberOfIterations)
&& not !allZeroesSyndrome do
allZeroesSyndrome := true;
for row = 0 to numberOfParityChecks - 1 do
if not (processCheckNode channelLLR
oneColumnInRow.(row)
jOneColumnInRow.(row)
syndrome.(row)
messagesSymbolNodeIn
messageCheckNodeIn
messageCheckNodeOut) then
allZeroesSyndrome := false;
done;
incr numberOfIterations;
done;
let numberOfBitErrors = ref 0 in
for column = 0 to codewordLength - 1 do
decodedCodeword.(column)
<- hardBitOut channelLLR.(column) messagesSymbolNodeIn.(column);
if (decodedCodeword.(column) <> transmittedCodeword.(column)) then
incr numberOfBitErrors;
done;
!numberOfBitErrors
let () =
(* Simulation parameters *)
let ebOverNoDbLow = 3.4 in
let ebOverNoDbHigh = 3.4 in
let ebOverNoDbStep = 0.1 in
let numberOfCodewordsToSimulate = 30000 in
let ebOverNoDb = ref ebOverNoDbLow in
let totalNumberOfFrameErrors = ref 1 in
while !ebOverNoDb <= ebOverNoDbHigh && !totalNumberOfFrameErrors > 0 do
printf "\n\nEbOverNoDb = %g, time = %g seconds!\n%!"
!ebOverNoDb (Sys.time());
let ebOverNo = 10. ** (!ebOverNoDb /. 10.) in
let sigmaNoise = 1. /. sqrt (2. *. ebOverNo *. codeRate) in
let totalNumberOfChannelBitErrors = ref 0. in
let totalNumberOfBits = ref 0. in
let totalNumberOfDecodedBitErrors = ref 0. in
let totalNumberOfFrames = ref 0 in
let totalNumberOfIterations = ref 0. in
let transmittedCodeword = Array.create codewordLength 0 in
let decodedCodeword = Array.create codewordLength 0 in
let syndrome = Array.create numberOfParityChecks 0 in
let channelLLR = Array.create codewordLength 0. in
let messagesSymbolNodeIn = Array.create_matrix codewordLength ldpc_j 0. in
totalNumberOfFrameErrors := 0;
while (!totalNumberOfFrames < numberOfCodewordsToSimulate)
&& ((!totalNumberOfDecodedBitErrors < 500.)
|| (!totalNumberOfFrameErrors < 200)) do
makeRandomBit transmittedCodeword;
for r = 0 to numberOfParityChecks-1 do
syndrome.(r) <- 0;
for k = 0 to ldpc_k-1 do
syndrome.(r)
<- syndrome.(r)+transmittedCodeword.(oneColumnInRow.(r).(k))
done;
syndrome.(r) <- syndrome.(r) mod 2
done;
for i = 0 to codewordLength - 1 do
channelLLR.(i) <- float_of_int (1 - 2*transmittedCodeword.(i))
+. randomGaussian sigmaNoise
done;
let numberOfChannelBitErrors = ref 0 in (
for i = 0 to codewordLength-1 do
if transmittedCodeword.(i) <> hardBit channelLLR.(i) then
incr numberOfChannelBitErrors
done;
totalNumberOfChannelBitErrors := !totalNumberOfChannelBitErrors
+. float !numberOfChannelBitErrors;
totalNumberOfBits := !totalNumberOfBits +. float codewordLength;
);
let numberOfIterations = ref 0 in
let numberOfDecodedBitErrors = ldpcDecode transmittedCodeword
syndrome
channelLLR
decodedCodeword
oneColumnInRow
jOneColumnInRow
messagesSymbolNodeIn
numberOfIterations in (
totalNumberOfDecodedBitErrors := !totalNumberOfDecodedBitErrors
+. float numberOfDecodedBitErrors;
totalNumberOfIterations := !totalNumberOfIterations
+. float !numberOfIterations;
if numberOfDecodedBitErrors > 0 then incr totalNumberOfFrameErrors;
incr totalNumberOfFrames;
);
if (!totalNumberOfFrames mod 100) = 0 then (
let avgChannelBER =
!totalNumberOfChannelBitErrors /. !totalNumberOfBits in
let avgDecodedBER =
!totalNumberOfDecodedBitErrors /. !totalNumberOfBits in
let avgFER =
float !totalNumberOfFrameErrors /. float !totalNumberOfFrames in
let avgIter =
!totalNumberOfIterations /. float !totalNumberOfFrames in
printf
"%i\tChannelBER = %f\tDecodedBER = %f\tFER = %f\tavgIter = %f\n%!"
!totalNumberOfFrames avgChannelBER avgDecodedBER avgFER avgIter;
)
done;
ebOverNoDb := !ebOverNoDb +. ebOverNoDbStep;
done;
printf "\nSuccesful completion in %f seconds!\n" (Sys.time())
[-- Attachment #3: Main2.ml --]
[-- Type: application/octet-stream, Size: 13448 bytes --]
open Printf (** Allows direct access to Printf functions *)
(* Parameters of the LDPC code, left and right degree*)
let ldpc_j = 3
let ldpc_k = 30
let codeRate = 1. -. float_of_int ldpc_j /. float_of_int ldpc_k
let targetCodewordLength = (ldpc_k * 300)
let blockSize = targetCodewordLength / ldpc_k
let codewordLength = blockSize * ldpc_k
let numberOfParityChecks = blockSize * ldpc_j
let fudgeFactor = 0.75
let maxNumberOfIterations = 30
(** With « let () = » you make sure your code returns a unit type --
thus you did not forgot an argument. Moreover, the variables
declared with a "let" binding are clearly local to the computation
instead of being global (see (+)). Printf is useful! *)
let () =
printf "LDPC decoder for rate 1 - %i/%i Gallager code with length %i\n%!"
ldpc_j ldpc_k codewordLength
(* Construction of a sparse parity check matrix*)
type vec = int array
type matrix = vec array (** An easy way to apply Gerd Stolpmann "trick" *)
(** Imperative style may be more readable and slightly faster (?) *)
let applyToMatrix f (a:matrix) =
let w = Array.length a.(0) - 1 in
for i = 0 to Array.length a - 1 do
let ai = a.(i) in
for j = 0 to w do
ai.(j) <- f i j ai.(j)
done
done
let order = Array.create_matrix ldpc_j codewordLength 0
let permutation = Array.create_matrix ldpc_j codewordLength 0
let () =
applyToMatrix (fun _ _ _ -> (Random.int (100 * codewordLength))) order;
applyToMatrix (fun _ i _ -> i) permutation
let swapElement i j (a:vec) = let x = a.(i) in a.(i) <- a.(j); a.(j) <- x
let qsort op (a:vec) =
let l = Array.length a in
let rec qsortAux i j op a =
if (i < j) then (
let x = a.((i + j) / 2) in
let i1 = ref i in
let j1 = ref j in
while (!i1 <= !j1) do
while (op a.(!i1) x) do incr i1 done;
while (op x a.(!j1)) do decr j1 done;
if (i1 <= j1) then (
swapElement !i1 !j1 a;
incr i1;
decr j1
)
done;
qsortAux !i1 j op a;
qsortAux i !j1 op a;
)
in qsortAux 0 (l - 1) op a;;
let qsort2 op (a:vec) (b:vec) =
let l = Array.length a in
let rec qsort2Aux i j op (a:vec) (b:vec) =
if (i < j) then (
let x = a.((i + j) / 2) in
let i1 = ref i in
let j1 = ref j in
while (!i1 <= !j1) do
while (op a.(!i1) x) do incr i1 done;
while (op x a.(!j1)) do decr j1 done;
if (i1 <= j1) then (
swapElement !i1 !j1 a;
swapElement !i1 !j1 b;
incr i1;
decr j1
)
done;
qsort2Aux !i1 j op a b;
qsort2Aux i !j1 op a b;
)
in qsort2Aux 0 (l - 1) op a b;;
let sortRows a =
for i = 0 to Array.length a - 1 do qsort (<) a.(i) done
let sortRows2 a b =
for i = 0 to Array.length a - 1 do qsort2 (<) a.(i) b.(i) done
(** This is a matter of taste but this is mostly imperative, thus a
for loop is appropriate (still it is good to have the choice).
For float arrays (not the case here), for loops are slightly
better as references to floats are "unboxed" (not as arguments of
a recursive function). *)
let oneColumnInRow = Array.create_matrix numberOfParityChecks ldpc_k 0
let jOneColumnInRow = Array.create_matrix numberOfParityChecks ldpc_k 0
let () =
sortRows2 order permutation;
(** (+) The following variable is local to the current computation.
This is clear (and devoid of error) now that it is inside a
« let () = ... ». *)
let inversePermutation = Array.create_matrix ldpc_j codewordLength 0 in
for j = 0 to ldpc_j - 1 do
for i = 0 to codewordLength - 1 do
inversePermutation.(j).(permutation.(j).(i)) <- i
done
done;
(* Here, i indexes parity checks *)
applyToMatrix (fun i k x ->
(** Local variable to help readability *)
let cwl = k + ldpc_k * (i mod blockSize) in
inversePermutation.(i / blockSize).(cwl)) oneColumnInRow;
applyToMatrix (fun i _ _ -> i / blockSize) jOneColumnInRow;
sortRows2 oneColumnInRow jOneColumnInRow
let () =
let oneRowInColumn = Array.create_matrix codewordLength ldpc_j 0
and kOneRowInColumn = Array.create_matrix codewordLength ldpc_j 0 in
for j = 0 to ldpc_j-1 do
for r = 0 to blockSize -1 do
let i = r + j * blockSize in
for k = 0 to ldpc_k-1 do
oneRowInColumn .(oneColumnInRow.(i).(k)).(j) <- i;
kOneRowInColumn.(oneColumnInRow.(i).(k)).(j) <- k;
done
done
done;
sortRows2 oneRowInColumn kOneRowInColumn;
let numberOfColumnPairsWithOverlap = ref 0 in
for i1 = 0 to codewordLength - 1 do
let oi1 = oneRowInColumn.(i1) in
for i2 = 0 to i1-1 do
let oi2 = oneRowInColumn.(i2) in
let overlap = ref 0 in
for j1 = 0 to ldpc_j-1 do
let oi1j1 = oi1.(j1) in
for j2 = 0 to ldpc_j-1 do
if oi1j1 = oi2.(j2) then
incr overlap;
done
done;
if !overlap > 1 then incr numberOfColumnPairsWithOverlap;
done
done;
printf "numberOfColumnPairsWithOverlap = %i\n%!"
!numberOfColumnPairsWithOverlap;;
(* Random bit vector *) (* First I used 'app)' of array : App)ing to an array is very inefficient *)
let bitArrayFromInt n iii =
let iii = ref iii in
Array.init n (fun _ -> let v = !iii mod 2 in iii := !iii / 2; v)
let makeRandomBit a =
let i = ref(Array.length a - 1) in
while !i > 0 do
let b = bitArrayFromInt 8 (Random.int 256) in
let count = ref 0 in
while !i >= 0 && !count < 8 do
a.(!i) <- b.(!count);
decr i;
incr count;
done;
done
(* Normal distribution *)
let pi = 2. *. asin 1.
let randomGaussianInStore = ref false
let storedRandomGaussian = ref 0.
let randomGaussian sigma =
if !randomGaussianInStore then (
randomGaussianInStore := false;
(sigma *. !storedRandomGaussian)
) else
let u1 = Random.float 1. in
let u2 = Random.float 1. in
(* log(.) is a slow computation *)
let r = sqrt (-2. *. log (u1)) in
let phi = 2. *. pi *. u2 in (
randomGaussianInStore := true;
storedRandomGaussian := r *. sin (phi);
sigma *. r *. cos (phi)
)
(* Message passing functions *)
let sum a =
let r = ref 0. in
for i = 0 to Array.length a - 1 do
r := !r +. a.(i)
done;
!r
let hardBit softBit = if (softBit >= 0.) then 0 else 1
let reliability softBit = if (softBit >= 0.) then softBit else (-. softBit)
(* LDPC decode *)
let softBitOut channelLLR1 messageSymbolNodeIn1 =
channelLLR1 +. sum messageSymbolNodeIn1
let hardBitOut channelLLR1 messageSymbolNodeIn1 =
hardBit (softBitOut channelLLR1 messageSymbolNodeIn1)
let processCheckNode channelLLR
oneColumnInRow
jOneColumnInRow
syndrome
messagesSymbolNodeIn
messageCheckNodeIn =
let messageSymbolNodeOut channelLLR1 messageSymbolNodeIn1 j =
(softBitOut channelLLR1 messageSymbolNodeIn1) -. messageSymbolNodeIn1.(j) in
(* Get the input messages to the current check node *)
for k = 0 to ldpc_k-1 do
messageCheckNodeIn.(k) <- fudgeFactor *.
(messageSymbolNodeOut channelLLR.(oneColumnInRow.(k))
messagesSymbolNodeIn.(oneColumnInRow.(k))
jOneColumnInRow.(k))
done;
(* Compute the total EXOR, the min. and second min. reliability and
the index *)
let k1 = ref 0 in
let minReliability1 = ref 1e10 in
let minReliability2 = ref 1e10 in
let totalExor = ref syndrome in (** ( : unnecessary and cause
extra identation *)
for k = 0 to ldpc_k-1 do
totalExor := !totalExor + hardBit messageCheckNodeIn.(k);
let h = reliability messageCheckNodeIn.(k) in
if !minReliability1 >= h then (
minReliability2 := !minReliability1;
minReliability1 := h;
k1 := k;
) else (
if !minReliability2 >= h then minReliability2 := h
)
done;
totalExor := !totalExor mod 2;
(* Compute the output messages *)
for k = 0 to ldpc_k-1 do
let h = if k = !k1 then !minReliability2 else !minReliability1 in
let ok = oneColumnInRow.(k) in
let jk = jOneColumnInRow.(k) in
if ((!totalExor + (hardBit messageCheckNodeIn.(k))) mod 2) = 0 then
messagesSymbolNodeIn.(ok).(jk) <- h
else
messagesSymbolNodeIn.(ok).(jk) <- (-. h)
done;
!totalExor = 0
let ldpcDecode transmittedCodeword
syndrome
channelLLR
decodedCodeword
oneColumnInRow
jOneColumnInRow
messagesSymbolNodeIn
numberOfIterations =
let messageCheckNodeIn = Array.create ldpc_k 0. in
let allZeroesSyndrome = ref false in
numberOfIterations := 0;
while (!numberOfIterations < maxNumberOfIterations)
&& not !allZeroesSyndrome do
allZeroesSyndrome := true;
for row = 0 to numberOfParityChecks - 1 do
if not (processCheckNode channelLLR
oneColumnInRow.(row)
jOneColumnInRow.(row)
syndrome.(row)
messagesSymbolNodeIn
messageCheckNodeIn) then
allZeroesSyndrome := false;
done;
incr numberOfIterations;
done;
let numberOfBitErrors = ref 0 in
for column = 0 to codewordLength - 1 do
decodedCodeword.(column)
<- hardBitOut channelLLR.(column) messagesSymbolNodeIn.(column);
if (decodedCodeword.(column) <> transmittedCodeword.(column)) then
incr numberOfBitErrors;
done;
!numberOfBitErrors
let () =
(* Simulation parameters *)
let ebOverNoDbLow = 3.4 in
let ebOverNoDbHigh = 3.4 in
let ebOverNoDbStep = 0.1 in
let numberOfCodewordsToSimulate = 30000 in
let ebOverNoDb = ref ebOverNoDbLow in
let totalNumberOfFrameErrors = ref 1 in
while !ebOverNoDb <= ebOverNoDbHigh && !totalNumberOfFrameErrors > 0 do
printf "\n\nEbOverNoDb = %g, time = %g seconds!\n%!"
!ebOverNoDb (Sys.time());
let ebOverNo = 10. ** (!ebOverNoDb /. 10.) in
let sigmaNoise = 1. /. sqrt (2. *. ebOverNo *. codeRate) in
let totalNumberOfChannelBitErrors = ref 0. in
let totalNumberOfBits = ref 0. in
let totalNumberOfDecodedBitErrors = ref 0. in
let totalNumberOfFrames = ref 0 in
let totalNumberOfIterations = ref 0. in
let transmittedCodeword = Array.create codewordLength 0 in
let decodedCodeword = Array.create codewordLength 0 in
let syndrome = Array.create numberOfParityChecks 0 in
let channelLLR = Array.create codewordLength 0. in
let messagesSymbolNodeIn = Array.create_matrix codewordLength ldpc_j 0. in
totalNumberOfFrameErrors := 0;
while (!totalNumberOfFrames < numberOfCodewordsToSimulate)
&& ((!totalNumberOfDecodedBitErrors < 500.)
|| (!totalNumberOfFrameErrors < 200)) do
makeRandomBit transmittedCodeword;
for r = 0 to numberOfParityChecks-1 do
syndrome.(r) <- 0;
for k = 0 to ldpc_k-1 do
syndrome.(r)
<- syndrome.(r)+transmittedCodeword.(oneColumnInRow.(r).(k))
done;
syndrome.(r) <- syndrome.(r) mod 2
done;
for i = 0 to codewordLength - 1 do
channelLLR.(i) <- float_of_int (1 - 2*transmittedCodeword.(i))
+. randomGaussian sigmaNoise
done;
let numberOfChannelBitErrors = ref 0 in (
for i = 0 to codewordLength-1 do
if transmittedCodeword.(i) <> hardBit channelLLR.(i) then
incr numberOfChannelBitErrors
done;
totalNumberOfChannelBitErrors := !totalNumberOfChannelBitErrors
+. float !numberOfChannelBitErrors;
totalNumberOfBits := !totalNumberOfBits +. float codewordLength;
);
let numberOfIterations = ref 0 in
let numberOfDecodedBitErrors = ldpcDecode transmittedCodeword
syndrome
channelLLR
decodedCodeword
oneColumnInRow
jOneColumnInRow
messagesSymbolNodeIn
numberOfIterations in (
totalNumberOfDecodedBitErrors := !totalNumberOfDecodedBitErrors
+. float numberOfDecodedBitErrors;
totalNumberOfIterations := !totalNumberOfIterations
+. float !numberOfIterations;
if numberOfDecodedBitErrors > 0 then incr totalNumberOfFrameErrors;
incr totalNumberOfFrames;
);
if (!totalNumberOfFrames mod 100) = 0 then (
let avgChannelBER =
!totalNumberOfChannelBitErrors /. !totalNumberOfBits in
let avgDecodedBER =
!totalNumberOfDecodedBitErrors /. !totalNumberOfBits in
let avgFER =
float !totalNumberOfFrameErrors /. float !totalNumberOfFrames in
let avgIter =
!totalNumberOfIterations /. float !totalNumberOfFrames in
printf
"%i\tChannelBER = %f\tDecodedBER = %f\tFER = %f\tavgIter = %f\n%!"
!totalNumberOfFrames avgChannelBER avgDecodedBER avgFER avgIter;
)
done;
ebOverNoDb := !ebOverNoDb +. ebOverNoDbStep;
done;
printf "\nSuccesful completion in %f seconds!\n" (Sys.time())
^ permalink raw reply [flat|nested] 2+ messages in thread
* Re: [Caml-list] Funny effect computation time
2005-12-18 14:17 Funny effect computation time Andries Hekstra
@ 2005-12-18 15:04 ` malc
0 siblings, 0 replies; 2+ messages in thread
From: malc @ 2005-12-18 15:04 UTC (permalink / raw)
To: Andries Hekstra; +Cc: caml-list
On Sun, 18 Dec 2005, Andries Hekstra wrote:
> Dear OCaml-list,
>
> I would like to thank several people of the ocaml list for constructive
> comments to the LDPC decoder program I posted earlier. I attached the
> revised Main.ml to this email.
>
> I noticed that when compiling Main.ml with -inline 100 -noassert -unsafe
> on a Pentium M, and I removing the unused argument messageCheckNodeOut of
> the function processCheckNode (resulting in Main2.ml), the computation
> time increases ca. 1.5-fold. I wonder whether this has to do with
> (un)boxing?
It's hard to say why OCaml somtimes produces the code the way it does,
it's even harder to predict how it will be intrepreted by different
"modern" CPUs. Here's a list of potentially useful links.
http://groups.google.com/group/fa.caml/browse_frm/thread/9c152a39b30c88d5/080df2a83aa2d878?lnk=st&q=group%3Afa.caml+author%3Amalc+faster&rnum=3#080df2a83aa2d878
http://groups.google.com/group/fa.caml/browse_thread/thread/624ed8213174b634/f5371b023927cc8e?lnk=st&q=group%3Afa.caml+author%3Amalc&rnum=15#f5371b023927cc8e
http://groups.google.com/group/fa.caml/browse_frm/thread/c9e40e633120a79/5f14ce70f4a35373?tvc=1&q=group%3Afa.caml+author%3Amalc#5f14ce70f4a35373
http://groups.google.com/group/comp.lang.functional/browse_frm/thread/4a6159c5beef192a/689140a6436f558f?lnk=st&q=malc+stand+corrected+ocaml&rnum=1#689140a6436f558f
--
mailto:malc@pulsesoft.com
^ permalink raw reply [flat|nested] 2+ messages in thread
end of thread, other threads:[~2005-12-18 15:03 UTC | newest]
Thread overview: 2+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2005-12-18 14:17 Funny effect computation time Andries Hekstra
2005-12-18 15:04 ` [Caml-list] " malc
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for NNTP newsgroup(s).