diff --git a/ocaml/Id.ml b/ocaml/Id.ml index 660c3452..3e616922 100644 --- a/ocaml/Id.ml +++ b/ocaml/Id.ml @@ -1,26 +1,22 @@ -open Core.Std - -module Id : sig - type t - val of_int : int -> t - val to_int : t -> int - val of_string : string -> t - val to_string : t -> string - val increment : t -> t - val decrement : t -> t -end -= struct +module Id = struct type t = int + let of_int x = assert (x>0); x + let to_int x = x + let of_string x = - Int.of_string x + int_of_string x |> of_int + let to_string x = - Int.to_string x + string_of_int x + let increment x = x + 1 let decrement x = x - 1 + + let compare = compare end module Task = struct diff --git a/ocaml/Id.mli b/ocaml/Id.mli new file mode 100644 index 00000000..02d1efca --- /dev/null +++ b/ocaml/Id.mli @@ -0,0 +1,23 @@ +module Id : + sig + type t + val of_int : int -> t + val to_int : t -> int + val of_string : string -> t + val to_string : t -> string + val increment : t -> t + val decrement : t -> t + val compare : t -> t -> int + end + + +module Task : + sig + include (module type of Id) + end + + +module Client : + sig + include (module type of Id) + end diff --git a/ocaml/Progress_bar.ml b/ocaml/Progress_bar.ml index 2ca8bd00..b8e97a59 100644 --- a/ocaml/Progress_bar.ml +++ b/ocaml/Progress_bar.ml @@ -14,13 +14,13 @@ type t = let init ?(bar_length=20) ?(start_value=0.) ?(end_value=1.) ~title = { title ; start_value ; end_value ; bar_length ; cur_value=start_value ; - init_time= Time.now () ; dirty = true ; next = Time.now () } + init_time= Time.now () ; dirty = false ; next = Time.now () } let update ~cur_value bar = { bar with cur_value ; dirty=true } let increment_end bar = - { bar with end_value=(bar.end_value +. 1.) ; dirty=true } + { bar with end_value=(bar.end_value +. 1.) ; dirty=false } let increment_cur bar = { bar with cur_value=(bar.cur_value +. 1.) ; dirty=true } diff --git a/ocaml/Queuing_system.ml b/ocaml/Queuing_system.ml index 29a60538..0c668e16 100644 --- a/ocaml/Queuing_system.ml +++ b/ocaml/Queuing_system.ml @@ -1,27 +1,35 @@ -open Core.Std -open Qptypes - +module RunningMap = Map.Make (Id.Task) +module TasksMap = Map.Make (Id.Task) +module ClientsSet = Set.Make (Id.Client) type t = -{ queued : Id.Task.t list ; - running : (Id.Task.t, Id.Client.t) Map.Poly.t ; - tasks : (Id.Task.t, string) Map.Poly.t; - clients : Id.Client.t Set.Poly.t; +{ queued_front : Id.Task.t list ; + queued_back : Id.Task.t list ; + running : Id.Client.t RunningMap.t; + tasks : string TasksMap.t; + clients : ClientsSet.t; next_client_id : Id.Client.t; next_task_id : Id.Task.t; - number_of_queued : int; + number_of_queued : int; + number_of_running : int; + number_of_tasks : int; + number_of_clients : int; } let create () = - { queued = [] ; - running = Map.Poly.empty ; - tasks = Map.Poly.empty; - clients = Set.Poly.empty; + { queued_front = [] ; + queued_back = [] ; + running = RunningMap.empty ; + tasks = TasksMap.empty; + clients = ClientsSet.empty; next_client_id = Id.Client.of_int 1; next_task_id = Id.Task.of_int 1; - number_of_queued = 0; + number_of_queued = 0; + number_of_running = 0; + number_of_tasks = 0; + number_of_clients = 0; } @@ -32,10 +40,11 @@ let add_task ~task q = q.next_task_id in { q with - queued = task_id :: q.queued ; - tasks = Map.add q.tasks ~key:task_id ~data:task ; + queued_front = task_id :: q.queued_front ; + tasks = TasksMap.add task_id task q.tasks; next_task_id = Id.Task.increment task_id ; number_of_queued = q.number_of_queued + 1; + number_of_tasks = q.number_of_tasks + 1; } @@ -46,56 +55,73 @@ let add_client q = q.next_client_id in { q with - clients = Set.add q.clients client_id; + clients = ClientsSet.add client_id q.clients; next_client_id = Id.Client.increment client_id; + number_of_clients = q.number_of_clients + 1; }, client_id let pop_task ~client_id q = - let { queued ; running ; _ } = + let { queued_front ; queued_back ; running ; _ } = q in - assert (Set.mem q.clients client_id); - match queued with + assert (ClientsSet.mem client_id q.clients); + let queued_front', queued_back' = + match queued_front, queued_back with + | (l, []) -> ( [], List.rev l) + | t -> t + in + match queued_back' with | task_id :: new_queue -> let new_q = { q with - queued = new_queue ; - running = Map.add running ~key:task_id ~data:client_id ; - number_of_queued = q.number_of_queued - 1; + queued_front= queued_front' ; + queued_back = new_queue ; + running = RunningMap.add task_id client_id running; + number_of_queued = q.number_of_queued - 1; + number_of_running = q.number_of_running + 1; } - in new_q, Some task_id, (Map.find q.tasks task_id) + and found = + try Some (TasksMap.find task_id q.tasks) + with Not_found -> None + in new_q, Some task_id, found | [] -> q, None, None let del_client ~client_id q = - assert (Set.mem q.clients client_id); + assert (ClientsSet.mem client_id q.clients); { q with - clients = Set.remove q.clients client_id } + clients = ClientsSet.remove client_id q.clients; + number_of_clients = q.number_of_clients - 1 + } let end_task ~task_id ~client_id q = let { running ; tasks ; _ } = q in - assert (Set.mem q.clients client_id); - let () = - match Map.Poly.find running task_id with - | None -> failwith "Task already finished" - | Some client_id_check -> assert (client_id_check = client_id) + assert (ClientsSet.mem client_id q.clients); + let () = + let client_id_check = + try RunningMap.find task_id running with + Not_found -> failwith "Task already finished" + in + assert (client_id_check = client_id) in { q with - running = Map.remove running task_id ; + running = RunningMap.remove task_id running ; + number_of_running = q.number_of_running - 1 } - + let del_task ~task_id q = let { tasks ; _ } = q in - if (Map.mem tasks task_id) then + if (TasksMap.mem task_id tasks) then { q with - tasks = Map.remove tasks task_id ; + tasks = TasksMap.remove task_id tasks; + number_of_tasks = q.number_of_tasks - 1; } else Printf.sprintf "Task %d is already deleted" (Id.Task.to_int task_id) @@ -103,36 +129,81 @@ let del_task ~task_id q = -let number q = - Map.length q.tasks +let number_of_tasks q = + assert (q.number_of_tasks >= 0); + q.number_of_tasks let number_of_queued q = + assert (q.number_of_queued >= 0); q.number_of_queued let number_of_running q = - Map.length q.running + assert (q.number_of_running >= 0); + q.number_of_running + +let number_of_clients q = + assert (q.number_of_clients >= 0); + q.number_of_clients -let to_string { queued ; running ; tasks ; _ } = +let to_string qs = + let { queued_back ; queued_front ; running ; tasks ; _ } = qs in let q = - List.map ~f:Id.Task.to_string queued - |> String.concat ~sep:" ; " + (List.map Id.Task.to_string queued_front) @ + (List.map Id.Task.to_string @@ List.rev queued_back) + |> String.concat " ; " and r = - Map.Poly.to_alist running - |> List.map ~f:(fun (t,c) -> "("^(Id.Task.to_string t)^", " + RunningMap.bindings running + |> List.map (fun (t,c) -> "("^(Id.Task.to_string t)^", " ^(Id.Client.to_string c)^")") - |> String.concat ~sep:" ; " + |> String.concat " ; " and t = - Map.Poly.to_alist tasks - |> List.map ~f:(fun (t,c) -> "("^(Id.Task.to_string t)^", \"" + TasksMap.bindings tasks + |> List.map (fun (t,c) -> "("^(Id.Task.to_string t)^", \"" ^c^"\")") - |> String.concat ~sep:" ; " + |> String.concat " ; " in Printf.sprintf "{ +Tasks : %d Queued : %d Running : %d Clients : %d queued : { %s } running : { %s } tasks : [ %s ] -}" q r t +}" +(number_of_tasks qs) (number_of_queued qs) (number_of_running qs) (number_of_clients qs) +q r t + +let test () = + let q = + create () + |> add_task ~task:"First Task" + |> add_task ~task:"Second Task" + in + let q, client_id = + add_client q + in + let q, task_id, task_content = + match pop_task ~client_id q with + | q, Some x, Some y -> q, Id.Task.to_int x, y + | _ -> assert false + in + Printf.printf "Task_id : %d \t\t Task : %s\n" task_id task_content; + to_string q |> print_endline ; + let q, task_id, task_content = + match pop_task ~client_id q with + | q, Some x, Some y -> q, Id.Task.to_int x, y + | _ -> assert false + in + Printf.printf "Task_id : %d \t\t Task : %s\n" task_id task_content; + let q, task_id, task_content = + match pop_task ~client_id q with + | q, None, None -> q, 0, "None" + | _ -> assert false + in + Printf.printf "Task_id : %d \t\t Task : %s\n" task_id task_content; + q + |> to_string + |> print_endline + diff --git a/ocaml/Queuing_system.mli b/ocaml/Queuing_system.mli new file mode 100644 index 00000000..dc6836d2 --- /dev/null +++ b/ocaml/Queuing_system.mli @@ -0,0 +1,63 @@ +module RunningMap : Map.S with type key = Id.Task.t +module TasksMap : Map.S with type key = Id.Task.t +module ClientsSet : Set.S with type elt = Id.Client.t + +type t = { + queued_front : Id.Task.t list ; + queued_back : Id.Task.t list ; + running : Id.Client.t RunningMap.t ; + tasks : string TasksMap.t ; + clients : ClientsSet.t ; + next_client_id : Id.Client.t ; + next_task_id : Id.Task.t ; + number_of_queued : int ; + number_of_running : int ; + number_of_tasks : int ; + number_of_clients : int ; +} + +(** Creates a new queuing system. Returns the new queue. *) +val create : unit -> t + +(** Add a new task represented as a string. Returns the queue with the added task. *) +val add_task : task:string -> t -> t + +(** Add a new client. Returns the queue and a new client_id. *) +val add_client : t -> t * Id.Client.t + +(** Pops a task from the queue. The task is set as running on client client_id. + Returns the queue, a task_id and the content of the task. If the queue contains + no task, the task_id and the task content are None. *) +val pop_task : + client_id:ClientsSet.elt -> t -> t * Id.Task.t option * string option + +(** Deletes a client from the queuing system *) +val del_client : client_id:ClientsSet.elt -> t -> t + +(** Deletes a client from the queuing system. The client is assumed to be a member + of the set of clients. Returns the queue without the removed client. *) +val end_task : task_id:RunningMap.key -> client_id:ClientsSet.elt -> t -> t + +(** Deletes a task from the queuing system. The task is assumed to be a member + of the map of tasks. Returns the queue without the removed task. *) +val del_task : task_id:TasksMap.key -> t -> t + +(** Returns the number of tasks, assumed >= 0 *) +val number_of_tasks : t -> int + +(** Returns the number of queued tasks, assumed >= 0 *) +val number_of_queued : t -> int + +(** Returns the number of running tasks, assumed >= 0 *) +val number_of_running : t -> int + +(** Returns the number of connected clients, assumed >= 0 *) +val number_of_clients : t -> int + +(** Prints the content of the queue *) +val to_string : t -> string + +(** Test function for debug *) +val test : unit -> unit + + diff --git a/ocaml/TaskServer.ml b/ocaml/TaskServer.ml index cfc22cfc..9a1797f8 100644 --- a/ocaml/TaskServer.ml +++ b/ocaml/TaskServer.ml @@ -306,7 +306,7 @@ let del_task msg program_state rep_socket = } in let more = - (Queuing_system.number new_program_state.queue > 0) + (Queuing_system.number_of_tasks new_program_state.queue > 0) in Message.DelTaskReply (Message.DelTaskReply_msg.create ~task_id ~more) |> Message.to_string @@ -678,9 +678,9 @@ let run ~port = (** Debug input *) Printf.sprintf "q:%d r:%d n:%d : %s\n%!" - (Queuing_system.number_of_queued program_state.queue) + (Queuing_system.number_of_queued program_state.queue) (Queuing_system.number_of_running program_state.queue) - (Queuing_system.number program_state.queue) + (Queuing_system.number_of_tasks program_state.queue) (Message.to_string message) |> debug; diff --git a/plugins/Full_CI_ZMQ/fci_zmq.irp.f b/plugins/Full_CI_ZMQ/fci_zmq.irp.f index c81b1266..a5dd8dcf 100644 --- a/plugins/Full_CI_ZMQ/fci_zmq.irp.f +++ b/plugins/Full_CI_ZMQ/fci_zmq.irp.f @@ -134,7 +134,7 @@ subroutine ZMQ_selection(N_in, pt2) step = int(5000000.d0 / dble(N_int * N_states * elec_num * elec_num * mo_tot_num * mo_tot_num )) step = max(1,step) - do i= N_det_generators, 1, -step + do i= 1,N_det_generators, step i_generator_start = max(i-step+1,1) i_generator_max = i write(task,*) i_generator_start, i_generator_max, 1, N diff --git a/src/Davidson/diagonalization_hs2.irp.f b/src/Davidson/diagonalization_hs2.irp.f index c8ac3733..2db6b4cd 100644 --- a/src/Davidson/diagonalization_hs2.irp.f +++ b/src/Davidson/diagonalization_hs2.irp.f @@ -95,7 +95,6 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s double precision :: u_dot_v, u_dot_u - integer, allocatable :: kl_pairs(:,:) integer :: k_pairs, kl integer :: iter2 @@ -107,12 +106,14 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s character*(16384) :: write_buffer double precision :: to_print(3,N_st) double precision :: cpu, wall - integer :: shift, shift2 + integer :: shift, shift2, itermax include 'constants.include.F' !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, R, S, y, h, lambda - if (N_st_diag > sze) then - stop 'error in Davidson : N_st_diag > sze' + if (N_st_diag*3 > sze) then + print *, 'error in Davidson :' + print *, 'Increase n_det_max_jacobi to ', N_st_diag*3 + stop -1 endif PROVIDE nuclear_repulsion @@ -147,26 +148,26 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s integer, external :: align_double sze_8 = align_double(sze) + itermax = min(davidson_sze_max, sze/N_st_diag) allocate( & - kl_pairs(2,N_st_diag*(N_st_diag+1)/2), & - W(sze_8,N_st_diag*davidson_sze_max), & - U(sze_8,N_st_diag*davidson_sze_max), & + W(sze_8,N_st_diag*itermax), & + U(sze_8,N_st_diag*itermax), & R(sze_8,N_st_diag), & - S(sze_8,N_st_diag*davidson_sze_max), & - h(N_st_diag*davidson_sze_max,N_st_diag*davidson_sze_max), & - y(N_st_diag*davidson_sze_max,N_st_diag*davidson_sze_max), & - s_(N_st_diag*davidson_sze_max,N_st_diag*davidson_sze_max), & - s_tmp(N_st_diag*davidson_sze_max,N_st_diag*davidson_sze_max), & + S(sze_8,N_st_diag*itermax), & + h(N_st_diag*itermax,N_st_diag*itermax), & + y(N_st_diag*itermax,N_st_diag*itermax), & + s_(N_st_diag*itermax,N_st_diag*itermax), & + s_tmp(N_st_diag*itermax,N_st_diag*itermax), & residual_norm(N_st_diag), & - c(N_st_diag*davidson_sze_max), & - s2(N_st_diag*davidson_sze_max), & - lambda(N_st_diag*davidson_sze_max)) + c(N_st_diag*itermax), & + s2(N_st_diag*itermax), & + lambda(N_st_diag*itermax)) h = 0.d0 s_ = 0.d0 s_tmp = 0.d0 - c = 0.d0 U = 0.d0 + W = 0.d0 S = 0.d0 R = 0.d0 y = 0.d0 @@ -183,10 +184,6 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s converged = .False. - do k=1,N_st - call normalize(u_in(1,k),sze) - enddo - do k=N_st+1,N_st_diag do i=1,sze double precision :: r1, r2 @@ -194,14 +191,6 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s call random_number(r2) u_in(i,k) = dsqrt(-2.d0*dlog(r1))*dcos(dtwo_pi*r2) enddo - - ! Gram-Schmidt - ! ------------ - call dgemv('T',sze,k-1,1.d0,u_in,size(u_in,1), & - u_in(1,k),1,0.d0,c,1) - call dgemv('N',sze,k-1,-1.d0,u_in,size(u_in,1), & - c,1,1.d0,u_in(1,k),1) - call normalize(u_in(1,k),sze) enddo @@ -213,11 +202,12 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s enddo enddo - do iter=1,davidson_sze_max-1 + do iter=1,itermax-1 shift = N_st_diag*(iter-1) shift2 = N_st_diag*iter + call ortho_qr(U,size(U,1),sze,shift2) ! Compute |W_k> = \sum_i |i> ! ----------------------------------------- @@ -229,20 +219,6 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s ! Compute h_kl = = ! ------------------------------------------- - -! do l=1,N_st_diag -! do k=1,N_st_diag -! do iter2=1,iter-1 -! h(k,iter2,l,iter) = u_dot_v(U(1,k,iter2),W(1,l,iter),sze) -! h(k,iter,l,iter2) = h(k,iter2,l,iter) -! enddo -! enddo -! do k=1,l -! h(k,iter,l,iter) = u_dot_v(U(1,k,iter),W(1,l,iter),sze) -! h(l,iter,k,iter) = h(k,iter,l,iter) -! enddo -! enddo - call dgemm('T','N', shift2, N_st_diag, sze, & 1.d0, U, size(U,1), W(1,shift+1), size(W,1), & 0.d0, h(1,shift+1), size(h,1)) @@ -295,22 +271,6 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s ! Express eigenvectors of h in the determinant basis ! -------------------------------------------------- -! do k=1,N_st_diag -! do i=1,sze -! U(i,shift2+k) = 0.d0 -! W(i,shift2+k) = 0.d0 -! S(i,shift2+k) = 0.d0 -! enddo -! do l=1,N_st_diag*iter -! do i=1,sze -! U(i,shift2+k) = U(i,shift2+k) + U(i,l)*y(l,k) -! W(i,shift2+k) = W(i,shift2+k) + W(i,l)*y(l,k) -! S(i,shift2+k) = S(i,shift2+k) + S(i,l)*y(l,k) -! enddo -! enddo -! enddo -! -! call dgemm('N','N', sze, N_st_diag, shift2, & 1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1)) call dgemm('N','N', sze, N_st_diag, shift2, & @@ -321,13 +281,6 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s ! Compute residual vector ! ----------------------- -! do k=1,N_st_diag -! print *, s2(k) -! s2(k) = u_dot_v(U(1,shift2+k), S(1,shift2+k), sze) + S_z2_Sz -! print *, s2(k) -! print *, '' -! pause -! enddo do k=1,N_st_diag do i=1,sze R(i,k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & @@ -338,14 +291,17 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s to_print(1,k) = lambda(k) + nuclear_repulsion to_print(2,k) = s2(k) to_print(3,k) = residual_norm(k) - if (residual_norm(k) > 1.e9) then - stop 'Davidson failed' - endif endif enddo - write(iunit,'(X,I3,X,100(X,F16.10,X,F11.6,X,E11.3))') iter, to_print(:,1:N_st) + write(iunit,'(X,I3,X,100(X,F16.10,X,F11.6,X,E11.3,A20))') iter, to_print(:,1:N_st), '' call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged) + do k=1,N_st + if (residual_norm(k) > 1.e9) then + print *, '' + stop 'Davidson failed' + endif + enddo if (converged) then exit endif @@ -359,42 +315,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s enddo enddo - ! Gram-Schmidt - ! ------------ - - do k=1,N_st_diag - -! do l=1,N_st_diag*iter -! c(1) = u_dot_v(U(1,shift2+k),U(1,l),sze) -! do i=1,sze -! U(i,k,iter+1) = U(i,shift2+k) - c(1) * U(i,l) -! enddo -! enddo -! - call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1), & - U(1,shift2+k),1,0.d0,c,1) - call dgemv('N',sze,N_st_diag*iter,-1.d0,U,size(U,1), & - c,1,1.d0,U(1,shift2+k),1) -! -! do l=1,k-1 -! c(1) = u_dot_v(U(1,shift2+k),U(1,shift2+l),sze) -! do i=1,sze -! U(i,k,iter+1) = U(i,shift2+k) - c(1) * U(i,shift2+l) -! enddo -! enddo -! - call dgemv('T',sze,k-1,1.d0,U(1,shift2+1),size(U,1), & - U(1,shift2+k),1,0.d0,c,1) - call dgemv('N',sze,k-1,-1.d0,U(1,shift2+1),size(U,1), & - c,1,1.d0,U(1,shift2+k),1) - - call normalize( U(1,shift2+k), sze ) - enddo - enddo if (.not.converged) then - iter = davidson_sze_max-1 + iter = itermax-1 endif ! Re-contract to u_in @@ -404,20 +328,14 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s energies(k) = lambda(k) enddo -! do k=1,N_st_diag -! do i=1,sze -! do l=1,iter*N_st_diag -! u_in(i,k) += U(i,l)*y(l,k) -! enddo -! enddo -! enddo -! enddo - call dgemm('N','N', sze, N_st_diag, N_st_diag*iter, 1.d0, & U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) enddo + do k=1,N_st_diag + S2_jj(k) = s2(k) + enddo write_buffer = '===== ' do i=1,N_st write_buffer = trim(write_buffer)//' ================ =========== ===========' @@ -427,7 +345,6 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s call write_time(iunit) deallocate ( & - kl_pairs, & W, residual_norm, & U, & R, c, S, & diff --git a/src/Davidson/u0Hu0.irp.f b/src/Davidson/u0Hu0.irp.f index a1a72100..3787370a 100644 --- a/src/Davidson/u0Hu0.irp.f +++ b/src/Davidson/u0Hu0.irp.f @@ -252,7 +252,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8) ave_workload = ave_workload/dble(shortcut(0,1)) - do sh=shortcut(0,1),1,-1 + do sh=1,shortcut(0,1),1 workload = shortcut(0,1)+dble(shortcut(sh+1,1) - shortcut(sh,1))**2 do i=sh, shortcut(0,2), shortcut(0,1) do j=i, min(i, shortcut(0,2)) diff --git a/src/Determinants/H_apply_zmq.template.f b/src/Determinants/H_apply_zmq.template.f index d59f2994..59544b79 100644 --- a/src/Determinants/H_apply_zmq.template.f +++ b/src/Determinants/H_apply_zmq.template.f @@ -35,7 +35,7 @@ subroutine $subroutine($params_main) call zmq_put_psi(zmq_to_qp_run_socket,1,energy,size(energy)) - do i_generator=N_det_generators,1,-1 + do i_generator=1,N_det_generators $skip write(task,*) i_generator call add_task_to_taskserver(zmq_to_qp_run_socket,task) diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index 2ebb402e..9eadbf35 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -368,7 +368,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] call new_parallel_job(zmq_to_qp_run_socket,'ao_integrals') - do l=1,ao_num + do l=ao_num,1,-1 write(task,*) "triangle ", l call add_task_to_taskserver(zmq_to_qp_run_socket,task) enddo diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 00f61101..e44e8c2c 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -11,9 +11,9 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) integer, intent(in) :: LDA, LDU, LDVt, m, n double precision, intent(in) :: A(LDA,n) - double precision, intent(out) :: U(LDU,n) + double precision, intent(out) :: U(LDU,m) double precision,intent(out) :: Vt(LDVt,n) - double precision,intent(out) :: D(n) + double precision,intent(out) :: D(min(m,n)) double precision,allocatable :: work(:) integer :: info, lwork, i, j, k @@ -24,13 +24,13 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) ! Find optimal size for temp arrays allocate(work(1)) lwork = -1 - call dgesvd('A','A', n, n, A_tmp, LDA, & + call dgesvd('A','A', m, n, A_tmp, LDA, & D, U, LDU, Vt, LDVt, work, lwork, info) lwork = work(1) deallocate(work) allocate(work(lwork)) - call dgesvd('A','A', n, n, A_tmp, LDA, & + call dgesvd('A','A', m, n, A_tmp, LDA, & D, U, LDU, Vt, LDVt, work, lwork, info) deallocate(work,A_tmp) @@ -125,6 +125,40 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) end +subroutine ortho_qr(A,LDA,m,n) + implicit none + BEGIN_DOC + ! Orthogonalization using Q.R factorization + ! + ! A : matrix to orthogonalize + ! + ! LDA : leftmost dimension of A + ! + ! n : Number of rows of A + ! + ! m : Number of columns of A + ! + END_DOC + integer, intent(in) :: m,n, LDA + double precision, intent(inout) :: A(LDA,n) + + integer :: lwork, info + integer, allocatable :: jpvt(:) + double precision, allocatable :: tau(:), work(:) + + allocate (jpvt(n), tau(n), work(1)) + LWORK=-1 +! call dgeqp3(m, n, A, LDA, jpvt, tau, WORK, LWORK, INFO) + call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO ) + LWORK=WORK(1) + deallocate(WORK) + allocate(WORK(LWORK)) +! call dgeqp3(m, n, A, LDA, jpvt, tau, WORK, LWORK, INFO) + call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO ) + call dorgqr(m, n, n, A, LDA, tau, WORK, LWORK, INFO) + deallocate(WORK,jpvt,tau) +end + subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) implicit none BEGIN_DOC @@ -161,7 +195,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) allocate(U(ldc,n),Vt(lda,n),S_half(lda,n),D(n)) - call svd(overlap,lda,U,ldc,D,Vt,lda,m,n) + call svd(overlap,lda,U,ldc,D,Vt,lda,n,n) !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(S_half,U,D,Vt,n,C,m) &