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.)
See Question&Answers more detail:
os