diff --git a/Utils/Davidson.ml b/Utils/Davidson.ml index bfca313..1450fcf 100644 --- a/Utils/Davidson.ml +++ b/Utils/Davidson.ml @@ -6,7 +6,7 @@ let make ?guess ?(n_states=8) ?(n_iter=10) - ?(threshold=1.e-8) + ?(threshold=1.e-6) diagonal matrix_prod = @@ -24,16 +24,16 @@ let make let random_vectors = let random_vector k = Vec.init n (fun i -> - let r1 = Random.float 1. - and r2 = Random.float 1. - in - let a = sqrt (-2. *. log r1) - and b = Constants.two_pi *. r2 - in if ik then - a *. cos b - else 100.0 + else if i=k then 1.e5 + else + let r1 = Random.float 1. + and r2 = Random.float 1. + in + let a = sqrt (-2. *. log r1) + and b = Constants.two_pi *. r2 in + let c = a *. cos b in + if abs_float c > 1.e-1 then c else 0. ) |> Util.normalize in @@ -73,7 +73,7 @@ let make matrix_prod ( u_new_ortho |> Mat.of_col_vecs_list - |> Matrix.dense_of_mat ) + |> Matrix.sparse_of_mat ) |> Matrix.to_mat |> Mat.to_col_vecs_list in diff --git a/Utils/Vector.ml b/Utils/Vector.ml index 2c5421e..0eacd52 100644 --- a/Utils/Vector.ml +++ b/Utils/Vector.ml @@ -2,10 +2,16 @@ open Lacaml.D let epsilon = Constants.epsilon +type index_value = + { + index: int; + value: float + } + type sparse_vector = { n: int; - v: (int*float) list + v: index_value list } type t = @@ -23,13 +29,22 @@ let is_dense = function | Dense _ -> true +exception Found of float + let get = function | Dense v -> (fun i -> v.{i}) | Sparse { n ; v } -> (fun i -> if i < 1 || i > n then invalid_arg "index out of bounds"; - match List.assoc_opt i v with - | Some x -> x - | None -> 0. ) + try + List.iter (fun {index ; value} -> + if index=i then + raise (Found value)) v; + raise Not_found + with + | Not_found -> 0. + | Found x -> x + ) + let dim = function @@ -47,14 +62,14 @@ let sparse_of_dense ?(threshold=epsilon) = function if abs_float x < threshold then aux accu (i-1) else - aux ((i, x)::accu) (i-1) + aux ({index=i ; value=x}::accu) (i-1) in let n = Vec.dim v in Sparse { n ; v=aux [] n } let rec to_assoc_list ?(threshold=epsilon) = function - | Sparse {n ; v} -> v + | Sparse {n ; v} -> List.map (fun {index ; value} -> (index, value)) v | Dense v -> to_assoc_list @@ sparse_of_dense ~threshold (Dense v) @@ -62,7 +77,7 @@ let dense_of_sparse = function | Dense _ -> invalid_arg "Expected a sparse vector" | Sparse {n ; v} -> let v' = Vec.make0 n in - List.iter (fun (i, x) -> v'.{i} <- x) v; + List.iter (fun {index ; value} -> v'.{index} <- value) v; Dense v' @@ -74,7 +89,10 @@ let sparse_of_vec ?(threshold=epsilon) v = |> sparse_of_dense ~threshold -let sparse_of_assoc_list n v = Sparse { n ; v } +let sparse_of_assoc_list n v = + Sparse { n ; + v = List.map (fun (index, value) -> {index ; value}) v + } let rec to_vec = function @@ -85,67 +103,22 @@ let rec to_vec = function let scale ?(threshold=epsilon) x = function | Dense v -> let v' = copy v in (scal x v'; Dense v') - | Sparse {n ; v} -> Sparse {n ; v=List.map (fun (i,y) -> let z = x *. y in - if abs_float z > threshold then Some (i, z) else None ) v |> Util.list_some } + | Sparse {n ; v} -> + Sparse {n ; v=List.map (fun {index ; value} -> + let z = x *. value in + if abs_float z > threshold then + Some {index ; value=z} + else + None + ) v |> Util.list_some } let rec neg = function | Dense v -> Dense (Vec.neg v) - | Sparse {n ; v} -> Sparse {n ; v=List.map (fun (i,y) -> (i, -. y)) v} + | Sparse {n ; v} -> + Sparse {n ; v=List.map (fun {index ; value} -> {index ; value = -. value}) v} -(* -let rec add ?(threshold=epsilon) x y = - if dim x <> dim y then - invalid_arg "Inconsistent dimensions"; - - match x, y with - | Dense x , Dense y -> Dense (Vec.add x y) - | Sparse {n ; v}, Dense y -> - let v' = copy y in - List.iter (fun (i, x) -> v'.{i} <- v'.{i} +. x) v; - sparse_of_vec ~threshold v' - | Sparse {n ; v}, Sparse {n=n' ; v=v'} -> - begin - let rec aux accu v1 v2 = - match v1, v2 with - | [], [] -> {n ; v=List.rev accu} - | ((i, x)::v1), [] -> - aux ((i, x)::accu) v1 [] - | [], ((j, y)::v2) -> - aux ((j, y)::accu) [] v2 - | ((i, x)::v1), ((j, y)::v2) -> - if i = j then - begin - let z = x +. y in - if abs_float z > threshold then - aux ((i, (x +. y))::accu) v1 v2 - else - aux accu v1 v2 - end - else if i < j then - begin - if abs_float x > threshold then - aux ((i, x)::accu) v1 ((j, y)::v2) - else - aux accu v1 ((j, y)::v2) - end - else - begin - if abs_float y > threshold then - aux ((j, y)::accu) ((i, x)::v1) v2 - else - aux accu ((i, x)::v1) v2 - end - in - Sparse (aux [] v v') - end - | x, y -> add ~threshold y x - - -let sub ?(threshold=epsilon) x y = add ~threshold x (neg y) -*) - let axpy ?(threshold=epsilon) ?(alpha=1.) x y = if dim x <> dim y then @@ -156,47 +129,51 @@ let axpy ?(threshold=epsilon) ?(alpha=1.) x y = | Sparse {n ; v}, Dense y -> begin let v' = copy y in - List.iter (fun (i, x) -> v'.{i} <- v'.{i} +. alpha *. x) v; + List.iter (fun {index ; value} -> v'.{index} <- v'.{index} +. alpha *. value) v; sparse_of_vec ~threshold v' end | Dense x , Sparse {n ; v} -> begin let v' = copy x in scal alpha v'; - List.iter (fun (i, y) -> v'.{i} <- v'.{i} +. y) v; + List.iter (fun {index ; value} -> v'.{index} <- v'.{index} +. value) v; sparse_of_vec ~threshold v' end | Sparse {n ; v}, Sparse {n=n' ; v=v'} -> begin let rec aux accu v1 v2 = match v1, v2 with - | [] , [] -> {n ; v=List.rev accu} - | ((i, x)::v1), [] -> aux ((i, x)::accu) v1 [] - | [] , ((j, y)::v2) -> aux ((j, y)::accu) [] v2 - | ((i, x)::v1), ((j, y)::v2) -> - if i = j then - begin - let z = alpha *. x +. y in - if abs_float z > threshold then - aux ((i, z)::accu) v1 v2 - else - aux accu v1 v2 - end - else if i < j then + | ({index=i ; value=x}::r1), ({index=j ; value=y}::r2) -> + begin + match compare i j with + | -1 -> let z = alpha *. x in - begin + let new_accu = if abs_float z > threshold then - aux ((i, z)::accu) v1 ((j, y)::v2) + {index=i ; value=z} :: accu else - aux accu v1 ((j, y)::v2) - end - else - begin + accu + in aux new_accu r1 v2 + | 1 -> + let new_accu = if abs_float y > threshold then - aux ((j, y)::accu) ((i, x)::v1) v2 + {index=j ; value=y} :: accu else - aux accu ((i, x)::v1) v2 - end + accu + in aux new_accu v1 r2 + | 0 -> + let z = alpha *. x +. y in + let new_accu = + if abs_float z > threshold then + {index=i ; value=z} :: accu + else + accu + in aux new_accu r1 r2 + | _ -> assert false + end + | ({index=i ; value=x}::r1), [] -> aux ({index=i ; value=x}::accu) r1 [] + | [] , ({index=j ; value=y}::r2) -> aux ({index=j ; value=y}::accu) [] r2 + | [] , [] -> {n ; v=List.rev accu} in Sparse (aux [] v v') end @@ -210,7 +187,7 @@ let pp_vector ppf = function | Sparse {n ; v} -> begin Format.fprintf ppf "@[[ %d | " n; - List.iter (fun (i,x) -> Format.fprintf ppf "@[(%d, %f); @]" i x) v; + List.iter (fun {index ; value} -> Format.fprintf ppf "@[(%d, %f); @]" index value) v; Format.fprintf ppf "]@]" end @@ -225,18 +202,26 @@ let dot v v' = let d_sp v' {n ; v} = if n <> Vec.dim v' then invalid_arg "Inconsistent dimensions"; - List.fold_left (fun accu (i, v_i) -> accu +. v_i *. v'.{i}) 0. v + List.fold_left (fun accu {index ; value} -> accu +. value *. v'.{index}) 0. v in let sp_sp {n ; v} {n=n' ; v=v'} = if n <> n' then invalid_arg "Inconsistent dimensions"; - List.fold_left (fun accu (i, v_i) -> - match List.assoc_opt i v' with - | Some w_i -> accu +. v_i *. w_i - | None -> accu - ) 0. v + let rec aux accu = function + | (({index=i ; value=v1} :: r1) as s1), (({index=j ; value=v2}::r2) as s2)-> + begin + match compare i j with + | -1 -> aux accu (r1, s2) + | 1 -> aux accu (s1, r2) + | 0 -> aux (accu +. v1 *. v2) (r1, r2) + | _ -> assert false + end + | ([], _ ) + | (_ , []) -> accu + in + aux 0. (v, v') in match v, v' with diff --git a/Utils/Vector.mli b/Utils/Vector.mli index f7b99c4..7a7b40a 100644 --- a/Utils/Vector.mli +++ b/Utils/Vector.mli @@ -60,7 +60,7 @@ val sub : ?threshold:float -> t -> t -> t (** Subtract two vectors *) val axpy : ?threshold:float -> ?alpha:float -> t -> t -> t -(** $a \mathfb{x} + \mathfb{y}$ *) +(** {% $a \mathbf{x} + \mathbf{y}$ %} *) val dot : t -> t -> float (** Dot product. *)