1 subroutine matvecReal(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 ================================================================== 8 implicit none 9#include "errquit.fh" 10#include "inp.fh" 11#include "rtdb.fh" 12#include "stdio.fh" 13#include "nwc_const.fh" 14#include "mafdecls.fh" 15#include "global.fh" 16#include "testutil.fh" 17#include "dimqm_constants.fh" 18#include "dimqm.fh" 19c 20c Input Variables 21 integer rtdb 22 integer nlen 23 integer g_vec 24 integer g_rhs 25 logical ldbg 26 double precision xyz(3,nDIM) 27c 28c Variables from dimqm common blocks 29c integer nDIM, nDIMTypes, nOrder, g_DIM_diag 30c 31c Local Variables 32 integer ntype, mtype 33 integer m, n, ms, me, ns, ne, id 34 integer mLo, mUp, mOnNode, nodes 35 double precision time 36 double precision sPol(nDIMTypes) 37 double precision vec(nOrder) 38 integer tx(nDIM) ! Type Index 39 double precision diagS(nOrder) 40 double precision fld(3) 41 logical stat 42 integer lo, hi 43 call ga_sync() 44c 45c Determine node ID 46 id = ga_nodeid() 47 if(ldbg .and. id.eq.0) write(LuOut,*) 48 $ "Start MatvecReal routine" 49 time = util_timer() 50c 51c Zero right-hand side 52 call ga_zero(g_rhs) 53c 54c Determine chunk to be done on this node 55 call pphilo(id, nDIM, mLo, mUp, mOnNode) 56c 57c Pull values from RTDB and GAs 58 stat = rtdb_get(rtdb, 'dimpar:sPol', mt_dbl, nDIMTypes, sPol) 59 stat = rtdb_get(rtdb, 'dimpar:typeindex', mt_int, nDIM, tx) 60 call ga_get(g_DIM_diag, 1, nOrder, 1, 1, diagS, 1) 61 call ga_get(g_vec, 1, nOrder, 1, 1, vec, 1) 62 call ga_init_fence() 63c 64c Loop over atoms on this node 65 do m = mLo, mUp 66 mtype = tx(m) 67 fld = ZERO 68 ms = 3 * (m - 1) + 1 69 me = 3 * (m - 1) + 3 70 do n = 1, nDIM 71 ntype = tx(n) 72 ns = 3 * (n - 1) + 1 73 ne = 3 * (n - 1) + 3 74 fld(1:3) = fld(1:3) + atomPairField(m, n, 75 $ xyz(:,m), ! Coordinates to atom m 76 $ xyz(:,n), ! Coordinates to atom n 77 $ sPol(mtype), ! Static polarizability of atom m 78 $ sPol(ntype), ! Static polarizability of atom n 79 $ diagS(ns:ne), ! Diagonal elements 80 $ vec(ns:ne), ! Vector x in A*x = b 81 $ ldbg) ! debug 82 end do 83 call ga_acc(g_rhs, ms, me, 1, 1, fld, 1, ONE) 84 end do 85 call ga_fence() 86 time = util_timer() - time 87 if(ldbg .and. id.eq.0) write(LuOut,*) 88 $ "End MatvecReal Routine" 89 90 contains 91 function atomPairField(m, n, mxyz, nxyz, mpol, npol, 92 $ diag, vec, debug) result (fld) 93 implicit none 94#include "stdio.fh" 95#include "dimqm_constants.fh" 96 integer, intent(in) :: m, n 97 double precision, intent(in) :: mxyz(3), nxyz(3) ! Coordinates of the atoms 98 double precision, intent(in) :: mpol, npol ! Static polarizability of the atoms 99 double precision, intent(in) :: diag(3) ! Diagonal elements 100 double precision, intent(in) :: vec(3) ! Vector being multiplied 101 logical, intent(in) :: debug 102c 103c Output 104 double precision fld(3) 105c 106c Local variables 107 double precision r(3) 108 double precision t2(3,3) 109 double precision dist, invdist 110 interface 111 function t2r(r, dist, invdist, a1, a2, pol1, pol2) 112 integer, intent(in) :: a1, a2 113 double precision, intent(in) :: dist, invdist 114 double precision, intent(in) :: pol1, pol2 115 double precision t2r(3,3) 116 double precision, intent(in) :: r(3) 117 end function t2r 118 end interface 119c 120c Diagonal element 121 if(m .eq. n) then 122 fld(1:3) = vec(1:3) * diag(1:3) 123 else 124C Off diagonal 125 r(1:3) = nxyz(1:3) - mxyz(1:3) 126 dist = SQRT(DOT_PRODUCT(r, r)) 127 invdist = ONE / dist 128 t2 = t2r(r, dist, invdist, m, n, mpol, npol) 129 fld(1:3) = - MATMUL(t2, vec) 130 end if 131c write(LuOut,*) "Fld:", fld 132 end function atomPairField 133 end subroutine matvecReal 134 135