1c $Id$
2      subroutine util_test_eig
3      implicit none
4#include "mafdecls.fh"
5#include "global.fh"
6c     integer heap, stack
7c
8c
9c***  Intitialize a message passing library
10c
11#ifdef MPI
12c     integer ierr
13c     call mpi_init(ierr)
14#else
15c     call pbeginf
16#endif
17c
18c***  initialize PEIGS
19#ifdef PAR_DIAG
20c      call mxinit()   ! PEIGS needs mxinit
21#endif
22c
23c     Intitialize the GA package
24c
25c     call ga_initialize()
26c
27c     Initialize the MA package
28c
29c     heap = 190000
30c     stack= 190000
31c     if (.not. ma_init(MT_DBL, heap, stack))
32c    $    call ga_error("ma init failed",heap+stack)
33c
34      call testit5()
35c     call ga_terminate()
36c
37#ifdef MPI
38c     call mpi_finalize(ierr)
39#else
40c     call pend()
41#endif
42      end
43
44
45c-----------------
46
47      subroutine testit5()
48      implicit none
49#include "mafdecls.fh"
50#include "global.fh"
51c
52      integer n
53      parameter (n = 100)
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
61c
62      nproc = ga_nnodes()
63      me    = ga_nodeid()
64c
65c***  a() is a local copy of what the global array should start as
66c
67      do j = 1, n
68         do i = 1, n
69            a(i,j) = 1d0 * (i+j)
70            b(i,j) = DSIN(1d0* (i+j))
71	    if(i.eq.j) then
72               b(i,j) = 2d0 *n
73               a(i,j) = i
74            endif
75	    if(i.le.j)then
76               c(i,j) = a(i,j)
77            else
78               c(i,j) = 0d0
79            endif
80         enddo
81      enddo
82c
83c***  Create global arrays
84c
85      if (.not. ga_create(MT_DBL, n, n, 'A', 1, 1, g_a))
86     $     call ga_error(' ga_create failed ',2)
87      if (.not. ga_create(MT_DBL, n, n, 'B', 1, 1, g_b))
88     $     call ga_error(' ga_create failed ',2)
89      if (.not. ga_create(MT_DBL, n, n, 'C', 1, 1, g_c))
90     $     call ga_error(' ga_create failed ',2)
91      if (.not. ga_create(MT_DBL, n, n, 'D', 1, 1, g_d))
92     $     call ga_error(' ga_create failed ',2)
93c
94c***  Fill in arrays A & B
95c
96      if (me .eq. 0) then
97c        write(6,21)
98c21      format(/' filling in ... ')
99c        call ffflush(6)
100         do j = 1, n
101	    call ga_put(g_a, 1,n, j,j, a(1,j),n)
102	    call ga_put(g_b, 1,n, j,j, b(1,j),n)
103	    call ga_put(g_c, 1,n, j,j, c(1,j),n)
104         enddo
105c	print *,'A'
106        do j = 1, n
107       	  call GA_GET(g_a, 1,n, j,j, eva,1)
108c         write(6,'(10e8.2)')(eva(i),i=1,n)
109        enddo
110      endif
111c
112c***  check symmetrization
113c
114c     if (me .eq. 0) then
115c       print *,' '
116c	print *,'>checking ga_symmetrize '
117c       print *,' '
118c       call ffflush(6)
119c     endif
120      call ga_symmetrize(g_c)
121c
122      call GA_GET(g_c,  1,n, 1,n,c,n)
123      do j = ga_nodeid()+1, n, ga_nnodes()
124         do i = j+1, n
125            if(c(i,j).ne.c(j,i))then
126                 print *, me, ' symmetrize ',i,j,c(i,j),c(j,i)
127                 call ffflush(6)
128                 call ga_error('exiting',-1)
129            endif
130         enddo
131      enddo
132      call ga_sync()
133      if (me .eq. 0) then
134        write(6,'(A)') ' ga_symmetrize .................... OK'
135c       print *,' '
136c       print *,' ga_symmetrize is OK'
137c       print *,' '
138        call ffflush(6)
139      endif
140
141c
142c***  check symmetrization
143c
144c     if (me .eq. 0) then
145c	print *,' '
146c       print *,'>checking ga_transpose '
147c	print *,' '
148c       call ffflush(6)
149c     endif
150c
151      call ga_transpose(g_c,g_d)
152*     call ga_print(g_c)
153*     call ga_print(g_d)
154      call GA_GET(g_d,  1,n, 1,n,a,n)
155      do j = ga_nodeid()+1, n, ga_nnodes()
156         call GA_GET(g_d,  1,n, j,j, a,n)
157         do i = 1, n
158            if(a(i,1).ne.c(j,i))then
159                 print *, me, ' transpose ',i,j,a(i,1),c(j,i)
160                 call ffflush(6)
161                 call ga_error('exiting',-1)
162            endif
163         enddo
164      enddo
165      call ga_sync()
166      if (me .eq. 0) then
167        write(6,'(A)') ' ga_transpose ..................... OK'
168c	print *,' '
169c	print *,' ga_transpose is OK'
170c	print *,' '
171        call ffflush(6)
172      endif
173c
174c
175c***  solve the eigenproblem
176c     if (me .eq. 0)then
177c	print *,' '
178c       write(6,*) '>checking the generalized eigensolver ... '
179c	print *,' '
180c       call ffflush(6)
181c     endif
182      call ga_sync()
183#ifndef PAR_DIAG
184      call ga_diag_seq(g_a,g_b,g_c,evals)
185#else
186      call ga_diag(g_a,g_b,g_c,evals)
187#endif
188c     if (me .eq. 0) then
189c	print *,' '
190c	print *,' checking multiplication'
191c	print *,' '
192c       call ffflush(6)
193c     endif
194c
195      call ga_sync()
196      call ga_dgemm('t','n',n,n,n, 1d0, g_c, g_a, 0d0, g_d)
197      call ga_dgemm('n','n',n,n,n, 1d0, g_d, g_c, 0d0, g_a)
198c
199      call ga_sync()
200      if (me .eq. 0) then
201         do j = 1, n
202	    call GA_GET(g_a, j,j, j,j, eva(j),1)
203         enddo
204      endif
205
206      call ga_sync()
207      call ga_dgemm('t','n',n,n,n, 1d0, g_c, g_b, 0d0, g_d)
208      call ga_sync()
209      call ga_dgemm('n','n',n,n,n, 1d0, g_d, g_c, 0d0, g_b)
210c
211      call ga_sync()
212      if (me .eq. 0) then
213         do j = 1, n
214	    call GA_GET(g_b, j,j, j,j, evb(j),1)
215         enddo
216c        write(6,*)'  j   lambda      eva      evb      eva/evb'
217c        write(6,*)'  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
218         call ffflush(6)
219         do j = 1, n
220            if(ABS(evals(j) - eva(j)/evb(j)).gt.1d-5)
221     &         write(6,'(i4,1h_,6(e10.3,1h,))')
222     @            j,evals(j), eva(j), evb(j),eva(j)/evb(j)
223         enddo
224c        write(6,*)'       OK         OK       OK        OK'
225        write(6,'(A)') ' ga_dgemm ......................... OK'
226        write(6,'(A)') ' ga_diag .......................... OK'
227         call ffflush(6)
228      endif
229      if (me .eq. 0) then
230c	print *,' '
231c	print *,'  eigensolver & multiplication are OK'
232c	print *,' '
233c	print *,' '
234        call ffflush(6)
235      endif
236c
237c..................................................................
238c
239c***  solve the std eigenproblem
240c     if (me .eq. 0)then
241c       print *,' '
242c       write(6,*) '>checking the standard eigensolver ... '
243c       print *,' '
244c       call ffflush(6)
245c     endif
246      do j =1,n
247         index(j)=j
248         ind(j)=j
249      enddo
250
251      call ga_sync()
252#ifdef PAR_DIAG
253      call ga_diag_std(g_a,g_c,evals)
254#else
255      call ga_diag_std_seq(g_a,g_c,evals)
256#endif
257c
258      call ga_sync()
259
260      call ga_zero(g_b)
261
262      call ga_dgemm('n','n',n,n,n, 1d0, g_a, g_c, 0d0, g_d) ! d := a*c
263c
264c
265      if (me .eq. 0) call ga_scatter(g_b, evals, index, ind, n)
266
267      call ga_sync()
268      call ga_dgemm('n','n',n,n,n, 1d0, g_c, g_b, 0d0, g_a) ! a := c*b
269
270      call ga_sync()
271      call ga_add(1d0, g_d, -1d0, g_a, g_c)
272      sum = ga_ddot(g_c,g_c)
273      if (me .eq. 0) then
274c       print *,' '
275        if(dsqrt(sum)/n.lt.1d-11)then
276        write(6,'(A)') ' ga_diag_std ...................... OK'
277c          print *, ' std eigensolver is OK '
278        else
279           print *, ' test failed: norm = ', dsqrt(sum)/n
280        endif
281c       print *,' '
282        call ffflush(6)
283      endif
284c     status =  MA_summarize_allocated_blocks()
285      status =  ga_destroy(g_d)
286      status =  ga_destroy(g_c)
287      status =  ga_destroy(g_b)
288      status =  ga_destroy(g_a)
289      end
290
291