Commit 9c40bb50 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

Split ordered constraints into disjoint blocks

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@166 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent c51622e7
......@@ -205,7 +205,6 @@ SUBROUTINE clbits(in,jreqpe,jhispe,jsngpe,jcmprs,jextnd,idimb,iencdb,ispc)
nencdm=ishft(1,nencdb)-1
nthrd=1
!$ NTHRD=OMP_GET_MAX_THREADS()
print *, ' <<< inbits ', idimb,iencdb,ispc
RETURN
END SUBROUTINE clbits
......
......@@ -119,6 +119,9 @@ MODULE mpmod
INTEGER(mpi) :: nfgb !< number of fit parameters
INTEGER(mpi) :: ncgb !< number of constraints
INTEGER(mpi) :: ncgbe !< number of empty constraints (no variable parameters)
INTEGER(mpi) :: ncblck !< number of (disjoint) constraint blocks
INTEGER(mpi) :: mszcon !< (integrated block) matrix size for constraint matrix
INTEGER(mpi) :: mszprd !< (integrated block) matrix size for (constraint) product matrix
INTEGER(mpi), DIMENSION(2) :: nprecond !< number of constraints, matrix size for preconditioner
INTEGER(mpi) :: nagbn !< max number of global paramters per record
INTEGER(mpi) :: nalcn !< max number of local paramters per record
......@@ -198,6 +201,10 @@ MODULE mpmod
REAL(mpd), DIMENSION(:), ALLOCATABLE :: matConsProduct !< product matrix of constraints
REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecConsResiduals !< residuals of constraints
REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecConsSolution !< solution for constraint elimination
! constraint sorting, blocks
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: vecConsStart !< start of constraint in listConstraints (unsorted input)
INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: matConsSort !< keys and index for sorting
INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: matConsBlocks !< start of constraint blocks, parameter range
! monitoring of input residuals
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: measIndex !< mapping of 1. global label to measurement index
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: measHists !< measurement histograms (100 bins per thread)
......
......@@ -67,6 +67,7 @@
!! HEAPF heap sort reals direct
!! SORT1K sort 1-dim key-array (CHK)
!! SORT2K sort 2-dim key-array
!! SORT2I sort 2-dim key-array with index (CHK)
!!
!----------------------------------------------------------------------
......@@ -1689,6 +1690,105 @@ SUBROUTINE sort2k(a,n)
GO TO 10
END SUBROUTINE sort2k
!> Quick sort 2 with index.
!!
!! Quick sort of A(3,N) integer.
!!
!! \param[in,out] a vector (pair) of integers, sorted at return and an index
!! \param[in] n size of vector
SUBROUTINE sort2i(a,n)
USE mpdef
IMPLICIT NONE
INTEGER(mpi) :: nlev ! stack size
PARAMETER (nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
INTEGER(mpi) :: i
INTEGER(mpi) ::j
INTEGER(mpi) ::l
INTEGER(mpi) ::r
INTEGER(mpi) ::lev
INTEGER(mpi) ::lr(nlev)
INTEGER(mpi) ::lrh
INTEGER(mpi) ::maxlev
INTEGER(mpi) ::a1 ! pivot key
INTEGER(mpi) ::a2 ! pivot key
INTEGER(mpi) ::at ! pivot key
INTEGER(mpi), INTENT(IN OUT) :: a(3,*)
INTEGER(mpi), INTENT(IN) :: n
! ...
maxlev=0
lev=0
l=1
r=n
10 IF(r-l == 1) THEN ! sort two elements L and R
IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r))) THEN
at=a(1,l) ! exchange L <-> R
a(1,l)=a(1,r)
a(1,r)=at
at=a(2,l)
a(2,l)=a(2,r)
a(2,r)=at
at=a(3,l)
a(3,l)=a(3,r)
a(3,r)=at
END IF
r=l
END IF
IF(r == l) THEN
IF(lev <= 0) THEN
WRITE(*,*) 'SORT2I (quicksort): maxlevel used/available =', maxlev,'/64'
RETURN
END IF
lev=lev-2
l=lr(lev+1)
r=lr(lev+2)
ELSE
! LRH=(L+R)/2
lrh=(l/2)+(r/2) ! avoid bit overflow
IF(MOD(l,2) == 1.AND.MOD(r,2) == 1) lrh=lrh+1
a1=a(1,lrh) ! middle
a2=a(2,lrh)
i=l-1 ! find limits [J,I] with [L,R]
j=r+1
20 i=i+1
IF(a(1,i) < a1) GO TO 20
IF(a(1,i) == a1.AND.a(2,i) < a2) GO TO 20
30 j=j-1
IF(a(1,j) > a1) GO TO 30
IF(a(1,j) == a1.AND.a(2,j) > a2) GO TO 30
IF(i <= j) THEN
at=a(1,i) ! exchange I <-> J
a(1,i)=a(1,j)
a(1,j)=at
at=a(2,i)
a(2,i)=a(2,j)
a(2,j)=at
at=a(3,i)
a(3,i)=a(3,j)
a(3,j)=at
GO TO 20
END IF
IF(lev+2 > nlev) THEN
CALL peend(33,'Aborted, stack overflow in quicksort')
STOP 'SORT2I (quicksort): stack overflow'
END IF
IF(r-i < j-l) THEN
lr(lev+1)=l
lr(lev+2)=j
l=i
ELSE
lr(lev+1)=i
lr(lev+2)=r
r=j
END IF
lev=lev+2
maxlev=MAX(maxlev,lev)
END IF
GO TO 10
END SUBROUTINE sort2i
!> Chi2/ndf cuts.
!!
!! Return limit in Chi^2/ndf for N sigmas (N=1, 2 or 3).
......
......@@ -32,6 +32,8 @@ MODULE mpqldec
INTEGER(mpi) :: ncon !< number of constraints
REAL(mpd), DIMENSION(:), ALLOCATABLE :: matV !< unit normals (v_i) of Householder reflectors
REAL(mpd), DIMENSION(:), ALLOCATABLE :: matL !< lower diagonal matrix L
REAL(mpd), DIMENSION(:), ALLOCATABLE :: vecN !< normal vector
END MODULE mpqldec
!> Initialize QL decomposition.
......@@ -56,6 +58,8 @@ SUBROUTINE qlini(n,m)
CALL mpalloc(matV,length,'QLDEC: V')
length=ncon*ncon
CALL mpalloc(matL,length,'QLDEC: L')
length=npar
CALL mpalloc(vecN,length,'QLDEC: v')
END SUBROUTINE qlini
! 141217 C. Kleinwort, DESY-FH1
......@@ -95,8 +99,6 @@ SUBROUTINE qldec(a)
REAL(mpd), INTENT(IN) :: a(*)
REAL(mpd) :: v(npar)
! prepare
length=npar*ncon
matV=a(1:length)
......@@ -108,35 +110,166 @@ SUBROUTINE qldec(a)
! column offset
ioff1=(k-1)*npar
! get column
v(1:kn)=matV(ioff1+1:ioff1+kn)
nrm = SQRT(dot_product(v(1:kn),v(1:kn)))
vecN(1:kn)=matV(ioff1+1:ioff1+kn)
nrm = SQRT(dot_product(vecN(1:kn),vecN(1:kn)))
IF (nrm == 0.0_mpd) CYCLE
!
IF (v(kn) >= 0.0_mpd) THEN
v(kn)=v(kn)+nrm
IF (vecN(kn) >= 0.0_mpd) THEN
vecN(kn)=vecN(kn)+nrm
ELSE
v(kn)=v(kn)-nrm
vecN(kn)=vecN(kn)-nrm
END IF
! create normal vector
nrm = SQRT(dot_product(v(1:kn),v(1:kn)))
v(1:kn)=v(1:kn)/nrm
nrm = SQRT(dot_product(vecN(1:kn),vecN(1:kn)))
vecN(1:kn)=vecN(1:kn)/nrm
! transformation
ioff2=0
DO i=1,k
sp=dot_product(v(1:kn),matV(ioff2+1:ioff2+kn))
matV(ioff2+1:ioff2+kn)=matV(ioff2+1:ioff2+kn)-2.0_mpd*v(1:kn)*sp
sp=dot_product(vecN(1:kn),matV(ioff2+1:ioff2+kn))
matV(ioff2+1:ioff2+kn)=matV(ioff2+1:ioff2+kn)-2.0_mpd*vecN(1:kn)*sp
ioff2=ioff2+npar
END DO
! store column of L
ioff3=(k-1)*ncon
matL(ioff3+k:ioff3+ncon)=matV(ioff1+kn:ioff1+npar)
! store normal vector
matV(ioff1+1:ioff1+kn)=v(1:kn)
matV(ioff1+1:ioff1+kn)=vecN(1:kn)
matV(ioff1+kn+1:ioff1+npar)=0.0_mpd
END DO
END SUBROUTINE qldec
! 190312 C. Kleinwort, DESY-BELLE
!> QL decomposition (for disjoint block matrix).
!!
!! QL decomposition with Householder transformations.
!! Decompose N-By-M matrix A into orthogonal N-by-N matrix Q and a
!! N-by-M matrix containing zeros except for a lower triangular
!! M-by-M matrix L (at the bottom):
!!
!! | 0 |
!! A = Q * | |
!! | L |
!!
!! The decomposition is stored in a N-by-M matrix matV containing the unit
!! normal vectors v_i of the hyperplanes (Householder reflectors) defining Q.
!! The lower triangular matrix L is stored in the M-by-M matrix matL.
!!
!! \param [in] a block compressed Npar-by-Ncon matrix
!! \param [in] nb number of blocks
!! \param [in] b 3-by-Ncon+1 matrix (with block definition)
!!
SUBROUTINE qldecb(a,nb,b)
USE mpqldec
USE mpdalc
! cost[dot ops] ~= Npar*Ncon*Ncon
IMPLICIT NONE
INTEGER(mpi) :: i
INTEGER(mpi) :: ib
INTEGER(mpi) :: ifirst
INTEGER(mpi) :: ilast
INTEGER(mpl) :: ioff1
INTEGER(mpl) :: ioff2
INTEGER(mpl) :: ioff3
INTEGER(mpi) :: k
INTEGER(mpi) :: k1
INTEGER(mpi) :: kn
INTEGER(mpi) :: ncb
INTEGER(mpi) :: npb
INTEGER(mpl) :: length
REAL(mpd) :: nrm
REAL(mpd) :: sp
REAL(mpd), INTENT(IN) :: a(*)
INTEGER(mpi), INTENT(IN) :: nb
INTEGER(mpi), INTENT(IN) :: b(3,*)
! prepare
length=npar*ncon
matV=0.0_mpd
matL=0.0_mpd
! expand a into matV
ioff1=0
ioff2=0
DO ib=1,nb
ncb=b(1,ib+1)-b(1,ib) ! number of constraints in block
npb=b(3,ib)+1-b(2,ib) ! number of parameters in block
ifirst=b(2,ib)
ilast=b(3,ib)
DO i=1,ncb
matV(ioff1+ifirst:ioff1+ilast)=a(ioff2+1:ioff2+npb)
ioff1=ioff1+npar
ioff2=ioff2+npb
END DO
END DO
ib=nb ! start with last block
k1=b(1,ib) ! first constraint in block
! Householder procedure
DO k=ncon,1,-1
kn=npar+k-ncon
! different block?
IF (k < k1) THEN
ib=ib-1
k1=b(1,ib)
END IF
! index if first non-zero element
ifirst=b(2,ib)
IF (ifirst > kn) CYCLE
! index of last element
ilast=min(b(3,ib),kn)
! column offsets
ioff1=(k-1)*npar
ioff2=(k1-1)*npar
! get column
vecN(kn)=0.0_mpd
vecN(ifirst:ilast)=matV(ioff1+ifirst:ioff1+ilast)
nrm = SQRT(dot_product(vecN(ifirst:ilast),vecN(ifirst:ilast)))
IF (nrm == 0.0_mpd) CYCLE
!
IF (vecN(kn) >= 0.0_mpd) THEN
vecN(kn)=vecN(kn)+nrm
ELSE
vecN(kn)=vecN(kn)-nrm
END IF
IF (ilast < kn) THEN
! create normal vector
nrm = SQRT(dot_product(vecN(ifirst:ilast),vecN(ifirst:ilast))+vecN(kn)*vecN(kn))
vecN(ifirst:ilast)=vecN(ifirst:ilast)/nrm
vecN(kn)=vecN(kn)/nrm
! transformation
DO i=k1,k
sp=dot_product(vecN(ifirst:ilast),matV(ioff2+ifirst:ioff2+ilast))
matV(ioff2+ifirst:ioff2+ilast)=matV(ioff2+ifirst:ioff2+ilast)-2.0_mpd*vecN(ifirst:ilast)*sp
matV(ioff2+kn)=-2.0_mpd*vecN(kn)*sp
ioff2=ioff2+npar
END DO
ELSE
! create normal vector
nrm = SQRT(dot_product(vecN(ifirst:ilast),vecN(ifirst:ilast)))
vecN(ifirst:ilast)=vecN(ifirst:ilast)/nrm
! transformation
DO i=k1,k
sp=dot_product(vecN(ifirst:ilast),matV(ioff2+ifirst:ioff2+ilast))
matV(ioff2+ifirst:ioff2+ilast)=matV(ioff2+ifirst:ioff2+ilast)-2.0_mpd*vecN(ifirst:ilast)*sp
ioff2=ioff2+npar
END DO
END IF
! store column of L
ioff3=(k-1)*ncon
matL(ioff3+k:ioff3+ncon)=matV(ioff1+kn:ioff1+npar)
! store normal vector
matV(ioff1+1:ioff1+npar)=0.0_mpd
matV(ioff1+ifirst:ioff1+ilast)=vecN(ifirst:ilast)
matV(ioff1+kn)=vecN(kn)
END DO
END SUBROUTINE qldecb
!> Multiply left by Q(t).
!!
......
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment