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