Fortran Do loop bug

[guruduttc]

Hello,

I am using Intel Visual Fortran 11 (version 075) for writing subroutines to implement them in Abaqus software in a Windows OS . I have come across a problem with the compiler over the past few days. Here is a small DO loop from a big subroutine. .


case 1 DO i = 1,ndofel
rhs(i,1)=0.d0
DO j = 1, ndofel
AMATRX(i,j) = 0.d0

END DO
ENDDO

where ndofel has a size of 54, RHS(54,1),AMATRX(54,54). So,basically I am just initializing them to zero.

Here is the same DO loop written in a different way,

case 2 DO i = 1,ndofel
DO j = 1, ndofel
AMATRX(i,j) = 0.d0
rhs(i,1)=0.d0

END DO
ENDDO

I am not an expert in Fortran,but I can say that Case 1 and Case 2 are the same.

However,when case 1 is used as a loop in my subroutine,I get alarmingly bad results. The results which I get when case 2 is used seems more realistic.

My question is ,shouldn't the final result be the same irrespective of using Case 1 or Case2? They both look the same to me. Or am I missing a point somewhere?!!

[mecej4]

Your second version sets rhs(i,1) to 0.d0 ndofel times, whereas doing that once within the DO I loop is sufficient. Neither version is good in terms of memory access, since the access stride is ndofel, which amounts to over 400 bytes, a rather large value.

You say you get "alarmingly bad results". I have no idea what you mean. Since 0.d0 is zero at all times, if your variables are dimensioned as stated there should be no difference in the values of the results. The running time can be different, as I stated in the first paragraph.

[TimP]

If the arrays don't overlap, the compiler should be able to change the order so as to zero out the data by setting blocks of memory to zero, e.g. by calling __intel_fast_memset.

[Les_Neilson]

As mecej4 says "alarmingly bad results" tells us nothing about what you mean, or get, in reality.
But ...
why don't you just use array assignment ?

rhs = 0.0d0
amatrx = 0.0d0

Les

[guruduttc]

Hello,

I have tried using a simple array assignment as recommended by Les, but the compiler gives me an error as follows :
error #6364: The upper bound shall not be omitted in the last dimension of a reference to an assumed size array. [RHS]
rhs=0.d0


Basically here is what happens,
I am writing a subroutine called UEL to implement in Abaqus

SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,&
PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME,&
KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF,&
LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD)

DIMENSION RHS(mlvarx,*),AMATRX(NDOFEL,NDOFEL),PROPS(*),&
SVARS(*),ENERGY(8),COORDS(MCRD,NNODE),U(NDOFEL),&
DU(MLVARX,*),V(NDOFEL),A(NDOFEL),TIME(2),PARAMS(*),&
JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),&
PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*)

*Note: mlvarx has the same dimension as NDOFEL


The dimensions of the RHS array is (54x1) . I need to initialize all the RHS and AMATRX(54x54) to zero.

Keeping all the other parts UEL subroutine same, the results I obtain when I use case 1 of array initialization is different from Case 2 of array initialization. Case 3 generates the error#6364 as mentioned above.
I must also mention that the results I get when I use case 1 are totally different from case 2.

case 1 DO i = 1,ndofel
rhs(i,1)=0.d0
DO j = 1, ndofel
AMATRX(i,j) = 0.d0

END DO
ENDDO

case 2 DO i = 1,ndofel
DO j = 1, ndofel
AMATRX(i,j) = 0.d0
rhs(i,1)=0.d0

END DO
ENDDO

case 3 rhs =0.d0
AMATRX= 0.d0
error #6364: The upper bound shall not be omitted in the last dimension of a reference to an assumed size array. [RHS]
rhs=0.d0

[Arjen_Markus]

Assumed-size arrays - arrays declared with an * as the last dimension - do not propagate
the actual size, so there is no way for the compiler to know how many elements they have.
Instead, use assumed-shape arrays: dimension(:,:,:) for instance.

(These require that the interface be known when the compiler encounters a call, and that
is easiest when you put the code in a module)

Regards,

Arjen

[TimP]

According to your usage, it seems you could have declared DIMENSION RHS(mlvarx,ndofel),...
so as to permit whole array assignment. Your report of problems raises the suspicion that something is wrong in the caller; possibly the arrays aren't originally declared with sufficient size, or you have violated the original Fortran rule about not passing the same array in multiple arguments.
ifort -gen-interfaces catches some possible violations in calling sequence. As Arjen said, Fortran features of the last 2 decades give you a portable way of preventing those violations.
ifort -assume dummy_aliases attempts to get around code breakage due to violations of the dummy array aliasing rule.

[guruduttc]

ok....Here is a modification I made to the code after I read ur reply...The problem persists even now,

SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,&
PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME,&
KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF,&
LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD)

DIMENSION RHS(mlvarx,nrhs),AMATRX(NDOFEL,NDOFEL),PROPS(*),&
SVARS(*),ENERGY(8),COORDS(MCRD,NNODE),U(NDOFEL),&
DU(MLVARX,nrhs),V(NDOFEL),A(NDOFEL),TIME(2),PARAMS(*),&
JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),&
PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*)

where, nrhs has a value 1

So,I no longer have assumed size arrays when it comes to RHS. But the problem with the case 1 and case 2 still persists . The error with case 3 no loner exists as there is no assumed size arrays when it comes to RHS.


There is also no declaration error with RHS as TIM pointed it out. Even if there is a declaration error how can the final results be different for case 1 and case 2 ? They are both the same ....

[jimdempseyatthecove]

[cpp]DIMENSION RHS(mlvarx,nrhs),AMATRX(NDOFEL,NDOFEL),...
               ^--------------------^

DO i = 1,ndofel
 rhs(i,1)=0.d0 ! ?? does nodfel == mlvarx ??
 DO j = 1, ndofel
   AMATRX(i,j) = 0.d0
        
[/cpp]

Jim Dempsey

[mecej4]

You have a complex problem. Given that the subroutine you wrote is being called by a closed-source package (Abaqus), the best that one can do in this forum is

(i) rely upon you to check that the subroutine adheres to the specifications laid down by Abaqus

(ii) check the source code for the subroutine to spot errors and inconsistencies in declarations or executable statements. This is possible only if you show all the relevant declarations and statements.

(iii) help you to ensure that, at subroutine exit, those variables that have INTENT(OUT) or INTENT(IN OUT) have correct values -- regardless of whether you have declared INTENTs, which you would not have done in Fortran-77 code.

So far, you have only shown disjointed pieces of code, and made an unsubstantiated statement about the results being wrong.

I don't dispute that the results are wrong, but I claim that we have no way of agreeing or disagreeing with your claim based on the facts that you have presented. We also have the handicap of not knowing the specifications for the argument list of the subroutine.

It does not help when you make confusing statements such as this one in #6

"Even if there is a declaration error how can the final results be different for case 1 and case 2 ? They are both the same"

If the declaration error caused failure to compile, there is is no question of results being different. There are no results in this case.

If the declaration error generates a warning, it is up to you to note it, perhaps report the warning here, and take appropriate action.