1
2! Copyright (C) 2011 J. K. Dewhurst, S. Sharma and E. K. U. Gross.
3! This file is distributed under the terms of the GNU Lesser General Public
4! License. See the file COPYING for license details.
5
6subroutine mixbroyden(iscl,n,msd,alpha,w0,nu,mu,f,df,u,a,d)
7use modomp
8implicit none
9! arguments
10integer, intent(in) :: iscl,n,msd
11real(8), intent(in) :: alpha,w0
12real(8), intent(inout) :: nu(n),mu(n,2)
13real(8), intent(inout) :: f(n,2),df(n,msd)
14real(8), intent(inout) :: u(n,msd)
15real(8), intent(inout) :: a(msd,msd)
16real(8), intent(out) :: d
17! local variables
18integer jc,kp,kc
19integer k,l,m
20integer info,nthd
21real(8) t1
22! automatic arrays
23integer ipiv(msd)
24real(8) c(msd),beta(msd,msd),gamma(msd),work(msd)
25! external functions
26real(8), external :: dnrm2
27if (n.lt.1) then
28  write(*,*)
29  write(*,'("Error(mixbroyden): n < 1 : ",I8)') n
30  write(*,*)
31  stop
32end if
33if (msd.lt.2) then
34  write(*,*)
35  write(*,'("Error(mixbroyden): msd < 2 : ",I8)') msd
36  write(*,*)
37  stop
38end if
39! initialise mixer
40if (iscl.le.0) then
41  call dcopy(n,nu,1,mu,1)
42  call dcopy(n,nu,1,mu(:,2),1)
43  f(:,1)=0.d0
44  df(:,1)=0.d0
45  u(:,1)=0.d0
46  a(:,:)=0.d0
47  d=1.d0
48  return
49end if
50! current subspace dimension
51m=min(iscl+1,msd)
52! current index modulo m
53jc=mod(iscl,m)+1
54! previous index modulo 2
55kp=mod(iscl-1,2)+1
56! current index modulo 2
57kc=mod(iscl,2)+1
58f(:,kc)=nu(:)-mu(:,kp)
59d=sum(f(:,kc)**2)
60d=sqrt(d/dble(n))
61df(:,jc)=f(:,kc)-f(:,kp)
62t1=dnrm2(n,df(:,jc),1)
63if (t1.gt.1.d-8) t1=1.d0/t1
64call dscal(n,t1,df(:,jc),1)
65u(:,jc)=alpha*df(:,jc)+t1*(mu(:,kp)-mu(:,kc))
66call holdthd(m,nthd)
67!$OMP PARALLEL DEFAULT(SHARED) &
68!$OMP NUM_THREADS(nthd)
69!$OMP DO
70do k=1,m
71  c(k)=dot_product(df(:,k),f(:,kc))
72end do
73!$OMP END DO NOWAIT
74!$OMP DO
75do k=1,m
76  a(k,jc)=dot_product(df(:,jc),df(:,k))
77  a(jc,k)=a(k,jc)
78end do
79!$OMP END DO
80!$OMP END PARALLEL
81call freethd(nthd)
82beta(:,:)=a(:,:)
83do k=1,m
84  beta(k,k)=beta(k,k)+w0**2
85end do
86! invert beta
87call dgetrf(m,m,beta,msd,ipiv,info)
88if (info.eq.0) call dgetri(m,beta,msd,ipiv,work,m,info)
89if (info.ne.0) then
90  write(*,*)
91  write(*,'("Error(mixbroyden): could not invert matrix")')
92  write(*,*)
93  stop
94end if
95do l=1,m
96  gamma(l)=0.d0
97  do k=1,m
98    gamma(l)=gamma(l)+c(k)*beta(k,l)
99  end do
100end do
101nu(:)=mu(:,kp)+alpha*f(:,kc)
102do l=1,m
103  call daxpy(n,-gamma(l),u(:,l),1,nu,1)
104end do
105call dcopy(n,nu,1,mu(:,kc),1)
106end subroutine
107
108