diff --git a/CI/CI.ml b/CI/CI.ml index 7752661..2fae19a 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -107,7 +107,7 @@ let make det_space = |> Farm.run ~ordered:false ~f:(fun i -> let ki = det.(i) in - Printf.eprintf "%d / %d\n%!" i ndet; + Printf.eprintf "%8d / %8d\r%!" i ndet; List.init ndet (fun j -> let kj = det.(j) in let x = h_ij mo_basis ki kj in @@ -130,9 +130,15 @@ let make det_space = in let eigensystem = lazy ( let m_H = Lazy.force m_H in + (* Parallel.broadcast @@ lazy (Util.diagonalize_symm m_H) + *) + let diagonal = Vec.init (Mat.dim1 m_H) (fun i -> m_H.{i,i}) in + let matrix_vector psi = symv m_H psi in + Davidson.make diagonal matrix_vector ) in { det_space ; m_H ; m_S2 ; eigensystem } + diff --git a/Utils/Constants.ml b/Utils/Constants.ml index a468873..a25723d 100644 --- a/Utils/Constants.ml +++ b/Utils/Constants.ml @@ -3,6 +3,7 @@ let epsilon = 1.e-20 (** Constants *) let pi = acos (-1.) +let two_pi = 2. *. pi let sq_pi = sqrt pi let sq_pi_over_two = sq_pi *. 0.5 let pi_inv = 1. /. pi diff --git a/Utils/Constants.mli b/Utils/Constants.mli index eff7450..5997ef3 100644 --- a/Utils/Constants.mli +++ b/Utils/Constants.mli @@ -14,6 +14,9 @@ val integrals_cutoff : float val pi : float (** {% $\pi$ %} = 3.141_592_653_589_793_12 *) +val two_pi : float +(** {% $2\pi$ %} *) + val sq_pi : float (** {% $\sqrt{\pi}$ %} *) diff --git a/Utils/Davidson.ml b/Utils/Davidson.ml new file mode 100644 index 0000000..0decc04 --- /dev/null +++ b/Utils/Davidson.ml @@ -0,0 +1,133 @@ +open Lacaml.D + +type t + +let make + ?guess + ?(n_states=8) + ?(n_iter=10) + ?(threshold=1.e-10) + diagonal + matrix_vector + = + + let n = Vec.dim diagonal in (* Size of the matrix to diagonalize *) + + + let m = (* Number of requested states *) + match guess with + | Some vectors -> (Mat.dim2 vectors) * n_states + | None -> n_states + in + + (* Create guess vectors u, with randomly initialized unknown vectors. *) + 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 + ) + |> Util.normalize + in + List.init m (fun i -> random_vector i) + in + + let pick_new u = + Mat.to_col_vecs_list u + |> Util.list_pack m + |> List.rev + |> List.hd + in + + let u_new = + match guess with + | Some vectors -> Mat.to_col_vecs_list vectors + | None -> random_vectors + in + + let rec iteration u u_new w iter = + (* u is a list of orthonormal vectors, on which the operator has + been applied : w = op.u + u_new is a list of vector which will increase the size of the + space. + *) + + (* Orthonormalize input vectors u_new *) + let u_new_ortho = + List.concat [u ; u_new] + |> Mat.of_col_vecs_list + |> Util.qr_ortho + |> pick_new + in + + (* Apply the operator the m last vectors *) + let w_new = + List.map matrix_vector u_new_ortho + in + + (* Data for the next iteration *) + let u_next = + List.concat [ u ; u_new_ortho ] + and w_next = + List.concat [ w ; w_new ] + in + + (* Build the small matrix h = *) + let m_U = + Mat.of_col_vecs_list u_next + and m_W = + Mat.of_col_vecs_list w_next + in + let m_h = + gemm ~transa:`T m_U m_W + in + + (* Diagonalize h *) + let y, lambda = + Util.diagonalize_symm m_h + in + + (* Express m lowest eigenvectors of h in the large basis *) + let m_new_U = + gemm ~n:m m_U y + and m_new_W = + gemm ~n:m m_W y + in + + (* Compute the residual as proposed new vectors *) + let u_proposed = + Mat.init_cols n m (fun i k -> (lambda.{k} *. m_new_U.{i,k} -. m_new_W.{i,k}) /. + (max (diagonal.{i} -. lambda.{k}) 0.01) ) + |> Mat.to_col_vecs_list + in + + let residual_norms = List.map nrm2 u_proposed in + let residual_norm = List.fold_left (fun accu i -> max accu i) 0. residual_norms in + + Printf.printf "%3d %16.10f %16.8e%!\n" iter lambda.{1} residual_norm; + + if residual_norm > threshold then + let u_next, w_next, iter = + if iter = n_iter then + m_new_U |> pick_new, + m_new_W |> pick_new, + 0 + else + u_next, w_next, iter + in + iteration u_next u_proposed w_next (iter+1) + else + (Mat.of_col_vecs_list u_next |> pick_new |> Mat.of_col_vecs_list), lambda + + in + iteration [] u_new [] 1 + + diff --git a/Utils/Util.ml b/Utils/Util.ml index 02f18a9..6c7f2f9 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -178,6 +178,29 @@ let list_some l = |> List.map (function Some x -> x | _ -> assert false) +let list_range first last = + if last < first then [] else + let rec aux accu = function + | 0 -> first :: accu + | i -> aux ( (first+i)::accu ) (i-1) + in + aux [] (last-first) + + +let list_pack n l = + let rec aux i accu1 accu2 = function + | [] -> if accu1 = [] then + List.rev accu2 + else + List.rev ((List.rev accu1) :: accu2) + | a :: rest -> + match i with + | 0 -> aux (n-1) [] ((List.rev (a::accu1)) :: accu2) rest + | _ -> aux (i-1) (a::accu1) accu2 rest + in + aux (n-1) [] [] l + + (** {2 Stream functions} *) let stream_range first last = @@ -195,16 +218,6 @@ let stream_to_list stream = in aux [] -let list_range first last = - if last < first then [] else - let rec aux accu = function - | 0 -> first :: accu - | i -> aux ( (first+i)::accu ) (i-1) - in - aux [] (last-first) - - - (** {2 Linear algebra} *) @@ -250,6 +263,27 @@ let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c = gemm c u +let qr_ortho m = + (** Performed twice for precision *) + let result = lacpy m in + let tau = geqrf result in + orgqr ~tau result; + let tau = geqrf result in + orgqr ~tau result; + result + + +let normalize v = + let result = copy v in + scal (1. /. (nrm2 v)) result; + result + + +let normalize_mat m = + Mat.to_col_vecs m + |> Array.map (fun v -> normalize v) + |> Mat.of_col_vecs + (** {2 Bitstring functions} *) let bit_permtutations m n = diff --git a/Utils/Util.mli b/Utils/Util.mli index b2ff95e..47c59eb 100644 --- a/Utils/Util.mli +++ b/Utils/Util.mli @@ -1,3 +1,5 @@ +open Lacaml.D + (** All utilities which should be included in all source files are defined here *) (** {2 Functions from libm} *) @@ -69,6 +71,15 @@ val list_some : 'a option list -> 'a list val list_range : int -> int -> int list (** [list_range first last] returns a list [first; first+1 ; ... ; last-1 ; last ]. *) +val list_pack : int -> 'a list -> 'a list list +(** Example: +{[ + list_pack 3 [ 1; 2; 3; ...; 18; 19; 20 ] = + [[1; 2; 3]; [4; 5; 6]; [7; 8; 9]; [10; 11; 12]; [13; 14; 15]; + [16; 17; 18]; [19; 20]] +]} +*) + (** {2 Useful streams} *) val stream_range : int -> int -> int Stream.t (** [stream_range first last] returns a stream . *) @@ -78,16 +89,25 @@ val stream_to_list : 'a Stream.t -> 'a list (** {2 Linear algebra } *) -val diagonalize_symm : Lacaml.D.mat -> Lacaml.D.mat * Lacaml.D.vec +val normalize : Vec.t -> Vec.t +(** Returns a the vector normalized *) + +val normalize_mat : Mat.t -> Mat.t +(** Returns the matrix with all the column vectors normalized *) + +val diagonalize_symm : Mat.t -> Mat.t * Vec.t (** Diagonalize a symmetric matrix. Returns the eigenvectors and the eigenvalues. *) -val xt_o_x : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat +val xt_o_x : o:Mat.t -> x:Mat.t -> Mat.t (** Computes {% $\mathbf{X^\dag\, O\, X}$ %} *) -val x_o_xt : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat +val x_o_xt : o:Mat.t -> x:Mat.t -> Mat.t (** Computes {% $\mathbf{X\, O\, X^\dag}$ %} *) -val canonical_ortho: ?thresh:float -> overlap:Lacaml.D.mat -> Lacaml.D.mat -> Lacaml.D.mat +val qr_ortho : Mat.t -> Mat.t +(** QR orthogonalization of the input matrix *) + +val canonical_ortho: ?thresh:float -> overlap:Mat.t -> Mat.t -> Mat.t (** Canonical orthogonalization. [overlap] is the overlap matrix {% $\mathbf{S}$ %}, and the last argument contains the vectors {% $\mathbf{C}$ %} to orthogonalize. @@ -100,14 +120,14 @@ $$ *) -val debug_matrix: string -> Lacaml.D.mat -> unit +val debug_matrix: string -> Mat.t -> unit (** Prints a matrix in stdout for debug *) -val matrix_of_file : string -> Lacaml.D.mat +val matrix_of_file : string -> Mat.t (** Reads a matrix from a file with format "%d %d %f" corresponding to [i, j, A.{i,j}]. *) -val sym_matrix_of_file : string -> Lacaml.D.mat +val sym_matrix_of_file : string -> Mat.t (** Reads a symmetric matrix from a file with format "%d %d %f" corresponding to [i, j, A.{i,j}]. *) @@ -164,5 +184,5 @@ val pp_float_2darray : Format.formatter -> float array array -> unit val pp_bitstring : int -> Format.formatter -> Z.t -> unit (** Example: [ pp_bitstring 14 ppf z -> +++++------+-- ] *) -val pp_matrix : Format.formatter -> Lacaml.D.mat -> unit +val pp_matrix : Format.formatter -> Mat.t -> unit