Laurence Marks píše v Po 14. 10. 2013 v 08:46 -0500:
> I was travelling so did not have access to the source. You are correct
> and it is technically a bug which apparently neither -ftrapuv nor
> gfortran picked up. I have added
>
> ROKMIX = 0.D0 ; ROKSC = 0.D0 ; ROKVL = 0.D0 ; ROKOLD = 0.D0
This is probably just a typo on your side, however ROKMIX = 0.D0 will
initialize just the real part... ROKMIX = (0.D0, 0.D0) is needed (at
lest on gfortran and xlf90)
> KVECVL = 0 ; KVECSC = 0 ; KVCOLD = 0
>
> after the allocation, and at some stage will add comparable to other
> allocations.
>
> That said, it should have no effect on the performance of mixer, and
> if it does then there is a serious bug somewhere else. The test on
> line 1504 of your version is just for some output of how the PW's are
> changing for monitoring as Peter has noticed that converging these can
> sometimes be tricky. There are a few other similar traps just for
> cleaning things up.
I didn't saw any actual problems except the crash in f7splt, I was just
being extra cautious, since I imagine the AIX 5.3 build doesn't see much
testing.
>
> Can you please send the lines for LimitDMIX8n.F in your version since
> it is probably different from mine and I do not see anything relevant
> at line 135.
LimitDMIX8n.F:135 PM1=min(1.D0,PM1)
attached is the stock LimitDMIX8n.F from WIEN2k_13.1 version
> N.B., concerning valgrind and mixer (for instance) not deallocating
> everything. There is a routine called "cleaner.f" which is in
> SRC_mixer for this purpose. However, with a multiply branching code it
> is very, very, very hard to trap everything and I have ended up hoping
> that the compiler does what it is supposed to. It is also very hard to
> protect for all possible systems and compilers particularly if I have
> no access to them. And.....mpi is much worse.
It wasn't actually leaking that much at all, just the small leak in
lapw1 (not counting the allocated memory still accessible on exit)
Basically even though my platform of interest is AIX 5.3, the testing
with valgrind was done on linux with gfortran, because valgrind isn't
working on AIX. As I've said before, my main motivation for this vagrind
exercise was just that I was really concerned that clearly invalid code
wasn't crashing at all neither with ifort (because it is initialized to
0 defaultly) nor with gfortran (it was probably optimized away, because
with -O0 it crash in the same way) and I was wondering if there are more
similar cases that are not severe enough to result in crash but rather
some unwanted behavior.
Also all the testing was done just with the TiC sample, so I probably
haven't tested that many codepaths at all...
Anyway thanks for looking into this and let me know if you need any more
info.
> On Thu, Oct 10, 2013 at 10:47 AM, Pavel Ondračka
> <pavel.ondra...@email.cz> wrote:
> > Laurence Marks píše v Čt 10. 10. 2013 v 09:40 -0500:
> >> Sorry, but I agree with Peter & I am 99.9% certain that this is not a
> >> bug in the mixer. The way to test this is (with ifort) to use
> >> -ftrapuv ;
> > That is not true. According to ifort docs:
> > -ftrapuv
> > The option sets any uninitialized local variables that are allocated on
> > the _stack_ to a value that is typically interpreted as a very large
> > integer or an invalid address.
> >
> > However the ROKMIX is allocated on the _heap_
> >
> > Just try the test I've described earlier, set the imaginary part by hand
> > right after allocation and check the value right before the if check
> > either with debugger or just write it to stdout and you will see the
> > value doesn't change.
> >> the arrays are not set in mixer.F but elsewhere within some
> >> complicated subroutines. This flag forces a fault if an undefined
> >> variable is used.
> >
> >
> >> ---------------------------
> >> Professor Laurence Marks
> >> Department of Materials Science and Engineering
> >> Northwestern University
> >> www.numis.northwestern.edu 1-847-491-3996
> >> "Research is to see what everybody else has seen, and to think what
> >> nobody else has thought"
> >> Albert Szent-Gyorgi
> >>
> >>
> >> On Oct 10, 2013 7:23 AM, "Pavel Ondračka" <pavel.ondra...@email.cz>
> >> wrote:
> >> Peter Blaha píše v Čt 10. 10. 2013 v 12:22 +0200:
> >> > I'm not very familiar with valgrid, but all those reports
> >> seem not
> >> > relevant to me.
> >> > It seems to complain about all allocated arrays, which are
> >> not set to
> >> > zero globally. But this does NOT mean that one uses
> >> uninitialized
> >> > variables or assumes that the compiler sets them to zero.
> >> Valgrind doesn't complain about uninitialized variables in
> >> general. Just
> >> when they are used in such way that they can alter the flow of
> >> the
> >> program or generally alter the externally-visible behavior.
> >>
> >> I'll take the problem from mixer as an example, because
> >> contrary to the
> >> rest of the problems, here I'm 99% sure it is a real bug.
> >>
> >> The complex array ROKMIX in mixer.F, is only partially
> >> initialized. The
> >> real part of the complex number is initialized somewhere
> >> (don't know
> >> exactly where, since I haven't tracked it that much), however
> >> the
> >> imaginary part is not. You can check this by setting the
> >> imaginary part
> >> to some specific value right after the allocation, and then
> >> set a
> >> breakpoint or print its value at line
> >> mixer.F:1504 if(abs(ROKMIX(J,I)).gt.1D-12)then
> >> and you will see, it has the same value as set before, so it
> >> wasn't
> >> written anywhere.
> >> And now basically if real part of ROKMIX(J,I) is lower than
> >> 1D-12 we
> >> would randomly fail or pass the if test, depending on what
> >> stuff is
> >> written in the corresponding uninitialized memory.
> >>
> >> >
> >> > PS: I don't have a variable "WS" in dergl.f
> >> > WS=0.0D0 in dergl.f
> >> My mistake, this should have been RW
> >> >
> >> > We can verify this also with ifort by setting -C in the
> >> flags, and in
> >> > fact doing so would result in an error in f7splt due to the
> >> mistake you
> >> > mentioned before.
> >> >
> >> > So the only problem was in f7splt.f and that should be fixed
> >> as
> >> > previously mentioned.
> >> Well, this is not what valgrind does. Valgrind doesn't do
> >> static code
> >> analysis as compilers do during compilation, but it replaces
> >> calls to
> >> malloc and similar at runtime with his own versions and tracks
> >> the used
> >> memory at runtime.
> >>
> >> To be more specific, I'll try to explain on my previous
> >> example why I
> >> think it can't be discovered by the ifort checking.
> >> Although ifort probably can track, that imaginary part of
> >> ROKMIX(J,I) is
> >> uninitialized at that point, it has no way of knowing if
> >> return value of
> >> abs(ROKMIX(J,I)) actually depends on the uninitialized value
> >> or not,
> >> because abs function has not been linked at that point.
> >>
> >> >
> >> > (PS: I was wrong with my first reply that you need ISPLIT=15
> >> (this
> >> > relates to an older version). The present if statement is
> >> intentional,
> >> > but still, no result from f7splt will be used (except for
> >> ISPLIT=15).
> >> Thanks, this explains it.
> >>
> >> I'm very grateful and I appreciate it very much, that you are
> >> looking
> >> into this, even though this is not a problem at your platform.
> >> Thanks again.
> >>
> >> _______________________________________________
> >> Wien mailing list
> >> Wien@zeus.theochem.tuwien.ac.at
> >> http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
> >> SEARCH the MAILING-LIST at:
> >>
> >> http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html
> >> _______________________________________________
> >> Wien mailing list
> >> Wien@zeus.theochem.tuwien.ac.at
> >> http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
> >> SEARCH the MAILING-LIST at:
> >> http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html
> >
> >
> > _______________________________________________
> > Wien mailing list
> > Wien@zeus.theochem.tuwien.ac.at
> > http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
> > SEARCH the MAILING-LIST at:
> > http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html
>
>
>
Subroutine LimitDMIX8n(Y,S,YY,SY,F,MAXMIX,MEMALL,MEMORY, &
qmx_input, dmix_last, dmixM, &
PM1, qtot, &
RedPred, RedOld, DIAG,STEP,YPROJ,ISKIP)
!
! Bound the step based just upon prior history and NFL bound
! Returns the projected part as well as part of the step
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program 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 General Public License for more details.
!
! Copyright (C) L. D. Marks, July 2010, 2011
!
use varyatoms
use TrustRegion
use MatrixNorms
implicit real*8 (a-h,o-z)
dimension Y(MAXMIX,MEMORY),YY(MEMORY,MEMORY),F(MAXMIX),S(MAXMIX,MEMORY)
dimension STEP(MAXMIX),YPROJ(MEMORY), SY(MEMORY,MEMORY)
real*8 LowerBound, ChangeStep
parameter (LowerBound=0.025 )
!
! Some Debug
dimension A(MEMORY,MEMORY),WR(MEMORY),WI(MEMORY),VL(1,MEMORY),VR(1,MEMORY)
DIMENSION WORK(4*MEMORY), WR2(MEMORY),WI2(MEMORY), WR3(MEMORY),WR4(MEMORY)
logical debug
#ifdef DEBUG
parameter (debug=.true.)
#else
parameter (debug=.false.)
#endif
! logical IsOK, Ismore
!
LWORK=4*MEMORY
A=SY
call DGEEV('N', 'N', MEMORY, A, MEMORY, WR, WI, VL, 1, VR, &
1, WORK, LWORK, INFO )
A=YY
call DGEEV('N', 'N', MEMORY, A, MEMORY, WR2, WI2, VL, 1, VR, &
1, WORK, LWORK, INFO )
A=SS
call DGEEV('N', 'N', MEMORY, A, MEMORY, WR4, WI2, VL, 1, VR, &
1, WORK, LWORK, INFO )
A=YY*Ylambda+SY*Slambda
call DGEEV('N', 'N', MEMORY, A, MEMORY, WR3, WI2, VL, 1, VR, &
1, WORK, LWORK, INFO )
!
write(21,222)
222 format(/'Eigenvalues, unscaled except for SY+YY'/ &
' # SY Real SY Imag SS YY SY+YY Real SY+YY Imag')
DO J=1,MEMORY
write(21,221)J,WR(J),WI(J),WR4(J),WR2(J),WR3(J),WI2(J)
ENDDO
!
write(21,*)
221 format(i4,8ES15.6)
!
!--------------------------------------------------------------------
! Limit based upon improvement
AllBetter= FHIST(MEMALL)/FHIST(MEMORY)
!
! Should this be a weighted average -- Average Change(M)*Pratt(M) ?
!
#ifndef NoMem
ChangeStep=1.5D0
HistPred(IHST)=RedOld
HistRed (IHST)=min(AllBetter,9.999D0)
JJ=1
J1=2
J2=IHST
UseBetter = min(ChangeStep,max(HistRed(JJ),1.D0/ChangeStep))
if(Rdebug)write(21,*)'Check UseBetter Loop'
if(Rdebug)write(21,*)JJ,UseBetter
DO M=J1,J2
JJ=JJ+1
Change=HistRed(JJ) ! FHIST(M+1)/FHIST(M)
Change=min(ChangeStep,max(Change,1.D0/ChangeStep)) ! Limit Range
UseBetter=UseBetter*RunX2+Change*RunX1
if(Rdebug)write(21,*)JJ,UseBetter,HistRed(JJ)
enddo
#else
HistRed(MEMORY)=AllBetter
ChangeStep=5.D0/3.D0
if(Safer)ChangeStep=4.D0/3.D0
UseBetter = min(ChangeStep,max(HistRed(1),1.D0/ChangeStep))
DO M=2,MEMoRY
Change=HistRed(M) ! FHIST(M+1)/FHIST(M)
Change=min(ChangeStep,max(Change,1.D0/ChangeStep)) ! Limit Range
UseBetter=UseBetter*RunX2+Change*RunX1
enddo
#endif
!
! ***************** Moved DMIX part here ***********************
if(memory .eq. 1.and.(.not.Restart))then
dmix_last=max(dmix_last,0.015D0)
! Experiment here, min of Shanno-Phua and <S>/<Y>
red=min(abs(SY(1,1))/YY(1,1),1.D0/limitscl)
if(red .gt. dmix_last)then
dmix_last=min(2.D0*dmix_last,0.05D0,red)
! else if(red .lt. dmix_last)then
! dmix_last=max(dmix_last*0.5,0.01D0,red)
endif
endif
! Limit based upon improvement over previous step
dmix_tmp=dmix_last/UseBetter
!
if(UseBetter .lt. 1.0D0)then
qlimit1=min(ChangeStep*dmix_last,dmix_tmp)
else
qlimit1=max(dmix_last/ChangeStep,dmix_tmp)
endif
!
if(memory .eq. 1)qlimit1=max(qlimit1,0.015d0)
if(abs(qmx_input-qlimit1).lt.0.01D0)qlimit1=qmx_input
! Why use qlimit2 at all ??????
qlimit2=qmx_input
if(.not.SAFER)qlimit2=min(1.D0,qmx_input*max(1.D0,1.D0/limitscl))
!
dmix =qlimit1
dmixM=qlimit1
21 format(a,6D11.3)
PM1=min(1.D0,PM1)
!
dmixM=dmix
dmixM=max(dmixM,dbase)
dmix=dmixM
! Multisecant Shanmo-Phua
UpShanno=0.25D0 ! Was 1/3
! We use Zlimit to help prevent stagnation
Zlimitscl = max(ZLimitscl,dbase*2.0D0,0.03D0)
if(dmix .lt. Zlimitscl*UpShanno)then !/3.D0)then
dmix = min(dmix*ChangeStep,Zlimitscl*UpShanno) !/3.D0)
write(21,21)':INFO : Increased using Multisecant Shanno-Phua'
else if (dmix .gt. Zlimitscl)then
dmix = max(dmix/ChangeStep,Zlimitscl)
write(21,21)':INFO : Decreased using Multisecant Shanno-Phua'
endif
! Obey user limits, enhanced
dmix=min(qlimit2,dmix)
! Limit Change to a factor of 2.5 (not ChangeStep)
if(dmix .gt. dmix_last)then
dmix=min(dmix,dmix_last*2.5D0) !ChangeStep)
else if(dmix .lt. dmix_last)then
dmix=max(dmix,dmix_last/2.5D0) !ChangeStep)
endif
dmixm=dmix
! ****************** End DMIX part here ***********************
!
! RedOld is what we predicted we would get as a reduction
! RedGot is what we in fact achieved
!
if(Regularize)then
if(MSec_Mode .ge. 0)then
call RegM2(Y,F,YY,MAXMIX,MEMORY,DIAG,ISKIP,IFAIL)
if(IFAIL .ne. 0)DIAG=3.5d-4
else
call RegGCV (Y,S,F,YY,SY,MAXMIX,MEMORY,DIAG,ISKIP,IFAIL)
if(IFAIL .ne. 0)then
MSec_Mode = 1
DIAG=3.5D-4
endif
endif
else if(MSec_Mode .lt. 0)then
call RegFixed (Y,S,F,YY,SY,MAXMIX,MEMORY,DIAG,ISKIP,IFAIL)
endif
!
! Look at projected fractions for No Free Lunch Control for different modes
if(UseTrust)then
allocate(YYRaw(MEMORY,MEMORY), YProj_MSec(MEMORY),MSEC_STEP(MAXMIX) )
YYRAW=YY
endif
if(MSec_Mode .lt. 0)then
PM1= FprojmemR1m(Y,S,F,YY,SY,MAXMIX,MEMORY,DIAG,STEP,YPROJ,CondN)
if(UseTrust)then
! Get the MSEC step
YY=YYRAW
DIAGC = 3.5D-4
PMM=Fprojmsec (Y,F,YY,MAXMIX,MEMORY,DIAGC,MSEC_STEP,YPROJ_MSec,CondN2)
!
! Next, get the MSEC Cauchy Step (creating SS in the process, needed later)
YY=YYRAW
DIAGC = 1D-6
call CauchyMSEC(S,Y,YY,SY,MAXMIX,MEMORY,DIAGC,YPROJ_MSEC )
! Last, the MSR1 Cauchy Step (with the SS from above)
YY=YYRAW
DIAGC = 1D-6
call CauchyMSR1(S,Y,YY,SY,MAXMIX,MEMORY,DIAGC,YPROJ)
YY=YYRAW
!
! Set Flags
HMSR1Cauchy=.true.
HMSR1 =.true.
HMSEC =.true.
HMSECCauchy=.true.
! Some Debug Information, let's probe the angles between the Pratt and Cauchy Steps
APMSEC=ddot(MAXMIX,MSEC_Cauchy,1,F,1)
APMSR1=ddot(MAXMIX,MSR1_Cauchy,1,F,1)
APMM =ddot(MAXMIX,MSR1_Cauchy,1,MSEC_Cauchy,1)
FF =dnrm2(MAXMIX,F,1)
DMSEC =dnrm2(MAXMIX,MSEC_Cauchy,1)
DMSR1 =dnrm2(MAXMIX,MSR1_Cauchy,1)
APMSEC=APMSEC/FF/DMSEC
APMSR1=APMSR1/FF/DMSR1
APMM =APMM /DMSR1/DMSEC
if(MEMORY .eq. 1)APMM=1.D0
if(APMM .gt. 0.999999D0)then
APMSEC=0.D0
else
APMSEC=acos(APMSEC)*180.D0/acos(-1.D0)
endif
if(APMSR1 .gt. 0.999999D0)then
APMSR1=0.D0
else
APMSR1=acos(APMSR1)*180/acos(-1.D0)
endif
endif
else
PM1= Fprojmsec (Y,F,YY,MAXMIX,MEMORY,DIAG,STEP,YPROJ,CondN)
! MSEC Cauchy Step
if(UseTrust)then
YY=YYRAW
DIAGC = DIAG
call CauchyMSEC(S,Y,YY,SY,MAXMIX,MEMORY,DIAGC, YPROJ )
HMSEC =.true.
HMSECCauchy=.true.
YY=YYRAW
endif
endif
CondN = DIAG/CondN/MEMORY
if(MEMORY .eq. 1)CondN=1
RedPred = PM1
!
MSECINFO(3)=RedOld
MSECINFO(4)=RedPred
MSECINFO(2)=min(AllBetter,9.999D0)
MSECINFO(6)=CondN
! Limit based upon improvement
#ifndef NoMem
HISTPred(IHST+1)=RedPred
#else
HistPred(MEMORY)=MSECINFO(3)
HistRed (MEMORY)=MSECINFO(2)
#endif
!
qlimit1=qmx_input
qlimit2=qmx_input
!
! Experiment: we have limitscl=<Y>/<S>, use it
!
! For Trust Region
if(UseTrust)then
if(MSec_Mode .lt. 0)then
!
! Add Unpredicted step to MSEC Cauchy
call daxpy(MAXMIX,DMIXM,MSEC_STEP,1,MSEC_Cauchy,1)
!
! And to MSR1 cauchy
call daxpy(MAXMIX,DMIXM,STEP,1,MSR1_Cauchy,1)
!
! Now generate the full MSR1 step
MSEC_Old = MSec_Mode
MSEC_Mode = 1
call MSEC8(S,MSEC_STEP,YPROJ_MSec,MAXMIX,MEMORY,DMIXM)
MSEC_Mode = MSEC_Old
FF =dnrm2(MAXMIX,F,1)
#ifdef DEBUG
write(21,1002)'Angle/Mag Pratt to MSEC Cauchy',APMSEC,FF/DMSEC*DMIXM
write(21,1002)'Angle/Mag Pratt to MSR1 Cauchy',APMSR1,FF/DMSR1*DMIXM
write(21,1002)'Angle MSEC to MSR1 Cauchy',APMM
1002 format(':INFOA: ',a,F9.2,10ES11.3)
#endif
else
call dcopy(MAXMIX,STEP,1,MSEC_STEP,1)
!
! Add Unpredicted step to Cauchy
call daxpy(MAXMIX,DMIXM,MSEC_STEP,1,MSEC_Cauchy,1)
!
! And the MSEC step (not really needed)
call MSEC8(S,MSEC_STEP,YPROJ,MAXMIX,MEMORY,DMIXM)
endif
endif
!
#ifdef DEBUG
write(21,21)':INFO : Bounds ',qlimit1,qlimit2,dmixM,dmix_last
#endif
if(allocated (YYRAW)) deallocate(YYRAW)
if(allocated (YProj_MSEC)) deallocate(YProj_MSEC)
return
end
_______________________________________________
Wien mailing list
Wien@zeus.theochem.tuwien.ac.at
http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
SEARCH the MAILING-LIST at:
http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html
