File: mixbroyden.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (99 lines) | stat: -rw-r--r-- 2,169 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99

! Copyright (C) 2011 J. K. Dewhurst, S. Sharma and E. K. U. Gross.
! This file is distributed under the terms of the GNU Lesser General Public
! License. See the file COPYING for license details.

subroutine mixbroyden(iscl,n,msd,alpha,w0,nu,mu,f,df,u,a,d)
implicit none
! arguments
integer, intent(in) :: iscl,n,msd
real(8), intent(in) :: alpha,w0
real(8), intent(inout) :: nu(n),mu(n,2)
real(8), intent(inout) :: f(n,2),df(n,msd)
real(8), intent(inout) :: u(n,msd)
real(8), intent(inout) :: a(msd,msd)
real(8), intent(out) :: d
! local variables
integer jc,kp,kc
integer k,l,m,info
real(8) t1
! automatic arrays
integer ipiv(msd)
real(8) c(msd),beta(msd,msd),gamma(msd)
real(8) work(msd)
! external functions
real(8) dnrm2
external dnrm2
if (n.lt.1) then
  write(*,*)
  write(*,'("Error(mixbroyden): n < 1 : ",I8)') n
  write(*,*)
  stop
end if
if (msd.lt.2) then
  write(*,*)
  write(*,'("Error(mixbroyden): msd < 2 : ",I8)') msd
  write(*,*)
  stop
end if
! initialise mixer
if (iscl.le.0) then
  call dcopy(n,nu,1,mu(:,1),1)
  call dcopy(n,nu,1,mu(:,2),1)
  f(:,1)=0.d0
  df(:,1)=0.d0
  u(:,1)=0.d0
  a(:,:)=0.d0
  d=1.d0
  return
end if
! current subspace dimension
m=min(iscl+1,msd)
! current index modulo m
jc=mod(iscl,m)+1
! previous index modulo 2
kp=mod(iscl-1,2)+1
! current index modulo 2
kc=mod(iscl,2)+1
f(:,kc)=nu(:)-mu(:,kp)
d=sum(f(:,kc)**2)
d=sqrt(d/dble(n))
df(:,jc)=f(:,kc)-f(:,kp)
t1=dnrm2(n,df(:,jc),1)
if (t1.gt.1.d-8) t1=1.d0/t1
call dscal(n,t1,df(:,jc),1)
u(:,jc)=alpha*df(:,jc)+t1*(mu(:,kp)-mu(:,kc))
do k=1,m
  c(k)=dot_product(df(:,k),f(:,kc))
end do
do k=1,m
  a(k,jc)=dot_product(df(:,jc),df(:,k))
  a(jc,k)=a(k,jc)
end do
beta(:,:)=a(:,:)
do k=1,m
  beta(k,k)=beta(k,k)+w0**2
end do
! invert beta
call dgetrf(m,m,beta,msd,ipiv,info)
if (info.eq.0) call dgetri(m,beta,msd,ipiv,work,m,info)
if (info.ne.0) then
  write(*,*)
  write(*,'("Error(mixbroyden): could not invert matrix")')
  write(*,*)
  stop
end if
do l=1,m
  gamma(l)=0.d0
  do k=1,m
    gamma(l)=gamma(l)+c(k)*beta(k,l)
  end do
end do
nu(:)=mu(:,kp)+alpha*f(:,kc)
do l=1,m
  call daxpy(n,-gamma(l),u(:,l),1,nu,1)
end do
call dcopy(n,nu,1,mu(:,kc),1)
return
end subroutine