1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
|
#if HAVE_CONFIG_H
# include "config.fh"
#endif
#define xgemm TEST_DGEMM
#define ygemm TEST_ZGEMM
program ga_test
c $Id: pg2testmatmult.F,v 1.3.6.1 2007-04-25 21:49:59 manoj Exp $
c
c test ga_dgemm
c compile with make FLD_REN="optimized blas" testmatmult.x
c Note: - change nummax for large arrays
c - turn off "verify" for large arrays due to memory
c limitations, as verify=TRUE for large arrays produces
c segfault, dumps core,or any crap.
c
c
implicit none
integer num_m,num_n,num_k,nummax,howmany,ntrans
parameter (nummax=1024,howmany=2,ntrans=4)
integer num1
parameter(num1=nummax)
integer i,ii
real*8 h0(num1*num1)
integer g_c,g_b,g_a
real*8 a
real*8 t1,mf,avg_t(ntrans),avg_mf(ntrans)
integer itime,ntimes,me
integer nums_m(howmany),nums_n(howmany),nums_k(howmany)
character*1 transa(ntrans),transb(ntrans),ta,tb
real *8 tmpa(nummax,nummax), tmpb(nummax,nummax)
real *8 tmpc(nummax,nummax)
logical verify, status
integer nproc, nprocs, proc_group(0:100), proclist(100)
integer grp_id, grp_me, grp_nproc, ndim, dims(2)
data transa/'n','t','n','t'/
data transb/'n','n','t','t'/
data nums_m/512,1024/
data nums_n/512,1024/
data nums_k/512,1024/
#include "mafdecls.fh"
#include "global.fh"
#include "testutil.fh"
#include "mp3.fh"
if (.not.ma_init(MT_DBL,1,20000000)) then
call ga_error('failed: ma_init(MT_DBL,1,20000000)',10)
endif
call ga_initialize()
me = ga_nodeid()
nproc = ga_nnodes()
if(nproc .ne. 8) then
write(*,*) 'nproc =', nproc
call ga_error('failed: need exactly 8 processors',10)
endif
call ga_sync()
proclist(1) = 0
proclist(2) = 1
proclist(3) = 2
proclist(4) = 3
nprocs = 4
proc_group(0) = ga_pgroup_create(proclist, nprocs)
if (me.lt.4) then
c get new proc ids for your default group
grp_id = proc_group(0)
grp_me = ga_pgroup_nodeid(grp_id)
grp_nproc = ga_pgroup_nnodes(grp_id)
write(*,*) 'GRP INFO =',me,grp_me,grp_nproc
ii = 0
do i = 1, num1*num1
ii = ii + 1
if(ii.gt.num1) then
ii = 0
endif
h0(i) = ii
enddo
c
c Compute times assuming 500 mflops and 5 second target time
c
c ntimes = max(3.0d0,5.0d0/(4.0d-9*num**3))
ntimes = 5
verify = .TRUE. ! to verify ga_dgemm
c verify = .FALSE. ! to verify ga_dgemm
do ii=1,howmany
num_m = nums_m(ii)
num_n = nums_n(ii)
num_k = nums_k(ii)
a = 0.5d0/(num_m*num_n)
if(num_m.gt.nummax .OR. num_n.gt.nummax .OR.
& num_k.gt.nummax) then
call ga_error('Insufficient memory: check nummax', 0)
endif
c
c
ndim = 2
c Creating C array
g_c = ga_create_handle()
dims(1) = num_m
dims(2) = num_n
call ga_set_data(g_c,ndim,dims,MT_DBL)
call ga_set_array_name(g_c,'c')
call ga_set_pgroup(g_c,grp_id)
status = ga_allocate(g_c)
if (.not. status) then
write(6,*) ' g_c ga_create failed'
call ffflush(6)
call ga_error('... exiting ',0)
endif
c Creating B array
g_b = ga_create_handle()
dims(1) = num_k
dims(2) = num_n
call ga_set_data(g_b,ndim,dims,MT_DBL)
call ga_set_array_name(g_b,'b')
call ga_set_pgroup(g_b,grp_id)
status = ga_allocate(g_b)
if (.not. status) then
write(6,*) ' g_b ga_create failed'
call ffflush(6)
call ga_error('... exiting ',0)
endif
c Creating A array
g_a = ga_create_handle()
dims(1) = num_m
dims(2) = num_k
call ga_set_data(g_a,ndim,dims,MT_DBL)
call ga_set_array_name(g_a,'a')
call ga_set_pgroup(g_a,grp_id)
status = ga_allocate(g_a)
if (.not. status) then
write(6,*) ' g_a ga_create failed'
call ffflush(6)
call ga_error('... exiting ',0)
endif
c
c Initialize matrices A and B
c call load_ga_from_square(g_a,num_m,h0,num1)
c call load_ga_from_square(g_b,num_m,h0,num1)
c
if(grp_me.eq.0) then
call load_ga(g_a, h0, num_m, num_k)
call load_ga(g_b, h0, num_k, num_n)
endif
call ga_zero(g_c)
call ga_pgroup_sync(grp_id)
if (ga_pgroup_nodeid(grp_id).eq.0) then
write(*,*) ! for new line
write(*,*) 'Matrix Multiplication C = A[', num_m, ',',
. num_k, '] x B[', num_k, ',', num_n, ']'
write(*,*) ! for new line
call ffflush(6)
endif
do i = 1, ntrans
avg_t(i) = 0.0d0
avg_mf(i) = 0.0d0
enddo
do itime = 1, ntimes
do i = 1, ntrans
call ga_pgroup_sync(grp_id)
ta = transa(i)
tb = transb(i)
t1 = util_timer()
call ga_dgemm(ta,tb,num_m,num_n,num_k,
& 1.0d0, g_a, g_b, 0.0d0, g_c)
t1 = util_timer() - t1
if (grp_me.eq.0) then
mf = 2d0*num_m*num_n*num_k/t1*1d-6/grp_nproc
avg_t(i) = avg_t(i)+t1
avg_mf(i) = avg_mf(i) + mf
write(6,200) ' Run#', itime, t1, mf, ta, tb
call ffflush(6)
if (verify .AND. itime.eq.1) then
call verify_ga_dgemm(ta, tb,num_m, num_n,num_k,
& 1.0d0, g_a,g_b,0.0d0,g_c,tmpa,tmpb,tmpc)
endif
endif
enddo
enddo
if (ga_pgroup_nodeid(grp_id).eq.0) then
write(*,*) ! for new line
do i = 1, ntrans
write(6,200) 'Average:',0,avg_t(i)/ntimes,
. avg_mf(i)/ntimes,transa(i),transb(i)
enddo
if(verify) write(*,*) 'All ga_dgemms are verified...O.K.'
write(*,*) ! for new line
call ffflush(6)
endif
c
c call ga_print(g_a)
c call ga_print(g_b)
c call ga_print(g_c)
c
if (.not.ga_destroy(g_c)) then
call ga_error('failed: destroy g_c',20)
endif
if (.not.ga_destroy(g_b)) then
call ga_error('failed: destroy g_b',30)
endif
if (.not.ga_destroy(g_a)) then
call ga_error('failed: destroy g_a',40)
endif
enddo
endif
200 format(a15, i2, ': ', e12.4, ' seconds ',f12.1,
. ' mflops/proc ', 3a2)
call ga_sync()
call ga_terminate
call MP_FINALIZE()
c
end
c
c-----------------------------------------------------------------------
c Verify for correctness. Process 0 computes BLAS dgemm
c locally. For larger arrays, disbale this test as memory
c might not be sufficient
c
subroutine verify_ga_dgemm(xt1, xt2, num_m, num_n, num_k,
$ alpha, g_a, g_b, beta, g_c, tmpa, tmpb, tmpc)
implicit none
character *1 xt1, xt2
integer num_m, num_n, num_k, g_a, g_b, g_c
double precision alpha, beta
real *8 tmpa(num_m,num_k), tmpb(num_k,num_n), tmpc(num_m,num_n)
integer i,j,type,dim1, dim2
real *8 abs_value
#include "mafdecls.fh"
#include "testutil.fh"
do i = 1,num_n
do j = 1, num_m
tmpc(j,i) = -1.0
tmpa(j,i) = -2.0
enddo
enddo
call ga_inquire(g_a, type, dim1, dim2)
call ga_get(g_a, 1, dim1, 1, dim2, tmpa, dim1)
call ga_inquire(g_b, type, dim1, dim2)
call ga_get(g_b, 1, dim1, 1, dim2, tmpb, dim1)
c compute dgemm sequentially
call xgemm(xt1, xt2, num_m, num_n, num_k, alpha,
& tmpa, num_m, tmpb, num_k, beta, tmpc, num_m)
c after computing c locally, verify it with the values in g_c
call ga_inquire(g_c, type, dim1, dim2)
call ga_get(g_c, 1, dim1, 1, dim2, tmpa, dim1)
do i = 1,num_n
do j = 1, num_m
abs_value = abs(tmpc(j,i)-tmpa(j,i))
if(abs_value .gt. 1.0 .OR. abs_value .lt. -1.0) then
write(*,*) 'Values are = ', tmpc(j,i), tmpa(j,i)
write(*,*) 'Values are = ', abs(tmpc(j,i)-tmpa(j,i)),
. abs_value
call ffflush(6)
call ga_error('verify ga_dgemm failed', 0)
endif
enddo
enddo
return
end
c
c-----------------------------------------------------------------------
c called by process '0' (or your master process )
c
subroutine load_ga(handle,f, dim1,dim2)
implicit none
integer handle
integer dim1,dim2,i
real*8 f(dim1,dim2)
#include "mafdecls.fh"
if(dim1.le.0)return
if(dim2.le.0)return
call ga_put(handle, 1, dim1, 1, dim2, f, dim1)
return
end
c
c-----------------------------------------------------------------------
c must be called by all processors, if you need to fillup the
c entire array
c
subroutine load_ga_from_square(handle,num,f,ndim)
implicit none
integer handle, memhandle
integer num,ndim
real*8 f(ndim,ndim)
integer ilo, ihi, jlo, jhi, nx, ny, ibuff
integer ga_nodeid, i1, i2, i, j, ix, jx
#include "mafdecls.fh"
call ga_distribution(handle, ga_nodeid(), ilo, ihi, jlo, jhi)
if(ihi.le.0)return
if(jhi.le.0)return
c nx = ihi - ilo + 1
c ny = jhi - jlo + 1
do i = ilo,ihi,ndim
do j = jlo,jhi,ndim
call ga_put(handle,i,min(ihi,i+ndim),j,min(jhi,j+ndim),
& f,ndim)
enddo
enddo
return
end
c
c-----------------------------------------------------------------------
c must be called by all processors, if you need to fillup the
c entire array
c
subroutine load_ga_from_triangle(handle,f,ndim)
implicit none
integer handle, memhandle
real*8 f(*)
integer ndim
integer ilo, ihi, jlo, jhi, nx, ny, ibuff
integer ga_nodeid, i1, i2, i, j, ix, jx
#include "mafdecls.fh"
call ga_distribution(handle, ga_nodeid(), ilo, ihi, jlo, jhi)
if(ihi.le.0)return
if(jhi.le.0)return
nx = ihi - ilo + 1
ny = jhi - jlo + 1
if (.not.ma_alloc_get(MT_DBL,nx*ny,'flap',memhandle,ibuff)) then
call ga_error('failed: allocate triangle',100)
endif
do i = 1,nx
do j = 1,ny
ix = i + ilo - 1
jx = j + jlo - 1
i1 = min(ix,jx)
i2 = max(ix,jx)
dbl_mb(ibuff + nx*(j-1) + (i-1) ) = f(i2*(i2-1)/2 + i1)
enddo
enddo
call ga_put(handle,ilo,ihi,jlo,jhi,
& dbl_mb(ibuff),nx)
if (.not.ma_free_heap(memhandle)) then
call ga_error('failed: free triangle',100)
endif
return
end
|
