← Back to branch summary

~ubuntu-branches/ubuntu/saucy/nwchem/saucy

« back to all changes in this revision

Viewing changes to src/nwpw/nwpwlib/utilities/nwpw_interp.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2012-02-09 20:02:41 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20120209200241-jgk03qfsphal4ug2
Tags: 6.1-1
http://bugs.debian.org/632611
https://launchpad.net/bugs/791308
* New upstream release.

[ Michael Banck ]
* debian/patches/02_makefile_flags.patch: Updated.
* debian/patches/02_makefile_flags.patch: Use internal blas and lapack code.
* debian/patches/02_makefile_flags.patch: Define GCC4 for LINUX and LINUX64
  (Closes: #632611 and LP: #791308).
* debian/control (Build-Depends): Added openssh-client.
* debian/rules (USE_SCALAPACK, SCALAPACK): Removed variables (Closes:
  #654658).
* debian/rules (LIBDIR, USE_MPIF4, ARMCI_NETWORK): New variables.
* debian/TODO: New file.
* debian/control (Build-Depends): Removed libblas-dev, liblapack-dev and
  libscalapack-mpi-dev.
* debian/patches/04_show_testsuite_diff_output.patch: New patch, shows the
  diff output for failed tests.
* debian/patches/series: Adjusted.
* debian/testsuite: Optionally run all tests if "all" is passed as option.
* debian/rules: Run debian/testsuite with "all" if DEB_BUILD_OPTIONS
  contains "checkall".

[ Daniel Leidert ]
* debian/control: Used wrap-and-sort. Added Vcs-Svn and Vcs-Browser fields.
  (Priority): Moved to extra according to policy section 2.5.
  (Standards-Version): Bumped to 3.9.2.
  (Description): Fixed a typo.
* debian/watch: Added.
* debian/patches/03_hurd-i386_define_path_max.patch: Added.
  - Define MAX_PATH if not defines to fix FTBFS on hurd.
* debian/patches/series: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/nwpw/nwpwlib/utilities/nwpw_interp.F
 
1
*
 
2
*     $Id: nwpw_scratch.F 19707 2010-10-29 17:59:36Z d3y133 $
 
3
*
 
4
 
 
5
*     *********************************
 
6
*     *                               *
 
7
*     *      nwpw_interp_init         *
 
8
*     *                               *
 
9
*     *********************************
 
10
      subroutine nwpw_interp_init(ndim0,norder0)
 
11
      implicit none
 
12
      integer ndim0,norder0
 
13
 
 
14
#include "mafdecls.fh"
 
15
#include "errquit.fh"
 
16
 
 
17
*     **** nwpw_interp_common block ****
 
18
      integer ndim_max,norder_max
 
19
      integer root(2),coeff(2),pbasis(2)
 
20
      common /nwpw_interp_common/ root,coeff,pbasis,ndim_max,norder_max
 
21
 
 
22
*     **** local variable ****
 
23
      logical value
 
24
      integer i,k
 
25
      real*8  tmp
 
26
 
 
27
      ndim_max   = ndim0
 
28
      norder_max = norder0
 
29
      value = MA_alloc_get(mt_dbl,norder_max,'root',root(2),root(1))
 
30
      value = value.and.
 
31
     >        MA_alloc_get(mt_dbl,norder_max,'coef',coeff(2),coeff(1))
 
32
      value = value.and.
 
33
     >        MA_alloc_get(mt_dbl,norder_max*ndim_max*2,
 
34
     >                     'pbasis',pbasis(2),pbasis(1))
 
35
      if (.not.value) 
 
36
     >   call errquit('nwpw_interp_init:out of heap',0,MA_ERR)
 
37
 
 
38
      do i=0,norder_max-1
 
39
         dbl_mb(root(1)+i) = -1.0d0 + i*2.0d0/dble(norder_max-1)
 
40
      end do
 
41
 
 
42
      do i=0,norder_max-1
 
43
         tmp = 1.0d0
 
44
         do k=i+1,(norder_max+i-1)
 
45
            tmp = tmp
 
46
     >           *(dbl_mb(root(1)+i)-dbl_mb(root(1)+mod(k,norder_max)))
 
47
         end do
 
48
         dbl_mb(coeff(1)+i) = 1.0d0/tmp
 
49
      end do
 
50
 
 
51
      return
 
52
      end
 
53
 
 
54
*     *********************************
 
55
*     *                               *
 
56
*     *      nwpw_interp_end          *
 
57
*     *                               *
 
58
*     *********************************
 
59
      subroutine nwpw_interp_end()
 
60
      implicit none
 
61
 
 
62
#include "mafdecls.fh"
 
63
#include "errquit.fh"
 
64
 
 
65
*     **** nwpw_interp_common block ****
 
66
      integer ndim_max,norder_max
 
67
      integer root(2),coeff(2),pbasis(2)
 
68
      common /nwpw_interp_common/ root,coeff,pbasis,ndim_max,norder_max
 
69
 
 
70
*     **** local variables ****
 
71
      logical value
 
72
 
 
73
      value =           MA_free_heap(pbasis(2))
 
74
      value = value.and.MA_free_heap(coeff(2))
 
75
      value = value.and.MA_free_heap(root(2))
 
76
      if (.not.value)
 
77
     >   call errquit('nwpw_interp_end:freeing heap',0,MA_ERR)
 
78
      return
 
79
      end
 
80
 
 
81
 
 
82
*     *********************************
 
83
*     *                               *
 
84
*     *      nwpw_interp_basis        *
 
85
*     *                               *
 
86
*     *********************************
 
87
      real*8 function nwpw_interp_basis(i,x)
 
88
      implicit none
 
89
      integer i
 
90
      real*8 x
 
91
 
 
92
#include "mafdecls.fh"
 
93
#include "errquit.fh"
 
94
 
 
95
*     **** nwpw_interp_common block ****
 
96
      integer ndim_max,norder_max
 
97
      integer root(2),coeff(2),pbasis(2)
 
98
      common /nwpw_interp_common/ root,coeff,pbasis,ndim_max,norder_max
 
99
 
 
100
*     **** local variables ****
 
101
      integer k
 
102
      real*8  f
 
103
 
 
104
      f = 1.0d0
 
105
      do k=i+1,(i+norder_max-1)
 
106
         f = f*(x-dbl_mb(root(1)+mod(k,norder_max)))
 
107
      end do
 
108
      f = f*dbl_mb(coeff(1)+i)
 
109
 
 
110
      nwpw_interp_basis = f
 
111
      return
 
112
      end
 
113
 
 
114
 
 
115
*     *********************************
 
116
*     *                               *
 
117
*     *      nwpw_interp_dbasis       *
 
118
*     *                               *
 
119
*     *********************************
 
120
      real*8 function nwpw_interp_dbasis(i,x)
 
121
      implicit none
 
122
      integer i
 
123
      real*8 x
 
124
 
 
125
#include "mafdecls.fh"
 
126
#include "errquit.fh"
 
127
 
 
128
*     **** nwpw_interp_common block ****
 
129
      integer ndim_max,norder_max
 
130
      integer root(2),coeff(2),pbasis(2)
 
131
      common /nwpw_interp_common/ root,coeff,pbasis,ndim_max,norder_max
 
132
 
 
133
*     **** local variables ****
 
134
      integer k,kk
 
135
      real*8  tmp,df
 
136
 
 
137
      df = 0.0d0
 
138
      do k=i+1,(i+norder_max-1)
 
139
         tmp = 1.0d0
 
140
         do kk=i+1,(i+norder_max-1)
 
141
            if (kk.ne.k) tmp=tmp*(x-dbl_mb(root(1)+mod(kk,norder_max)))
 
142
         end do
 
143
         df = df + tmp
 
144
      end do
 
145
      df = df*dbl_mb(coeff(1)+i)
 
146
 
 
147
      nwpw_interp_dbasis = df
 
148
      return
 
149
      end
 
150
 
 
151
 
 
152
*     *********************************
 
153
*     *                               *
 
154
*     *          nwpw_interp          *
 
155
*     *                               *
 
156
*     *********************************
 
157
      real*8 function nwpw_interp(ndim,N,a,b,mesh,x)
 
158
      implicit none
 
159
      integer ndim,N(*)
 
160
      real*8  a(*),b(*),mesh(*),x(*)
 
161
 
 
162
#include "mafdecls.fh"
 
163
#include "errquit.fh"
 
164
 
 
165
*     **** nwpw_interp_common block ****
 
166
      integer ndim_max,norder_max
 
167
      integer root(2),coeff(2),pbasis(2)
 
168
      common /nwpw_interp_common/ root,coeff,pbasis,ndim_max,norder_max
 
169
 
 
170
*     **** local variables ****
 
171
      logical failed
 
172
      integer d,n1,nn1,i,ii,im,ip,ishift,jm(10),j(10)
 
173
      real*8  xm,xp,xtilde,f,fi
 
174
 
 
175
*     **** external functions ****
 
176
      real*8   nwpw_interp_basis
 
177
      external nwpw_interp_basis
 
178
 
 
179
      failed = .false.
 
180
      do d=1,ndim
 
181
         failed = failed.or.(x(d).lt.a(d))
 
182
         failed = failed.or.(x(d).gt.b(d))
 
183
      end do
 
184
      if (failed) then
 
185
         nwpw_interp = 0.0d0
 
186
         return
 
187
      end if
 
188
 
 
189
 
 
190
      ishift = norder_max/2 - 1 + mod(norder_max,2)
 
191
 
 
192
      do d=1,ndim
 
193
         im = (N(d)-1)*(x(d)-a(d))/(b(d)-a(d)) - ishift
 
194
         ip = im+norder_max-1
 
195
         if (im<0) then
 
196
            im=0 
 
197
            ip=im+norder_max-1
 
198
         end if
 
199
         if (ip>(N(d)-1)) then
 
200
            ip=(N(d)-1)
 
201
            im=ip-norder_max+1
 
202
         end if
 
203
         xm = a(d) + im*(b(d)-a(d))/dble(N(d)-1)
 
204
         xp = a(d) + ip*(b(d)-a(d))/dble(N(d)-1)
 
205
         xtilde = (2.0d0*x(d) - xp - xm)/(xp-xm)
 
206
         jm(d) = im
 
207
         do i=0,norder_max-1
 
208
            dbl_mb(pbasis(1)+i+(d-1)*norder_max) 
 
209
     >         = nwpw_interp_basis(i,xtilde)
 
210
         end do
 
211
      end do
 
212
 
 
213
      do d=1,ndim
 
214
         j(d) = 0
 
215
      end do
 
216
      n1=1
 
217
      do d=1,ndim
 
218
         n1 = n1*norder_max
 
219
      end do
 
220
 
 
221
      f = 0.0d0
 
222
      do i=1,n1
 
223
         nn1 = N(1)
 
224
         ii = jm(1)+j(1)
 
225
         do d=2,ndim
 
226
            ii = ii + (jm(d)+j(d))*nn1
 
227
            nn1 = nn1*N(d)
 
228
         end do
 
229
         fi=1.0d0
 
230
         do d=1,ndim
 
231
            fi = fi*dbl_mb(pbasis(1)+j(d)+(d-1)*norder_max)
 
232
         end do
 
233
         f = f + mesh(ii+1)*fi
 
234
 
 
235
         j(1) = j(1) + 1
 
236
         do d=1,ndim-1
 
237
            if (j(d).ge.norder_max) then
 
238
               j(d) = 0
 
239
               j(d+1) = j(d+1)+1
 
240
            end if
 
241
         end do
 
242
      end do
 
243
 
 
244
      nwpw_interp = f
 
245
      return
 
246
      end
 
247
 
 
248
 
 
249
*     *********************************
 
250
*     *                               *
 
251
*     *          nwpw_dinterp         *
 
252
*     *                               *
 
253
*     *********************************
 
254
      subroutine nwpw_dinterp(ndim,N,a,b,mesh,x,df)
 
255
      implicit none
 
256
      integer ndim,N(*)
 
257
      real*8  a(*),b(*),mesh(*),x(*),df(*)
 
258
 
 
259
#include "mafdecls.fh"
 
260
#include "errquit.fh"
 
261
 
 
262
*     **** nwpw_interp_common block ****
 
263
      integer ndim_max,norder_max
 
264
      integer root(2),coeff(2),pbasis(2)
 
265
      common /nwpw_interp_common/ root,coeff,pbasis,ndim_max,norder_max
 
266
 
 
267
*     **** local variables ****
 
268
      logical failed
 
269
      integer d,dd,n1,nn1,i,ii,im,ip,ishift,k,jm(10),j(10)
 
270
      real*8  xm,xp,xtilde,dx,fi
 
271
 
 
272
*     **** external functions ****
 
273
      real*8   nwpw_interp_basis,nwpw_interp_dbasis
 
274
      external nwpw_interp_basis,nwpw_interp_dbasis
 
275
 
 
276
      failed = .false.
 
277
      do d=1,ndim
 
278
         failed = failed.or.(x(d).lt.a(d))
 
279
         failed = failed.or.(x(d).gt.b(d))
 
280
      end do
 
281
      if (failed) then
 
282
         return
 
283
      end if
 
284
 
 
285
      ishift = norder_max/2 - 1 + mod(norder_max,2)
 
286
 
 
287
      do d=1,ndim
 
288
         im = (N(d)-1)*(x(d)-a(d))/(b(d)-a(d)) - ishift
 
289
         ip = im+norder_max-1
 
290
         if (im<0) then
 
291
            im=0
 
292
            ip=im+norder_max-1
 
293
         end if
 
294
         if (ip>(N(d)-1)) then
 
295
            ip=(N(d)-1)
 
296
            im=ip-norder_max+1
 
297
         end if
 
298
         xm = a(d) + im*(b(d)-a(d))/dble(N(d)-1)
 
299
         xp = a(d) + ip*(b(d)-a(d))/dble(N(d)-1)
 
300
         xtilde = (2.0d0*x(d) - xp - xm)/(xp-xm)
 
301
         dx = 2.0d0/(xp-xm)
 
302
         jm(d) = im
 
303
         do i=0,norder_max-1
 
304
            dbl_mb(pbasis(1)+2*i+  (d-1)*2*norder_max)
 
305
     >         = nwpw_interp_basis(i,xtilde)
 
306
            dbl_mb(pbasis(1)+2*i+1+(d-1)*2*norder_max)
 
307
     >         = dx*nwpw_interp_dbasis(i,xtilde)
 
308
         end do
 
309
      end do
 
310
 
 
311
      do d=1,ndim
 
312
         j(d) = 0
 
313
      end do
 
314
      n1=1
 
315
      do d=1,ndim
 
316
         n1 = n1*norder_max
 
317
      end do
 
318
      do d=1,ndim
 
319
         df(d) = 0.0d0
 
320
      end do
 
321
 
 
322
      do i=1,n1
 
323
         nn1 = N(1)
 
324
         ii = jm(1)+j(1)
 
325
         do d=2,ndim
 
326
            ii = ii + (jm(d)+j(d))*nn1
 
327
            nn1 = nn1*N(d)
 
328
         end do
 
329
         do d=0,ndim-1
 
330
            fi = dbl_mb(pbasis(1)+2*j(d+1)+1+d*2*norder_max)
 
331
            do dd=d+1,(d+ndim-1)
 
332
               k = mod(dd,ndim)
 
333
               fi = fi*dbl_mb(pbasis(1)+2*j(k+1)+k*2*norder_max)
 
334
            end do
 
335
            df(d+1) = df(d+1) + mesh(ii+1)*fi
 
336
         end do
 
337
 
 
338
 
 
339
         j(1) = j(1) + 1
 
340
         do d=1,ndim-1
 
341
            if (j(d).ge.norder_max) then
 
342
               j(d) = 0
 
343
               j(d+1) = j(d+1)+1
 
344
            end if
 
345
         end do
 
346
      end do
 
347
 
 
348
 
 
349
      return
 
350
      end