TABLE OF CONTENTS
ABINIT/sincos [ 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