src/mixadapt.f90

   1
   2 ! Copyright (C) 2002-2008 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
   3 ! This file is distributed under the terms of the GNU Lesser General Public
   4 ! License. See the file COPYING for license details.
   5
   6 !BOP
   7 ! !ROUTINE: mixadapt
   8 ! !INTERFACE:
   9 subroutine mixadapt(iscl,beta0,betainc,betadec,n,nu,mu,beta,f,d)
  10 ! !INPUT/OUTPUT PARAMETERS:
  11 !   iscl    : self-consistent loop number (in,integer)
  12 !   beta0   : initial value for mixing parameter (in,real)
  13 !   betainc : mixing parameter increase (in,real)
  14 !   betadec : mixing parameter decrease (in,real)
  15 !   n       : vector length (in,integer)
  16 !   nu      : current output vector as well as the next input vector of the
  17 !             self-consistent loop (inout,real(n))
  18 !   mu      : used internally (inout,real(n))
  19 !   beta    : used internally (inout,real(n))
  20 !   f       : used internally (inout,real(n))
  21 !   d       : RMS difference between old and new output vectors (out,real)
  22 ! !DESCRIPTION:
  23 !   Given the input vector $\mu^i$ and output vector $\nu^i$ of the $i$th
  24 !   self-consistent loop, this routine generates the next input vector to the
  25 !   loop using an adaptive mixing scheme. The $j$th component of the output
  26 !   vector is mixed with a fraction of the same component of the input vector:
  27 !   $$ \mu^{i+1}_j=\beta^i_j\nu^i_j+(1-\beta^i_j)\mu^i_j, $$
  28 !   where $\beta^i_j$ is set to $\beta_0$ at initialisation and increased by
  29 !   scaling with $\beta_{\rm inc}$ ($>1$) if $f^i_j\equiv\nu^i_j-\mu^i_j$ does
  30 !   not change sign between loops. If $f^i_j$ does change sign, then $\beta^i_j$
  31 !   is scaled by $\beta_{\rm dec}$ ($<1$). Note that the array {\tt nu} serves
  32 !   for both input and output, and the arrays {\tt mu}, {\tt beta} and {\tt f}
  33 !   are used internally and should not be changed between calls. The routine is
  34 !   initialised at the first iteration and is thread-safe so long as each thread
  35 !   has its own independent work array. Complex arrays may be passed as real
  36 !   arrays with $n$ doubled.
  37 !
  38 ! !REVISION HISTORY:
  39 !   Created March 2003 (JKD)
  40 !   Modified, September 2008 (JKD)
  41 !EOP
  42 !BOC
  43 implicit none
  44 ! arguments
  45 integer, intent(in) :: iscl
  46 real(8), intent(in) :: beta0
  47 real(8), intent(in) :: betainc
  48 real(8), intent(in) :: betadec
  49 integer, intent(in) :: n
  50 real(8), intent(inout) :: nu(n)
  51 real(8), intent(inout) :: mu(n)
  52 real(8), intent(inout) :: beta(n)
  53 real(8), intent(inout) :: f(n)
  54 real(8), intent(out) :: d
  55 ! local variables
  56 integer i
  57 real(8) t1
  58 if (iscl.le.1) then
  59   mu(:)=nu(:)
  60   f(:)=0.d0
  61   beta(:)=beta0
  62   d=1.d0
  63   return
  64 end if
  65 do i=1,n
  66   t1=nu(i)-mu(i)
  67   if (t1*f(i).gt.0.d0) then
  68     beta(i)=beta(i)*betainc
  69     if (beta(i).gt.1.d0) beta(i)=1.d0
  70   else
  71     beta(i)=beta(i)*betadec
  72   end if
  73   f(i)=t1
  74 end do
  75 nu(:)=beta(:)*nu(:)+(1.d0-beta(:))*mu(:)
  76 d=0.d0
  77 do i=1,n
  78   d=d+f(i)**2
  79 end do
  80 d=sqrt(d/dble(n))
  81 mu(:)=nu(:)
  82 return
  83 end subroutine
  84 !EOC
  85