Source file `vector_complex_flat.ml`

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
(* gsl-ocaml - OCaml interface to GSL                       *)
(* Copyright (©) 2002-2012 - Olivier Andrieu                *)
(* Distributed under the terms of the GPL version 3         *)

let () = Error.init ()

type complex_vector_flat = 
    { data   : float array ;
      off    : int ;
      len    : int ;
      stride : int ; }

type vector = complex_vector_flat

let create ?(init=Complex.zero) len = 
  let arr = { 
    data = Array.make (2*len) init.Complex.re ; 
    off = 0; 
    len = len; 
    stride = 1; } in
  if init.Complex.im <> init.Complex.re
  then for i=0 to pred len do
    arr.data.(2*i+1) <- init.Complex.im
  done ;
  arr

let of_array arr =
  let carr = Gsl_complex.pack arr in
  { data = carr; off = 0; 
    len = Array.length arr; stride = 1; }

let length { len = len } = 
  len 

let get v i = 
  Gsl_complex.get v.data (v.off + i*v.stride)

let set v i d =
  Gsl_complex.set v.data (v.off + i*v.stride) d

let set_all v d = 
  for i=0 to pred v.len do
    set v i d
  done

let set_zero v = 
  set_all v Complex.zero

let set_basis v i = 
  set_zero v ;
  set v i Complex.one

let to_array v =
  Array.init v.len (get v)

let of_complex_array carr =
  { data = Array.copy carr; off = 0; 
    len = (Array.length carr)/2; stride = 1; }

let to_complex_array arr =
  let carr = Array.make (2*arr.len) 0. in
  for i=0 to pred arr.len do
    Gsl_complex.set carr i (get arr i)
  done ;
  carr

let real carr =
  Vector_flat.view_array
    ~stride:(2 * carr.stride)
    ~off:(2 * carr.off)
    ~len:carr.len
    carr.data

let imag carr =
  Vector_flat.view_array
    ~stride:(2 * carr.stride)
    ~off:(2 * carr.off + 1)
    ~len:carr.len
    carr.data

let subvector ?(stride=1) v ~off ~len =
  { v with 
    off = off * v.stride + v.off ;
    len = len ;
    stride = stride * v.stride ; }

let view_complex_array ?(stride=1) ?(off=0) ?len arr =
  let alen = Array.length arr in
  if alen mod 2 <> 0
  then invalid_arg "complex_array dim" ;
  let len = match len with
  | None -> alen / 2
  | Some l -> l in
  { data = arr ; off = off ;
    stride = stride ; len = len }

let memcpy v w = 
  if v.len <> w.len
  then invalid_arg "Vector.memcpy" ;
  for i=0 to pred v.len do
    set w i (get v i)
  done

let copy v = 
  { v with data = Array.copy v.data }

let swap_element v i j =
  let d  = get v i in
  let d' = get v j in
  set v j d ;
  set v i d'

let reverse v = 
  for i=0 to pred (v.len/2) do
    swap_element v i (pred v.len - i)
  done