TABLE OF CONTENTS


ABINIT/sincos [ Functions ]

[ Top ] [ Functions ]

NAME

 sincos

FUNCTION

 Update the sine and cosine values, needed inside the
 pseudopotential routines psp1lo and psp1nl.

COPYRIGHT

 Copyright (C) 1998-2017 ABINIT group (XG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  iq  = number of current wavevector q
  irmax = number of values  of r on the radial grid to be computed
  mmax = dimension of pspwk and rad
  pspwk(:,1,1) and pspwk(:,2,1) : sine and cosine of 2$\pi$ dq * rad
  pspwk(:,1,2) and pspwk(:,2,2) : sine and cosine of 2$\pi$ previous q * rad
  rad(mmax) radial grid
  tpiq = 2 $\pi$ * current wavevector q

OUTPUT

  pspwk(*,1,2) and pspwk(*,2,2) : sine and cosine of 2$\pi$ current q * rad

NOTES

 The speed was a special concern, so iterative computation
 based on addition formula is possible. Interestingly,
 this algorithm places strong constraints on accuracy,
 so this routine is machine-dependent.

PARENTS

      psp1lo,psp1nl

CHILDREN

SOURCE

 42 #if defined HAVE_CONFIG_H
 43 #include "config.h"
 44 #endif
 45 
 46 #include "abi_common.h"
 47 
 48 
 49 subroutine sincos(iq,irmax,mmax,pspwk,rad,tpiq)
 50 
 51  use defs_basis
 52  use m_profiling_abi
 53 
 54 !This section has been created automatically by the script Abilint (TD).
 55 !Do not modify the following lines by hand.
 56 #undef ABI_FUNC
 57 #define ABI_FUNC 'sincos'
 58 !End of the abilint section
 59 
 60  implicit none
 61 
 62 !Arguments ------------------------------------
 63 !scalars
 64  integer,intent(in) :: iq,irmax,mmax
 65  real(dp),intent(in) :: tpiq
 66 !arrays
 67  real(dp),intent(in) :: rad(mmax)
 68  real(dp),intent(inout) :: pspwk(mmax,2,2)
 69 
 70 !Local variables-------------------------------
 71 !scalars
 72  integer :: ir,nstep
 73  real(dp) :: prevcos,prevsin
 74  logical :: testmipspro
 75 #if defined HAVE_LINALG_MLIB
 76  real(dp) :: halfpi
 77 #endif
 78 
 79 
 80 ! *************************************************************************
 81 
 82 #if defined HAVE_LINALG_MLIB
 83  halfpi=asin(1.0d0)
 84 #endif
 85 
 86  if(iq==2)then
 87 
 88 !  Here set up the sin and cos at iq=2
 89    do ir=2,irmax
 90      pspwk(ir,1,2)=pspwk(ir,1,1)
 91      pspwk(ir,2,2)=pspwk(ir,2,1)
 92    end do
 93 
 94  else
 95 !  
 96 !  The sensitivity of the algorithm to changes of nstep
 97 !  has been tested : for all the machines except SGI - R10000 ,
 98 !  either using only the hard way, or
 99 !  using up to nstep=40 causes changes at the level
100 !  of 1.0d-16 in the total energy. Larger values of
101 !  nstep might be possible, but the associated residual
102 !  is already very small ! The accelerated computation of
103 !  sine and cosine is essential for a good speed on IBM, but,
104 !  fortunately, on the SGI - R10000 the normal computation is fast enough.
105 
106    testmipspro=.false.
107 #ifdef FC_MIPSPRO
108    testmipspro=.true.
109 #endif
110    nstep=40
111    if(iq-(iq/nstep)*nstep == 0 .or. testmipspro)then
112 
113 !    Every nstep steps, uses the hard way
114      do ir=2,irmax
115 #if defined HAVE_LINALG_MLIB
116 !      There is a bug in the hp library !! Sine is slightly inaccurate !
117        pspwk(ir,1,2)=cos(tpiq*rad(ir)-halfpi)
118 #else
119        pspwk(ir,1,2)=sin(tpiq*rad(ir))
120 #endif
121        pspwk(ir,2,2)=cos(tpiq*rad(ir))
122      end do
123 
124    else
125 
126 !    Here the fastest way, iteratively
127      do ir=2,irmax
128        prevsin=pspwk(ir,1,2)
129        prevcos=pspwk(ir,2,2)
130        pspwk(ir,1,2)=prevsin*pspwk(ir,2,1)+prevcos*pspwk(ir,1,1)
131        pspwk(ir,2,2)=prevcos*pspwk(ir,2,1)-prevsin*pspwk(ir,1,1)
132      end do
133 
134    end if 
135 
136  end if ! iq==2
137 
138 end subroutine sincos