← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

Viewing changes to src/tools/ga-5-1/global/testing/testeig.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/tools/ga-5-1/global/testing/testeig.F
1
 
#if HAVE_CONFIG_H
2
 
#   include "config.fh"
3
 
#endif
4
 
c $Id: testeig.F,v 1.10 2004-11-13 02:49:18 edo Exp $
5
 
      program test
6
 
      implicit none
7
 
#include "mafdecls.fh"
8
 
#include "global.fh"
9
 
      integer heap, stack
10
 
c#define BLOCK_CYCLIC 0
11
 
c#define PAR_DIAG 1
12
 
c     
13
 
c
14
 
c***  Intitialize a message passing library
15
 
c
16
 
#include "mp3.fh"
17
 
c
18
 
c***  initialize PEIGS
19
 
#ifdef PAR_DIAG
20
 
#ifdef HAVE_SCALAPACK
21
 
c     something here?
22
 
#else
23
 
       call mxinit()   ! PEIGS needs mxinit
24
 
#endif
25
 
#endif
26
 
c 
27
 
c     Intitialize the GA package
28
 
c
29
 
      call ga_initialize()
30
 
c
31
 
c     Initialize the MA package
32
 
c
33
 
      heap = 190000
34
 
      stack= 190000
35
 
      if (.not. ma_init(MT_DBL, heap, stack))
36
 
     $    call ga_error("ma init failed",heap+stack)
37
 
c
38
 
      call testit()
39
 
      call ga_terminate()
40
 
c
41
 
      call MP_FINALIZE() 
42
 
      end
43
 
 
44
 
 
45
 
c-----------------
46
 
 
47
 
      subroutine testit()
48
 
      implicit none
49
 
#include "mafdecls.fh"
50
 
#include "global.fh"
51
 
c     
52
 
      integer n
53
 
      parameter (n = 10)
54
 
      double precision a(n,n), b(n,n), c(n,n), evals(n)
55
 
      double precision eva(n), evb(n)
56
 
      integer g_a,g_b,g_c,g_d
57
 
      integer  i, j, index(n),ind(n)
58
 
      integer nproc, me
59
 
      double precision dsin, sum
60
 
      logical status
61
 
#if BLOCK_CYCLIC
62
 
      integer g_aa, g_bb, g_cc, g_dd
63
 
      integer dims(2), proc_grid(2), block(2)
64
 
      integer g1, g2
65
 
#endif
66
 
c     
67
 
      nproc = ga_nnodes()
68
 
      me    = ga_nodeid()
69
 
c     
70
 
c***  a() is a local copy of what the global array should start as
71
 
c
72
 
      do j = 1, n
73
 
         do i = 1, n
74
 
            a(i,j) = 1d0 * (i+j)  
75
 
            b(i,j) = DSIN(1d0* (i+j))   
76
 
        if(i.eq.j) then
77
 
               b(i,j) = 2d0 *n
78
 
               a(i,j) = i
79
 
            endif
80
 
        if(i.le.j)then
81
 
               c(i,j) = a(i,j)
82
 
            else
83
 
               c(i,j) = 0d0
84
 
            endif
85
 
         enddo
86
 
      enddo
87
 
c
88
 
c***  Create global arrays
89
 
c
90
 
#if BLOCK_CYCLIC
91
 
      if (me.eq.0) then
92
 
        write(6,*) ' '
93
 
        write(6,*) '> Creating Block-Cyclic Arrays'
94
 
        write(6,*) ' '
95
 
      endif
96
 
      dims(1) = n
97
 
      dims(2) = n
98
 
      block(1) = 16
99
 
      block(2) = 16
100
 
      call factor(nproc,g1,g2)
101
 
      proc_grid(1) = g1
102
 
      proc_grid(2) = g2
103
 
      if (me.eq.0) then
104
 
        write(6,'(a,i3,a,i3)') '  Proc grid dimensions: ',g1,' x ',g2
105
 
      endif
106
 
 
107
 
      g_a = ga_create_handle()
108
 
      call ga_set_data(g_a,2,dims,MT_DBL)
109
 
      call ga_set_block_cyclic_proc_grid(g_a,block,proc_grid)
110
 
      if (.not.ga_allocate(g_a))
111
 
     &     call ga_error(' ga_create a failed ', 2)
112
 
 
113
 
      g_b = ga_create_handle()
114
 
      call ga_set_data(g_b,2,dims,MT_DBL)
115
 
      call ga_set_block_cyclic_proc_grid(g_b,block,proc_grid)
116
 
      if (.not.ga_allocate(g_b))
117
 
     &     call ga_error(' ga_create b failed ', 2)
118
 
 
119
 
      g_c = ga_create_handle()
120
 
      call ga_set_data(g_c,2,dims,MT_DBL)
