1 subroutine matvecComplex(rtdb, nlen, g_vec, g_rhs, xyz, ldbg) 2c ===================================================================== 3c Performs the operation A*x = b, where x and b are real. 4c 5c In : rtdb - Handle to the RTDB 6c nlen - Length of the work array 7c g_vec - Global array handle for the vector to be multiplied 8c xyz - Cartesian coordinates to DIM atoms 9c ldebug - Debug flag 10c Out : g_rhs - Right-hand side of matrix-vector multiply (solution) 11c ===================================================================== 12 implicit none 13#include "errquit.fh" 14#include "inp.fh" 15#include "rtdb.fh" 16#include "stdio.fh" 17#include "nwc_const.fh" 18#include "mafdecls.fh" 19#include "global.fh" 20#include "testutil.fh" 21#include "dimqm_constants.fh" 22#include "dimqm.fh" 23c 24c Input Variables 25 integer rtdb 26 integer nlen 27 integer g_vec 28 integer g_rhs 29 logical ldbg 30 double precision xyz(3,nDIM) 31c 32c Local Variables 33 integer ntype, mtype 34 integer m, n, ms, me, ns, ne, id 35 integer mLo, mUp, mOnNode, nodes 36 double precision time 37 double precision sPol(nDIMTypes) 38 double complex vec(nOrder) 39 integer tx(nDIM) ! Type Index 40 double complex diagFD(nOrder) 41 double complex fld(3) 42 logical stat 43 integer lo, hi 44c 45c Common variables used from dimqm.fh 46c 47c integer nDIM, nDIMTypes, nOrder, d_DIM_diag 48c 49c Determine node ID 50 id = ga_nodeid() 51 if(id.eq.0 .and. ldbg) write(LuOut,*) 52 $ "Start MatvecComplex routine" 53 time = util_timer() 54c 55c Initialize right hand side 56 call ga_zero(g_rhs) 57c 58c Determine chunk to be done on this node 59 call pphilo(id, nDIM, mLo, mUp, mOnNode) 60c 61c Fill arrays 62 stat = rtdb_get(rtdb, 'dimpar:sPol', mt_dbl, nDIMTypes, sPol) 63 stat = rtdb_get(rtdb, 'dimpar:typeindex', mt_int, nDIM, tx) 64 call ga_get(g_DIM_diag, 1, nOrder, 1, 1, diagFD, 1) 65 call ga_get(g_vec, 1, nOrder, 1, 1, vec, 1) 66 call ga_sync() 67 68 call ga_init_fence() 69c 70c Loop over atoms on this node 71 do m = mLo, mUp 72 mtype = tx(m) 73 fld = ZERO_C 74 ms = 3 * (m - 1) + 1 75 me = 3 * (m - 1) + 3 76 do n = 1, nDIM 77 ntype = tx(n) 78 ns = 3 * (n - 1) + 1 79 ne = 3 * (n - 1) + 3 80 fld(1:3) = fld(1:3) + atomPairField(m, n, 81 $ xyz(:,m), ! xyz1 82 $ xyz(:,n), ! xyz2 83 $ sPol(mtype), ! static pol1 (for screening) 84 $ sPol(ntype), ! static pol2 (for screening) 85 $ diagFD(ns:ne), ! ainv 86 $ vec(ns:ne), ! vec 87 $ ldbg) ! debug 88 end do 89 call ga_acc(g_rhs, ms, me, 1, 1, fld, 1, ONE_C) 90 end do 91 call ga_fence() 92c 93c Clean up arrays 94 time = util_timer() - time 95 if(id.eq.0.and.ldbg) write(LuOut,*) 96 $ "End MatvecComplex Routine" 97 98 contains 99 function atomPairField(m, n, mxyz, nxyz, mpol, npol, 100 $ diag, vec, debug) result (fld) 101c 102c ============================================ 103c Calculates the field at atom m due to atom n 104c ============================================ 105c 106 107 implicit none 108#include "stdio.fh" 109#include "dimqm_constants.fh" 110 integer, intent(in) :: m, n 111 double precision, intent(in) :: mxyz(3), nxyz(3) 112 double precision, intent(in) :: mpol, npol 113 double complex, intent(in) :: diag(3) 114 double complex, intent(in) :: vec(3) 115 logical, intent(in) :: debug 116c 117c Output 118 double complex fld(3) 119c 120c Local variables 121 double precision r(3) 122 double precision t2(3,3) 123 double precision dist, invdist 124 interface 125 function t2r(r, dist, invdist, a1, a2, pol1, pol2) 126 integer, intent(in) :: a1, a2 127 double precision, intent(in) :: dist, invdist 128 double precision, intent(in) :: pol1, pol2 129 double precision t2r(3,3) 130 double precision, intent(in) :: r(3) 131 end function t2r 132 end interface 133c 134c Diagonal element 135 if(m .eq. n) then 136 fld(1:3) = vec(1:3) * diag(1:3) 137 else 138c Off diagonal 139 r(1:3) = nxyz(1:3) - mxyz(1:3) 140 dist = SQRT(DOT_PRODUCT(r, r)) 141 invdist = ONE / dist 142 t2 = t2r(r, dist, invdist, m, n, mpol, npol) 143 fld(1:3) = - MATMUL(t2, vec) 144 end if 145 end function atomPairField 146 end subroutine matvecComplex 147 148