LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-09-27 20:11:25 +02:00
Code Issues Releases Wiki Activity

Merge amazing pull request #193 from scemama/master

Browse Source 
Improved Stochastic PT2
This commit is contained in:
Thomas Applencourt 2017-05-09 21:39:58 -05:00 committed by GitHub
commit 5e2309246c
30 changed files with 1176 additions and 499 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
ocaml/Charge.ml
View File

@ -1,6 +1,6 @@
open Core.Std;;
open Core.Std
type t = float with sexp;;
type t = float with sexp
let of_float x = x
let of_int i = Float.of_int i
@ -14,5 +14,4 @@ let to_string x =
Printf.sprintf "+%f" x
else
Printf.sprintf "%f" x
;;

145
ocaml/Message.ml
View File

@ -455,6 +455,122 @@ end = struct
end
(** GetVector : get the current vector (Davidson) *)
module GetVector_msg : sig
type t =
{ client_id: Id.Client.t ;
}
val create : client_id:int -> t
val to_string : t -> string
end = struct
type t =
{ client_id: Id.Client.t ;
}
let create ~client_id =
{ client_id = Id.Client.of_int client_id }
let to_string x =
Printf.sprintf "get_vector %d"
(Id.Client.to_int x.client_id)
end
module Vector : sig
type t =
{
size : Strictly_positive_int.t;
data : string;
}
val create : size:Strictly_positive_int.t -> data:string -> t
end = struct
type t =
{
size : Strictly_positive_int.t;
data : string;
}
let create ~size ~data =
{ size ; data }
end
(** GetVectorReply_msg : Reply to the GetVector message *)
module GetVectorReply_msg : sig
type t =
{ client_id : Id.Client.t ;
vector : Vector.t }
val create : client_id:Id.Client.t -> vector:Vector.t -> t
val to_string : t -> string
val to_string_list : t -> string list
end = struct
type t =
{ client_id : Id.Client.t ;
vector : Vector.t }
let create ~client_id ~vector =
{ client_id ; vector }
let to_string x =
Printf.sprintf "get_vector_reply %d %d"
(Id.Client.to_int x.client_id)
(Strictly_positive_int.to_int x.vector.Vector.size)
let to_string_list x =
[ to_string x ; x.vector.Vector.data ]
end
(** PutVector : put the current variational wave function *)
module PutVector_msg : sig
type t =
{ client_id : Id.Client.t ;
size : Strictly_positive_int.t ;
vector : Vector.t option;
}
val create :
client_id:int -> size:int -> data:string option -> t
val to_string_list : t -> string list
val to_string : t -> string
end = struct
type t =
{ client_id : Id.Client.t ;
size : Strictly_positive_int.t ;
vector : Vector.t option;
}
let create ~client_id ~size ~data =
let size =
Strictly_positive_int.of_int size
in
let vector =
match data with
| None -> None
| Some s -> Some (Vector.create ~size ~data:s)
in
{ client_id = Id.Client.of_int client_id ;
vector ; size
}
let to_string x =
Printf.sprintf "put_vector %d %d"
(Id.Client.to_int x.client_id)
(Strictly_positive_int.to_int x.size)
let to_string_list x =
match x.vector with
| Some v -> [ to_string x ; v.Vector.data ]
| None -> failwith "Empty vector"
end
(** PutVectorReply_msg : Reply to the PutVector message *)
module PutVectorReply_msg : sig
type t
val create : client_id:Id.Client.t -> t
val to_string : t -> string
end = struct
type t =
{ client_id : Id.Client.t ;
}
let create ~client_id =
{ client_id; }
let to_string x =
Printf.sprintf "put_vector_reply %d"
(Id.Client.to_int x.client_id)
end
(** TaskDone : Inform the server that a task is finished *)
module TaskDone_msg : sig
type t =
@ -494,6 +610,17 @@ end = struct
let to_string x = "terminate"
end
(** Abort *)
module Abort_msg : sig
type t
val create : t
val to_string : t -> string
end = struct
type t = Abort
let create = Abort
let to_string x = "abort"
end
(** OK *)
module Ok_msg : sig
type t
@ -526,6 +653,10 @@ type t =
| PutPsi of PutPsi_msg.t
| GetPsiReply of GetPsiReply_msg.t
| PutPsiReply of PutPsiReply_msg.t
| GetVector of GetVector_msg.t
| PutVector of PutVector_msg.t
| GetVectorReply of GetVectorReply_msg.t
| PutVectorReply of PutVectorReply_msg.t
| Newjob of Newjob_msg.t
| Endjob of Endjob_msg.t
| Connect of Connect_msg.t
@ -540,6 +671,7 @@ type t =
| AddTaskReply of AddTaskReply_msg.t
| TaskDone of TaskDone_msg.t
| Terminate of Terminate_msg.t
| Abort of Abort_msg.t
| Ok of Ok_msg.t
| Error of Error_msg.t
| SetStopped
@ -580,7 +712,12 @@ let of_string s =
~n_det_generators:None ~n_det_selectors:None
~psi_det:None ~psi_coef:None ~energy:None )
end
| GetVector_ client_id ->
GetVector (GetVector_msg.create ~client_id)
| PutVector_ { client_id ; size } ->
PutVector (PutVector_msg.create ~client_id ~size ~data:None )
| Terminate_ -> Terminate (Terminate_msg.create )
| Abort_ -> Abort (Abort_msg.create )
| SetWaiting_ -> SetWaiting
| SetStopped_ -> SetStopped
| SetRunning_ -> SetRunning
@ -592,6 +729,8 @@ let of_string s =
let to_string = function
| GetPsi x -> GetPsi_msg.to_string x
| PutPsiReply x -> PutPsiReply_msg.to_string x
| GetVector x -> GetVector_msg.to_string x
| PutVectorReply x -> PutVectorReply_msg.to_string x
| Newjob x -> Newjob_msg.to_string x
| Endjob x -> Endjob_msg.to_string x
| Connect x -> Connect_msg.to_string x
@ -606,10 +745,13 @@ let to_string = function
| AddTaskReply x -> AddTaskReply_msg.to_string x
| TaskDone x -> TaskDone_msg.to_string x
| Terminate x -> Terminate_msg.to_string x
| Abort x -> Abort_msg.to_string x
| Ok x -> Ok_msg.to_string x
| Error x -> Error_msg.to_string x
| PutPsi x -> PutPsi_msg.to_string x
| GetPsiReply x -> GetPsiReply_msg.to_string x
| PutVector x -> PutVector_msg.to_string x
| GetVectorReply x -> GetVectorReply_msg.to_string x
| SetStopped -> "set_stopped"
| SetRunning -> "set_running"
| SetWaiting -> "set_waiting"
@ -618,4 +760,7 @@ let to_string = function
let to_string_list = function
| PutPsi x -> PutPsi_msg.to_string_list x
| GetPsiReply x -> GetPsiReply_msg.to_string_list x
| PutVector x -> PutVector_msg.to_string_list x
| GetVectorReply x -> GetVectorReply_msg.to_string_list x
| _ -> assert false

27
ocaml/Message_lexer.mll
View File

@ -15,8 +15,11 @@ type kw_type =
| NEW_JOB
| END_JOB
| TERMINATE
| ABORT
| GET_PSI
| PUT_PSI
| GET_VECTOR
| PUT_VECTOR
| OK
| ERROR
| SET_STOPPED
@ -29,7 +32,8 @@ type state_taskids_clientid = { state : string ; task_ids : int list ;
type state_clientid = { state : string ; client_id : int ; }
type state_tcp_inproc = { state : string ; push_address_tcp : string ; push_address_inproc : string ; }
type psi = { client_id: int ; n_state: int ; n_det: int ; psi_det_size: int ;
n_det_generators: int option ; n_det_selectors: int option }
n_det_generators: int option ; n_det_selectors: int option ; }
type vector = { client_id: int ; size: int }
type msg =
| AddTask_ of state_tasks
@ -41,8 +45,11 @@ type msg =
| NewJob_ of state_tcp_inproc
| EndJob_ of string
| Terminate_
| Abort_
| GetPsi_ of int
| PutPsi_ of psi
| GetVector_ of int
| PutVector_ of vector
| Ok_
| Error_ of string
| SetStopped_
@ -83,8 +90,11 @@ and kw = parse
| "new_job" { NEW_JOB }
| "end_job" { END_JOB }
| "terminate" { TERMINATE }
| "abort" { ABORT }
| "get_psi" { GET_PSI }
| "put_psi" { PUT_PSI }
| "get_vector" { GET_PSI }
| "put_vector" { PUT_PSI }
| "ok" { OK }
| "error" { ERROR }
| "set_stopped" { SET_STOPPED }
@ -179,6 +189,15 @@ and kw = parse
in
PutPsi_ { client_id ; n_state ; n_det ; psi_det_size ; n_det_generators ; n_det_selectors }
| GET_VECTOR ->
let client_id = read_int lexbuf in
GetVector_ client_id
| PUT_VECTOR ->
let client_id = read_int lexbuf in
let size = read_int lexbuf in
PutVector_ { client_id ; size }
| CONNECT ->
let socket = read_word lexbuf in
Connect_ socket
@ -202,6 +221,7 @@ and kw = parse
| SET_RUNNING -> SetRunning_
| SET_STOPPED -> SetStopped_
| TERMINATE -> Terminate_
| ABORT -> Abort_
| NONE -> parse_rec lexbuf
| _ -> failwith "Error in MessageLexer"
@ -226,6 +246,7 @@ and kw = parse
"new_job state_pouet tcp://test.com:12345 ipc:///dev/shm/x.socket";
"end_job state_pouet";
"terminate" ;
"abort" ;
"set_running" ;
"set_stopped" ;
"set_waiting" ;
@ -253,7 +274,11 @@ and kw = parse
| Some s, Some g -> Printf.sprintf "PUT_PSI client_id:%d n_state:%d n_det:%d psi_det_size:%d n_det_generators:%d n_det_selectors:%d" client_id n_state n_det psi_det_size g s
| _ -> Printf.sprintf "PUT_PSI client_id:%d n_state:%d n_det:%d psi_det_size:%d" client_id n_state n_det psi_det_size
end
| GetVector_ client_id -> Printf.sprintf "GET_VECTOR client_id:%d" client_id
| PutVector_ { client_id ; size } ->
Printf.sprintf "PUT_VECTOR client_id:%d size:%d" client_id size
| Terminate_ -> "TERMINATE"
| Abort_ -> "ABORT"
| SetWaiting_ -> "SET_WAITING"
| SetStopped_ -> "SET_STOPPED"
| SetRunning_ -> "SET_RUNNING"