121
 
      call ga_set_block_cyclic_proc_grid(g_c,block,proc_grid)
122
 
      if (.not.ga_allocate(g_c))
123
 
     &     call ga_error(' ga_create c failed ', 2)
124
 
 
125
 
      g_d = ga_create_handle()
126
 
      call ga_set_data(g_d,2,dims,MT_DBL)
127
 
      call ga_set_block_cyclic_proc_grid(g_d,block,proc_grid)
128
 
      if (.not.ga_allocate(g_d))
129
 
     &     call ga_error(' ga_create d failed ', 2)
130
 
 
131
 
      if (.not. ga_create(MT_DBL, n, n, 'AA', 1, 1, g_aa))
132
 
     $     call ga_error(' ga_create aa failed ',2)
133
 
      if (.not. ga_create(MT_DBL, n, n, 'BB', 1, 1, g_bb))
134
 
     $     call ga_error(' ga_create bb failed ',2)
135
 
      if (.not. ga_create(MT_DBL, n, n, 'CC', 1, 1, g_cc))
136
 
     $     call ga_error(' ga_create cc failed ',2)
137
 
      if (.not. ga_create(MT_DBL, n, n, 'DD', 1, 1, g_dd))
138
 
     $     call ga_error(' ga_create dd failed ',2)
139
 
#else
140
 
      if (.not. ga_create(MT_DBL, n, n, 'A', 1, 1, g_a))
141
 
     $     call ga_error(' ga_create a failed ',2)
142
 
      if (.not. ga_create(MT_DBL, n, n, 'B', 1, 1, g_b))
143
 
     $     call ga_error(' ga_create b failed ',2)
144
 
      if (.not. ga_create(MT_DBL, n, n, 'C', 1, 1, g_c))
145
 
     $     call ga_error(' ga_create c failed ',2)
146
 
      if (.not. ga_create(MT_DBL, n, n, 'D', 1, 1, g_d))
147
 
     $     call ga_error(' ga_create d failed ',2)
148
 
#endif
149
 
c     
150
 
c***  Fill in arrays A & B
151
 
c
152
 
      if (me .eq. 0) then
153
 
c        write(6,21) 
154
 
 21      format(/' filling in ... ')
155
 
c        call ffflush(6)
156
 
         do j = 1, n
157
 
        call ga_put(g_a, 1,n, j,j, a(1,j),n)
158
 
        call ga_put(g_b, 1,n, j,j, b(1,j),n)
159
 
        call ga_put(g_c, 1,n, j,j, c(1,j),n)
160
 
         enddo
161
 
c    print *,'A'
162
 
        do j = 1, n
163
 
c            call GA_GET(g_a, 1,n, j,j, eva,1)
164
 
c         write(6,'(10e8.2)')(eva(i),i=1,n)
165
 
        enddo
166
 
      endif
167
 
c
168
 
c***  check symmetrization
169
 
c
170
 
      if (me .eq. 0) then
171
 
        print *,' '
172
 
        print *,'>checking ga_symmetrize '
173
 
        print *,' '
174
 
        call ffflush(6)
175
 
      endif
176
 
      call ga_symmetrize(g_c)
177
 
c
178
 
      call GA_GET(g_c,  1,n, 1,n,c,n)
179
 
      do j = ga_nodeid()+1, n, ga_nnodes()
180
 
         do i = j+1, n
181
 
            if(c(i,j).ne.c(j,i))then
182
 
                 print *, me, ' symmetrize ',i,j,c(i,j),c(j,i)
183
 
                 call ffflush(6)
184
 
                 call ga_error('exiting',-1)
185
 
            endif
186
 
         enddo
187
 
      enddo
188
 
      call ga_sync()
189
 
      if (me .eq. 0) then
190
 
        print *,' '
191
 
        print *,' ga_symmetrize is OK'
192
 
        print *,' '
193
 
        call ffflush(6)
194
 
      endif
195
 
 
196
 
c
197
 
c***  check symmetrization 
198
 
c
199
 
      if (me .eq. 0) then
200
 
        print *,' '
201
 
        print *,'>checking ga_transpose '
202
 
        print *,' '
203
 
        call ffflush(6)
204
 
      endif
205
 
c
206
 
      call ga_transpose(g_c,g_d)
207
 
*     call ga_print(g_c)
208
 
*     call ga_print(g_d)
209
 
      call GA_GET(g_d,  1,n, 1,n,a,n)
210
 
      do j = ga_nodeid()+1, n, ga_nnodes()
211
 
         call GA_GET(g_d,  1,n, j,j, a,n)
212
 
         do i = 1, n
213
 
            if(a(i,1).ne.c(j,i))then
214
 
                 print *, me, ' transpose ',i,j,a(i,1),c(j,i) 
