SND@LHC Software
Loading...
Searching...
No Matches
Modules | Functions/Subroutines | Variables
mpbits.f90 File Reference

Bit field counters. More...

Go to the source code of this file.

Modules

module  mpbits
 Bit field data.
 

Functions/Subroutines

subroutine inbits (im, jm, inc)
 Fill bit fields (counters).
 
subroutine clbits (in, jreqpe, jhispe, jsngpe, jcmprs, jextnd, idimb, iencdb, ispc)
 Calculate bit (field) array size, encoding.
 
subroutine ndbits (ndims, ncmprs, nsparr, ihst)
 Analyze bit fields.
 
subroutine ckbits (ndims)
 Check sparsity of matrix.
 
subroutine spbits (nsparr, nsparc, ncmprs)
 Create sparsity information.
 
subroutine clbmap (in)
 Clear (additional) bit map.
 
subroutine inbmap (im, jm)
 Fill bit map.
 
subroutine gpbmap (npair)
 Get pairs (statistic) from map.
 

Variables

integer(mpl) mpbits::ndimb
 dimension for bit (field) array
 
integer(mpl) mpbits::ndimb2
 dimension for bit map
 
integer(mpi) mpbits::n
 matrix size (counters)
 
integer(mpi) mpbits::n2
 matrix size (map)
 
integer(mpi) mpbits::ibfw
 bit field width
 
integer(mpi) mpbits::ireqpe
 min number of pair entries
 
integer(mpi) mpbits::isngpe
 upper bound for pair entry single precision storage
 
integer(mpi) mpbits::icmprs
 compression flag for sparsity (column indices)
 
integer(mpi) mpbits::iextnd
 flag for extended storage (both 'halves' of sym. mat. for improved access patterns)
 
integer(mpi) mpbits::nspc
 number of precision for sparse global matrix (1=D, 2=D+f)
 
integer(mpi) mpbits::mxcnt
 max value for bit field counters
 
integer(mpi) mpbits::nencdm
 max value for column counter
 
integer(mpi) mpbits::nencdb
 number of bits for encoding column counter
 
integer(mpi) mpbits::nthrd
 number of threads
 
integer(mpi), dimension(:), allocatable mpbits::bitfieldcounters
 fit field counters for global parameters pairs (tracks)
 
integer(mpi), dimension(:), allocatable mpbits::bitmap
 fit field map for global parameters pairs (measurements)
 
integer(mpi), parameter mpbits::bs = BIT_SIZE(1_mpi)
 number of bits in INTEGER(mpi)
 

Detailed Description

Bit field counters.

Author
Volker Blobel, University Hamburg, 2005-2009 (initial Fortran77 version)
Claus Kleinwort, DESY (maintenance and developement)
Copyright
Copyright (c) 2009 - 2015 Deutsches Elektronen-Synchroton, Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY

This library is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more details.

You should have received a copy of the GNU Library General Public License along with this program (see the file COPYING.LIB for more details); if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

Count pairs of global parameters for sparse storage of global matrix, apply pair entries cut and build (compressed) sparsity structure (row offsets, column lists).

In sparse storage mode for each row the list of column indices with non zero elements (and those elements) are stored. With compression this list is represented by the first column and their number for continous regions (encoded in single INTEGER(mpi) words). Rare elements may be stored in single precision.

An additional bit map is used to monitor the parameter pairs for measurements (or 'equations').

Definition in file mpbits.f90.

Function/Subroutine Documentation

◆ ckbits()

subroutine ckbits ( integer(mpl), dimension(4), intent(out)  ndims)

Check sparsity of matrix.

Parameters
[out]ndims(1): (reduced) size of bit array; (2): size of column lists; (3/4): number of (double/single precision) off diagonal elements;

Definition at line 486 of file mpbits.f90.

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
mpbits
Bit field data.
Definition mpbits.f90:35
mpbits::ireqpe
integer(mpi) ireqpe
min number of pair entries
Definition mpbits.f90:44
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::isngpe
integer(mpi) isngpe
upper bound for pair entry single precision storage
Definition mpbits.f90:45
mpbits::icmprs
integer(mpi) icmprs
compression flag for sparsity (column indices)
Definition mpbits.f90:46
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

