Commit e2837442 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

Upgraded indices for constraint matrices to 64bit

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@208 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent 9ee817fc
...@@ -131,8 +131,8 @@ MODULE mpmod ...@@ -131,8 +131,8 @@ MODULE mpmod
INTEGER(mpi) :: napgrp !< number of all parameter groups (variable + Lagrange mult.) 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) :: npblck !< number of (disjoint) parameter blocks (>1: block diagonal storage)
INTEGER(mpi) :: ncblck !< number of (disjoint) constraint blocks INTEGER(mpi) :: ncblck !< number of (disjoint) constraint blocks
INTEGER(mpi) :: mszcon !< (integrated block) matrix size for constraint matrix INTEGER(mpl) :: mszcon !< (integrated block) matrix size for constraint matrix
INTEGER(mpi) :: mszprd !< (integrated block) matrix size for (constraint) product matrix INTEGER(mpl) :: mszprd !< (integrated block) matrix size for (constraint) product matrix
INTEGER(mpi), DIMENSION(2) :: nprecond !< number of constraints, matrix size for preconditioner INTEGER(mpi), DIMENSION(2) :: nprecond !< number of constraints, matrix size for preconditioner
INTEGER(mpi) :: nagbn !< max number of global paramters per record INTEGER(mpi) :: nagbn !< max number of global paramters per record
INTEGER(mpi) :: nalcn !< max number of local paramters per record INTEGER(mpi) :: nalcn !< max number of local paramters per record
......
...@@ -69,10 +69,10 @@ SUBROUTINE qlini(n,m,l,k) ...@@ -69,10 +69,10 @@ SUBROUTINE qlini(n,m,l,k)
! allocate ! allocate
length=5*ncon length=5*ncon
CALL mpalloc(sparseV,length,'QLDEC: sparsity structure of V') CALL mpalloc(sparseV,length,'QLDEC: sparsity structure of V')
length=npar*ncon length=INT(npar,mpl)*INT(ncon,mpl)
CALL mpalloc(matV,length,'QLDEC: V') CALL mpalloc(matV,length,'QLDEC: V')
matV=0. matV=0.
length=ncon*ncon length=INT(ncon,mpl)*INT(ncon,mpl)
CALL mpalloc(matL,length,'QLDEC: L') CALL mpalloc(matL,length,'QLDEC: L')
matL=0. matL=0.
length=npar length=npar
...@@ -123,7 +123,7 @@ SUBROUTINE qldec(a) ...@@ -123,7 +123,7 @@ SUBROUTINE qldec(a)
REAL(mpd), INTENT(IN) :: a(*) REAL(mpd), INTENT(IN) :: a(*)
! prepare ! prepare
length=npar*ncon length=INT(npar,mpl)*INT(ncon,mpl)
matV=a(1:length) matV=a(1:length)
matL=0.0_mpd matL=0.0_mpd
! implemented as single block ! implemented as single block
...@@ -137,7 +137,7 @@ SUBROUTINE qldec(a) ...@@ -137,7 +137,7 @@ SUBROUTINE qldec(a)
IF(monpg>0) CALL monpgs(ncon+1-k) IF(monpg>0) CALL monpgs(ncon+1-k)
kn=npar+k-ncon kn=npar+k-ncon
! column offset ! column offset
ioff1=(k-1)*npar ioff1=INT(k-1,mpl)*INT(npar,mpl)
! get column ! get column
vecN(1:kn)=matV(ioff1+1:ioff1+kn) vecN(1:kn)=matV(ioff1+1:ioff1+kn)
nrm = SQRT(dot_product(vecN(1:kn),vecN(1:kn))) nrm = SQRT(dot_product(vecN(1:kn),vecN(1:kn)))
...@@ -159,7 +159,7 @@ SUBROUTINE qldec(a) ...@@ -159,7 +159,7 @@ SUBROUTINE qldec(a)
ioff2=ioff2+npar ioff2=ioff2+npar
END DO END DO
! store column of L ! store column of L
ioff3=(k-1)*ncon ioff3=INT(k-1,mpl)*INT(ncon,mpl)
matL(ioff3+k:ioff3+ncon)=matV(ioff1+kn:ioff1+npar) matL(ioff3+k:ioff3+ncon)=matV(ioff1+kn:ioff1+npar)
! store normal vector ! store normal vector
matV(ioff1+1:ioff1+kn)=vecN(1:kn) matV(ioff1+1:ioff1+kn)=vecN(1:kn)
...@@ -217,7 +217,6 @@ SUBROUTINE qldecb(a,bpar,bcon) ...@@ -217,7 +217,6 @@ SUBROUTINE qldecb(a,bpar,bcon)
INTEGER(mpi) :: kn INTEGER(mpi) :: kn
INTEGER(mpi) :: ncb INTEGER(mpi) :: ncb
INTEGER(mpi) :: npb INTEGER(mpi) :: npb
INTEGER(mpl) :: length
REAL(mpd) :: nrm REAL(mpd) :: nrm
REAL(mpd) :: sp REAL(mpd) :: sp
...@@ -228,7 +227,6 @@ SUBROUTINE qldecb(a,bpar,bcon) ...@@ -228,7 +227,6 @@ SUBROUTINE qldecb(a,bpar,bcon)
!$POMP INST BEGIN(qldecb) !$POMP INST BEGIN(qldecb)
! prepare ! prepare
length=npar*ncon
matV=0.0_mpd matV=0.0_mpd
matL=0.0_mpd matL=0.0_mpd
...@@ -275,8 +273,8 @@ SUBROUTINE qldecb(a,bpar,bcon) ...@@ -275,8 +273,8 @@ SUBROUTINE qldecb(a,bpar,bcon)
! index of last element ! index of last element
ilast=min(bcon(3,ibcon),kn) ilast=min(bcon(3,ibcon),kn)
! column offsets ! column offsets
ioff1=(k-1)*npar ioff1=INT(k-1,mpl)*INT(npar,mpl)
ioff2=(k1-1)*npar ioff2=INT(k1-1,mpl)*INT(npar,mpl)
! get column ! get column
vecN(kn)=0.0_mpd vecN(kn)=0.0_mpd
vecN(ifirst:ilast)=matV(ioff1+ifirst:ioff1+ilast) vecN(ifirst:ilast)=matV(ioff1+ifirst:ioff1+ilast)
...@@ -314,7 +312,7 @@ SUBROUTINE qldecb(a,bpar,bcon) ...@@ -314,7 +312,7 @@ SUBROUTINE qldecb(a,bpar,bcon)
END IF END IF
! store column of L ! store column of L
ioff3=(k-1)*ncon ioff3=INT(k-1,mpl)*INT(ncon,mpl)
matL(ioff3+k-icoff:ioff3+iclast-icoff)=matV(ioff1+kn:ioff1+iplast) matL(ioff3+k-icoff:ioff3+iclast-icoff)=matV(ioff1+kn:ioff1+iplast)
! store normal vector ! store normal vector
matV(ioff1+1:ioff1+npar)=0.0_mpd matV(ioff1+1:ioff1+npar)=0.0_mpd
...@@ -374,7 +372,7 @@ SUBROUTINE qlmlq(x,m,t) ...@@ -374,7 +372,7 @@ SUBROUTINE qlmlq(x,m,t)
IF (t) k=nconb+1-j IF (t) k=nconb+1-j
kn=nparb+k-nconb kn=nparb+k-nconb
! column offset ! column offset
ioff1=(k-1+icoff)*npar ioff1=INT(k-1+icoff,mpl)*INT(npar,mpl)
! transformation ! transformation
ioff2=0 ioff2=0
DO i=1,m DO i=1,m
...@@ -418,7 +416,7 @@ SUBROUTINE qlmrq(x,m,t) ...@@ -418,7 +416,7 @@ SUBROUTINE qlmrq(x,m,t)
IF (.not.t) k=ncon+1-j IF (.not.t) k=ncon+1-j
kn=npar+k-ncon kn=npar+k-ncon
! column offset ! column offset
ioff1=(k-1)*npar ioff1=INT(k-1,mpl)*INT(npar,mpl)
! transformation ! transformation
iend=m*kn iend=m*kn
DO i=1,npar DO i=1,npar
...@@ -462,9 +460,9 @@ SUBROUTINE qlsmq(x,t) ...@@ -462,9 +460,9 @@ SUBROUTINE qlsmq(x,t)
IF (t) k=ncon+1-j IF (t) k=ncon+1-j
kn=npar+k-ncon kn=npar+k-ncon
! column offset ! column offset
ioff1=(k-1)*npar ioff1=INT(k-1,mpl)*INT(npar,mpl)
! transformation ! transformation
iend=npar*kn iend=INT(npar,mpl)*INT(kn,mpl)
DO i=1,npar DO i=1,npar
sp=dot_product(matV(ioff1+1:ioff1+kn),x(i:iend:npar)) sp=dot_product(matV(ioff1+1:ioff1+kn),x(i:iend:npar))
x(i:iend:npar)=x(i:iend:npar)-2.0_mpd*matV(ioff1+1:ioff1+kn)*sp x(i:iend:npar)=x(i:iend:npar)-2.0_mpd*matV(ioff1+1:ioff1+kn)*sp
...@@ -544,7 +542,7 @@ SUBROUTINE qlssq(aprod,A,roff,t) ...@@ -544,7 +542,7 @@ SUBROUTINE qlssq(aprod,A,roff,t)
IF (t) k=nconb+1-j IF (t) k=nconb+1-j
kn=nparb+k-nconb kn=nparb+k-nconb
! column offset ! column offset
ioff1=(k-1+icoff)*npar ioff1=INT(k-1+icoff,mpl)*INT(npar,mpl)
! A*v ! A*v
CALL aprod(nparb,ioffp,matV(ioff1+1:ioff1+nparb),Av(1:nparb)) CALL aprod(nparb,ioffp,matV(ioff1+1:ioff1+nparb),Av(1:nparb))
! transformation ! transformation
...@@ -633,8 +631,7 @@ SUBROUTINE qlpssq(aprod,B,m,t) ...@@ -633,8 +631,7 @@ SUBROUTINE qlpssq(aprod,B,m,t)
END INTERFACE END INTERFACE
!$POMP INST BEGIN(qlpssq) !$POMP INST BEGIN(qlpssq)
length=npar length=INT(npar,mpl)*INT(ncon,mpl)
length=npar*ncon
CALL mpalloc(Av,length,'qlpssq: Av') CALL mpalloc(Av,length,'qlpssq: Av')
mbnd=max(0,m-1) ! band width without diagonal mbnd=max(0,m-1) ! band width without diagonal
...@@ -652,7 +649,7 @@ SUBROUTINE qlpssq(aprod,B,m,t) ...@@ -652,7 +649,7 @@ SUBROUTINE qlpssq(aprod,B,m,t)
IF (t) k=ncon+1-j IF (t) k=ncon+1-j
kn=npar+k-ncon kn=npar+k-ncon
! column offset ! column offset
ioff1=(k-1)*npar ioff1=INT(k-1,mpl)*INT(npar,mpl)
! redion offset ! redion offset
ioffr=(k-1)*5 ioffr=(k-1)*5
! transformation (diagonal block) ! transformation (diagonal block)
...@@ -680,7 +677,7 @@ SUBROUTINE qlpssq(aprod,B,m,t) ...@@ -680,7 +677,7 @@ SUBROUTINE qlpssq(aprod,B,m,t)
DO j2=j+1,ncon DO j2=j+1,ncon
k2=j2 k2=j2
IF (t) k2=ncon+1-j2 IF (t) k2=ncon+1-j2
ioff2=(k2-1)*npar ioff2=INT(k2-1,mpl)*INT(npar,mpl)
! loop over non-zero regions ! loop over non-zero regions
vtvp=0._mpd vtvp=0._mpd
vtAvp=0._mpd vtAvp=0._mpd
...@@ -737,7 +734,7 @@ SUBROUTINE qlgete(emin,emax) ...@@ -737,7 +734,7 @@ SUBROUTINE qlgete(emin,emax)
DO ibpar=1,nblock ! parameter block DO ibpar=1,nblock ! parameter block
icoff=ioffBlock(ibpar) ! constraint offset in parameter block icoff=ioffBlock(ibpar) ! constraint offset in parameter block
iclast=ioffBlock(ibpar+1) ! last constraint in parameter block iclast=ioffBlock(ibpar+1) ! last constraint in parameter block
idiag=ncon*icoff+1 idiag=INT(ncon,mpl)*INT(icoff,mpl)+1
DO i=icoff+1,iclast DO i=icoff+1,iclast
IF (ABS(emax) < ABS(matL(idiag))) emax=matL(idiag) IF (ABS(emax) < ABS(matL(idiag))) emax=matL(idiag)
IF (ABS(emin) > ABS(matL(idiag))) emin=matL(idiag) IF (ABS(emin) > ABS(matL(idiag))) emin=matL(idiag)
...@@ -772,7 +769,7 @@ SUBROUTINE qlbsub(d,y) ...@@ -772,7 +769,7 @@ SUBROUTINE qlbsub(d,y)
icoff=ioffBlock(iblock) ! constraint offset in parameter block icoff=ioffBlock(iblock) ! constraint offset in parameter block
iclast=ioffBlock(iblock+1) ! last constraint in parameter block iclast=ioffBlock(iblock+1) ! last constraint in parameter block
nconb=iclast-icoff ! number of constraints in block nconb=iclast-icoff ! number of constraints in block
idiag=ncon*(iclast-1)+nconb idiag=INT(ncon,mpl)*INT(iclast-1,mpl)+nconb
DO k=nconb,1,-1 DO k=nconb,1,-1
y(k)=(d(k)-dot_product(matL(idiag+1:idiag+nconb-k),y(k+1:nconb)))/matL(idiag) y(k)=(d(k)-dot_product(matL(idiag+1:idiag+nconb-k),y(k+1:nconb)))/matL(idiag)
idiag=idiag-ncon-1 idiag=idiag-ncon-1
......
...@@ -51,7 +51,7 @@ ...@@ -51,7 +51,7 @@
!! 1. Download the software package from the DESY \c svn server to !! 1. Download the software package from the DESY \c svn server to
!! \a target directory, e.g.: !! \a target directory, e.g.:
!! !!
!! svn checkout http://svnsrv.desy.de/public/MillepedeII/tags/V04-09-00 target !! svn checkout http://svnsrv.desy.de/public/MillepedeII/tags/V04-09-01 target
!! !!
!! 2. Create **Pede** executable (in \a target directory): !! 2. Create **Pede** executable (in \a target directory):
!! !!
...@@ -256,7 +256,7 @@ ...@@ -256,7 +256,7 @@
!! !!
!! \subsection ch-lapack LAPACK !! \subsection ch-lapack LAPACK
!! For these optional methods the solution is obtained by matrix factorization from an external LAPACK library. !! For these optional methods the solution is obtained by matrix factorization from an external LAPACK library.
!! There exist (open or properiatary) implementations heavily optimized for specific hardware !! There exist (open or proprietary) implementations heavily optimized for specific hardware
!! (and partially multi-threaded) !! (and partially multi-threaded)
!! which could easily be an order of magnitude faster than e.g. the custom code used for !! which could easily be an order of magnitude faster than e.g. the custom code used for
!! \ref ch-mchdec "decomposition". !! \ref ch-mchdec "decomposition".
...@@ -1492,8 +1492,8 @@ SUBROUTINE prpcon ...@@ -1492,8 +1492,8 @@ SUBROUTINE prpcon
npar=ilast+1-ifrst npar=ilast+1-ifrst
nconmx=max(nconmx,ncon) nconmx=max(nconmx,ncon)
nparmx=max(nparmx,npar) nparmx=max(nparmx,npar)
mszcon=mszcon+ncon*npar ! (sum of) block size for constraint matrix mszcon=mszcon+INT(ncon,mpl)*INT(npar,mpl) ! (sum of) block size for constraint matrix
mszprd=mszprd+(ncon*ncon+ncon)/2 ! (sum of) block size for product matrix mszprd=mszprd+INT(ncon,mpl)*INT(ncon+1,mpl)/2 ! (sum of) block size for product matrix
IF (icheck > 0) THEN IF (icheck > 0) THEN
IF (ilast > ifrst) THEN IF (ilast > ifrst) THEN
labelf=globalParLabelIndex(1,globalParVarToTotal(ifrst)) labelf=globalParLabelIndex(1,globalParVarToTotal(ifrst))
...@@ -1535,18 +1535,18 @@ SUBROUTINE feasma ...@@ -1535,18 +1535,18 @@ SUBROUTINE feasma
INTEGER(mpi) :: i INTEGER(mpi) :: i
INTEGER(mpi) :: iblck INTEGER(mpi) :: iblck
INTEGER(mpi) :: icgb INTEGER(mpi) :: icgb
INTEGER(mpi) :: ij INTEGER(mpl) :: ij
INTEGER(mpi) :: ifirst INTEGER(mpi) :: ifirst
INTEGER(mpi) :: ilast INTEGER(mpi) :: ilast
INTEGER(mpi) :: ioffc INTEGER(mpl) :: ioffc
INTEGER(mpi) :: ioffp INTEGER(mpl) :: ioffp
INTEGER(mpi) :: irank INTEGER(mpi) :: irank
INTEGER(mpi) :: ipar0 INTEGER(mpi) :: ipar0
INTEGER(mpi) :: itgbi INTEGER(mpi) :: itgbi
INTEGER(mpi) :: ivgb INTEGER(mpi) :: ivgb
INTEGER(mpi) :: j INTEGER(mpi) :: j
INTEGER(mpi) :: jcgb INTEGER(mpi) :: jcgb
INTEGER(mpi) :: l INTEGER(mpl) :: ll
INTEGER(mpi) :: label INTEGER(mpi) :: label
INTEGER(mpi) :: ncon INTEGER(mpi) :: ncon
INTEGER(mpi) :: npar INTEGER(mpi) :: npar
...@@ -1601,7 +1601,7 @@ SUBROUTINE feasma ...@@ -1601,7 +1601,7 @@ SUBROUTINE feasma
factr=listConstraints(j)%value factr=listConstraints(j)%value
itgbi=inone(label) ! -> ITGBI= index of parameter label itgbi=inone(label) ! -> ITGBI= index of parameter label
ivgb =globalParLabelIndex(2,itgbi) ! -> variable-parameter index ivgb =globalParLabelIndex(2,itgbi) ! -> variable-parameter index
IF(ivgb > 0) matConstraintsT(ivgb-ipar0+ioffc+(jcgb-ifirst)*npar)=factr ! matrix element IF(ivgb > 0) matConstraintsT(INT(jcgb-ifirst,mpl)*INT(npar,mpl)+ivgb-ipar0+ioffc)=factr ! matrix element
globalParCons(itgbi)=globalParCons(itgbi)+1 globalParCons(itgbi)=globalParCons(itgbi)+1
rhs=rhs-factr*globalParameter(itgbi) ! reduce residuum rhs=rhs-factr*globalParameter(itgbi) ! reduce residuum
END DO END DO
...@@ -1609,12 +1609,14 @@ SUBROUTINE feasma ...@@ -1609,12 +1609,14 @@ SUBROUTINE feasma
END DO END DO
! get rank of blocks ! get rank of blocks
DO l=ioffc+1,ioffc+npar DO ll=ioffc+1,ioffc+npar
ij=ioffp ij=ioffp
DO i=1,ncon DO i=1,ncon
DO j=1,i DO j=1,i
ij=ij+1 ij=ij+1
matConsProduct(ij)=matConsProduct(ij)+matConstraintsT((i-1)*npar+l)*matConstraintsT((j-1)*npar+l) matConsProduct(ij)=matConsProduct(ij)+ &
matConstraintsT(INT(i-1,mpl)*INT(npar,mpl)+ll)* &
matConstraintsT(INT(j-1,mpl)*INT(npar,mpl)+ll)
END DO END DO
END DO END DO
END DO END DO
...@@ -1622,7 +1624,7 @@ SUBROUTINE feasma ...@@ -1622,7 +1624,7 @@ SUBROUTINE feasma
CALL sqminv(matConsProduct(ioffp+1:ij),vecConsResiduals(ifirst:ilast),ncon,irank, auxVectorD, auxVectorI) CALL sqminv(matConsProduct(ioffp+1:ij),vecConsResiduals(ifirst:ilast),ncon,irank, auxVectorD, auxVectorI)
IF (icheck > 1) WRITE(*,*) ' Constraint block rank', iblck, ncon, irank IF (icheck > 1) WRITE(*,*) ' Constraint block rank', iblck, ncon, irank
nrank=nrank+irank nrank=nrank+irank
ioffc=ioffc+npar*ncon ioffc=ioffc+INT(npar,mpl)*INT(ncon,mpl)
ioffp=ij ioffp=ij
END DO END DO
...@@ -7064,7 +7066,7 @@ SUBROUTINE loop2 ...@@ -7064,7 +7066,7 @@ SUBROUTINE loop2
WRITE(lu,101) 'NFGB',nfgb,'number of fit parameters' WRITE(lu,101) 'NFGB',nfgb,'number of fit parameters'
IF(metsol >= 4.AND. metsol <7) THEN ! band matrix as MINRES preconditioner IF(metsol >= 4.AND. metsol <7) THEN ! band matrix as MINRES preconditioner
WRITE(lu,101) 'MBANDW',mbandw,'band width of preconditioner matrix' WRITE(lu,101) 'MBANDW',mbandw,'band width of preconditioner matrix'
WRITE(lu,102) '(if =0, no preconditioner matrix)' WRITE(lu,102) '(if <0, no preconditioner matrix)'
END IF END IF
IF (nagb >= 65536) THEN IF (nagb >= 65536) THEN
WRITE(lu,101) 'NOFF/K',noff,'max number of off-diagonal elements' WRITE(lu,101) 'NOFF/K',noff,'max number of off-diagonal elements'
...@@ -7213,7 +7215,7 @@ SUBROUTINE loop2 ...@@ -7213,7 +7215,7 @@ SUBROUTINE loop2
CALL feasma ! prepare constraint matrices CALL feasma ! prepare constraint matrices
CALL vmprep(matsiz) ! prepare matrix and gradient storage IF (icheck <= 0) CALL vmprep(matsiz) ! prepare matrix and gradient storage
WRITE(*,*) ' ' WRITE(*,*) ' '
IF (matwords < 250000) THEN IF (matwords < 250000) THEN
WRITE(*,*) 'Size of global matrix: < 1 MB' WRITE(*,*) 'Size of global matrix: < 1 MB'
......
Supports Markdown
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