215
 
                 call ffflush(6)
216
 
                 call ga_error('exiting',-1)
217
 
            endif
218
 
         enddo
219
 
      enddo
220
 
      call ga_sync()
221
 
      if (me .eq. 0) then
222
 
        print *,' '
223
 
        print *,' ga_transpose is OK'
224
 
        print *,' '
225
 
        call ffflush(6)
226
 
      endif
227
 
c
228
 
c
229
 
c***  solve the eigenproblem
230
 
      if (me .eq. 0)then
231
 
        print *,' '
232
 
        write(6,*) '>checking the generalized eigensolver ... '
233
 
        print *,' '
234
 
        call ffflush(6)
235
 
      endif
236
 
      call ga_sync()
237
 
#ifndef PAR_DIAG
238
 
      call ga_diag_seq(g_a,g_b,g_c,evals)
239
 
#else
240
 
#ifdef HAVE_SCALAPACK
241
 
#if BLOCK_CYCLIC
242
 
      call ga_copy(g_a, g_aa)
243
 
      call ga_copy(g_b, g_bb)
244
 
#endif
245
 
      call ga_pdsygv( g_a,  g_b, g_c, evals)
246
 
#if BLOCK_CYCLIC
247
 
      call ga_copy(g_aa, g_a)
248
 
      call ga_copy(g_bb, g_b)
249
 
#endif
250
 
#else
251
 
      call ga_diag(g_a,g_b,g_c,evals)
252
 
#endif
253
 
#endif
254
 
      if (me .eq. 0) then
255
 
        print *,' '
256
 
        print *,' checking multiplication'
257
 
        print *,' '
258
 
        call ffflush(6)
259
 
      endif
260
 
c
261
 
      call ga_sync()
262
 
#if BLOCK_CYCLIC
263
 
      call ga_copy(g_a, g_aa)
264
 
      call ga_copy(g_c, g_cc)
265
 
      call ga_dgemm('t','n',n,n,n, 1d0, g_cc, g_aa, 0d0, g_dd)
266
 
      call ga_dgemm('n','n',n,n,n, 1d0, g_dd, g_cc, 0d0, g_aa)
267
 
      call ga_copy(g_aa, g_a)
268
 
#else
269
 
      call ga_dgemm('t','n',n,n,n, 1d0, g_c, g_a, 0d0, g_d)
270
 
      call ga_dgemm('n','n',n,n,n, 1d0, g_d, g_c, 0d0, g_a)
271
 
#endif
272
 
c
273
 
      call ga_sync()
274
 
      if (me .eq. 0) then
275
 
         do j = 1, n
276
 
        call GA_GET(g_a, j,j, j,j, eva(j),1)
277
 
         enddo
278
 
      endif
279
 
 
280
 
      call ga_sync()
281
 
#if BLOCK_CYCLIC
282
 
      call ga_copy(g_b, g_bb)
283
 
      call ga_dgemm('t','n',n,n,n, 1d0, g_cc, g_bb, 0d0, g_dd)
284
 
      call ga_sync()
285
 
      call ga_dgemm('n','n',n,n,n, 1d0, g_dd, g_cc, 0d0, g_bb)
286
 
      call ga_copy(g_bb, g_b)
287
 
#else
288
 
      call ga_dgemm('t','n',n,n,n, 1d0, g_c, g_b, 0d0, g_d)
289
 
      call ga_sync()
290
 
      call ga_dgemm('n','n',n,n,n, 1d0, g_d, g_c, 0d0, g_b)
291
 
#endif
292
 
c
293
 
      call ga_sync()
294
 
      if (me .eq. 0) then
295
 
         do j = 1, n
296
 
        call GA_GET(g_b, j,j, j,j, evb(j),1)
297
 
         enddo
298
 
         write(6,*)'  j   lambda      eva      evb      eva/evb' 
299
 
         write(6,*)'  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^' 
300
 
         call ffflush(6)
301
 
         do j = 1, n
302
 
            if(ABS(evals(j) - eva(j)/evb(j)).gt.1d-5)
303
 
     &         write(6,'(i4,1h_,6(e10.3,1h,))')
304
 
     @            j,evals(j), eva(j), evb(j),eva(j)/evb(j)
305
 
         enddo
306
 
         write(6,*)'       OK         OK       OK        OK' 
307
 
         call ffflush(6)
308
 
      endif
309
 
      if (me .eq. 0) then
310
 
        print *,' '
311
 
        print *,'  eigensolver & multiplication are OK'
312
 
        print *,' '
313
 
        print *,' '
314
 
        call ffflush(6)
315
 
      endif
316
 
c
317
 
c..................................................................
318
 
c
319
 
c***  solve the std eigenproblem
320
 
      if (me .eq. 0)then
321
 
        print *,' '
