SND@LHC Software
Loading...
Searching...
No Matches
mpbits.f90
Go to the documentation of this file.
1
33
35MODULE mpbits
36 USE mpdef
37 IMPLICIT NONE
38
39 INTEGER(mpl) :: ndimb
40 INTEGER(mpl) :: ndimb2
41 INTEGER(mpi) :: n
42 INTEGER(mpi) :: n2
43 INTEGER(mpi) :: ibfw
44 INTEGER(mpi) :: ireqpe
45 INTEGER(mpi) :: isngpe
46 INTEGER(mpi) :: icmprs
47 INTEGER(mpi) :: iextnd
48 INTEGER(mpi) :: nspc
49 INTEGER(mpi) :: mxcnt
50 INTEGER(mpi) :: nencdm
51 INTEGER(mpi) :: nencdb
52 INTEGER(mpi) :: nthrd
53 INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: bitfieldcounters
54 INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: bitmap
55 INTEGER(mpi), PARAMETER :: bs = bit_size(1_mpi)
56
57END MODULE mpbits
58
65SUBROUTINE inbits(im,jm,inc) ! include element (I,J)
66 USE mpbits
67
68 INTEGER(mpi), INTENT(IN) :: im
69 INTEGER(mpi), INTENT(IN) :: jm
70 INTEGER(mpi), INTENT(IN) :: inc
71
72 INTEGER(mpl) :: l
73 INTEGER(mpl) :: ll
74 INTEGER(mpi) :: i
75 INTEGER(mpi) :: j
76 INTEGER(mpi) :: noffj
77 INTEGER(mpi) :: m
78 INTEGER(mpi) :: mm
79 INTEGER(mpi) :: icount
80 INTEGER(mpi) :: ib
81 INTEGER(mpi) :: jcount
82 INTEGER(mpl) :: noffi
83 LOGICAL :: btest
84
85 IF(im == jm) RETURN ! diagonal
86 j=min(im,jm)
87 i=max(im,jm)
88 IF(j <= 0) RETURN ! out low
89 IF(i > n) RETURN ! out high
90 noffi=int(i-1,mpl)*int(i-2,mpl)*int(ibfw,mpl)/2 ! for J=1
91 noffj=(j-1)*ibfw
92 l=noffi/bs+i+noffj/bs ! row offset + column offset
93 ! add I instead of 1 to keep bit maps of different rows in different words (openMP !)
94 m=mod(noffj,bs)
95 IF (ibfw <= 1) THEN
96 bitfieldcounters(l)=ibset(bitfieldcounters(l),m)
97 ELSE
98 ! get counter from bit field
99 ll=l
100 mm=m
101 icount=0
102 DO ib=0,ibfw-1
103 IF (btest(bitfieldcounters(ll),mm)) icount=ibset(icount,ib)
104 mm=mm+1
105 IF (mm >= bs) THEN
106 ll=ll+1
107 mm=mm-bs
108 END IF
109 END DO
110 ! increment
111 jcount=icount
112 icount=min(icount+inc,mxcnt)
113 ! store counter into bit field
114 IF (icount /= jcount) THEN
115 ll=l
116 mm=m
117 DO ib=0,ibfw-1
118 IF (btest(icount,ib)) THEN
119 bitfieldcounters(ll)=ibset(bitfieldcounters(ll),mm)
120 ELSE
121 bitfieldcounters(ll)=ibclr(bitfieldcounters(ll),mm)
122 END IF
123 mm=mm+1
124 IF (mm >= bs) THEN
125 ll=ll+1
126 mm=mm-bs
127 END IF
128 END DO
129 END IF
130 END IF
131 RETURN
132
133END SUBROUTINE inbits
134
147SUBROUTINE clbits(in,jreqpe,jhispe,jsngpe,jcmprs,jextnd,idimb,iencdb,ispc)
148 USE mpbits
149 USE mpdalc
150
151 INTEGER(mpi), INTENT(IN) :: in
152 INTEGER(mpi), INTENT(IN) :: jreqpe
153 INTEGER(mpi), INTENT(IN) :: jhispe
154 INTEGER(mpi), INTENT(IN) :: jsngpe
155 INTEGER(mpi), INTENT(IN) :: jcmprs
156 INTEGER(mpi), INTENT(IN) :: jextnd
157 INTEGER(mpl), INTENT(OUT) :: idimb
158 INTEGER(mpi), INTENT(OUT) :: iencdb
159 INTEGER(mpi), INTENT(OUT) :: ispc
160
161 INTEGER(mpl) :: noffd
162 INTEGER(mpi) :: i
163 INTEGER(mpi) :: icount
164 INTEGER(mpi) :: mb
165 INTEGER(mpi) :: nbcol
166 !$ INTEGER(mpi) :: OMP_GET_MAX_THREADS
167 ! save input parameter
168 n=in
169 ireqpe=jreqpe
170 isngpe=jsngpe
171 icmprs=jcmprs+jextnd ! enforce compression for extended storage
172 iextnd=jextnd
173 ! number of precision types (D, F)
174 ispc=1
175 if (jsngpe>0) ispc=2
176 nspc = ispc
177 ! bit field size
178 icount=max(jsngpe+1,jhispe)
179 icount=max(jreqpe,icount)
180 ibfw=1 ! number of bits needed to count up to ICOUNT
181 mxcnt=1
182 DO i=1,30
183 IF (icount > mxcnt) THEN
184 ibfw=ibfw+1
185 mxcnt=mxcnt*2+1
186 END IF
187 END DO
188 ! bit field array size
189 noffd=int(n,mpl)*int(n-1,mpl)*int(ibfw,mpl)/2
190 ndimb=noffd/bs+n
191 idimb=ndimb
192 mb=int(4.0e-6*real(ndimb,mps),mpi)
193 WRITE(*,*) ' '
194 WRITE(*,*) 'CLBITS: symmetric matrix of dimension',n
195 WRITE(*,*) 'CLBITS: off-diagonal elements',noffd
196 IF (mb > 0) THEN
197 WRITE(*,*) 'CLBITS: dimension of bit-array',ndimb , '(',mb,'MB)'
198 ELSE
199 WRITE(*,*) 'CLBITS: dimension of bit-array',ndimb , '(< 1 MB)'
200 END IF
201 CALL mpalloc(bitfieldcounters,ndimb,'INBITS: bit storage')
202 bitfieldcounters=0
203 ! encoding for compression
204 nbcol=bs/2 ! one half of the bits for column number, other for column counter
205 DO i=bs/2,bs-2
206 IF (btest(n,i)) nbcol=i+1 ! more bits for column number
207 END DO
208 nencdb=bs-nbcol
209 iencdb=nencdb
210 nencdm=ishft(1,nencdb)-1
211 nthrd=1
212 !$ NTHRD=OMP_GET_MAX_THREADS()
213 RETURN
214END SUBROUTINE clbits
215
224SUBROUTINE ndbits(ndims,ncmprs,nsparr,ihst)
225 USE mpbits
226
227 INTEGER(mpl), DIMENSION(4), INTENT(OUT) :: ndims
228 INTEGER(mpi), DIMENSION(:), INTENT(OUT) :: ncmprs
229 INTEGER(mpl), DIMENSION(:,:), INTENT(OUT) :: nsparr
230 INTEGER(mpi), INTENT(IN) :: ihst
231
232 INTEGER(mpi) :: nwcp(0:1)
233 INTEGER(mpi) :: irgn(2)
234 INTEGER(mpi) :: inr(2)
235 INTEGER(mpi) :: ichunk
236 INTEGER(mpi) :: i
237 INTEGER(mpi) :: j
238 INTEGER(mpi) :: m
239 INTEGER(mpi) :: last
240 INTEGER(mpi) :: lrgn
241 INTEGER(mpi) :: next
242 INTEGER(mpi) :: icp
243 INTEGER(mpi) :: mm
244 INTEGER(mpi) :: jp
245 INTEGER(mpi) :: nj
246 INTEGER(mpi) :: ib
247 INTEGER(mpi) :: ir
248 INTEGER(mpi) :: icount
249 INTEGER(mpi) :: iproc
250 INTEGER(mpi) :: iword
251 INTEGER(mpi) :: k
252 INTEGER(mpi) :: mb
253 INTEGER(mpi) :: n1
254 INTEGER(mpl) :: ll
255 INTEGER(mpl) :: lb
256 INTEGER(mpl) :: nin
257 INTEGER(mpl) :: ntot
258 INTEGER(mpl) :: noffi
259 REAL(mps) :: cpr
260 REAL(mps) :: fracu
261 REAL(mps) :: fracz
262 LOGICAL :: btest
263 !$ INTEGER(mpi) :: OMP_GET_THREAD_NUM
264
265 ndims(1)=ndimb
266 ndims(2)=0
267 ndims(3)=0
268 ndims(4)=0
269 ntot=0
270 ll=0
271 lb=0
272 ichunk=min((n+nthrd-1)/nthrd/32+1,256)
273 IF (ibfw > 1.OR.icmprs > 0) THEN
274 ! reduce bit field counters to (precision type) bits, analyze precision type bit fields ('1st half' (j<i))
275
276 ! parallelize row loop
277 ! private copy of NTOT for each thread, combined at end, init with 0.
278 !$OMP PARALLEL DO &
279 !$OMP PRIVATE(I,NOFFI,LL,MM,LB,MB,IWORD,IPROC,J,ICOUNT,IB,INR,IRGN,LAST,LRGN,NEXT,JP,IR) &
280 !$OMP REDUCTION(+:NTOT) &
281 !$OMP SCHEDULE(DYNAMIC,ICHUNK)
282 DO i=1,n
283 noffi=int(i-1,mpl)*int(i-2,mpl)*int(ibfw,mpl)/2
284 ll=noffi/bs+i
285 mm=0
286 lb=ll
287 mb=0
288 iword=0 ! temporary bit fields
289 iproc=0
290 !$ IPROC=OMP_GET_THREAD_NUM() ! thread number
291 inr(1)=0
292 inr(2)=0
293 irgn(1)=0
294 irgn(2)=0
295 last=0
296 lrgn=0
297
298 DO j=1,i-1
299 ! get (pair) counter
300 icount=0
301 next=0
302 DO ib=0,ibfw-1
303 IF (btest(bitfieldcounters(ll),mm)) icount=ibset(icount,ib)
304 mm=mm+1
305 IF (mm >= bs) THEN
306 ll=ll+1
307 mm=mm-bs
308 END IF
309 END DO
310
311 IF (icount > 0) THEN
312 ntot=ntot+1
313 IF (iproc == 0.AND.ihst > 0) CALL hmpent(ihst,real(icount,mps))
314 END IF
315
316 ! keep pair ?
317 IF (icount >= ireqpe) THEN
318 next=1 ! double
319 IF (icount <= isngpe) next=2 ! single
320 iword=ibset(iword,mb+next-1)
321 inr(next)=inr(next)+1
322 IF (next /= last.OR.lrgn >= nencdm) THEN
323 irgn(next)=irgn(next)+1
324 lrgn=0
325 END IF
326 lrgn=lrgn+1
327 END IF
328 last=next
329
330 mb=mb+nspc
331 IF (mb >= bs) THEN
332 bitfieldcounters(lb)=iword ! store
333 iword=0
334 lb=lb+1
335 mb=mb-bs
336 END IF
337 END DO
338 bitfieldcounters(lb)=iword ! store
339
340 ! save row statistics
341 ir=i+1
342 DO jp=1,nspc
343 nsparr(1,ir)=irgn(jp) ! number of regions per row and precision
344 nsparr(2,ir)=inr(jp) ! number of columns per row and precision
345 ir=ir+n+1
346 END DO
347 END DO
348
349 ! analyze precision type bit fields for extended storage, check for row compression
350
351 ! parallelize row loop
352 ! private copy of NDIMS for each thread, combined at end, init with 0.
353 !$OMP PARALLEL DO &
354 !$OMP PRIVATE(I,NOFFI,NOFFJ,LL,MM,INR,IRGN,LAST,LRGN,J,NEXT,ICP,NWCP,JP,IR,IB) &
355 !$OMP REDUCTION(+:NDIMS) &
356 !$OMP SCHEDULE(DYNAMIC,ICHUNK)
357 DO i=1,n
358 ! restore row statistics
359 irgn(1)=int(nsparr(1,i+1),mpi)
360 irgn(2)=int(nsparr(1,i+n+2),mpi)
361 inr(1)=int(nsparr(2,i+1),mpi)
362 inr(2)=int(nsparr(2,i+n+2),mpi)
363
364 ! analyze precision type bit fields for extended storage ('2nd half' (j>i) too) ?
365 IF (iextnd > 0) THEN
366
367 noffj=(i-1)*nspc
368 mm=mod(noffj,bs)
369
370 last=0
371 lrgn=0
372
373 ! remaining columns
374 DO j=i+1, n
375 ! index for pair (J,I)
376 noffi=int(j-1,mpl)*int(j-2,mpl)*int(ibfw,mpl)/2 ! for I=1
377 ll=noffi/bs+j+noffj/bs ! row offset + column offset
378
379 ! get precision type
380 next=0
381 DO ib=0,nspc-1
382 IF (btest(bitfieldcounters(ll),mm+ib)) next=ibset(next,ib)
383 END DO
384
385 ! keep pair ?
386 IF (next > 0) THEN
387 inr(next)=inr(next)+1
388 IF (next /= last.OR.lrgn >= nencdm) THEN
389 irgn(next)=irgn(next)+1
390 lrgn=0
391 END IF
392 lrgn=lrgn+1
393 END IF
394 last=next
395 END DO
396 END IF
397
398 ! row statistics, compression
399 ir=i+1
400 DO jp=1,nspc
401 icp=0
402 nwcp(0)=inr(jp) ! list of column indices (default)
403 IF (inr(jp) > 0) THEN
404 nwcp(1)=irgn(jp)+(irgn(jp)+7)/8 ! list of regions of consecutive columns (and group offsets)
405 ! compress row ?
406 IF ((nwcp(1) < nwcp(0).AND.icmprs > 0).OR.iextnd > 0) THEN
407 icp=1
408 ncmprs(i+n*(jp-1))=irgn(jp) ! number of regions per row and precision
409 END IF
410 ! total space
411 ndims(2) =ndims(2) +nwcp(icp)
412 ndims(jp+2)=ndims(jp+2)+nwcp(0)
413 END IF
414 ! per row and precision
415 nsparr(1,ir)=nwcp(icp)
416 nsparr(2,ir)=nwcp(0)
417 ir=ir+n+1
418 END DO
419 END DO
420 !$OMP END PARALLEL DO
421
422 ! sum up, fill row offsets
423 lb=1
424 n1=0
425 ll=n+1
426 DO jp=1,nspc
427 DO i=1,n
428 n1=n1+1
429 nsparr(1,n1)=lb
430 nsparr(2,n1)=ll
431 lb=lb+nsparr(1,n1+1)
432 ll=ll+nsparr(2,n1+1)
433 END DO
434 n1=n1+1
435 nsparr(1,n1)=lb
436 nsparr(2,n1)=ll
437 ll=1
438 END DO
439
440 ELSE
441
442 nin=0
443 nsparr(1,1)=1
444 nsparr(2,1)=n+1
445 n1=1
446 DO i=1,n
447 noffi=int(i-1,mpl)*int(i-2,mpl)/2
448 ll=noffi/bs+i
449 nj=(i-1)/bs
450 DO k=0,nj
451 DO m=0,bs-1
452 IF(btest(bitfieldcounters(ll+k),m)) nin=nin+1
453 END DO
454 END DO
455 n1=n1+1
456 nsparr(1,n1)=nsparr(1,1)+nin
457 nsparr(2,n1)=nsparr(2,1)+nin
458 END DO
459 ndims(2)=nin
460 ndims(3)=nin
461 ntot=nin
462
463 END IF
464
465 nin=ndims(3)+ndims(4)
466 fracz=200.0*real(ntot,mps)/real(n,mps)/real(n-1,mps)
467 fracu=200.0*real(nin,mps)/real(n,mps)/real(n-1,mps)
468 WRITE(*,*) ' '
469 WRITE(*,*) 'NDBITS: number of diagonal elements',n
470 WRITE(*,*) 'NDBITS: number of used off-diagonal elements',nin
471 WRITE(*,1000) 'fraction of non-zero off-diagonal elements', fracz
472 WRITE(*,1000) 'fraction of used off-diagonal elements', fracu
473 IF (icmprs /= 0) THEN
474 cpr=100.0*real(mpi*ndims(2)+mpd*ndims(3)+mps*ndims(4),mps)/real((mpd+mpi)*nin,mps)
475 WRITE(*,1000) 'compression ratio for off-diagonal elements', cpr
476 END IF
4771000 FORMAT(' NDBITS: ',a,f6.2,' %')
478 RETURN
479END SUBROUTINE ndbits
480
486SUBROUTINE ckbits(ndims)
487 USE mpbits
488
489 INTEGER(mpl), DIMENSION(4), INTENT(OUT) :: ndims
490
491 INTEGER(mpi) :: nwcp(0:1)
492 INTEGER(mpi) :: irgn(2)
493 INTEGER(mpi) :: inr(2)
494 INTEGER(mpl) :: ll
495 INTEGER(mpl) :: noffi
496 INTEGER(mpi) :: i
497 INTEGER(mpi) :: j
498 INTEGER(mpi) :: last
499 INTEGER(mpi) :: lrgn
500 INTEGER(mpi) :: next
501 INTEGER(mpi) :: icp
502 INTEGER(mpi) :: ib
503 INTEGER(mpi) :: icount
504 INTEGER(mpi) :: kbfw
505 INTEGER(mpi) :: jp
506 INTEGER(mpi) :: mm
507 LOGICAL :: btest
508
509 DO i=1,4
510 ndims(i)=0
511 END DO
512 kbfw=1
513 IF (ibfw > 1.AND.icmprs > 0) kbfw=2
514 ll=0
515
516 DO i=1,n
517 noffi=int(i-1,mpl)*int(i-2,mpl)*int(ibfw,mpl)/2
518 ll=noffi/bs+i
519 mm=0
520 inr(1)=0
521 inr(2)=0
522 irgn(1)=0
523 irgn(2)=0
524 last=0
525 lrgn=0
526 DO j=1,i-1
527 icount=0
528 next=0
529 DO ib=0,ibfw-1
530 IF (btest(bitfieldcounters(ll),mm)) icount=ibset(icount,ib)
531 mm=mm+1
532 IF (mm >= bs) THEN
533 ll=ll+1
534 mm=mm-bs
535 END IF
536 END DO
537
538 IF (icount > 0) ndims(1)=ndims(1)+1
539 ! keep pair ?
540 IF (icount >= ireqpe) THEN
541 next=1 ! double
542 IF (icount <= isngpe) next=2 ! single
543 inr(next)=inr(next)+1
544 IF (next /= last.OR.lrgn >= nencdm) THEN
545 irgn(next)=irgn(next)+1
546 lrgn=0
547 END IF
548 lrgn=lrgn+1
549 END IF
550 last=next
551 END DO
552
553 IF (icmprs > 0) THEN
554 DO jp=1,kbfw
555 IF (inr(jp) > 0) THEN
556 icp=0
557 nwcp(0)=inr(jp) ! list of column indices (default)
558 nwcp(1)=irgn(jp)+(irgn(jp)+7)/8 ! list of regions of consecutive columns
559 ! compress row ?
560 IF (nwcp(1) < nwcp(0).OR. iextnd > 0) icp=1
561 ndims(2) =ndims(2) +nwcp(icp)
562 ndims(jp+2)=ndims(jp+2)+nwcp(0)
563 END IF
564 END DO
565 ELSE
566 ndims(2)=ndims(2)+inr(1)
567 ndims(3)=ndims(3)+inr(1)
568 END IF
569
570 END DO
571
572 RETURN
573END SUBROUTINE ckbits
574
581SUBROUTINE spbits(nsparr,nsparc,ncmprs) ! collect elements
582 USE mpbits
583 USE mpdalc
584
585 INTEGER(mpl), DIMENSION(:,:), INTENT(IN) :: nsparr
586 INTEGER(mpi), DIMENSION(:), INTENT(OUT) :: nsparc
587 INTEGER(mpi), DIMENSION(:), INTENT(INOUT) :: ncmprs
588
589 INTEGER(mpl) :: kl
590 INTEGER(mpl) :: l
591 INTEGER(mpl) :: ll
592 INTEGER(mpl) :: l1
593 INTEGER(mpl) :: k8
594 INTEGER(mpl) :: n1
595 INTEGER(mpl) :: noffi
596 INTEGER(mpi) :: i
597 INTEGER(mpi) :: j
598 INTEGER(mpi) :: j1
599 INTEGER(mpi) :: jb
600 INTEGER(mpi) :: jn
601 INTEGER(mpi) :: m
602 INTEGER(mpi) :: ichunk
603 INTEGER(mpi) :: next
604 INTEGER(mpi) :: last
605 INTEGER(mpi) :: lrgn
606 INTEGER(mpi) :: nrgn
607 INTEGER(mpi) :: nrgn8
608 LOGICAL :: btest
609
610 ichunk=min((n+nthrd-1)/nthrd/32+1,256)
611
612 DO jb=0,nspc-1
613 ! parallelize row loop
614 !$OMP PARALLEL DO &
615 !$OMP PRIVATE(I,N1,NOFFI,NOFFJ,L,M,KL,L1,NRGN,NRGN8,K8,LAST,LRGN,LL,J1,JN,J,NEXT) &
616 !$OMP SCHEDULE(DYNAMIC,ICHUNK)
617 DO i=1,n
618 n1=i+jb*(n+1)
619 noffi=int(i-1,mpl)*int(i-2,mpl)*int(ibfw,mpl)/2
620 l=noffi/bs+i
621 m=jb
622 kl=nsparr(1,n1)-1 ! pointer to row in NSPARC
623 l1=nsparr(2,n1) ! pointer to row in sparse matrix
624 nrgn=ncmprs(i+n*jb)! compression (number of consecutive regions)
625 nrgn8=(nrgn+7)/8 ! number of groups (1 offset per group)
626 k8=kl
627 kl=kl+nrgn8 ! reserve space of offsets
628 last=0
629 lrgn=0
630 ll=l1-1
631 j1=0
632 jn=0
633
634 DO j=1,i-1 ! loop for off-diagonal elements
635 next=0
636 IF(bitfieldcounters(l) /= 0) THEN
637 IF(btest(bitfieldcounters(l),m)) THEN
638 ll=ll+1
639 IF (nrgn <= 0) THEN
640 kl=kl+1
641 nsparc(kl)=j ! column index
642 ELSE
643 next=1
644 IF (last == 0.OR.jn >= nencdm) THEN
645 IF (mod(lrgn,8) == 0) THEN
646 k8=k8+1
647 nsparc(k8)=int(ll-l1,mpi)
648 END IF
649 lrgn=lrgn+1
650 kl=kl+1
651 j1=ishft(j,nencdb)
652 jn=0
653 END IF
654 jn=jn+1
655 nsparc(kl)=ior(j1,jn)
656 END IF
657 END IF
658 END IF
659 last=next
660 m=m+nspc
661 IF (m >= bs) THEN
662 m=m-bs
663 l=l+1
664 END IF
665 END DO
666
667 ! extended storage ('2nd half' too) ?
668 IF (iextnd > 0) THEN
669 noffj=(i-1)*nspc
670 m=mod(noffj,bs)+jb
671 last=0
672 ncmprs(i+n*jb)=lrgn ! remember number of regions in 1st half (j<i)
673 ! remaining columns
674 DO j=i+1, n
675 ! index for pair (J,I)
676 noffi=int(j-1,mpl)*int(j-2,mpl)*int(ibfw,mpl)/2 ! for I=1
677 l=noffi/bs+j+noffj/bs ! row offset + column offset
678 next=0
679 IF(btest(bitfieldcounters(l),m)) THEN
680 ll=ll+1
681 IF (nrgn <= 0) THEN
682 kl=kl+1
683 nsparc(kl)=j ! column index
684 ELSE
685 next=1
686 IF (last == 0.OR.jn >= nencdm) THEN
687 IF (mod(lrgn,8) == 0) THEN
688 k8=k8+1
689 nsparc(k8)=int(ll-l1,mpi)
690 END IF
691 lrgn=lrgn+1
692 kl=kl+1
693 j1=ishft(j,nencdb)
694 jn=0
695 END IF
696 jn=jn+1
697 nsparc(kl)=ior(j1,jn)
698 END IF
699 END IF
700 last=next
701
702 END DO
703 END IF
704
705 END DO
706 !$OMP END PARALLEL DO
707
708 END DO
709
710 n1=(n+1)*ibfw
711 WRITE(*,*) ' '
712 WRITE(*,*) 'SPBITS: sparse structure constructed ',nsparr(1,n1), ' words'
713 WRITE(*,*) 'SPBITS: dimension parameter of matrix',nsparr(2,1)-1
714 IF (ibfw <= 1) THEN
715 WRITE(*,*) 'SPBITS: index of last used location',nsparr(2,n1)-1
716 ELSE
717 WRITE(*,*) 'SPBITS: index of last used double',nsparr(2,n1/2)-1
718 WRITE(*,*) 'SPBITS: index of last used single',nsparr(2,n1)-1
719 END IF
720 CALL mpdealloc(bitfieldcounters)
721 RETURN
722END SUBROUTINE spbits
723
724
728!
729SUBROUTINE clbmap(in)
730 USE mpbits
731 USE mpdalc
732
733 INTEGER(mpi), INTENT(IN) :: in
734
735 INTEGER(mpl) :: noffd
736 INTEGER(mpi) :: mb
737
738 ! save input parameter
739 n2=in
740 ! bit field array size
741 noffd=int(n2,mpl)*int(n2-1,mpl)/2
742 ndimb2=noffd/bs+n2
743 mb=int(4.0e-6*real(ndimb2,mps),mpi)
744 WRITE(*,*) ' '
745 IF (mb > 0) THEN
746 WRITE(*,*) 'CLBMAP: dimension of bit-map',ndimb2 , '(',mb,'MB)'
747 ELSE
748 WRITE(*,*) 'CLBMAP: dimension of bit-map',ndimb2 , '(< 1 MB)'
749 END IF
750 CALL mpalloc(bitmap,ndimb2,'INBMAP: bit storage')
751 bitmap=0
752 RETURN
753END SUBROUTINE clbmap
754
760SUBROUTINE inbmap(im,jm) ! include element (I,J)
761 USE mpbits
762
763 INTEGER(mpi), INTENT(IN) :: im
764 INTEGER(mpi), INTENT(IN) :: jm
765
766 INTEGER(mpl) :: l
767 INTEGER(mpi) :: i
768 INTEGER(mpi) :: j
769 INTEGER(mpi) :: noffj
770 INTEGER(mpl) :: noffi
771 INTEGER(mpi) :: m
772
773 IF(im == jm) RETURN ! diagonal
774 j=min(im,jm)
775 i=max(im,jm)
776 IF(j <= 0) RETURN ! out low
777 IF(i > n2) RETURN ! out high
778 noffi=int(i-1,mpl)*int(i-2,mpl)/2 ! for J=1
779 noffj=(j-1)
780 l=noffi/bs+i+noffj/bs ! row offset + column offset
781 ! add I instead of 1 to keep bit maps of different rows in different words (openMP !)
782 m=mod(noffj,bs)
783 bitmap(l)=ibset(bitmap(l),m)
784 RETURN
785END SUBROUTINE inbmap
786
791SUBROUTINE gpbmap(npair)
792 USE mpbits
793
794 INTEGER(mpi), DIMENSION(:), INTENT(OUT) :: npair
795
796 INTEGER(mpl) :: l
797 INTEGER(mpl) :: noffi
798 INTEGER(mpi) :: i
799 INTEGER(mpi) :: j
800 INTEGER(mpi) :: m
801 LOGICAL :: btest
802
803 npair(1:n2)=0
804 l=0
805
806 DO i=1,n2
807 noffi=int(i-1,mpl)*int(i-2,mpl)/2
808 l=noffi/bs+i
809 m=0
810 DO j=1,i-1
811 IF (btest(bitmap(l),m)) THEN
812 npair(i)=npair(i)+1
813 npair(j)=npair(j)+1
814 END IF
815 m=m+1
816 IF (m >= bs) THEN
817 l=l+1
818 m=m-bs
819 END IF
820 END DO
821 END DO
822
823 RETURN
824END SUBROUTINE gpbmap
mpdalc::mpalloc
allocate array
Definition mpdalc.f90:36
mpdalc::mpdealloc
deallocate array
Definition mpdalc.f90:41
clbits
subroutine clbits(in, jreqpe, jhispe, jsngpe, jcmprs, jextnd, idimb, iencdb, ispc)
Calculate bit (field) array size, encoding.
Definition mpbits.f90:148
gpbmap
subroutine gpbmap(npair)
Get pairs (statistic) from map.
Definition mpbits.f90:792
ckbits
subroutine ckbits(ndims)
Check sparsity of matrix.
Definition mpbits.f90:487
spbits
subroutine spbits(nsparr, nsparc, ncmprs)
Create sparsity information.
Definition mpbits.f90:582
clbmap
subroutine clbmap(in)
Clear (additional) bit map.
Definition mpbits.f90:730
inbmap
subroutine inbmap(im, jm)
Fill bit map.
Definition mpbits.f90:761
inbits
subroutine inbits(im, jm, inc)
Fill bit fields (counters).
Definition mpbits.f90:66
ndbits
subroutine ndbits(ndims, ncmprs, nsparr, ihst)
Analyze bit fields.
Definition mpbits.f90:225
mpbits
Bit field data.
Definition mpbits.f90:35
mpbits::n2
integer(mpi) n2
matrix size (map)
Definition mpbits.f90:42
mpbits::ireqpe
integer(mpi) ireqpe
min number of pair entries
Definition mpbits.f90:44
mpbits::ndimb2
integer(mpl) ndimb2
dimension for bit map
Definition mpbits.f90:40
mpbits::ibfw
integer(mpi) ibfw
bit field width
Definition mpbits.f90:43
mpbits::iextnd
integer(mpi) iextnd
flag for extended storage (both 'halves' of sym. mat. for improved access patterns)
Definition mpbits.f90:47
mpbits::n
integer(mpi) n
matrix size (counters)
Definition mpbits.f90:41
mpbits::bs
integer(mpi), parameter bs
number of bits in INTEGER(mpi)
Definition mpbits.f90:55
mpbits::nspc
integer(mpi) nspc
number of precision for sparse global matrix (1=D, 2=D+f)
Definition mpbits.f90:48
mpbits::nencdb
integer(mpi) nencdb
number of bits for encoding column counter
Definition mpbits.f90:51
mpbits::isngpe
integer(mpi) isngpe
upper bound for pair entry single precision storage
Definition mpbits.f90:45
mpbits::bitmap
integer(mpi), dimension(:), allocatable bitmap
fit field map for global parameters pairs (measurements)
Definition mpbits.f90:54
mpbits::mxcnt
integer(mpi) mxcnt
max value for bit field counters
Definition mpbits.f90:49
mpbits::icmprs
integer(mpi) icmprs
compression flag for sparsity (column indices)
Definition mpbits.f90:46
mpbits::nthrd
integer(mpi) nthrd
number of threads
Definition mpbits.f90:52
mpbits::ndimb
integer(mpl) ndimb
dimension for bit (field) array
Definition mpbits.f90:39
mpbits::bitfieldcounters
integer(mpi), dimension(:), allocatable bitfieldcounters
fit field counters for global parameters pairs (tracks)
Definition mpbits.f90:53
mpbits::nencdm
integer(mpi) nencdm
max value for column counter
Definition mpbits.f90:50
mpdalc
(De)Allocate vectors and arrays.
Definition mpdalc.f90:24
mpdef
Definition of constants.
Definition mpdef.f90:24