10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-17 02:35:26 +02:00

Merge branch 'master' of github.com:scemama/quantum_package

This commit is contained in:
Anthony Scemama 2016-10-18 23:07:17 +02:00
commit 43aa7a3de4
12 changed files with 292 additions and 188 deletions

View File

@ -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

23
ocaml/Id.mli Normal file
View File

@ -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

View File

@ -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 }

View File

@ -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

63
ocaml/Queuing_system.mli Normal file
View File

@ -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

View File

@ -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;

View File

@ -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

View File

@ -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><i|H|u_k>
! -----------------------------------------
@ -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 = <u_k | W_l> = <u_k| H |u_l>
! -------------------------------------------
! 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, &

View File

@ -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))

View File

@ -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)

View File

@ -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

View File

@ -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) &