322
 
        write(6,*) '>checking the standard eigensolver ... '
323
 
        print *,' '
324
 
        call ffflush(6)
325
 
      endif
326
 
      do j =1,n
327
 
         index(j)=j
328
 
         ind(j)=j
329
 
      enddo
330
 
 
331
 
      call ga_sync()
332
 
#ifdef PAR_DIAG
333
 
#ifdef HAVE_SCALAPACK
334
 
#if BLOCK_CYCLIC
335
 
      call ga_copy(g_a, g_aa)
336
 
#endif
337
 
      call ga_pdsyev( g_a,  g_c, evals, 16)
338
 
#if BLOCK_CYCLIC
339
 
      call ga_copy(g_aa, g_a)
340
 
#endif
341
 
#else
342
 
      call ga_diag_std(g_a,g_c,evals)
343
 
#endif
344
 
#else
345
 
      call ga_diag_std_seq(g_a,g_c,evals)
346
 
#endif
347
 
c
348
 
 
349
 
      call ga_zero(g_b)
350
 
      
351
 
#if BLOCK_CYCLIC
352
 
      call ga_copy(g_a, g_aa)
353
 
      call ga_copy(g_c, g_cc)
354
 
      call ga_dgemm('n','n',n,n,n, 1d0, g_aa, g_cc, 0d0, g_dd) ! d := a*c
355
 
      call ga_copy(g_dd, g_d)
356
 
#else
357
 
      call ga_dgemm('n','n',n,n,n, 1d0, g_a, g_c, 0d0, g_d) ! d := a*c
358
 
#endif
359
 
c
360
 
c  copy eigenvalues to diagonal of g_b
361
 
c
362
 
      if (me .eq. 0) call ga_scatter(g_b, evals, index, ind, n)
363
 
 
364
 
      call ga_sync()
365
 
#if BLOCK_CYCLIC
366
 
      call ga_copy(g_b, g_bb)
367
 
      call ga_dgemm('n','n',n,n,n, 1d0, g_cc, g_bb, 0d0, g_aa) ! a := c*b
368
 
      call ga_copy(g_aa, g_a)
369
 
#else
370
 
      call ga_dgemm('n','n',n,n,n, 1d0, g_c, g_b, 0d0, g_a) ! a := c*b
371
 
#endif
372
 
 
373
 
      call ga_sync()
374
 
      call ga_add(1d0, g_d, -1d0, g_a, g_c)
375
 
      sum = ga_ddot(g_c,g_c)
376
 
      if (me .eq. 0) then
377
 
        if(dsqrt(sum)/n.lt.1d-11)then
378
 
           print *, ' std eigensolver is OK '
379
 
        else
380
 
           print *, ' test failed: norm = ', dsqrt(sum)/n
381
 
        endif
382
 
        print *,' '
383
 
        call ffflush(6)
384
 
      endif
385
 
c     status =  MA_summarize_allocated_blocks()
386
 
      status =  ga_destroy(g_d)
387
 
      status =  ga_destroy(g_c)
388
 
      status =  ga_destroy(g_b)
389
 
      status =  ga_destroy(g_a)
390
 
      end
391
 
 
392
 
#if BLOCK_CYCLIC
393
 
      subroutine factor(p,idx,idy)
394
 
      implicit none
395
 
      integer i,j,p,idx,idy,it
396
 
      integer ip,ifac,pmax,prime(1280)
397
 
      integer fac(1280)
398
 
c
399
 
      i = 1
400
 
      ip = p
401
 
c
402
 
c    factor p completely
403
 
c    first, find all prime numbers less than or equal to p
404
 
c
405
 
      pmax = 0
406
 
      do i = 2, p
407
 
        do j = 1, pmax
408
 
          if (mod(i,prime(j)).eq.0) go to 100
409
 
        end do
410
 
        pmax = pmax + 1
411
 
        prime(pmax) = i
412
 
  100   continue
413
 
      end do
414
 
c
415
 
c    find all prime factors of p
416
 
c
417
 
      ifac = 0
418
 
      do i = 1, pmax
419
 
  200   if (mod(ip,prime(i)).eq.0) then
420
 
          ifac = ifac + 1
421
 
          fac(ifac) = prime(i)
422
 
          ip = ip/prime(i)
423
 
          go to 200
424
 
        endif
425
 
      end do
426
 
c
427
 
c    determine two factors of p of approximately the
428
 
c    same size
429
 
c
430
 
      idx = 1
431
 
      idy = 1
432
 
      do i = ifac, 1, -1
433
 
        if (idx.le.idy) then
434
 
          idx = fac(i)*idx
435
 
        else
436
 
          idy = fac(i)*idy
437
 
        endif
438
 
      end do
439
 
      return
440
 
      end
441
 
#endif