9
ocaml/Progress_bar.ml
View File

@ -22,6 +22,11 @@ let update ~cur_value bar =
let increment_end bar =
{ bar with end_value=(bar.end_value +. 1.) ; dirty=false }
let clear bar =
Printf.eprintf " \r%!";
None
let increment_cur bar =
{ bar with cur_value=(bar.cur_value +. 1.) ; dirty=true }
@ -56,7 +61,7 @@ let display_tty bar =
else
Time.Span.of_float 0.
in
Printf.printf "%s : [%s] %4.1f%% | %10s, ~%10s left\r%!"
Printf.eprintf "%s : [%s] %4.1f%% | %10s, ~%10s left\r%!"
bar.title
hashes
percent
@ -83,7 +88,7 @@ let display_file bar =
else
Time.Span.of_float 0.
in
Printf.printf "%5.2f %% in %20s, ~%20s left\n%!"
Printf.eprintf "%5.2f %% in %20s, ~%20s left\n%!"
percent
(Time.Span.to_string running_time)
(Time.Span.to_string stop_time);

125
ocaml/TaskServer.ml
View File

@ -26,6 +26,7 @@ type t =
address_tcp : Address.Tcp.t option ;
address_inproc : Address.Inproc.t option ;
psi : Message.Psi.t option;
vector : Message.Vector.t option;
progress_bar : Progress_bar.t option ;
running : bool;
}
@ -48,12 +49,7 @@ let zmq_context =
ZMQ.Context.create ()
let () =
let nproc =
match Sys.getenv "OMP_NUM_THREADS" with
| Some m -> int_of_string m
| None -> 2
in
ZMQ.Context.set_io_threads zmq_context nproc
ZMQ.Context.set_io_threads zmq_context 2
let bind_socket ~socket_type ~socket ~port =
@ -69,16 +65,7 @@ let bind_socket ~socket_type ~socket ~port =
with
| Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_float 1. ; loop (i-1) )
| other_exception -> raise other_exception
in loop 60;
let filename =
Printf.sprintf "/tmp/qp_run:%d" port
in
begin
match Sys.file_exists filename with
| `Yes -> Sys.remove filename
| _ -> ()
end;
ZMQ.Socket.bind socket ("ipc://"^filename)
in loop 60
let hostname = lazy (
@ -132,7 +119,7 @@ let stop ~port =
let req_socket =
ZMQ.Socket.create zmq_context ZMQ.Socket.req
and address =
Printf.sprintf "ipc:///tmp/qp_run:%d" port
Printf.sprintf "tcp://localhost:%d" port
in
ZMQ.Socket.set_linger_period req_socket 1_000_000;
ZMQ.Socket.connect req_socket address;
@ -212,7 +199,7 @@ let end_job msg program_state rep_socket pair_socket =
reply_ok rep_socket;
{ program_state with
state = None ;
progress_bar = None ;
progress_bar = Progress_bar.clear ();
}
in
@ -389,7 +376,12 @@ let get_task msg program_state rep_socket pair_socket =
let new_queue, task_id, task =
Queuing_system.pop_task ~client_id program_state.queue
in
if (Queuing_system.number_of_queued new_queue = 0) then
let no_task =
Queuing_system.number_of_queued new_queue = 0
in
if no_task then
string_of_pub_state Waiting
|> ZMQ.Socket.send pair_socket
else
@ -518,16 +510,100 @@ let get_psi msg program_state rep_socket =
let put_vector msg rest_of_msg program_state rep_socket =
let vector_local =
match msg.Message.PutVector_msg.vector with
| Some x -> x
| None ->
begin
let data =
match rest_of_msg with
| [ x ] -> x
| _ -> failwith "Badly formed put_vector message"
in
Message.Vector.create
~size:msg.Message.PutVector_msg.size
~data
end
in
let new_program_state =
{ program_state with
vector = Some vector_local
}
and client_id =
msg.Message.PutVector_msg.client_id
in
Message.PutVectorReply (Message.PutVectorReply_msg.create ~client_id)
|> Message.to_string
|> ZMQ.Socket.send rep_socket;
new_program_state
let get_vector msg program_state rep_socket =
let client_id =
msg.Message.GetVector_msg.client_id
in
match program_state.vector with
| None -> failwith "No wave function saved in TaskServer"
| Some vector ->
Message.GetVectorReply (Message.GetVectorReply_msg.create ~client_id ~vector)
|> Message.to_string_list
|> ZMQ.Socket.send_all rep_socket;
program_state
let terminate program_state rep_socket =
reply_ok rep_socket;
{ program_state with
psi = None;
vector = None;
address_tcp = None;
address_inproc = None;
running = false
}
let abort program_state rep_socket =
let queue, client_id =
Queuing_system.add_client program_state.queue
in
let rec aux accu queue = function
| 0 -> (queue, accu)
| rest ->
let new_queue, task_id, _ =
Queuing_system.pop_task ~client_id queue
in
let new_accu =
match task_id with
| Some task_id -> task_id::accu
| None -> accu
in
Queuing_system.number_of_queued new_queue
|> aux new_accu new_queue
in
let queue, tasks =
aux [] queue 1
in
let queue =
List.fold ~f:(fun queue task_id ->
Queuing_system.end_task ~task_id ~client_id queue)
~init:queue tasks
in
let queue =
List.fold ~f:(fun queue task_id -> Queuing_system.del_task ~task_id queue)
~init:queue tasks
in
reply_ok rep_socket;
{ program_state with
queue
}
let error msg program_state rep_socket =
Message.Error (Message.Error_msg.create msg)
|> Message.to_string
@ -605,6 +681,7 @@ let run ~port =
{ queue = Queuing_system.create () ;
running = true ;
psi = None;
vector = None;
state = None;
address_tcp = None;
address_inproc = None;
@ -658,17 +735,25 @@ let run ~port =
in
(** Debug input *)
let () =
if debug_env then
begin
Printf.sprintf "q:%d r:%d n:%d : %s\n%!"
(Queuing_system.number_of_queued program_state.queue)
(Queuing_system.number_of_running program_state.queue)
(Queuing_system.number_of_tasks program_state.queue)
(Message.to_string message)
|> debug;
|> debug
end
in
let new_program_state =
try
match program_state.state, message with
| _ , Message.Terminate _ -> terminate program_state rep_socket
| _ , Message.Abort _ -> abort program_state rep_socket
| _ , Message.PutVector x -> put_vector x rest program_state rep_socket
| _ , Message.GetVector x -> get_vector x program_state rep_socket
| _ , Message.PutPsi x -> put_psi x rest program_state rep_socket
| _ , Message.GetPsi x -> get_psi x program_state rep_socket
| None , Message.Newjob x -> new_job x program_state rep_socket pair_socket

1
ocaml/TaskServer.mli
View File

@ -5,6 +5,7 @@ type t =
address_tcp : Address.Tcp.t option ;
address_inproc : Address.Inproc.t option ;
psi : Message.Psi.t option;
vector : Message.Vector.t option ;
progress_bar : Progress_bar.t option ;
running : bool;
}

52
plugins/CAS_SD_ZMQ/cassd_zmq.irp.f
View File

@ -1,16 +1,23 @@
program fci_zmq
program cassd_zmq
implicit none
integer :: i,j,k
logical, external :: detEq
double precision, allocatable :: pt2(:)
integer :: degree
integer :: n_det_before, to_select
double precision :: threshold_davidson_in
allocate (pt2(N_states))
double precision :: hf_energy_ref
logical :: has
pt2 = -huge(1.d0)
threshold_davidson_in = threshold_davidson
threshold_davidson = threshold_davidson_in * 100.d0
SOFT_TOUCH threshold_davidson
call diagonalize_CI
call save_wavefunction
call ezfio_has_hartree_fock_energy(has)
if (has) then
call ezfio_get_hartree_fock_energy(hf_energy_ref)
@ -18,15 +25,7 @@ program fci_zmq
hf_energy_ref = ref_bitmask_energy
endif
pt2 = 1.d0
threshold_davidson_in = threshold_davidson
threshold_davidson = threshold_davidson_in * 100.d0
SOFT_TOUCH threshold_davidson
if (N_det > N_det_max) then
call diagonalize_CI
call save_wavefunction
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
N_det = N_det_max
@ -46,22 +45,26 @@ program fci_zmq
double precision :: E_CI_before(N_states)
integer :: n_det_before, to_select
print*,'Beginning the selection ...'
if (.True.) then ! Avoid pre-calculation of CI_energy
E_CI_before(1:N_states) = CI_energy(1:N_states)
endif
n_det_before = 0
double precision :: correlation_energy_ratio
correlation_energy_ratio = E_CI_before(1) - hf_energy_ref
correlation_energy_ratio = correlation_energy_ratio / (correlation_energy_ratio + pt2(1))
correlation_energy_ratio = 0.d0
if (.True.) then ! Avoid pre-calculation of CI_energy
do while ( &
(N_det < N_det_max) .and. &
(maxval(abs(pt2(1:N_states))) > pt2_max) .and. &
(correlation_energy_ratio < correlation_energy_ratio_max) &
(correlation_energy_ratio <= correlation_energy_ratio_max) &
)
correlation_energy_ratio = E_CI_before(1) - hf_energy_ref
correlation_energy_ratio = correlation_energy_ratio / (correlation_energy_ratio + pt2(1))
correlation_energy_ratio = (CI_energy(1) - hf_energy_ref) / &
(E_CI_before(1) + pt2(1) - hf_energy_ref)
correlation_energy_ratio = min(1.d0,correlation_energy_ratio)
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
@ -73,6 +76,7 @@ program fci_zmq
print *, 'E = ', CI_energy(k)
print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k)
enddo
print *, '-----'
if(N_states.gt.1)then
print*,'Variational Energy difference'
@ -87,11 +91,10 @@ program fci_zmq
enddo
endif
E_CI_before(1:N_states) = CI_energy(1:N_states)
call ezfio_set_cas_sd_zmq_energy(CI_energy(1))
n_det_before = N_det
to_select = 2*N_det
to_select = max(64-to_select, to_select)
to_select = N_det
to_select = max(N_det, to_select)
to_select = min(to_select, N_det_max-n_det_before)
call ZMQ_selection(to_select, pt2)
@ -99,18 +102,17 @@ program fci_zmq
PROVIDE psi_det
PROVIDE psi_det_sorted
if (N_det == N_det_max) then
if (N_det >= N_det_max) then
threshold_davidson = threshold_davidson_in
SOFT_TOUCH threshold_davidson
end if
call diagonalize_CI
call save_wavefunction
call ezfio_set_cas_sd_zmq_energy(CI_energy(1))
enddo
endif
if (N_det < N_det_max) then
threshold_davidson = threshold_davidson_in
SOFT_TOUCH threshold_davidson
call diagonalize_CI
call save_wavefunction
call ezfio_set_cas_sd_zmq_energy(CI_energy(1))
@ -147,10 +149,8 @@ program fci_zmq
print *, 'E+PT2 = ', E_CI_before(k)+pt2(k)
print *, '-----'
enddo
call ezfio_set_cas_sd_zmq_energy_pt2(E_CI_before+pt2)
endif
call save_wavefunction
call ezfio_set_cas_sd_zmq_energy(CI_energy(1))
call ezfio_set_cas_sd_zmq_energy_pt2(E_CI_before(1)+pt2(1))
endif
end

4
plugins/CAS_SD_ZMQ/run_selection_slave.irp.f
View File

@ -111,7 +111,7 @@ subroutine push_selection_results(zmq_socket_push, pt2, b, task_id, ntask)
if(rc /= 4*ntask) stop "push"
! Activate is zmq_socket_push is a REQ
rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0)
! rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0)
end subroutine
@ -145,7 +145,7 @@ subroutine pull_selection_results(zmq_socket_pull, pt2, val, det, N, task_id, nt
if(rc /= 4*ntask) stop "pull"
! Activate is zmq_socket_pull is a REP
rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0)
! rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0)
end subroutine

29
plugins/Full_CI/jmmrpt2.irp.f Normal file
View File

@ -0,0 +1,29 @@
program pouet
implicit none
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, degree
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
character*(64) :: perturbation
double precision :: i_H_psi_array(N_states),diag_H_mat_elem,h,i_O1_psi_array(N_states)
double precision :: E_CI_before(N_states)
integer :: n_det_before
threshold_generators = threshold_generators_pt2
threshold_selectors = threshold_selectors_pt2
SOFT_TOUCH threshold_generators threshold_selectors
call H_apply_FCI_PT2_new(pt2, norm_pert, H_pert_diag, N_st)
print *, 'Final step'
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', CI_energy(1:N_states)
print *, 'E+PT2 = ', CI_energy(1:N_states)+pt2(1:N_states)
print *, '-----'
call ezfio_set_full_ci_energy_pt2(CI_energy(1)+pt2(1))
deallocate(pt2,norm_pert)
end

70
plugins/Full_CI_ZMQ/fci_zmq.irp.f
View File

@ -1,8 +1,6 @@
program fci_zmq
implicit none
integer :: i,j,k
logical, external :: detEq
double precision, allocatable :: pt2(:)
integer :: degree
integer :: n_det_before, to_select
@ -12,6 +10,14 @@ program fci_zmq
double precision :: hf_energy_ref
logical :: has
pt2 = -huge(1.d0)
threshold_davidson_in = threshold_davidson
threshold_davidson = threshold_davidson_in * 100.d0
SOFT_TOUCH threshold_davidson
call diagonalize_CI
call save_wavefunction
call ezfio_has_hartree_fock_energy(has)
if (has) then
call ezfio_get_hartree_fock_energy(hf_energy_ref)
@ -19,14 +25,7 @@ program fci_zmq
hf_energy_ref = ref_bitmask_energy
endif
pt2 = 1.d0
threshold_davidson_in = threshold_davidson
threshold_davidson = threshold_davidson_in * 100.d0
SOFT_TOUCH threshold_davidson
if (N_det > N_det_max) then
call diagonalize_CI
call save_wavefunction
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
N_det = N_det_max
@ -47,21 +46,25 @@ program fci_zmq
print*,'Beginning the selection ...'
if (.True.) then ! Avoid pre-calculation of CI_energy
E_CI_before(1:N_states) = CI_energy(1:N_states)
endif
n_det_before = 0
double precision :: correlation_energy_ratio
correlation_energy_ratio = E_CI_before(1) - hf_energy_ref
correlation_energy_ratio = correlation_energy_ratio / (correlation_energy_ratio + pt2(1))
correlation_energy_ratio = 0.d0
if (.True.) then ! Avoid pre-calculation of CI_energy
do while ( &
(N_det < N_det_max) .and. &
(maxval(abs(pt2(1:N_states))) > pt2_max) .and. &
(correlation_energy_ratio < correlation_energy_ratio_max) &
(correlation_energy_ratio <= correlation_energy_ratio_max) &
)
correlation_energy_ratio = E_CI_before(1) - hf_energy_ref
correlation_energy_ratio = correlation_energy_ratio / (correlation_energy_ratio + pt2(1))
correlation_energy_ratio = (CI_energy(1) - hf_energy_ref) / &
(E_CI_before(1) + pt2(1) - hf_energy_ref)
correlation_energy_ratio = min(1.d0,correlation_energy_ratio)
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
@ -88,7 +91,6 @@ program fci_zmq
enddo
endif
E_CI_before(1:N_states) = CI_energy(1:N_states)
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
n_det_before = N_det
to_select = N_det
@ -100,18 +102,17 @@ program fci_zmq
PROVIDE psi_det
PROVIDE psi_det_sorted
if (N_det == N_det_max) then
if (N_det >= N_det_max) then
threshold_davidson = threshold_davidson_in
SOFT_TOUCH threshold_davidson
end if
call diagonalize_CI
call save_wavefunction
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
enddo
endif
if (N_det < N_det_max) then
threshold_davidson = threshold_davidson_in
SOFT_TOUCH threshold_davidson
call diagonalize_CI
call save_wavefunction
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
@ -119,32 +120,35 @@ program fci_zmq
if(do_pt2_end)then
print*,'Last iteration only to compute the PT2'
!threshold_selectors = max(threshold_selectors,threshold_selectors_pt2)
!threshold_generators = max(threshold_generators,threshold_generators_pt2)
!TOUCH threshold_selectors threshold_generators
threshold_selectors = 1.d0
threshold_generators = 1d0
E_CI_before(1:N_states) = CI_energy(1:N_states)
double precision :: relative_error
relative_error=1.d-3
pt2 = 0.d0
call ZMQ_pt2(pt2,relative_error) ! Stochastic PT2
!call ZMQ_selection(0, pt2) ! Deterministic PT2
if (N_states == 1) then
threshold_selectors = 1.d0
threshold_generators = 1d0
SOFT_TOUCH threshold_selectors threshold_generators
print *, 'Stochastic PT2'
call ZMQ_pt2(E_CI_before(1), pt2,relative_error) ! Stochastic PT2
else
threshold_selectors = max(threshold_selectors,threshold_selectors_pt2)
threshold_generators = max(threshold_generators,threshold_generators_pt2)
SOFT_TOUCH threshold_selectors threshold_generators
print *, 'Deterministic PT2'
call ZMQ_selection(0, pt2) ! Deterministic PT2
endif
print *, 'Final step'
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
do k=1,N_states
print *, 'State', k
print *, 'PT2 = ', pt2
print *, 'E = ', E_CI_before
print *, 'E+PT2 = ', E_CI_before+pt2
print *, 'PT2 = ', pt2(k)
print *, 'E = ', E_CI_before(k)
print *, 'E+PT2 = ', E_CI_before(k)+pt2(k)
print *, '-----'
enddo
call ezfio_set_full_ci_zmq_energy(E_CI_before(1))
call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before(1)+pt2(1))
endif
call save_wavefunction
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before(1)+pt2(1))
end

3
plugins/Full_CI_ZMQ/pt2_stoch.irp.f
View File

@ -25,7 +25,8 @@ subroutine run
threshold_selectors = 1.d0
threshold_generators = 1d0
relative_error = 1.d-3
call ZMQ_pt2(pt2, relative_error)
! relative_error = 1.d-8
call ZMQ_pt2(E_CI_before, pt2, relative_error)
print *, 'Final step'
print *, 'N_det = ', N_det
print *, 'PT2 = ', pt2

130
plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
View File

@ -3,7 +3,7 @@ BEGIN_PROVIDER [ integer, fragment_first ]
fragment_first = first_det_of_teeth(1)
END_PROVIDER
subroutine ZMQ_pt2(pt2,relative_error)
subroutine ZMQ_pt2(E, pt2,relative_error)
use f77_zmq
use selection_types
@ -13,7 +13,7 @@ subroutine ZMQ_pt2(pt2,relative_error)
integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_to_qp_run_socket2
type(selection_buffer) :: b
integer, external :: omp_get_thread_num
double precision, intent(in) :: relative_error
double precision, intent(in) :: relative_error, E
double precision, intent(out) :: pt2(N_states)
@ -25,16 +25,14 @@ subroutine ZMQ_pt2(pt2,relative_error)
double precision :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
double precision, external :: omp_get_wtime
double precision :: time0, time
double precision :: time
allocate(pt2_detail(N_states, N_det_generators), comb(N_det_generators/2), computed(N_det_generators), tbc(0:size_tbc))
allocate(pt2_detail(N_states, N_det_generators), comb(N_det_generators), computed(N_det_generators), tbc(0:size_tbc))
sumabove = 0d0
sum2above = 0d0
Nabove = 0d0
provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral pt2_weight
!call random_seed()
provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral pt2_weight psi_selectors
computed = .false.
@ -45,26 +43,20 @@ subroutine ZMQ_pt2(pt2,relative_error)
end do
pt2_detail = 0d0
time0 = omp_get_wtime()
print *, "time - avg - err - n_combs"
generator_per_task = 1
do while(.true.)
print *, '========== ================= ================= ================='
print *, ' Samples Energy Stat. Error Seconds '
print *, '========== ================= ================= ================='
call write_time(6)
call new_parallel_job(zmq_to_qp_run_socket,"pt2")
call new_parallel_job(zmq_to_qp_run_socket,'pt2')
call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator))
call create_selection_buffer(1, 1*2, b)
Ncomb=size(comb)
call get_carlo_workbatch(computed, comb, Ncomb, tbc)
call write_time(6)
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
integer :: ipos
logical :: tasks
tasks = .False.
ipos=1
do i=1,tbc(0)
@ -72,55 +64,39 @@ subroutine ZMQ_pt2(pt2,relative_error)
write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') 0, tbc(i)
ipos += 20
if (ipos > 63980) then
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos-20)))
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos)))
ipos=1
tasks = .True.
endif
else
do j=1,fragment_count
write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') j, tbc(i)
ipos += 20
if (ipos > 63980) then
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos-20)))
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos)))
ipos=1
tasks = .True.
endif
end do
end if
end do
if (ipos > 1) then
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos-20)))
tasks = .True.
call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos)))
endif
if (tasks) then
call zmq_set_running(zmq_to_qp_run_socket)
!$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc+1) &
!$OMP PRIVATE(i)
i = omp_get_thread_num()
if (i==0) then
call pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, pt2)
call pt2_collector(E, b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, pt2)
else
call pt2_slave_inproc(i)
endif
!$OMP END PARALLEL
call delete_selection_buffer(b)
call end_parallel_job(zmq_to_qp_run_socket, 'pt2')
else
pt2 = 0.d0
do i=1,N_det_generators
do k=1,N_states
pt2(k) = pt2(k) + pt2_detail(k,i)
enddo
enddo
endif
tbc(0) = 0
if (pt2(1) /= 0.d0) then
exit
endif
end do
print *, '========== ================= ================= ================='
deallocate(pt2_detail, comb, computed, tbc)
@ -162,7 +138,7 @@ subroutine pt2_slave_inproc(i)
call run_pt2_slave(1,i,pt2_e0_denominator)
end
subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, pt2)
subroutine pt2_collector(E, b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, pt2)
use f77_zmq
use selection_types
use bitmasks
@ -171,7 +147,7 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
integer, intent(in) :: Ncomb
double precision, intent(inout) :: pt2_detail(N_states, N_det_generators)
double precision, intent(in) :: comb(Ncomb), relative_error
double precision, intent(in) :: comb(Ncomb), relative_error, E
logical, intent(inout) :: computed(N_det_generators)
integer, intent(in) :: tbc(0:size_tbc)
double precision, intent(inout) :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
@ -194,14 +170,17 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
integer :: done, Nindex
integer, allocatable :: index(:)
double precision, save :: time0 = -1.d0
double precision :: time, timeLast
double precision :: time, timeLast, Nabove_old
double precision, external :: omp_get_wtime
integer :: tooth, firstTBDcomb, orgTBDcomb
integer, allocatable :: parts_to_get(:)
logical, allocatable :: actually_computed(:)
character*(512) :: task
Nabove_old = -1.d0
allocate(actually_computed(N_det_generators), parts_to_get(N_det_generators), &
pt2_mwen(N_states, N_det_generators) )
pt2_mwen(1:N_states, 1:N_det_generators) =0.d0
do i=1,N_det_generators
actually_computed(i) = computed(i)
enddo
@ -225,12 +204,16 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det_generators), index(1))
more = 1
if (time0 < 0.d0) then
time0 = omp_get_wtime()
call wall_time(time0)
endif
timeLast = time0
print *, 'N_deterministic = ', first_det_of_teeth(1)-1
call get_first_tooth(actually_computed, tooth)
Nabove_old = Nabove(tooth)
! print *, 'N_deterministic = ', first_det_of_teeth(1)-1
pullLoop : do while (more == 1)
call pull_pt2_results(zmq_socket_pull, Nindex, index, pt2_mwen, task_id, ntask)
do i=1,Nindex
pt2_detail(1:N_states, index(i)) += pt2_mwen(1:N_states,i)
@ -253,11 +236,10 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
time = omp_get_wtime()
if(time - timeLast > 1d1 .or. more /= 1) then
if(time - timeLast > 10d0 .or. more /= 1) then
timeLast = time
do i=1, first_det_of_teeth(1)-1
if(.not.(actually_computed(i))) then
print *, "PT2 : deterministic part not finished"
cycle pullLoop
end if
end do
@ -265,7 +247,7 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
double precision :: E0, avg, eqt, prop
call do_carlo(tbc, Ncomb+1-firstTBDcomb, comb(firstTBDcomb), pt2_detail, actually_computed, sumabove, sum2above, Nabove)
firstTBDcomb = int(Nabove(1)) - orgTBDcomb + 1
if(Nabove(1) < 2d0) cycle
if(Nabove(1) < 5d0) cycle
call get_first_tooth(actually_computed, tooth)
done = 0
@ -279,16 +261,30 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
E0 += pt2_detail(1,first_det_of_teeth(tooth)) * prop
avg = E0 + (sumabove(tooth) / Nabove(tooth))
eqt = sqrt(1d0 / (Nabove(tooth)-1) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2))
time = omp_get_wtime()
call wall_time(time)
if (dabs(eqt/avg) < relative_error) then
! Termination
pt2(1) = avg
! exit pullLoop
print '(G10.3, 2X, F16.10, 2X, G16.3, 2X, F16.4, A20)', Nabove(tooth), avg+E, eqt, time-time0, ''
call zmq_abort(zmq_to_qp_run_socket)
else
print "(4(G22.13), 4(I9))", time - time0, avg, eqt, Nabove(tooth), tooth, first_det_of_teeth(tooth)-1, done, first_det_of_teeth(tooth+1)-first_det_of_teeth(tooth)
if (Nabove(tooth) > Nabove_old) then
print '(G10.3, 2X, F16.10, 2X, G16.3, 2X, F16.4, A20)', Nabove(tooth), avg+E, eqt, time-time0, ''
!print "(4(G23.13), 4(I9))", time - time0, avg, eqt, Nabove(tooth), tooth, first_det_of_teeth(tooth)-1, done, first_det_of_teeth(tooth+1)-first_det_of_teeth(tooth)
Nabove_old = Nabove(tooth)
endif
endif
!print "(4(G22.13), 4(I9))", time - time0, avg, eqt, Nabove(tooth), tooth, first_det_of_teeth(tooth)-1, done, first_det_of_teeth(tooth+1)-first_det_of_teeth(tooth)
end if
end do pullLoop
E0 = sum(pt2_detail(1,:first_det_of_teeth(tooth)-1))
prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - pt2_cweight(first_det_of_teeth(tooth)-1))
prop = prop * pt2_weight_inv(first_det_of_teeth(tooth))
E0 += pt2_detail(1,first_det_of_teeth(tooth)) * prop
pt2(1) = E0 + (sumabove(tooth) / Nabove(tooth))
eqt = sqrt(1d0 / (Nabove(tooth)-1) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2))
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_pull_socket(zmq_socket_pull)
call sort_selection_buffer(b)
@ -361,7 +357,7 @@ subroutine get_last_full_tooth(computed, last_tooth)
last_tooth = 0
combLoop : do i=comb_teeth, 1, -1
missing = 1+ ishft(first_det_of_teeth(i+1)-first_det_of_teeth(i),-12) ! /4096
missing = 1+ ishft(first_det_of_teeth(i+1)-first_det_of_teeth(i),-4) ! /16
do j=first_det_of_teeth(i), first_det_of_teeth(i+1)-1
if(.not.computed(j)) then
missing -= 1
@ -379,7 +375,7 @@ BEGIN_PROVIDER [ integer, size_tbc ]
BEGIN_DOC
! Size of the tbc array
END_DOC
size_tbc = N_det_generators + fragment_count*fragment_first
size_tbc = 2*N_det_generators + fragment_count*fragment_first
END_PROVIDER
subroutine get_carlo_workbatch(computed, comb, Ncomb, tbc)
@ -388,30 +384,25 @@ subroutine get_carlo_workbatch(computed, comb, Ncomb, tbc)
double precision, intent(out) :: comb(Ncomb)
integer, intent(inout) :: tbc(0:size_tbc)
logical, intent(inout) :: computed(N_det_generators)
integer :: i, j, last_full, dets(comb_teeth), tbc_save
integer :: i, j, last_full, dets(comb_teeth)
integer :: icount, n
n = tbc(0)
icount = 0
integer :: k, l
l=1
call RANDOM_NUMBER(comb)
do i=1,size(comb)
comb(i) = comb(i) * comb_step
tbc_save = tbc(0)
!DIR$ FORCEINLINE
call add_comb(comb(i), computed, tbc, size_tbc, comb_teeth)
if (tbc(0) < size(tbc)) then
Ncomb = i
else
tbc(0) = tbc_save
return
endif
icount = icount + tbc(0) - tbc_save
if ((i>1000).and.(icount > n)) then
call get_filling_teeth(computed, tbc)
icount = 0
n = ishft(tbc_save,-4)
endif
do while (computed(l))
l=l+1
if (l == N_det_generators+1) return
enddo
k=tbc(0)+1
tbc(k) = l
computed(l) = .True.
tbc(0) = k
enddo
call get_filling_teeth(computed, tbc)
end subroutine
@ -535,12 +526,13 @@ end subroutine
comb_step = 1d0/dfloat(comb_teeth)
first_det_of_comb = 1
do i=1,N_det_generators
if(pt2_weight(i)/norm_left < comb_step*.5d0) then
if(pt2_weight(i)/norm_left < .5d0*comb_step) then
first_det_of_comb = i
exit
end if
norm_left -= pt2_weight(i)
end do
call write_int(6, first_det_of_comb-1, 'Size of deterministic set')
comb_step = (1d0 - pt2_cweight(first_det_of_comb-1)) * comb_step

30
plugins/Full_CI_ZMQ/run_pt2_slave.irp.f
View File

@ -17,7 +17,7 @@ subroutine run_pt2_slave(thread,iproc,energy)
integer(ZMQ_PTR), external :: new_zmq_push_socket
integer(ZMQ_PTR) :: zmq_socket_push
type(selection_buffer) :: buf, buf2
type(selection_buffer) :: buf
logical :: done
double precision :: pt2(N_states)
@ -47,18 +47,12 @@ subroutine run_pt2_slave(thread,iproc,energy)
if (done) then
ctask = ctask - 1
else
integer :: i_generator, i_i_generator, N, subset
integer :: i_generator, i_i_generator, subset
read (task,*) subset, index
!!!!!
N=1
!!!!!
if(buf%N == 0) then
! Only first time
call create_selection_buffer(N, N*2, buf)
call create_selection_buffer(N, N*3, buf2)
else
if(N /= buf%N) stop "N changed... wtf man??"
call create_selection_buffer(1, 2, buf)
end if
do i_i_generator=1, Nindex
i_generator = index
@ -67,18 +61,13 @@ subroutine run_pt2_slave(thread,iproc,energy)
enddo
endif
if(done .or. ctask == size(task_id)) then
if(done .or. (ctask == size(task_id)) ) then
if(buf%N == 0 .and. ctask > 0) stop "uninitialized selection_buffer"
do i=1, ctask
call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i))
end do
if(ctask > 0) then
call push_pt2_results(zmq_socket_push, Nindex, index, pt2_detail, task_id(1), ctask)
do i=1,buf%cur
call add_to_selection_buffer(buf2, buf%det(1,1,i), buf%val(i))
enddo
call sort_selection_buffer(buf2)
buf%mini = buf2%mini
pt2 = 0d0
pt2_detail(:,:Nindex) = 0d0
buf%cur = 0
@ -92,6 +81,7 @@ subroutine run_pt2_slave(thread,iproc,energy)
call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id)
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_push_socket(zmq_socket_push,thread)
call delete_selection_buffer(buf)
end subroutine
@ -123,8 +113,8 @@ subroutine push_pt2_results(zmq_socket_push, N, index, pt2_detail, task_id, ntas
if(rc /= 4*ntask) stop "push"
! Activate is zmq_socket_push is a REQ
character*(2) :: ok
rc = f77_zmq_recv( zmq_socket_push, ok, 2, 0)
! character*(2) :: ok
! rc = f77_zmq_recv( zmq_socket_push, ok, 2, 0)
end subroutine
@ -154,11 +144,7 @@ subroutine pull_pt2_results(zmq_socket_pull, N, index, pt2_detail, task_id, ntas
if(rc /= 4*ntask) stop "pull"
! Activate is zmq_socket_pull is a REP
rc = f77_zmq_send( zmq_socket_pull, 'ok', 2, 0)
do i=N+1,N_det_generators
pt2_detail(1:N_states,i) = 0.d0
enddo
! rc = f77_zmq_send( zmq_socket_pull, 'ok', 2, 0)
end subroutine

18
plugins/Full_CI_ZMQ/run_selection_slave.irp.f
View File

@ -58,13 +58,9 @@ subroutine run_selection_slave(thread,iproc,energy)
call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i))
end do
if(ctask > 0) then
call sort_selection_buffer(buf)
call push_selection_results(zmq_socket_push, pt2, buf, task_id(1), ctask)
do i=1,buf%cur
call add_to_selection_buffer(buf2, buf%det(1,1,i), buf%val(i))
if (buf2%cur == buf2%N) then
call sort_selection_buffer(buf2)
endif
enddo
call merge_selection_buffers(buf,buf2)
buf%mini = buf2%mini
pt2 = 0d0
buf%cur = 0
@ -78,6 +74,10 @@ subroutine run_selection_slave(thread,iproc,energy)
call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id)
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_push_socket(zmq_socket_push,thread)
if (buf%N > 0) then
call delete_selection_buffer(buf)
call delete_selection_buffer(buf2)
endif
end subroutine
@ -92,8 +92,6 @@ subroutine push_selection_results(zmq_socket_push, pt2, b, task_id, ntask)
integer, intent(in) :: ntask, task_id(*)
integer :: rc
call sort_selection_buffer(b)
rc = f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE)
if(rc /= 4) stop "push"
rc = f77_zmq_send( zmq_socket_push, pt2, 8*N_states, ZMQ_SNDMORE)
@ -112,7 +110,7 @@ subroutine push_selection_results(zmq_socket_push, pt2, b, task_id, ntask)
if(rc /= 4*ntask) stop "push"
! Activate is zmq_socket_push is a REQ
rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0)
! rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0)
end subroutine
@ -146,7 +144,7 @@ subroutine pull_selection_results(zmq_socket_pull, pt2, val, det, N, task_id, nt
if(rc /= 4*ntask) stop "pull"
! Activate is zmq_socket_pull is a REP
rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0)
! rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0)
end subroutine

10
plugins/Full_CI_ZMQ/selection.irp.f
View File

@ -509,7 +509,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
logical :: ok
integer :: s1, s2, p1, p2, ib, j, istate
integer(bit_kind) :: mask(N_int, 2), det(N_int, 2)
double precision :: e_pert, delta_E, val, Hii, max_e_pert,tmp
double precision :: e_pert, delta_E, val, Hii, min_e_pert,tmp
double precision, external :: diag_H_mat_elem_fock
logical, external :: detEq
@ -536,7 +536,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
max_e_pert = 0d0
min_e_pert = 0d0
do istate=1,N_states
delta_E = E0(istate) - Hii
@ -547,12 +547,12 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
endif
e_pert = 0.5d0 * (tmp - delta_E)
pt2(istate) = pt2(istate) + e_pert
max_e_pert = min(e_pert,max_e_pert)
min_e_pert = min(e_pert,min_e_pert)
! ci(istate) = e_pert / mat(istate, p1, p2)
end do
if(dabs(max_e_pert) > buf%mini) then
call add_to_selection_buffer(buf, det, max_e_pert)
if(min_e_pert <= buf%mini) then
call add_to_selection_buffer(buf, det, min_e_pert)
end if
end do
end do

80
plugins/Full_CI_ZMQ/selection_buffer.irp.f
View File

@ -15,6 +15,18 @@ subroutine create_selection_buffer(N, siz, res)
res%cur = 0
end subroutine
subroutine delete_selection_buffer(b)
use selection_types
implicit none
type(selection_buffer), intent(inout) :: b
if (associated(b%det)) then
deallocate(b%det)
endif
if (associated(b%val)) then
deallocate(b%val)
endif
end
subroutine add_to_selection_buffer(b, det, val)
use selection_types
@ -25,7 +37,7 @@ subroutine add_to_selection_buffer(b, det, val)
double precision, intent(in) :: val
integer :: i
if(dabs(val) >= b%mini) then
if(val <= b%mini) then
b%cur += 1
b%det(1:N_int,1:2,b%cur) = det(1:N_int,1:2)
b%val(b%cur) = val
@ -35,39 +47,79 @@ subroutine add_to_selection_buffer(b, det, val)
end if
end subroutine
subroutine merge_selection_buffers(b1, b2)
use selection_types
implicit none
BEGIN_DOC
! Merges the selection buffers b1 and b2 into b2
END_DOC
type(selection_buffer), intent(in) :: b1
type(selection_buffer), intent(inout) :: b2
integer(bit_kind), pointer :: detmp(:,:,:)
double precision, pointer :: val(:)
integer :: i, i1, i2, k, nmwen
nmwen = min(b1%N, b1%cur+b2%cur)
allocate( val(size(b1%val)), detmp(N_int, 2, size(b1%det,3)) )
i1=1
i2=1
do i=1,nmwen
if ( (i1 > b1%cur).and.(i2 > b2%cur) ) then
exit
else if (i1 > b1%cur) then
val(i) = b2%val(i2)
detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2)
detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2)
i2=i2+1
else if (i2 > b2%cur) then
val(i) = b1%val(i1)
detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1)
detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1)
i1=i1+1
else
if (b1%val(i1) <= b2%val(i2)) then
val(i) = b1%val(i1)
detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1)
detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1)
i1=i1+1
else
val(i) = b2%val(i2)
detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2)
detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2)
i2=i2+1
endif
endif
enddo
deallocate(b2%det, b2%val)
b2%det => detmp
b2%val => val
b2%mini = min(b2%mini,b2%val(b2%N))
b2%cur = nmwen
end
subroutine sort_selection_buffer(b)
use selection_types
implicit none
type(selection_buffer), intent(inout) :: b
double precision, allocatable:: absval(:)
integer, allocatable :: iorder(:)
double precision, pointer :: vals(:)
integer(bit_kind), pointer :: detmp(:,:,:)
integer :: i, nmwen
logical, external :: detEq
nmwen = min(b%N, b%cur)
allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3)), absval(b%cur), vals(size(b%val)))
absval = -dabs(b%val(:b%cur))
allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3)))
do i=1,b%cur
iorder(i) = i
end do
call dsort(absval, iorder, b%cur)
call dsort(b%val, iorder, b%cur)
do i=1, nmwen
detmp(1:N_int,1,i) = b%det(1:N_int,1,iorder(i))
detmp(1:N_int,2,i) = b%det(1:N_int,2,iorder(i))
vals(i) = b%val(iorder(i))
end do
do i=nmwen+1, size(vals)
vals(i) = 0.d0
enddo
deallocate(b%det, b%val)
deallocate(b%det,iorder)
b%det => detmp
b%val => vals
b%mini = max(b%mini,dabs(b%val(b%N)))
b%mini = min(b%mini,b%val(b%N))
b%cur = nmwen
end subroutine

2
plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f
View File

@ -37,7 +37,7 @@ subroutine run_wf
do
call wait_for_states(states,zmq_state,2)
call wait_for_states(states,zmq_state,3)
if(trim(zmq_state) == 'Stopped') then

32
plugins/Full_CI_ZMQ/zmq_selection.irp.f
View File

@ -45,7 +45,7 @@ subroutine ZMQ_selection(N_in, pt2)
!$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc+1)
i = omp_get_thread_num()
if (i==0) then
call selection_collector(b, pt2)
call selection_collector(b, N, pt2)
else
call selection_slave_inproc(i)
endif
@ -59,6 +59,7 @@ subroutine ZMQ_selection(N_in, pt2)
endif
call save_wavefunction
endif
call delete_selection_buffer(b)
end subroutine
@ -70,7 +71,7 @@ subroutine selection_slave_inproc(i)
call run_selection_slave(1,i,pt2_e0_denominator)
end
subroutine selection_collector(b, pt2)
subroutine selection_collector(b, N, pt2)
use f77_zmq
use selection_types
use bitmasks
@ -78,6 +79,7 @@ subroutine selection_collector(b, pt2)
type(selection_buffer), intent(inout) :: b
integer, intent(in) :: N
double precision, intent(out) :: pt2(N_states)
double precision :: pt2_mwen(N_states)
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
@ -87,40 +89,38 @@ subroutine selection_collector(b, pt2)
integer(ZMQ_PTR) :: zmq_socket_pull
integer :: msg_size, rc, more
integer :: acc, i, j, robin, N, ntask
double precision, allocatable :: val(:)
integer(bit_kind), allocatable :: det(:,:,:)
integer :: acc, i, j, robin, ntask
double precision, pointer :: val(:)
integer(bit_kind), pointer :: det(:,:,:)
integer, allocatable :: task_id(:)
integer :: done
real :: time, time0
type(selection_buffer) :: b2
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
zmq_socket_pull = new_zmq_pull_socket()
allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det_generators))
done = 0
call create_selection_buffer(N, N*2, b2)
allocate(task_id(N_det_generators))
more = 1
pt2(:) = 0d0
call CPU_TIME(time0)
do while (more == 1)
call pull_selection_results(zmq_socket_pull, pt2_mwen, val(1), det(1,1,1), N, task_id, ntask)
pt2 += pt2_mwen
do i=1, N
call add_to_selection_buffer(b, det(1,1,i), val(i))
end do
call pull_selection_results(zmq_socket_pull, pt2_mwen, b2%val(1), b2%det(1,1,1), b2%cur, task_id, ntask)
pt2 += pt2_mwen
call merge_selection_buffers(b2,b)
do i=1, ntask
if(task_id(i) == 0) then
print *, "Error in collector"
endif
call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id(i),more)
end do
done += ntask
call CPU_TIME(time)
! print *, "DONE" , done, time - time0
end do
call delete_selection_buffer(b2)
call sort_selection_buffer(b)
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_pull_socket(zmq_socket_pull)
call sort_selection_buffer(b)
end subroutine

24
plugins/mrcepa0/dressing_slave.irp.f
View File

@ -316,12 +316,12 @@ subroutine push_mrsc2_results(zmq_socket_push, I_i, J, delta, delta_s2, task_id)
endif
! Activate is zmq_socket_push is a REQ
integer :: idummy
rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
if (rc /= 4) then
print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
stop 'error'
endif
! integer :: idummy
! rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
! if (rc /= 4) then
! print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
! stop 'error'
! endif
end
@ -390,12 +390,12 @@ subroutine pull_mrsc2_results(zmq_socket_pull, I_i, J, n, idx, delta, delta_s2,
! Activate is zmq_socket_pull is a REP
integer :: idummy
rc = f77_zmq_send( zmq_socket_pull, idummy, 4, 0)
if (rc /= 4) then
print *, irp_here, 'f77_zmq_send( zmq_socket_pull, idummy, 4, 0)'
stop 'error'
endif
! integer :: idummy
! rc = f77_zmq_send( zmq_socket_pull, idummy, 4, 0)
! if (rc /= 4) then
! print *, irp_here, 'f77_zmq_send( zmq_socket_pull, idummy, 4, 0)'
! stop 'error'
! endif
end

21
promela/collector.pml Normal file
View File

@ -0,0 +1,21 @@
proctype collector(byte state) {
byte task;
req_message msg;
rep_message reply;
bit loop = 1;
xr pull_socket;
do
:: (loop == 0) -> break
:: else ->
pull_socket ? task;
/* Handle result */
send_req(DELTASK, task);
assert (reply.m == OK);
loop = reply.value;
od;
}

48
promela/fortran.pml Normal file
View File

@ -0,0 +1,48 @@
active proctype fortran() {
req_message msg;
rep_message reply;
byte state;
byte count, wait;
/* New parallel job */
state=1;
send_req( NEWJOB, state );
assert (reply.m == OK);
send_req( PUTPSI, state );
assert (reply.m == PUTPSI_REPLY);
/* Add tasks */
count = 0;
do
:: (count == NTASKS) -> break;
:: else ->
count++;
send_req( ADDTASK, count );
assert (reply.m == OK);
od
wait = _nr_pr;
/* Run collector */
run collector(state);
/* Run slaves */
count = 0;
do
:: (count == NPROC) -> break;
:: else -> count++; run slave();
od
/* Wait for collector and slaves to finish */
(_nr_pr == wait);
send_req( ENDJOB, state );
assert (reply.m == OK);
state = reply.value;
send_req( TERMINATE, 0);
assert (reply.m == OK);
}

35
promela/model.pml Normal file
View File

@ -0,0 +1,35 @@
#define NPROC 3
#define BUFSIZE 2
#define NTASKS 4
mtype = { NONE, OK, WRONG_STATE, TERMINATE, GETPSI, PUTPSI, NEWJOB, ENDJOB, SETRUNNING,
SETWAITING, SETSTOPPED, CONNECT, DISCONNECT, ADDTASK, DELTASK, TASKDONE, GETTASK,
PSI, TASK, PUTPSI_REPLY, WAITING, RUNNING, STOPPED
}
#define send_req( MESSAGE, VALUE ) atomic { msg.m=MESSAGE ; msg.value=VALUE ; msg.state=state; } ; rep_socket ! msg; msg.reply ? reply
/* Request/Reply pattern */
typedef rep_message {
mtype m = NONE;
byte value = 0;
}
typedef req_message {
mtype m = NONE;
byte state = 0;
byte value = 0;
chan reply = [BUFSIZE] of { rep_message };
}
/* Channels */
chan rep_socket = [NPROC] of { req_message };
chan pull_socket = [NPROC] of { byte };
#include "task_server.pml"
#include "fortran.pml"
#include "collector.pml"
#include "slave.pml"

29
promela/slave.pml Normal file
View File

@ -0,0 +1,29 @@
proctype slave() {
req_message msg;
rep_message reply;
byte task;
byte state;
send_req(CONNECT, 0);
assert (reply.m == OK);
state = reply.value;
send_req(GETPSI, 0);
assert (reply.m == PSI);
task=255;
do
:: (task == 0) -> break;
:: else ->
send_req( GETTASK, 0);
if
:: (reply.m == NONE) -> task = 0;
:: (reply.m == TASK) ->
task = reply.value;
/* Compute task */
send_req( TASKDONE, task);
assert (reply.m == OK);
pull_socket ! task;
fi
od
}

138
promela/task_server.pml Normal file
View File

@ -0,0 +1,138 @@
/* State of the task server */
typedef state_t {
chan queue = [NTASKS+2] of { byte };
byte state = 0;
bit address_tcp = 0;
bit address_inproc = 0;
bit psi = 0;
bit running = 0;
byte ntasks;
byte nclients = 0;
}
active proctype task_server() {
xr rep_socket;
state_t state;
req_message msg;
rep_message reply;
byte task;
state.running = 1;
do
:: ( state.running + state.nclients == 0 ) -> break
:: else ->
rep_socket ? msg;
printf("req: "); printm(msg.m); printf("\t%d\n",msg.value);
if
:: ( msg.m == TERMINATE ) ->
atomic {
assert (state.state == 0);
assert (msg.state == state.state);
state.running = 0;
reply.m = OK;
}
:: ( msg.m == CONNECT ) ->
atomic {
assert (state.state != 0);
state.nclients++;
reply.m = OK;
reply.value = state.state;
}
/*
:: ( msg.m == DISCONNECT ) ->
atomic {
assert (state.state != 0);
assert (msg.state == state.state);
state.nclients--;
reply.m = OK;
}
*/
:: ( msg.m == PUTPSI ) ->
atomic {
assert (state.state != 0);
assert (msg.state == state.state);
assert (state.psi == 0);
state.psi = 1;
reply.m = PUTPSI_REPLY;
}
:: ( msg.m == GETPSI ) ->
atomic {
assert (state.state != 0);
assert (msg.state == state.state);
assert (state.psi == 1);
reply.m = PSI;
}
:: ( msg.m == NEWJOB ) ->
atomic {
assert (state.state == 0);
state.state = msg.value;
reply.m = OK;
reply.value = state.state;
}
:: ( msg.m == ENDJOB ) ->
atomic {
assert (state.state != 0);
assert (msg.state == state.state);
state.state = 0;
reply.m = OK;
}
:: ( msg.m == TASKDONE ) ->
atomic {
assert (state.state != 0);
assert (state.ntasks > 0);
assert (msg.state == state.state);
reply.m = OK;
}
:: ( msg.m == GETTASK ) ->
assert (state.state != 0);
assert (state.nclients > 0);
assert (msg.state == state.state);
if
:: ( state.queue ?[task] ) ->
reply.m = TASK;
state.queue ? reply.value
:: else ->
atomic {
reply.m = NONE;
reply.value = 0;
state.nclients--;
}
fi;
:: ( msg.m == DELTASK ) ->
assert (state.state != 0);
assert (msg.state == state.state);
state.ntasks--;
if
:: (state.ntasks > 0) -> reply.value = 1;
:: else -> reply.value = 0;
fi;
reply.m = OK;
:: ( msg.m == ADDTASK ) ->
assert (state.state != 0);
assert (msg.state == state.state);
atomic {
state.ntasks++;
reply.m = OK;
}
state.queue ! msg.value;
fi;
msg.reply ! reply;
printf("rep: "); printm(reply.m); printf("\t%d\n",reply.value);
od;
}

3
src/Bitmask/bitmasks.irp.f
View File

@ -315,10 +315,7 @@ BEGIN_PROVIDER [ integer(bit_kind), cas_bitmask, (N_int,2,N_cas_bitmask) ]
call ezfio_has_bitmasks_cas(exists)
if (exists) then
print*,'---------------------'
print*,'CAS BITMASK RESTART'
call ezfio_get_bitmasks_cas(cas_bitmask)
print*,'---------------------'
else
if(N_generators_bitmask == 1)then
do j=1, N_cas_bitmask

31
src/Davidson/davidson_parallel.irp.f
View File

@ -70,6 +70,11 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, N_st, sze,
! -----------------------
integer :: rc
integer :: N_states_read, N_det_read, psi_det_size_read
integer :: N_det_selectors_read, N_det_generators_read
double precision :: energy(N_st)
write(msg, *) 'get_psi ', worker_id
rc = f77_zmq_send(zmq_to_qp_run_socket,trim(msg),len(trim(msg)),0)
if (rc /= len(trim(msg))) then
@ -84,10 +89,6 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, N_st, sze,
stop 'error'
endif
integer :: N_states_read, N_det_read, psi_det_size_read
integer :: N_det_selectors_read, N_det_generators_read
double precision :: energy(N_st)
read(msg(14:rc),*) rc, N_states_read, N_det_read, psi_det_size_read, &
N_det_generators_read, N_det_selectors_read
@ -167,12 +168,12 @@ subroutine davidson_push_results(zmq_socket_push, v_0, s_0, task_id)
if(rc /= 4) stop "davidson_push_results failed to push task_id"
! Activate is zmq_socket_push is a REQ
integer :: idummy
rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
if (rc /= 4) then
print *, irp_here, ': f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
stop 'error'
endif
! integer :: idummy
! rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
! if (rc /= 4) then
! print *, irp_here, ': f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
! stop 'error'
! endif
end subroutine
@ -199,11 +200,11 @@ subroutine davidson_pull_results(zmq_socket_pull, v_0, s_0, task_id)
if(rc /= 4) stop "davidson_pull_results failed to pull task_id"
! Activate if zmq_socket_pull is a REP
rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
if (rc /= 4) then
print *, irp_here, ' : f77_zmq_send (zmq_socket_pull,...'
stop 'error'
endif
! rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
! if (rc /= 4) then
! print *, irp_here, ' : f77_zmq_send (zmq_socket_pull,...'
! stop 'error'
! endif
end subroutine

22
src/Determinants/H_apply.irp.f
View File

@ -362,12 +362,12 @@ subroutine push_pt2(zmq_socket_push,pt2,norm_pert,H_pert_diag,i_generator,N_st,t
endif
! Activate if zmq_socket_push is a REQ
integer :: idummy
rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
if (rc /= 4) then
print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
stop 'error'
endif
! integer :: idummy
! rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
! if (rc /= 4) then
! print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
! stop 'error'
! endif
end
subroutine pull_pt2(zmq_socket_pull,pt2,norm_pert,H_pert_diag,i_generator,N_st,n,task_id)
@ -433,11 +433,11 @@ subroutine pull_pt2(zmq_socket_pull,pt2,norm_pert,H_pert_diag,i_generator,N_st,n
endif
! Activate if zmq_socket_pull is a REP
rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
if (rc /= 4) then
print *, irp_here, 'f77_zmq_send( zmq_socket_pull, 0, 4, 0)'
stop 'error'
endif
! rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
! if (rc /= 4) then
! print *, irp_here, 'f77_zmq_send( zmq_socket_pull, 0, 4, 0)'
! stop 'error'
! endif
end

22
src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f
View File

@ -57,12 +57,12 @@ subroutine push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value,
endif
! Activate is zmq_socket_push is a REQ
integer :: idummy
rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
if (rc /= 4) then
print *, irp_here, ': f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
stop 'error'
endif
! integer :: idummy
! rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0)
! if (rc /= 4) then
! print *, irp_here, ': f77_zmq_send( zmq_socket_push, idummy, 4, 0)'
! stop 'error'
! endif
end
@ -187,11 +187,11 @@ subroutine ao_bielec_integrals_in_map_collector
rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)
! Activate if zmq_socket_pull is a REP
rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
if (rc /= 4) then
print *, irp_here, ' : f77_zmq_send (zmq_socket_pull,...'
stop 'error'
endif
! rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
! if (rc /= 4) then
! print *, irp_here, ' : f77_zmq_send (zmq_socket_pull,...'
! stop 'error'
! endif
call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)

67
src/Utils/sort.irp.f
View File

@ -27,6 +27,56 @@ BEGIN_TEMPLATE
enddo
end subroutine insertion_$Xsort
subroutine quick_$Xsort(x, iorder, isize)
implicit none
BEGIN_DOC
! Sort array x(isize) using the quicksort algorithm.
! iorder in input should be (1,2,3,...,isize), and in output
! contains the new order of the elements.
END_DOC
integer,intent(in) :: isize
$type,intent(inout) :: x(isize)
integer,intent(inout) :: iorder(isize)
call rec_$X_quicksort(x,iorder,isize,1,isize)
end
recursive subroutine rec_$X_quicksort(x, iorder, isize, first, last)
implicit none
integer, intent(in) :: isize, first, last
integer,intent(inout) :: iorder(isize)
$type, intent(inout) :: x(isize)
$type :: c, tmp
integer :: itmp
integer :: i, j
c = x( ishft(first+last,-1) )
i = first
j = last
do
do while (x(i) < c)
i=i+1
end do
do while (c < x(j))
j=j-1
end do
if (i >= j) exit
tmp = x(i)
x(i) = x(j)
x(j) = tmp
itmp = iorder(i)
iorder(i) = iorder(j)
iorder(j) = itmp
i=i+1
j=j-1
enddo
if (first < i-1) then
call rec_$X_quicksort(x, iorder, isize, first, i-1)
endif
if (j+1 < last) then
call rec_$X_quicksort(x, iorder, isize, j+1, last)
endif
end
subroutine heap_$Xsort(x,iorder,isize)
implicit none
BEGIN_DOC
@ -202,14 +252,15 @@ BEGIN_TEMPLATE
if (isize < 2) then
return
endif
call sorted_$Xnumber(x,isize,n)
if (isize == n) then
return
endif
if ( isize < 32+n) then
! call sorted_$Xnumber(x,isize,n)
! if (isize == n) then
! return
! endif
if ( isize < 32) then
call insertion_$Xsort(x,iorder,isize)
else
call heap_$Xsort(x,iorder,isize)
! call heap_$Xsort(x,iorder,isize)
call quick_$Xsort(x,iorder,isize)
endif
end subroutine $Xsort
@ -359,6 +410,10 @@ BEGIN_TEMPLATE
integer :: err
!DIR$ ATTRIBUTES ALIGN : 128 :: iorder1,iorder2, x2, x1
if (isize < 2) then
return
endif
if (iradix == -1) then ! Sort Positive and negative
allocate(x1(isize),iorder1(isize), x2(isize),iorder2(isize),stat=err)

141
src/ZMQ/utils.irp.f
View File

@ -140,15 +140,15 @@ function new_zmq_to_qp_run_socket()
stop 'Unable to create zmq req socket'
endif
rc = f77_zmq_setsockopt(new_zmq_to_qp_run_socket, ZMQ_SNDTIMEO, 120000, 4)
if (rc /= 0) then
stop 'Unable to set send timeout in new_zmq_to_qp_run_socket'
endif
rc = f77_zmq_setsockopt(new_zmq_to_qp_run_socket, ZMQ_RCVTIMEO, 120000, 4)
if (rc /= 0) then
stop 'Unable to set recv timeout in new_zmq_to_qp_run_socket'
endif
! rc = f77_zmq_setsockopt(new_zmq_to_qp_run_socket, ZMQ_SNDTIMEO, 120000, 4)
! if (rc /= 0) then
! stop 'Unable to set send timeout in new_zmq_to_qp_run_socket'
! endif
!
! rc = f77_zmq_setsockopt(new_zmq_to_qp_run_socket, ZMQ_RCVTIMEO, 120000, 4)
! if (rc /= 0) then
! stop 'Unable to set recv timeout in new_zmq_to_qp_run_socket'
! endif
rc = f77_zmq_connect(new_zmq_to_qp_run_socket, trim(qp_run_address)//':'//trim(zmq_port(0)))
if (rc /= 0) then
@ -180,14 +180,14 @@ function new_zmq_pair_socket(bind)
endif
rc = f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_SNDHWM, 1, 4)
rc = f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_SNDHWM, 4, 4)
if (rc /= 0) then
stop 'f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_SNDHWM, 1, 4)'
stop 'f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_SNDHWM, 4, 4)'
endif
rc = f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_RCVHWM, 1, 4)
rc = f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_RCVHWM, 4, 4)
if (rc /= 0) then
stop 'f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_RCVHWM, 1, 4)'
stop 'f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_RCVHWM, 4, 4)'
endif
rc = f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_IMMEDIATE, 1, 4)
@ -232,24 +232,24 @@ function new_zmq_pull_socket()
if (zmq_context == 0_ZMQ_PTR) then
stop 'zmq_context is uninitialized'
endif
! new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_PULL)
new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_REP)
new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_PULL)
! new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_REP)
call omp_unset_lock(zmq_lock)
if (new_zmq_pull_socket == 0_ZMQ_PTR) then
stop 'Unable to create zmq pull socket'
endif
rc = f77_zmq_setsockopt(new_zmq_pull_socket,ZMQ_LINGER,300000,4)
if (rc /= 0) then
stop 'Unable to set ZMQ_LINGER on pull socket'
endif
! rc = f77_zmq_setsockopt(new_zmq_pull_socket,ZMQ_LINGER,300000,4)
! if (rc /= 0) then
! stop 'Unable to set ZMQ_LINGER on pull socket'
! endif
!
rc = f77_zmq_setsockopt(new_zmq_pull_socket,ZMQ_RCVBUF,100000000,4)
if (rc /= 0) then
stop 'Unable to set ZMQ_RCVBUF on pull socket'
endif
rc = f77_zmq_setsockopt(new_zmq_pull_socket,ZMQ_RCVHWM,1,4)
rc = f77_zmq_setsockopt(new_zmq_pull_socket,ZMQ_RCVHWM,4,4)
if (rc /= 0) then
stop 'Unable to set ZMQ_RCVHWM on pull socket'
endif
@ -309,26 +309,26 @@ function new_zmq_push_socket(thread)
if (zmq_context == 0_ZMQ_PTR) then
stop 'zmq_context is uninitialized'
endif
! new_zmq_push_socket = f77_zmq_socket(zmq_context, ZMQ_PUSH)
new_zmq_push_socket = f77_zmq_socket(zmq_context, ZMQ_REQ)
new_zmq_push_socket = f77_zmq_socket(zmq_context, ZMQ_PUSH)
! new_zmq_push_socket = f77_zmq_socket(zmq_context, ZMQ_REQ)
call omp_unset_lock(zmq_lock)
if (new_zmq_push_socket == 0_ZMQ_PTR) then
stop 'Unable to create zmq push socket'
endif
rc = f77_zmq_setsockopt(new_zmq_push_socket,ZMQ_LINGER,300000,4)
if (rc /= 0) then
stop 'Unable to set ZMQ_LINGER on push socket'
endif
! rc = f77_zmq_setsockopt(new_zmq_push_socket,ZMQ_LINGER,300000,4)
! if (rc /= 0) then
! stop 'Unable to set ZMQ_LINGER on push socket'
! endif
rc = f77_zmq_setsockopt(new_zmq_push_socket,ZMQ_SNDHWM,1,4)
rc = f77_zmq_setsockopt(new_zmq_push_socket,ZMQ_SNDHWM,4,4)
if (rc /= 0) then
stop 'Unable to set ZMQ_SNDHWM on push socket'
endif
rc = f77_zmq_setsockopt(new_zmq_push_socket,ZMQ_SNDBUF,100000000,4)
if (rc /= 0) then
stop 'Unable to set ZMQ_RCVBUF on push socket'
stop 'Unable to set ZMQ_SNDBUF on push socket'
endif
rc = f77_zmq_setsockopt(new_zmq_push_socket,ZMQ_IMMEDIATE,1,4)
@ -336,10 +336,10 @@ function new_zmq_push_socket(thread)
stop 'Unable to set ZMQ_IMMEDIATE on push socket'
endif
rc = f77_zmq_setsockopt(new_zmq_push_socket, ZMQ_SNDTIMEO, 100000, 4)
if (rc /= 0) then
stop 'Unable to set send timout in new_zmq_push_socket'
endif
! rc = f77_zmq_setsockopt(new_zmq_push_socket, ZMQ_SNDTIMEO, 100000, 4)
! if (rc /= 0) then
! stop 'Unable to set send timout in new_zmq_push_socket'
! endif
if (thread == 1) then
rc = f77_zmq_connect(new_zmq_push_socket, zmq_socket_push_inproc_address)
@ -373,10 +373,10 @@ function new_zmq_sub_socket()
stop 'Unable to create zmq sub socket'
endif
rc = f77_zmq_setsockopt(new_zmq_sub_socket,ZMQ_RCVTIMEO,10000,4)
if (rc /= 0) then
stop 'Unable to set timeout in new_zmq_sub_socket'
endif
! rc = f77_zmq_setsockopt(new_zmq_sub_socket,ZMQ_RCVTIMEO,10000,4)
! if (rc /= 0) then
! stop 'Unable to set timeout in new_zmq_sub_socket'
! endif
rc = f77_zmq_setsockopt(new_zmq_sub_socket,ZMQ_CONFLATE,1,4)
if (rc /= 0) then
@ -445,10 +445,10 @@ subroutine end_zmq_pull_socket(zmq_socket_pull)
integer :: rc
character*(8), external :: zmq_port
rc = f77_zmq_setsockopt(zmq_socket_pull,ZMQ_LINGER,0,4)
if (rc /= 0) then
stop 'Unable to set ZMQ_LINGER on pull socket'
endif
! rc = f77_zmq_setsockopt(zmq_socket_pull,ZMQ_LINGER,0,4)
! if (rc /= 0) then
! stop 'Unable to set ZMQ_LINGER on pull socket'
! endif
call omp_set_lock(zmq_lock)
rc = f77_zmq_close(zmq_socket_pull)
@ -472,10 +472,10 @@ subroutine end_zmq_push_socket(zmq_socket_push,thread)
integer :: rc
character*(8), external :: zmq_port
rc = f77_zmq_setsockopt(zmq_socket_push,ZMQ_LINGER,300000,4)
if (rc /= 0) then
stop 'Unable to set ZMQ_LINGER on push socket'
endif
! rc = f77_zmq_setsockopt(zmq_socket_push,ZMQ_LINGER,300000,4)
! if (rc /= 0) then
! stop 'Unable to set ZMQ_LINGER on push socket'
! endif
call omp_set_lock(zmq_lock)
rc = f77_zmq_close(zmq_socket_push)
@ -535,7 +535,7 @@ subroutine new_parallel_job(zmq_to_qp_run_socket,name_in)
rc = f77_zmq_recv(zmq_to_qp_run_socket,message,510,0)
message = trim(message(1:rc))
if (message(1:2) /= 'ok') then
print *, message
print *, trim(message(1:rc))
print *, 'Unable to start parallel job : '//name
stop 1
endif
@ -565,6 +565,7 @@ subroutine zmq_set_running(zmq_to_qp_run_socket)
rc = f77_zmq_recv(zmq_to_qp_run_socket,message,510,0)
message = trim(message(1:rc))
if (message(1:2) /= 'ok') then
print *, trim(message(1:rc))
print *, 'Unable to set qp_run to Running'
stop 1
endif
@ -718,6 +719,7 @@ subroutine add_task_to_taskserver(zmq_to_qp_run_socket,task)
rc = f77_zmq_recv(zmq_to_qp_run_socket, message, sze-1, 0)
if (message(1:rc) /= 'ok') then
print *, trim(message(1:rc))
print *, trim(task)
print *, 'Unable to add the next task'
stop -1
@ -762,12 +764,40 @@ subroutine add_task_to_taskserver_recv(zmq_to_qp_run_socket)
character*(512) :: message
rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0)
if (message(1:rc) /= 'ok') then
print *, trim(message(1:rc))
print *, 'Unable to add the next task'
stop -1
endif
end
subroutine zmq_abort(zmq_to_qp_run_socket)
use f77_zmq
implicit none
BEGIN_DOC
! Aborts a running parallel computation
END_DOC
integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket
integer :: rc, sze
character*(512) :: message
write(message,*) 'abort '
sze = len(trim(message))
rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
if (rc /= sze) then
print *, irp_here, 'f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)'
stop 'error'
endif
rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0)
if (trim(message(1:rc)) /= 'ok') then
print *, trim(message(1:rc))
print *, 'Unable to send abort message'
stop -1
endif
end
subroutine task_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_id)
use f77_zmq
implicit none
@ -790,6 +820,7 @@ subroutine task_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_id)
rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0)
if (trim(message(1:rc)) /= 'ok') then
print *, trim(message(1:rc))
print *, 'Unable to send task_done message'
stop -1
endif
@ -859,10 +890,10 @@ subroutine end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
character*(8), external :: zmq_port
integer :: rc
rc = f77_zmq_setsockopt(zmq_to_qp_run_socket,ZMQ_LINGER,1000,4)
if (rc /= 0) then
stop 'Unable to set ZMQ_LINGER on zmq_to_qp_run_socket'
endif
! rc = f77_zmq_setsockopt(zmq_to_qp_run_socket,ZMQ_LINGER,1000,4)
! if (rc /= 0) then
! stop 'Unable to set ZMQ_LINGER on zmq_to_qp_run_socket'
! endif
rc = f77_zmq_close(zmq_to_qp_run_socket)
if (rc /= 0) then
@ -901,11 +932,11 @@ subroutine zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more)
more = 1
else if (reply(16:19) == 'done') then
more = 0
rc = f77_zmq_setsockopt(zmq_socket_pull, ZMQ_RCVTIMEO, 1000, 4)
if (rc /= 0) then
print *, 'f77_zmq_setsockopt(zmq_socket_pull, ZMQ_RCVTIMEO, 3000, 4)'
stop 'error'
endif
! rc = f77_zmq_setsockopt(zmq_socket_pull, ZMQ_RCVTIMEO, 1000, 4)
! if (rc /= 0) then
! print *, 'f77_zmq_setsockopt(zmq_socket_pull, ZMQ_RCVTIMEO, 3000, 4)'
! stop 'error'
! endif
else
print *, reply
print *, 'f77_zmq_recv(zmq_to_qp_run_socket,reply,64,0)'