Commit 4f288098 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

Constraint groups implemented, documention (check input, trouble shooting) added

parent 6784325a
# ignore file from 'make'
*.o
*.mod
# ignore pede executable
/pede
# ignore editor backups
*.bck
Documentation is maintained at www.wiki.terascale.de/index.php/Millepede_II .
Documentation is maintained at
https://gitlab.desy.de/claus.kleinwort/millepede-ii/-/wikis/home .
The previous location (www.wiki.terascale.de/index.php/Millepede_II)
is not maintained anymore and may disappear.
......@@ -730,13 +730,15 @@ END SUBROUTINE inbmap
!> Get pairs (statistic) from map.
!!
!! \param [in] ngroup number of parameter groups
!! \param [in] npgrp parameter groups
!! \param [out] npair number of paired parameters
!!
SUBROUTINE gpbmap(npgrp,npair)
SUBROUTINE gpbmap(ngroup,npgrp,npair)
USE mpbits
IMPLICIT NONE
INTEGER(mpi), INTENT(IN) :: ngroup
INTEGER(mpi), DIMENSION(:,:), INTENT(IN) :: npgrp
INTEGER(mpi), DIMENSION(:), INTENT(OUT) :: npair
......@@ -747,10 +749,10 @@ SUBROUTINE gpbmap(npgrp,npair)
INTEGER(mpi) :: m
LOGICAL :: btest
npair(1:n2)=0
npair(1:ngroup)=0
l=0
DO i=1,n2
DO i=1,ngroup
npair(i)=npair(i)+npgrp(2,i)-1 ! from own group
noffi=INT(i-1,mpl)*INT(i-2,mpl)/2
l=noffi/bs+i
......@@ -771,3 +773,51 @@ SUBROUTINE gpbmap(npgrp,npair)
RETURN
END SUBROUTINE gpbmap
!> Get paired (parameter) groups from map.
!!
!! \param [in] ipgrp parameter group
!! \param [out] npair number of paired parameters
!! \param [in] npgrp paired parameter groups (for ipgrp)
!!
SUBROUTINE ggbmap(ipgrp,npair,npgrp)
USE mpbits
IMPLICIT NONE
INTEGER(mpi), INTENT(IN) :: ipgrp
INTEGER(mpi), INTENT(OUT) :: npair
INTEGER(mpi), DIMENSION(:), INTENT(OUT) :: npgrp
INTEGER(mpl) :: l
INTEGER(mpl) :: noffi
INTEGER(mpi) :: noffj
INTEGER(mpi) :: i
INTEGER(mpi) :: j
LOGICAL :: btest
npair=0
i=ipgrp
noffi=INT(i-1,mpl)*INT(i-2,mpl)/2 ! for J=1
l=noffi/bs+i! row offset
! add I instead of 1 to keep bit maps of different rows in different words (openMP !)
DO j=1,ipgrp-1
noffj=j-1
IF (btest(bitMap(l+noffj/bs),MOD(noffj,bs))) THEN
npair=npair+1
npgrp(npair)=j
END IF
END DO
noffj=ipgrp-1
DO i=ipgrp+1,n2
noffi=INT(i-1,mpl)*INT(i-2,mpl)/2 ! for J=1
l=noffi/bs+i ! row offset
IF (btest(bitMap(l+noffj/bs),MOD(noffj,bs))) THEN
npair=npair+1
npgrp(npair)=i
END IF
END DO
RETURN
END SUBROUTINE ggbmap
......@@ -130,7 +130,8 @@ MODULE mpmod
INTEGER(mpi) :: nvpgrp !< number of variable parameter groups
INTEGER(mpi) :: napgrp !< number of all parameter groups (variable + Lagrange mult.)
INTEGER(mpi) :: npblck !< number of (disjoint) parameter blocks (>1: block diagonal storage)
INTEGER(mpi) :: ncblck !< number of (disjoint) constraint blocks
INTEGER(mpi) :: ncgrp !< number of (disjoint) constraint groups
INTEGER(mpi) :: ncblck !< number of (non overlapping) constraint blocks
INTEGER(mpl) :: mszcon !< (integrated block) matrix size for constraint matrix
INTEGER(mpl) :: mszprd !< (integrated block) matrix size for (constraint) product matrix
INTEGER(mpi), DIMENSION(3) :: nprecond !< number of constraints (blocks), matrix size for preconditioner
......@@ -230,7 +231,11 @@ MODULE mpmod
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 :: matConsRanges !< parameter ranges for constraints
INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: matConsSort !< keys and index for sorting
INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: matConsGroups !< start of constraint groups, parameter range
INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: vecConsGroupCounts !< counter for constraint groups
INTEGER(mpi), DIMENSION(:,:), ALLOCATABLE :: matConsGroupStats !< statistics for constraint groups
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
......
......@@ -51,7 +51,7 @@ END MODULE mpqldec
!! \param [in] n number of rows (parameters)
!! \param [in] m number of columns (constraints)
!! \param [in] l number of disjoint blocks
!! \param [in] s size of constrints matrix (compressed or n*m)
!! \param [in] s size of constraints matrix (compressed or n*m)
!! \param [in] k flag for progress monitoring
!!
SUBROUTINE qlini(n,m,l,s,k)
......@@ -210,9 +210,9 @@ END SUBROUTINE qldec
!! \param [in] a block compressed Npar-by-Ncon matrix
!! \param [in] bpar 2-by-NparBlock+1 matrix (with parameter block definition)
!! \param [in] bcon 3-by-Ncon(Block)s+1 matrix (with constraint block definition)
!! \param [in] cons 3-by-Ncons+1 matrix (with constraint definition)
!! \param [in] rcon 4-by-Ncons matrix (with constraint ranges)
!!
SUBROUTINE qldecb(a,bpar,bcon,cons)
SUBROUTINE qldecb(a,bpar,bcon,rcon)
USE mpqldec
USE mpdalc
......@@ -246,7 +246,7 @@ SUBROUTINE qldecb(a,bpar,bcon,cons)
REAL(mpd), INTENT(IN) :: a(matsize)
INTEGER(mpi), INTENT(IN) :: bpar(2,nblock+1)
INTEGER(mpi), INTENT(IN) :: bcon(3,ncon+1)
INTEGER(mpi), INTENT(IN) :: cons(3,ncon+1)
INTEGER(mpi), INTENT(IN) :: rcon(4,ncon)
!$POMP INST BEGIN(qldecb)
! prepare
......@@ -264,10 +264,12 @@ SUBROUTINE qldecb(a,bpar,bcon,cons)
ifirst=bcon(2,ibcon)
ilast=bcon(3,ibcon)
DO i=bcon(1,ibcon),bcon(1,ibcon+1)-1
ioffRow(i+1)=ioffRow(i)+npb ! row offset
ioffPar(i)=bcon(2,ibcon)-1 ! parameter offset
irangeParNZ(1,i)=cons(1,i)
irangeParNZ(2,i)=cons(2,i)
! non-zero range: first, last parameter
irangeParNZ(1,i)=rcon(1,i)
irangeParNZ(2,i)=rcon(2,i)
! storage: parameter, row offset
ioffPar(i)=rcon(3,i)-1
ioffRow(i+1)=ioffRow(i)+rcon(4,i)-ioffPar(i)
END DO
iclast=iclast+ncb
END DO
......@@ -344,10 +346,12 @@ SUBROUTINE qldecb(a,bpar,bcon,cons)
vecN(ifirst:ilast)=vecN(ifirst:ilast)/nrm
END IF
! update L too
ioff3=INT(k1-1,mpl)*INT(ncon,mpl)
ioff3=INT(k1-2,mpl)*INT(ncon,mpl)
! transformation
DO i=k1,k
ioff2=ioffRow(i)
ioff3=ioff3+ncon
IF (irangeParNZ(1,k) > irangeParNZ(2,i)) CYCLE ! no overlap
ioff2=ioffRow(i)+ioffPar(k)-ioffPar(i)
sp=dot_product(vecN(ifirst:ilast),matV(ioff2+1:ioff2+ncol))
IF (sp /= 0.0_mpd) THEN
! update matV
......@@ -359,7 +363,6 @@ SUBROUTINE qldecb(a,bpar,bcon,cons)
irangeParNZ(1,i)=min(irangeParNZ(1,i),irangeParNZ(1,k))
irangeParNZ(2,i)=max(irangeParNZ(2,i),irangeParNZ(2,k))
END IF
ioff3=ioff3+ncon
END DO
! store secondary diagonal
vecVk(icoff+k)=vecN(kn)
......@@ -553,10 +556,11 @@ END SUBROUTINE qlsmq
!! \param [in] aprod external procedure to calculate A*v
!! \param [in,out] A symmetric Npar-by-Npar matrix A in symmetric or unpacked storage mode
!! overwritten with Q*A*Q^t (t=false) or Q^t*A*Q (t=true)
!! \param [in] s size of A
!! \param [in] roff row offsets for A
!! \param [in] t use transposed of Q
!!
SUBROUTINE qlssq(aprod,A,roff,t)
SUBROUTINE qlssq(aprod,A,s,roff,t)
USE mpqldec
USE mpdalc
......@@ -582,7 +586,8 @@ SUBROUTINE qlssq(aprod,A,roff,t)
REAL(mpd) :: vtAv
REAL(mpd), DIMENSION(:), ALLOCATABLE :: Av
REAL(mpd), INTENT(IN OUT) :: A((INT(npar,mpl)*INT(npar,mpl)+INT(npar,mpl))/2)
REAL(mpd), INTENT(IN OUT) :: A(s)
INTEGER(mpl), INTENT(IN) :: s
INTEGER(mpl), INTENT(IN) :: roff(npar)
LOGICAL, INTENT(IN) :: t
......@@ -1014,6 +1019,8 @@ SUBROUTINE qldump()
INTEGER(mpi) :: kn
REAL(mpd) :: v1
REAL(mpd) :: v2
REAL(mpd) :: v3
REAL(mpd) :: v4
print *
ioff1=0
......@@ -1043,15 +1050,19 @@ SUBROUTINE qldump()
ioff1=ioff1+npar
DO j=1,ncon
IF (matL(ioff2+j) /= 0.0_mpd) THEN
IF (istat(6) == 0) istat(4)=j
v4=matL(ioff2+j)
IF (istat(6) == 0) THEN
istat(4)=j
v3=v4
END IF
istat(5)=j
istat(6)=istat(6)+1
END IF
END DO
ioff2=ioff2+ncon
print 100, i, istat, v1, v2, irangeParNZ(:,i)
print 100, i, istat, v1, v2, v3, v4, irangeParNZ(:,i)
END DO
print *
100 FORMAT(" qldump",7I8,2G13.5,2I8)
100 FORMAT(" qldump",7I8,4G13.5,2I8)
END SUBROUTINE qldump
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