1C> \ingroup nwdft_xc
2C> @{
3C>
4C> \file dft_utilmap.F
5C> Various routines that build or use mappings between atoms, matrix
6C> elements and what not.
7C>
8      subroutine build_maps(basis, cntoce, cntobfr, cetobfr, natoms,
9     &                      nshells)
10*
11* $Id$
12*
13      implicit none
14#include "errquit.fh"
15      integer basis, natoms, nshells
16      integer cntoce(nshells), cntobfr(2,nshells), cetobfr(2,natoms)
17c
18#include "bas.fh"
19c
20      integer ish, iat
21c
22c     Build maps (for speed).
23c
24      do ish = 1, nshells
25         if (.not. bas_cn2ce(basis, ish, cntoce(ish)))
26     &        call errquit('build_maps: bad basis', 0, BASIS_ERR)
27      end do
28c
29      do ish = 1, nshells
30         if (.not. bas_cn2bfr(basis,ish,cntobfr(1,ish),cntobfr(2,ish)))
31     &        call errquit('build_maps: bad basis', 0, BASIS_ERR)
32      end do
33c
34      do iat = 1, natoms
35         if (.not. bas_ce2bfr(basis,iat,cetobfr(1,iat),cetobfr(2,iat)))
36     &        call errquit('build_maps: bad basis', 0, BASIS_ERR)
37      end do
38c
39      return
40      end
41      subroutine util_ga_mat_reduce(n, nr, map, g_a, n_a, r, op,
42     &                              scr, mxmap,lfirst)
43      implicit none
44c
45#include "global.fh"
46#include "msgids.fh"
47      integer n                         ! Original size [input]
48      integer nr                        ! Reduced size  [input]
49      integer map(2,nr)                 ! map(1,*)=lo, map(2,*)=hi [input]
50      integer n_a                       ! number of GAs to reduce
51      integer g_a(n_a)                  ! Original GA handle(s) [input]
52      double precision r(nr,nr,n_a)     ! Reduced matrix [output]
53      integer mxmap                     ! max size of map vector
54      double precision scr(mxmap,*) ! scratch
55      character*(*) op                  ! Reduction operation
56      logical lfirst
57c
58c     R(i,j) <= R(i,j) op A(map(1,i):map(2,i),map(1,j):map(2,j))
59c
60c     where op is one of 'abssum', 'absmax', 'rms' (extend as necessary)
61c
62      integer ir, jr, k, ielems, jelems
63      double precision sum
64      integer me,nproc,itask
65      double precision util_matops
66      external util_matops
67      integer ityp, ndim, dims(3)
68      integer lo(3),hi(3),ld(2),ld12,offs
69      integer distr_lo(3),distr_hi(3)
70      integer get_lo(3),get_hi(3)
71      logical l3d,getdone
72c
73      call nga_inquire(g_a, ityp, ndim, dims)
74      l3d=.false.
75      if(ndim.gt.2) then
76#ifdef DEBUG
77         write(6,*) ' matreduce: ndim gt 2!, =',ndim
78         write(6,*) ' na ',n_a, ' dims1 ',dims(1),dims(1)/n_a
79         write(6,*) '  dims ',dims
80#endif
81         l3d=.true.
82         offs=dims(1)/n_a
83      endif
84      me=ga_nodeid()
85      nproc=ga_nnodes()
86      itask=0
87      if(lfirst)
88     A     call dfill(n_a*nr*nr, 0.0d0, r, 1)
89      call nga_distribution(g_a,
90     .           ga_nodeid(), distr_lo, distr_hi)
91      if (distr_lo(1).gt.0 .and. distr_lo(1).le.distr_hi(1)) then
92
93      do k = 1, n_a
94c
95c         write(6,*) ' util_mat_reduce: input matrix '
96c         call ga_print(g_a(k))
97c
98         do jr = 1, nr
99cedo            do ir = 1, nr
100            do ir = 1, jr
101c               itask=itask+1
102c               if(mod(itask,nproc).eq.me) then
103                  ielems = map(2,ir) - map(1,ir) + 1
104                  jelems = map(2,jr) - map(1,jr) + 1
105                  if (ielems.gt.0 .and. jelems.gt.0) then
106                     getdone=.false.
107                     if(l3d) then
108                        ielems=ielems*offs
109                        lo(1)=1+(k-1)*offs
110                        hi(1)=lo(1)+offs-1
111#if 0
112                        lo(2)=map(1,ir)
113                        hi(2)=map(2,ir)
114                        lo(3)=map(1,jr)
115                        hi(3)=map(2,jr)
116#else
117                        lo(2)=max(map(1,ir),distr_lo(1))
118                        hi(2)=min(map(2,ir),distr_hi(1))
119                        lo(3)=max(map(1,jr),distr_lo(2))
120                        hi(3)=min(map(2,jr),distr_hi(2))
121#endif
122                        ld(1)=offs
123                        ld(2)=mxmap
124                        getdone=(get_lo(2).le.get_hi(2)).and.
125     A                       (get_lo(3).le.get_hi(3))
126                        if(getdone)
127     G                   call nga_get(g_a, lo,hi,
128     .                       scr, ld)
129                        ld12=ld(1)*ld(2)
130                     else
131c                        call ga_get(g_a(k), map(1,ir), map(2,ir),
132c     .                       map(1,jr), map(2,jr), scr, mxmap)
133                        get_lo(1)=max(map(1,ir),distr_lo(1))
134                        get_hi(1)=min(map(2,ir),distr_hi(1))
135                        get_lo(2)=max(map(1,jr),distr_lo(2))
136                        get_hi(2)=min(map(2,jr),distr_hi(2))
137                        getdone=(get_lo(1).le.get_hi(1)).and.
138     A                       (get_lo(2).le.get_hi(2))
139                        if(getdone)
140     G                  call ga_get(g_a(k), get_lo(1), get_hi(1),
141     .                       get_lo(2), get_hi(2), scr, mxmap)
142                        ld12=mxmap
143                     endif
144                     if(getdone) then
145                     sum=util_matops(op,ielems,jelems,scr,ld12)
146                     r(ir,jr,k) = max(r(ir,jr,k),sum)
147c
148c     copy upper triangle of r(ij) to upper triangle
149c
150                     if(ir.ne.jr) r(jr,ir,k)=r(ir,jr,k)
151                     endif
152                  endif
153cc               endif
154            end do
155         end do
156c         write(6,*) ' util_mat_reduce: reduced matrix '
157c         call output(r(1,1,k), 1, nr, 1, nr, nr, nr, 1)
158      enddo
159      endif
160c
161c     global sum
162c
163c      write(6,*) ' length ',nr*nr*n_a
164      call ga_dgop(msg_gop_rdens,r,nr*nr*n_a, 'absmax')
165c
166      end
167      double precision function util_matops(op,ielems,jelems,scr,ld)
168      implicit none
169#include "errquit.fh"
170      character*(*) op
171      integer ielems,jelems
172      integer ld
173      double precision scr(ld,*)
174c
175      integer i,j
176      double precision sum
177      sum=0d0
178      if (op .eq. 'abssum') then
179         do j = 1, jelems
180            do i = 1, ielems
181               sum = sum + abs(scr(i,j))
182            end do
183         end do
184      else if (op .eq. 'absmax') then
185         do j = 1, jelems
186            do i = 1, ielems
187               sum = max(sum, abs(scr(i,j)))
188            end do
189         end do
190      else if (op .eq. 'rms') then
191         do j = 1, jelems
192            do i = 1, ielems
193               sum = sum + scr(i,j)*scr(i,j)
194            end do
195         enddo
196         sum = sqrt(sum)
197      else
198         call errquit('util_ga_mat_red: unknown op',0, UNKNOWN_ERR)
199      end if
200      util_matops=sum
201      return
202      end
203      subroutine util_irreg_mat_reduce(n_row, n_col, nr_row, nr_col,
204     &                                 row_map, col_map, a, r, op)
205      implicit none
206#include "errquit.fh"
207c
208      integer n_row                      ! Original row size [input]
209      integer n_col                      ! Original col size [input]
210      integer nr_row                     ! Reduced row size  [input]
211      integer nr_col                     ! Reduced col size  [input]
212      integer row_map(2,nr_row)          ! map(1,*)=lo, map(2,*)=hi [input]
213      integer col_map(2,nr_col)          ! map(1,*)=lo, map(2,*)=hi [input]
214      double precision a(n_row, n_col)   ! Original matrix [input]
215      double precision r(nr_row, nr_col) ! Reduced matrix [output]
216      character*(*) op                   ! Reduction operation
217c
218c     R(i,j) <= R(i,j) op A(map(1,i):map(2,i),map(1,j):map(2,j))
219c
220c     where op is one of 'abssum', 'absmax', 'rms' (extend as necessary)
221c
222      integer ir, jr, i, j
223      double precision sum
224c
225      do jr = 1, nr_col
226         do ir = 1, nr_row
227            sum = 0.0d0
228            if (op .eq. 'abssum') then
229               do j = col_map(1,jr), col_map(2,jr)
230                  do i = row_map(1,ir), row_map(2,ir)
231                     if (i.ne.0.and.j.ne.0)sum = sum + abs(a(i,j))
232                  end do
233               end do
234            else if (op .eq. 'absmax') then
235               do j = col_map(1,jr), col_map(2,jr)
236                  do i = row_map(1,ir), row_map(2,ir)
237                     if (i.ne.0.and.j.ne.0)sum = max(sum, abs(a(i,j)))
238                  end do
239               end do
240            else if (op .eq. 'rms') then
241               do j = col_map(1,jr), col_map(2,jr)
242                  do i = row_map(1,ir), row_map(2,ir)
243                     if (i.ne.0.and.j.ne.0)sum = sum + a(i,j)*a(i,j)
244                  end do
245               enddo
246               sum = sqrt(sum)
247            else
248               call errquit('util_irreg_mat_reduce: unknown op', 0,
249     &       UNKNOWN_ERR)
250            end if
251            r(ir,jr) = max(r(ir,jr),sum)
252         end do
253      end do
254c
255c      write(6,*) ' util_irreg_mat_reduce: input matrix '
256c      call output(a, 1, n_row, 1, n_col, n_row, n_col, 1)
257c
258c      write(6,*) ' Row map begin: ',(row_map(1,ir),ir = 1,nr_row)
259c      write(6,*) ' Row map end: ',(row_map(2,ir),ir = 1,nr_row)
260c      write(6,*) ' Col map begin: ',(col_map(1,ir),ir = 1,nr_col)
261c      write(6,*) ' Col map end: ',(col_map(2,ir),ir = 1,nr_col)
262c
263c      write(6,*) ' util_irreg_mat_reduce: reduced matrix '
264c      call output(r, 1, nr_row, 1, nr_col, nr_row, nr_col, 1)
265c
266      end
267C>
268C> \brief Compute the root-mean-square value of the basis functions
269C> or basis function gradients for each atom
270C>
271C> Compute the root-mean-square value of basis functions or basis
272C> function gradients for each atom. This data is usefull for screening
273C> purposes and is a measure for the size of the basis functions on
274C> a particular atom at the grid points.
275C>
276C> The exact quantity computed depends on the inputs for `chi` and
277C> `nqsmall`. If `chi` represents basis functions then `rchi` will
278C> be computed to be
279C> \f{eqnarray*}{
280C>   \mathrm{rchi}(i_{at})
281C>   &=& \sqrt{\sum_{j\in\{\chi(i_{at})\}}\sum_q\chi_j^2(r_q)}
282C> \f}
283C> If `chi` represents the gradients of the basis functions then `rchi`
284C> will be computed to be (the number of components in the gradient is
285C> in practice encoded in `nqsmall`)
286C> \f{eqnarray*}{
287C>   \mathrm{rchi}(i_{at})
288C>   &=& \sqrt{\sum_{j\in\{\chi(i_{at})\}}\sum_q\sum_{c\in\{x,y,z\}}
289C>             \nabla_c\chi_j^2(r_q)}
290C> \f}
291C> where \f$\{\chi(i_{at})\}\f$ designates the set of all basis
292C> functions on atom \f$i_{at}\f$.
293C>
294      subroutine util_rmsatbf(nqsmall, natoms,iniz,ifin,
295     .     chi,rchi)
296      implicit none
297      integer nqsmall !< [Input] The number of data points
298                      !< - the number of grid points, `nq`, for basis
299                      !< functions
300                      !< - `3*nq` for gradients of basis functions
301      integer natoms  !< [Input] The number of atoms
302      integer iniz(*) !< [Input] The initial basis function on each
303                      !< atom
304      integer ifin(*) !< [Input] The final basis function on each atom
305      double precision chi(nqsmall,*) !< [Input] The basis function
306                                      !< (gradient) amplitudes
307      double precision rchi(*) !< [Output] The root-mean-square of
308                               !< of all amplitudes on each atom
309c
310      integer q,iat,jf
311      double precision sum,mxsum
312c
313      mxsum=0d0
314      do iat=1,natoms
315        if(iniz(iat).eq.0) then
316          rchi(iat)=0d0
317        else
318           sum=0d0
319           do jf=iniz(iat),ifin(iat)
320              do q=1,nqsmall
321                 sum=sum+chi(q,jf)*chi(q,jf)
322              enddo
323           enddo
324           rchi(iat)=max(mxsum,sqrt(sum))
325        endif
326      enddo
327      return
328      end
329C>
330C> @}
331