◆ clbits()

subroutine clbits ( integer(mpi), intent(in)  in,
integer(mpi), intent(in)  jreqpe,
integer(mpi), intent(in)  jhispe,
integer(mpi), intent(in)  jsngpe,
integer(mpi), intent(in)  jcmprs,
integer(mpi), intent(in)  jextnd,
integer(mpl), intent(out)  idimb,
integer(mpi), intent(out)  iencdb,
integer(mpi), intent(out)  ispc 
)

Calculate bit (field) array size, encoding.

Parameters
[in]inmatrix size
[in]jreqpemin number of pair entries
[in]jhispemupper bound for pair entry histogrammimg
[in]jsngpeupper bound for pair entry single precision storage
[in]jcmprscompression flag for sparsity (column indices)
[in]jextndflag for extended storage
[out]idimbdimension for bit (field) array
[out]iencdbnumber of bits for encoding column counter
[out]ispcnumber of precision for sparse global matrix

Definition at line 147 of file mpbits.f90.

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
mpdalc::mpalloc
allocate array
Definition mpdalc.f90:36
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::mxcnt
integer(mpi) mxcnt
max value for bit field counters
Definition mpbits.f90:49
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
mpdalc
(De)Allocate vectors and arrays.
Definition mpdalc.f90:24

◆ clbmap()

subroutine clbmap ( integer(mpi), intent(in)  in)

Clear (additional) bit map.

Parameters
[in]inmatrix size

Definition at line 729 of file mpbits.f90.

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
mpbits::n2
integer(mpi) n2
matrix size (map)
Definition mpbits.f90:42
mpbits::ndimb2
integer(mpl) ndimb2
dimension for bit map
Definition mpbits.f90:40
mpbits::bitmap
integer(mpi), dimension(:), allocatable bitmap
fit field map for global parameters pairs (measurements)
Definition mpbits.f90:54

◆ gpbmap()

subroutine gpbmap ( integer(mpi), dimension(:), intent(out)  npair)

Get pairs (statistic) from map.

Parameters
[out]npairnumber of paired parameters

Definition at line 791 of file mpbits.f90.

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

◆ inbits()

subroutine inbits ( integer(mpi), intent(in)  im,
integer(mpi), intent(in)  jm,
integer(mpi), intent(in)  inc 
)

Fill bit fields (counters).

Parameters
[in]imfirst index
[in]jmsecond index
[in]incincrement (usually 1)

Definition at line 65 of file mpbits.f90.

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

◆ inbmap()

subroutine inbmap ( integer(mpi), intent(in)  im,
integer(mpi), intent(in)  jm 
)

Fill bit map.

Parameters
[in]imfirst index
[in]jmsecond index

Definition at line 760 of file mpbits.f90.

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

◆ ndbits()

subroutine ndbits ( integer(mpl), dimension(4), intent(out)  ndims,
integer(mpi), dimension(:), intent(out)  ncmprs,
integer(mpl), dimension(:,:), intent(out)  nsparr,
integer(mpi), intent(in)  ihst 
)

Analyze bit fields.

Parameters
[out]ndims(1): (reduced) size of bit array; (2): size of column lists; (3/4): number of (double/single precision) off diagonal elements;
[out]ncmprscompression info (per row)
[out]nsparrrow offsets
[in]ihst>0: histogram number

Definition at line 224 of file mpbits.f90.

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

◆ spbits()

subroutine spbits ( integer(mpl), dimension(:,:), intent(in)  nsparr,
integer(mpi), dimension(:), intent(out)  nsparc,
integer(mpi), dimension(:), intent(inout)  ncmprs 
)

Create sparsity information.

Parameters
[in]nsparrrow offsets
[out]nsparccolumn indices
[in,out]ncmprscompression info (per row, in: number of all regions, out: number of regions in 1st half (for accessing 2nd half))

Definition at line 581 of file mpbits.f90.

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
mpdalc::mpdealloc
deallocate array
Definition mpdalc.f90:41