checking for self-assignment in fortran overloaded assignment

Question

I am trying to implement a polynomial class with fortran 2003, with overloaded arithmetic operations and assignments. The derived type maintains allocatable list of term definitions and coefficients, like this

type polynomial
  private
    type(monomial),dimension(:),allocatable     ::  term
    double precision,dimension(:),allocatable   ::  coef
    integer                                     ::  nterms=0
  contains
   ...
end type polynomial

interface assignment(=)
   module procedure :: polynomial_assignment
end interface
   ...
contains
    elemental subroutine polyn_assignment(lhs,rhs)
      implicit none
      type(polynomial),intent(???)  :: lhs
      type(polynomial),intent(in)   :: rhs
 ...

I had to make it elemental because this is intended to be used as matrices of polynomials. That does work, for the most cases at least. However, I somehow got myself into concerns about self-assignment here. One can simply check the pointers to see if things are the same in C++, but it doesn't seem to be an option in Fortran. However the compiler do detect the self-assignment and gave me a warning. (gfortran 4.9.0)

When I have intent(out) for lhs, the allocatable entries for both lhs and rhs appeared to be deallocated on entry to the subroutine, which made sense since they were both p, and an intent(out) argument would first be finalized.

Then I tried to avoid the deallocation with an intent(inout), and check self-assignment by modifying one field in the lhs output

   elemental subroutine polyn_assignment(lhs,rhs)
      implicit none
      type(polynomial),intent(inout)  :: lhs
      type(polynomial),intent(in)   :: rhs
      lhs%nterms=rhs%nterms-5
      if(lhs%nterms==rhs%nterms)then
        lhs%nterms=rhs%nterms+5
        return
      end if
      lhs%nterms=rhs%nterms

Well, now this is what surprised me. When i do

p=p

It didn't make the test and proceeded, giving me a polynomial with 0 terms but no memory violations. Confused, I printed lhs%nterms and rhs%nterms inside the assignment, only to find that they are different!

What is even more confusing is that when I did the same thing with

call polyn_assignment(p,p)

It works perfectly and detected that both arguments are the same. I am puzzled how an interface of a subroutine can run differently from the subroutine itself.

Is there something special about assignment in Fortran 2003 that I've missed?

(First time to ask a question here. Please correct me if i didn't do it right.)

score 5 · Accepted Answer · answered Nov 19 '13 at 08:35

5

If you have a statement a = b that invokes defined assignment via a subroutine sub, the assignment statement is equivalent to call sub(a, (b)). Note the parentheses - the right hand side argument is the result of evaluating a parenthesised expression and is therefore not conceptually the same object as b. See F2008 12.4.3.4.3 for details.

Consequently, a = a is equivalent to call sub(a, (a)). The two arguments are not aliased. It is different from call sub(a,a), the latter may (depending on the specifics of the internals of sub, including dummy argument attributes) break Fortran's argument aliasing rules (e.g. in your example, a statement such as call polyn_subroutine(a,a) is illegal).

answered Nov 19 '13 at 08:35

IanH

21,026
2
37
59

Thank you very much! I have some difficulty understanding the reference of "parentheses enclosing" in the standard. unrelated to self-assignment, if we call `a=b`, then it says b is first copied to temporary c then `sub(a,c)` is called. However, to copy b to c, does the program invoke assignment? Otherwise, how is this copying done for a derived type with allocatable? It is fairly strange that this sounds like a iterative definition. Your last point was very important though! I guess I was indeed tempering with illegal argument. @IanH – Xiaolei Zhu Nov 19 '13 at 14:39
The conceptual "copy to create the temporary" is not considered assignment, so you don't have an issue with recursive definition. The temporary (assuming the compiler isn't smart enough to avoid one) "has the value" of the right hand side, or words like that. If a component in b was allocated, then the corresponding component in the temporary would be allocated, and the value of the allocated stuff would be copied across. – IanH Nov 19 '13 at 18:51
You're right, I did some more test and found that it generated another copy of the object by copying over every entry and the allocatable objects had their pointer copied. So the %nterms fields were distinct variables but the %term(:) field pointer to the same address and when lhs is deallocated rhs got deallocated as well. I changed lhs to intent(inout) and did deallocation/allocation only when the shape of lhs/rhs do not conform or lhs is unallocated, and now the problem seems to be solved. Thanks! – Xiaolei Zhu Nov 19 '13 at 19:30
That doesn't sound right (or perhaps I misunderstand what you are saying) - if `a = a` invokes defined assignment, the allocatable components of lhs and rhs should be completely independent. – IanH Nov 20 '13 at 02:03
i ran a short test program with ifort and gfortran. It appears that ifort do generate independent allocatable components for lhs and rhs but gfortran created lhs and rhs that shared the same allocatable object. Maybe its a bug in gfortran? gfortran: `allocation status of lhs and rhs: F F.\n allocate lhs%a.\n allocation status of lhs and rhs: T T`. with ifort: `allocation status of lhs and rhs: F F.\n allocate lhs%a.\n allocation status of lhs and rhs: T F` – Xiaolei Zhu Nov 20 '13 at 02:30
1

I agree - my tests also show erroneous behaviour in current gfortran. Reported as bug 59202. – IanH Nov 20 '13 at 06:04

checking for self-assignment in fortran overloaded assignment

1 Answers1

Linked

Related