Commit 7b06a922 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

Fixing possible integer overflows for elimination of constraints

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@200 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent 41b88125
......@@ -261,7 +261,7 @@ SUBROUTINE sqminl(v,b,n,nrank,diag,next,vk) !
next0=1
l=1
DO i=1,n
i8=int8(i)
i8=INT(i,mpl)
next(i)=i+1 ! set "next" pointer
diag(i)=ABS(v((i8*i8+i8)/2)) ! save abs of diagonal elements
END DO
......@@ -274,7 +274,7 @@ SUBROUTINE sqminl(v,b,n,nrank,diag,next,vk) !
j=next0
last=0
05 IF(j > 0) THEN
j8=int8(j)
j8=INT(j,mpl)
jj=(j8*j8+j8)/2
IF(ABS(v(jj)) > MAX(ABS(vkk),eps*diag(j))) THEN
vkk=v(jj)
......@@ -287,7 +287,7 @@ SUBROUTINE sqminl(v,b,n,nrank,diag,next,vk) !
END IF
IF(k /= 0) THEN ! pivot found
k8=int8(k)
k8=INT(k,mpl)
kk=(k8*k8+k8)/2
kkmk=kk-k8
IF(l == 0) THEN
......@@ -310,7 +310,7 @@ SUBROUTINE sqminl(v,b,n,nrank,diag,next,vk) !
IF(j < k) THEN
jk=jk+1
ELSE
jk=jk+int8(j)-1
jk=jk+INT(j,mpl)-1
END IF
v(jk)=v(jk)*vkk
vk(j)=v(jk)
......@@ -323,7 +323,7 @@ SUBROUTINE sqminl(v,b,n,nrank,diag,next,vk) !
!$OMP SCHEDULE(DYNAMIC,128)
DO j=n,1,-1
IF(j == k) CYCLE
j8=int8(j)
j8=INT(j,mpl)
jl=j8*(j8-1)/2
vjk =vk(j)/vkk
b(j) =b(j)-b(k)*vjk
......@@ -332,19 +332,19 @@ SUBROUTINE sqminl(v,b,n,nrank,diag,next,vk) !
!$OMP END PARALLEL DO
ELSE
DO k=1,n
k8=int8(k)
k8=INT(k,mpl)
kk=(k8*k8-k8)/2
IF(next(k) /= 0) THEN
b(k)=0.0_mpd ! clear vector element
DO j=1,k
IF(next(j) /= 0) v(kk+int8(j))=0.0_mpd ! clear matrix row/col
IF(next(j) /= 0) v(kk+INT(j,mpl))=0.0_mpd ! clear matrix row/col
END DO
END IF
END DO
GO TO 10
END IF
END DO ! end of loop
10 DO jj=1,(int8(n)*int8(n)+int8(n))/2
10 DO jj=1,(INT(n,mpl)*INT(n+1,mpl))/2
v(jj)=-v(jj) ! finally reverse sign of all matrix elements
END DO
END SUBROUTINE sqminl
......@@ -900,7 +900,7 @@ SUBROUTINE chdec2(g,n,nrank,evmax,evmin)
REAL(mpd), INTENT(OUT) :: evmax
nrank=0
ii=(INT8(n)*INT8(n+1))/2
ii=(INT(n,mpl)*INT(n+1,mpl))/2
DO i=n,1,-1
IF (g(ii) > 0.0_mpd) THEN
! update rank, min, max eigenvalue
......@@ -921,7 +921,7 @@ SUBROUTINE chdec2(g,n,nrank,evmax,evmin)
!$OMP PRIVATE(RATIO,JJ) &
!$OMP SCHEDULE(DYNAMIC,128)
DO j=1,i-1
jj=(INT8(j-1)*INT8(j))/2
jj=(INT(j-1,mpl)*INT(j,mpl))/2
ratio=g(ii+j)*g(ii+i) ! (I,J) (I,I)
g(jj+1:jj+j)=g(jj+1:jj+j)-g(ii+1:ii+j)*ratio ! (K,J) (K,I)
END DO ! J
......@@ -954,7 +954,7 @@ SUBROUTINE chslv2(g,x,n)
REAL(mpd), INTENT(IN OUT) :: x(n)
INTEGER(mpi), INTENT(IN) :: n
ii=(INT8(n)*INT8(n+1))/2
ii=(INT(n,mpl)*INT(n+1,mpl))/2
DO i=n,1,-1
dsum=x(i)
kk=ii
......@@ -1322,11 +1322,11 @@ SUBROUTINE dbsvxl(v,a,b,n) ! LARGE symm. matrix, vector
IF(j < i) THEN
ij=ij+1
ELSE
ij=ij+int8(j)
ij=ij+INT(j,mpl)
END IF
END DO
b(i)=dsum
ijs=ijs+int8(i)
ijs=ijs+INT(i,mpl)
END DO
END SUBROUTINE dbsvxl
......
......@@ -51,7 +51,7 @@
!! 1. Download the software package from the DESY \c svn server to
!! \a target directory, e.g.:
!!
!! svn checkout http://svnsrv.desy.de/public/MillepedeII/tags/V04-08-00 target
!! svn checkout http://svnsrv.desy.de/public/MillepedeII/tags/V04-08-01 target
!!
!! 2. Create **Pede** executable (in \a target directory):
!!
......@@ -583,6 +583,7 @@
!! + **25** Aborted, result vector contains NaNs
!! + **26** Aborted, too many rejects
!! + **27** Aborted, singular QL decomposition of constraints matrix
!! + **28** Aborted, no local parameters
!! + **30** Aborted, memory allocation failed
!! + **31** Aborted, memory deallocation failed
!! + **32** Aborted, iteration limit reached in diagonalization
......@@ -3335,7 +3336,7 @@ END SUBROUTINE mupdat
!! \param [in] j2 smaller index last group
!! \param [in] il subtrahends first row
!! \param [in] jl subtrahends first col
!! \param [in] sub subtrahends matrix
!! \param [in] sub subtrahends matrix ('small', size fits in 'mpi')
SUBROUTINE mgupdt(i,j1,j2,il,jl,sub)
USE mpmod
......@@ -3412,7 +3413,8 @@ SUBROUTINE mgupdt(i,j1,j2,il,jl,sub)
ib=ib-matParBlockOffsets(1,iblk)
ja=ja-matParBlockOffsets(1,iblk)
jb=jb-matParBlockOffsets(1,iblk)
ij=(ia*ia-ia)/2+vecParBlockOffsets(iblk) ! global ISYM offset
ij=ia
ij=(ij*ij-ij)/2+vecParBlockOffsets(iblk) ! global ISYM offset
k=il
ijl=(k*k-k)/2 ! ISYM index offset (subtrahends matrix)
DO ir=ia,ib
......@@ -6109,6 +6111,11 @@ SUBROUTINE loop1
! 111 FORMAT(I5,I10,F10.5,E12.4)
WRITE(*,101) 'NTGB',ntgb,'total number of parameters'
WRITE(*,101) 'NVGB',nvgb,'number of variable parameters'
! To avoid INT(mpi) overflows in diagonalization
IF (metsol == 2.AND.nvgb >= 46340) THEN
metsol=1
WRITE(*,101) 'Too many variable parameters for diagonalization, fallback is inversion'
END IF
! print overview over important numbers ----------------------------
......@@ -6460,7 +6467,7 @@ SUBROUTINE loop2
nprecond(1)=ncgb ! number of constraints for preconditioner
nprecond(2)=nvgb ! matrix size for preconditioner
ENDIF
noff8=int8(nagb)*int8(nagb-1)/2
noff8=INT(nagb,mpl)*INT(nagb-1,mpl)/2
! all (variable) parameter groups
length=napgrp+1
......@@ -7134,6 +7141,11 @@ SUBROUTINE loop2
END DO ! print loop
IF(nalcn == 0) THEN
CALL peend(28,'Aborted, no local parameters')
STOP 'LOOP2: stopping due to missing local parameters'
END IF
! Wolfe conditions
IF(0.0 < wolfc1.AND.wolfc1 < wolfc2.AND.wolfc2 < 1.0) GO TO 32
......@@ -7156,7 +7168,7 @@ SUBROUTINE loop2
! prepare matrix and gradient storage ------------------------------
!32 CONTINUE
32 matsiz(1)=int8(nagb)*int8(nagb+1)/2 ! number of words for double precision storage 'j'
32 matsiz(1)=INT(nagb,mpl)*INT(nagb+1,mpl)/2 ! number of words for double precision storage 'j'
matsiz(2)=0 ! number of words for single precision storage '.'
IF (matsto == 2) THEN ! sparse matrix
matsiz(1)=ndimsa(3)+nagb
......@@ -7486,7 +7498,7 @@ SUBROUTINE minver
CALL qlmlq(globalCorrections(ipoff+1:),1,.true.) ! Q^t*b
! correction from eliminated part
DO i=1,nfit
ioff1=((nfit+1)*nfit)/2+i+imoff
ioff1=(INT(nfit+1,mpl)*INT(nfit,mpl))/2+i+imoff
DO j=1,ncon
globalCorrections(i+ipoff)=globalCorrections(i+ipoff)-globalMatD(ioff1)*vecConsSolution(j)
ioff1=ioff1+nfit+j
......@@ -7601,7 +7613,7 @@ SUBROUTINE mchdec
CALL qlmlq(globalCorrections(ipoff+1:),1,.true.) ! Q^t*b
! correction from eliminated part
DO i=1,nfit
ioff1=((nfit+1)*nfit)/2+i+imoff
ioff1=(INT(nfit+1,mpl)*INT(nfit,mpl))/2+i+imoff
DO j=1,ncon
globalCorrections(i+ipoff)=globalCorrections(i+ipoff)-globalMatD(ioff1)*vecConsSolution(j)
ioff1=ioff1+nfit+j
......@@ -7690,7 +7702,7 @@ SUBROUTINE mdiags
CALL qlmlq(globalCorrections,1,.true.) ! Q^t*b
! correction from eliminated part
DO i=1,nfgb
ioff1=((nfgb+1)*nfgb)/2+i
ioff1=(INT(nfgb+1,mpl)*INT(nfgb,mpl))/2+i
DO j=1,ncgb
globalCorrections(i)=globalCorrections(i)-globalMatD(ioff1)*vecConsSolution(j)
ioff1=ioff1+nfgb+j
......@@ -8553,7 +8565,7 @@ SUBROUTINE xloopn !
npar=matParBlockOffsets(1,ib+1)-ipoff ! size of block (number of parameters)
ncon=matConsBlocks(1,icblst+1)-matConsBlocks(1,icboff+1) ! number of constraints in (parameter) block
DO i=npar-ncon+1,npar
ioff=((INT8(i)-1)*INT8(i))/2+imoff
ioff=(INT(i-1,mpl)*INT(i,mpl))/2+imoff
globalMatD(ioff+1:ioff+i)=0.0_mpd
END DO
END DO
......
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