Use of unlimited polymorphic object and performance question

[rudi-gaelzer]

I'm still learning how to do OO programing and have several doubts.  Although this could be a question better suited for a FORTRAN forum, I want to ask you about the performance of an executable created with polymorphic entities.

The attached code (tes_polymorphism.f90)  works.  It defines a upoly entity called "electron", which can be described (actually, the velocity distribution function) either by Maxwell-Boltzmann or by kappa (non-additive) statistics.

The caller program decides the definitive type of the electron.  I'm doing this using typed allocation, but it also works with either direct targeting or sourced allocation.  The advantage of typed allocation, IMHO, is that I don't need to explicitly declare variables of both types as would be necessary with sourced allocation.

Now, as I read in Metcalf's book , I cannot actually do anything with a upoly object.  It only serves as the selector in a "select type" construct.  In fact, I even have to define the pointer procedure component of the electron from inside the construct, otherwise I get the error:

error #6460: This is not a field name that is defined in the encompassing structure.   [PDF]

That's all right by me, if I only have to do this once in the program.  However, the function "disp_eq" needs to access the procedure component.  My intention was to use upoly entities so that I can reference the procedure by the same name, irrespective of the actual type of the electron. 

The only way I could make this work was by using select type also inside the function, as you can see in the code below.  However, it seems kinda silly to me, since I'm simply repeating twice the same line.

So, I have 2 questions for you:

1) (A fortran syntax question) Is this really the only way to do this?  Perhaps there's a trick I'm still not aware of...

2) (A performance question) The select type in the program is only accessed once, but the function disp_eq will be called several thousand times during run time.  Will the select type construct introduce an extra overhead, thereby reducing the code's performance?

[FortranFan]

Re: your question on performance, see a discussion at comp.lang.fortran on similar lines: link

I think you're on the right track with trying to pursue an object-oriented approach for such code.  See blog by Dr Fortran at Intel on recent books: https://software.intel.com/en-us/blogs/2013/12/30/doctor-fortran-in-its-a-modern-fortran-world and you may want to closely examine the one in the footnote by Rouson et al., Scientific Software Design: The Object-Oriented Way for ideas; the other books are helpful too.

I would suggest using type-bound procedures and have everything for "pre-wired", so the consumer the "class" (e.g., whoever is writing program tes_poly) doesn't have to do the right thing to ensure correct result such as pdf => zfn assignment which can be source of error, doing away with global variables in modules, using limited polymorphism as opposed to unlimited, consider pure and elemental procedures for your needs, and using allocatables instead of pointers as much as possible.  With such changes, the code can become simpler.  See a changed version below and you can decide for yourself:

module Electron_Poly

   implicit none

   type, abstract, public :: thermal
      real :: temp
   contains
      procedure, pass(this), public :: disp_eq
      procedure(Ipdf), pass(this), deferred, public :: pdf
   end type thermal
   
   type, extends(thermal), public :: part_maxwell
   contains
      procedure, pass(this), public :: pdf => zfn
   end type part_maxwell
   
   type, extends(thermal), public :: part_kappa_int
      real :: kappa
   contains
      procedure, pass(this), public :: pdf => zfnk
   end type part_kappa_int
   
   type, extends(part_kappa_int) :: part_kappa
   end type part_kappa
   
   abstract interface

      pure elemental function Ipdf(this, z) result(pdf_result)

         import :: thermal

         !.. argument list
         class(thermal), intent(in) :: this
         real, intent(in)           :: z
         !.. function result
         real :: pdf_result

      end function Ipdf

   end interface

contains

   pure elemental function zfn(this, z) result(pdf_result)

      class(part_maxwell), intent(in) :: this
      real, intent(in) :: z
      !.. function result
      real :: pdf_result

      !..
      pdf_result = erfc_scaled(z*z)

      return

   end function zfn

   pure elemental function zfnk(this, z) result(pdf_result)

      class(part_kappa_int), intent(in) :: this
      real, intent(in) :: z
      !.. function result
      real :: pdf_result

      !..
      pdf_result = this%kappa*erfc_scaled(z*z)

      return

   end function zfnk

   pure elemental function disp_eq(this, z) result(disp)

      class(thermal), intent(in) :: this
      real, intent(in) :: z
      !.. function result
      real :: disp

      !.. local variables
      real :: pdfe

      !..
      pdfe = this%pdf(z)

      disp = z*pdfe

      return

   end function disp_eq

end module Electron_Poly

program tes_poly

   use Electron_Poly, only : thermal, part_maxwell, part_kappa

   implicit none

   !.. local variables
   class(thermal), allocatable :: electron
   integer :: istat
   integer :: i
   real :: z, kappa
   character(len=80) :: err_alloc
   
   write(*, '(a)', advance= 'no')'Enter z: '
   read(*,*) z
   write(*, '(a)', advance= 'no')'Enter particle statistics (1/2): '
   read(*,*) i
   
   if (i == 1) then ! Electron is Maxwellian
      allocate(part_maxwell :: electron, stat=istat, errmsg=err_alloc)
   else           ! Electron is kappa.
      allocate(part_kappa :: electron, stat=istat, errmsg=err_alloc)
   end if
   if (istat /= 0) then
      write(*,*) " error allocating electron: stat = ", istat
      write(*,*) err_alloc
   end if
   
   select type (electron)
      class is (part_maxwell)
         print*, 'zfn(z)=', electron%pdf(z)
         print*, 'D.  eq=', electron%disp_eq(z)
      class is (part_kappa)
         write(*, '(a)', advance= 'no')'Enter kappa: '
         read(*,*) kappa
         electron%kappa= kappa
         print*, 'zfnk(kappa, z)=', electron%pdf(z)
         print*, 'Disp. eq=      ', electron%disp_eq(z)
      class default
         write(*,*) " electron is of unsupported type."
   end select

   deallocate(electron, stat=istat, errmsg=err_alloc)
   if (istat /= 0) then
      write(*,*) " error deallocating electron: stat = ", istat
      write(*,*) err_alloc
   end if

   stop

end program tes_poly

Upon execution,

Enter z: 2.0
Enter particle statistics (1/2): 1
 zfn(z)= 0.1369995
 D.  eq= 0.2739989
Press any key to continue . . .

[FortranFan]

Note the code in message #2 is just an example for you to consider; please do not see it as recommended design or the only way to write code.  There are many design possibilities and that's the beauty of OO approach.  OO analysis and design are arguably the most important aspects of OO approach and only you can do that best keeping in mind how the class will consumed, of course.  The analysis and design will inform you how to best setup your code.

[rudi-gaelzer]

FortranFan wrote:

Note the code in message #2 is just an example for you to consider; please do not see it as recommended design or the only way to write code.  There are many design possibilities and that's the beauty of OO approach.  OO analysis and design are arguably the most important aspects of OO approach and only you can do that best keeping in mind how the class will consumed, of course.  The analysis and design will inform you how to best setup your code.

Hi, FortranFan. Thank you for sending me the suggestions and the tips. I'll take a look at them.  Indeed, I want to invest more on OO techniques, as opposed to the traditional programming style commonly used with Fortran.  Alas, there's none around in my Department doing the same, so I have to slowly study by myself.

I still did not fully digest your tips, but a crucial point I noticed was that you bounded the disp_eq function to the electrons.  I don't think this works for me.  The code I sent is only a restricted example of what I'm attempting.  In real applications, I'll have other particles species (protons, alpha particles, etc), and the disp_ep has to handle all of them.  That's why I think it cannot be bounded to any particular particle species.

The problem is that each particle species can have its own properties: statistical distribution, kappa index (in some cases they can have 2 distinct kappas), distribution model, temperature ratio, etc.  There are several different combinations of particle properties that I want to run.  If I employed traditional programming to this problem, I either had to write several different version of the zfn* and disp_eq routines, and then possibly employ unbounded procedure pointers, or I could write a single routine and fill it with if's (which is ugly).

That is why I'm attempting to write the code employing unlimited polymorphic entities.  In this case it seems that I can define a single routine, but the price I pay is having to insert a select type construct for each particle species (either for a scalar particle class or an array of particles).  This means that I'll have to repeat the same lines, such as "pdfe= electron%pdf(z)" twice of more, as in the code I sent.  This seems redundant to me and raises the question of performance loss.

But this is just a first impression I got from your suggestions.  I'll take a closer look at them to see if they indeed can be applied to my problem.  But in my VERY limited knowledge of OO techniques, using limited polymorphism I'll be forced to define several different derived types.  If I could find a solution employing upoly objects, the code could be more compact (of course, with the possible downfall on the performance...).

Thanks again for the tips.

[FortranFan]

rudi-gaelzer wrote:

..

I still did not fully digest your tips, but a crucial point I noticed was that you bounded the disp_eq function to the electrons.  I don't think this works for me.  The code I sent is only a restricted example of what I'm attempting.  In real applications, I'll have other particles species (protons, alpha particles, etc), and the disp_ep has to handle all of them.  That's why I think it cannot be bounded to any particular particle species.

..

Clearly I have little understanding of the basic physics of your calculations as well as your needs.  Note, however, your starting point can be an abstract type called particle which may have certain basic properties such as mass, velocity/momentum, etc. and certain methods such as accelerate() and disp_eq that operate on it.  The other particles - electrons, protons, alpha particles, etc. - can all derive from this base type, each having all the properties of the base type but additional ones as well, depending on their characteristics such as charge.  What you may find advantageous over time is to have objects in your code (i.e., the actual particles in your computations) to derive from class(particle), rather than class(*).  This is what I'd suggest to you to keep in mind.

rudi-gaelzer wrote:

..

there's none around in my Department doing the same, so I have to slowly study by myself

..

Assume there are many C++ "dabblers" and "experts" in your department touting OO can only be done using their approach for compute-intensive scientific code (thus excluding Python, etc.)!?  If you study "modern" Fortran and employing the many feature sets from Fortran 2003, 2008, (and 2015), you can show to your physics world Fortran is fully first-class in terms of OO design patterns for scientific programming and that the resultant code is both high performing as well as easier to understand and maintain and improve for scientists.  This thread may be of some interest if you ever wonder about Fortran relative to C++: https://software.intel.com/en-us/forums/topic/509099; note I'm only bringing this up because you've indicated you're in academia and I've experienced first-hand some biases there and I know how Fortranners can suffer in their midst, with few people around to help them.

[rudi-gaelzer]

FortranFan wrote:

Note, however, your starting point can be an abstract type called particle which may have certain basic properties such as mass, velocity/momentum, etc. and certain methods such as accelerate() and disp_eq that operate on it.  The other particles - electrons, protons, alpha particles, etc. - can all derive from this base type, each having all the properties of the base type but additional ones as well, depending on their characteristics such as charge.  What you may find advantageous over time is to have objects in your code (i.e., the actual particles in your computations) to derive from class(particle), rather than class(*).

Assume there are many C++ "dabblers" and "experts" in your department touting OO can only be done using their approach for compute-intensive scientific code (thus excluding Python, etc.)!?  If you study "modern" Fortran and employing the many feature sets from Fortran 2003, 2008, (and 2015), you can show to your physics world Fortran is fully first-class in terms of OO design patterns for scientific programming and that the resultant code is both high performing as well as easier to understand and maintain and improve for scientists.  This thread may be of some interest if you ever wonder about Fortran relative to C++: https://software.intel.com/en-us/forums/topic/509099; note I'm only bringing this up because you've indicated you're in academia and I've experienced first-hand some biases there and I know how Fortranners can suffer in their midst, with few people around to help them.

Again, thanks for taking your time to help me out.  I was already aware of deferred procedures in abstract types, but didn't really know the difference between them and type-bound pointer procedures.  By your comments it seems that you would recommend the use of the former, right?

I was thinking how I could implement what I want with abstract interfaces and deferred procedures.  Let's see if I can explain it to you.   If I can build up an hierarchy of abstract types, with each level extending the previous one, and all of them with non-allocated deferred-type bound procedures, then I can create a final type "particle", which contains a bunch of non-allocated procedure components.  Then, when I create  particle-type variables called "electron" or "proton", etc, I can decide which one of the deferred procedures will point to a particular "zfn" function for each particle species.  If I could do that, I could take advantage of the property of extended types that references to parent components need not mention explicitly the parent type, exactly as the example in Metcalf's book:

type :: person
   character(len= 10) :: name
   ...
end type person
!
type, extends(person) :: employee
   ...
end type employee

type(employee) :: director

then director%name is the same as director%person%name.   In my case, the "zfn" function employed for each particle species will be the only one that was allocated.

However, the only hierarchy I was allowed to build was something like:

module temp
implicit none
type, abstract :: t1
   contains
      procedure(at1p), deferred :: t1p
end type t1
!
type, abstract, extends(t1) :: t2
   contains
      procedure(at2p), deferred :: t1p
end type t2
!
type, extends(t2) :: t3
   real :: a
   contains
      procedure t1p => ft1
end type t3
!
abstract interface
   function at1p(pars, x)
   import :: t1
   real :: at1p
   class(t1), intent(in) :: pars
   real, intent(in) :: x
   end function at1p
!***
   function at2p(pars, x)
   import :: t2
   real :: at2p
   class(t2), intent(in) :: pars
   real, intent(in) :: x
   end function at2p
end interface
 CONTAINS
   function ft1(pars, x)
   real :: ft1
   class(t3), intent(in) :: pars
   real, intent(in) :: x
   ft1= pars%a*x**2
   return
   end function ft1
end module temp

When I tried in the type t3 the instruction

procedure t1p => t2%ft1

I got a "wrong syntax" error message.  If it was possible to extend 2 abstract types then I guess I could implement what I want.  What do you think?

FortranFan wrote:

Assume there are many C++ "dabblers" and "experts" in your department touting OO can only be done using their approach for compute-intensive scientific code (thus excluding Python, etc.)!? 

I recently attended a PhD defense, where the candidate needed to solve numerically a bunch of nonlinear PDE's in 2+1 dimensions in order to get the results for the thesis.  To my surprise, one of the board members, once he learned that the candidate used fortran to write the code, suggested he'd use python instead... Well, of course, if the guy had 10 years to finish his PhD...

[FortranFan]

rudi-gaelzer wrote:

 ..  When I tried in the type t3 the instruction

procedure t1p => t2%ft1

I got a "wrong syntax" error message.  If it was possible to extend 2 abstract types then I guess I could implement what I want.  What do you think? ..

The expected syntax is along the lines of [fortran]procedure :: t1p => ft1[/fortran]I personally prefer to be explicit about describing the passed object such as pass(this) (it can also be pass(self), pass(me), pass(pars), etc.) and also about applying the "private" attribute as default to all derived type components and type bound-procedures and then giving public attributes to those methods need to be invoked by consumers of the derived type; this is as per the OO paradigm of "information hiding".

But coming back to your question on "If I can build up an hierarchy of abstract types, with each level extending the previous one, and all of them with .. deferred-type bound procedures..", yes you can do so in Fortran.  Once you  correct your syntax as suggested above, you can continue building your "class", and try it out.  By the way, I'm not sure what you mean by "non-allocated" deferred procedures though.

[Steven_L_Intel1]

I'll make a performance-related comment. When you pass a non-polymorphic item to a polymorphic dummy argument, the compiler has to create a "class descriptor" for each call. This is a lot of code and can have a performance impact.

[rudi-gaelzer]

FortranFan wrote:

The expected syntax is along the lines of

procedure :: t1p => ft1

I personally prefer to be explicit about describing the passed object such as pass(this) (it can also be pass(self), pass(me), pass(pars), etc.) and also about applying the "private" attribute as default to all derived type components and type bound-procedures and then giving public attributes to those methods need to be invoked by consumers of the derived type; this is as per the OO paradigm of "information hiding".

But coming back to your question on "If I can build up an hierarchy of abstract types, with each level extending the previous one, and all of them with .. deferred-type bound procedures..", yes you can do so in Fortran.  Once you  correct your syntax as suggested above, you can continue building your "class", and try it out.  By the way, I'm not sure what you mean by "non-allocated" deferred procedures though.

Ok, "non-allocated deferred procedures" is incorrect.  What I meant to say was, if I could have a hierarchy of types, each one with an unallocated procedure pointer component (all referenced with the same name t1p) and then only determine which specific routine the name t1p would execute at the top of the hierarchy, then I could choose which link in the hierarchy would point to the specific routine with an instruction similar to

procedure :: t2%t1p => ft1

However, since this is not permitted, I guess that the instruction

procedure :: t1p => ft1

means that the same routine (ft1) will be referenced in all levels of the hierarchy, right?

[FortranFan]

rudi-gaelzer wrote:

...  What I meant to say was, if I could have a hierarchy of types, each one with an unallocated procedure pointer component (all referenced with the same name t1p) and then only determine which specific routine the name t1p would execute at the top of the hierarchy, then I could choose which link in the hierarchy would point to the specific routine ..

means that the same routine (ft1) will be referenced in all levels of the hierarchy, right?

Perhaps it makes sense for you to take some time and go through the review material.  You mentioned Metcalf's book in the original post - I assume it's "Modern Fortran Explained."  If so, take another look at chapter 14, especially section 14.6 on type-bound procedures and subsections in there on specific and generic bindings, type extension and overriding procedure bindings.  I think your interest is in having flexibility with the bindings to different types within a hierarchy and this section will show you the various options that are available.  If I understood what you're saying correctly, your requirements should be quite feasible to achieve - see a sample below that you can try out.

I'll now step back.   But note the crucial aspect is for you to perform a thorough object-oriented analysis and design (OOAD) of your "classes" and only then start to write code which follows the design; this while knowing there is quite a bit of flexibility available at your disposal - it's up to you to design your classes well; it is not about right or wrong but it is more about functionality, utility, efficiency, elegance, aesthetics. etc. - think of architecture which is usually mentioned in connection with OO design.  Look into UML diagrams, if you're not using them already - these can be of great help with development, implementation, maintenance, and communication of OO classes.

module temp_m

   implicit none

   private

   type, abstract :: t1
      private
   contains
      private
      procedure(It1p), pass(pars), deferred, public :: t1p
   end type t1

   type, extends(t1), abstract :: t2
      private
   contains
      private
      procedure(It2p), pass(pars), deferred, public :: t2p
   end type t2

   type, extends(t2), public :: t3   !.. t3 extends abstract t2 while providing a concrete implementation
      private
      real :: m_a
   contains
      private
      procedure, pass(pars), public :: t1p => t1p_t3
      procedure, pass(pars), public :: t2p => t2p_t3
      procedure, pass(pars), public :: set_a => set_a_t3
   end type t3

   type, extends(t3), public :: t4
      private
      real :: m_b
   contains
      private
      procedure, pass(pars), public :: t1p => t1p_t4  !.. note how this overrides the binding in t3
      procedure, pass(pars), public :: set_b => set_b_t4
   end type t4

   abstract interface

      function It1p(pars, x) result(y)

         import :: t1

         !.. argument list
         class(t1), intent(in) :: pars
         real, intent(in)      :: x
         !.. function result
         real :: y

      end function It1p

      function It2p(pars, x) result(y)

         import :: t2

         !.. argument list
         class(t2), intent(in) :: pars
         real, intent(in)      :: x
         !.. function result
         real :: y

      end function It2p

   end interface

contains

   function t1p_t3(pars, x) result(y)

      !.. argument list
      class(t3), intent(in) :: pars
      real, intent(in)      :: x
      !.. function result
      real :: y

      y = pars%m_a * x**2

      return

   end function t1p_t3

   function t2p_t3(pars, x) result(y)

      !.. argument list
      class(t3), intent(in) :: pars
      real, intent(in)      :: x
      !.. function result
      real :: y

      y = pars%m_a * x**3

      return

   end function t2p_t3

   subroutine set_a_t3(pars, a)

      class(t3), intent(inout) :: pars
      real, intent(in)         :: a

      pars%m_a = a

      return

   end subroutine set_a_t3

   function t1p_t4(pars, x) result(y)

      !.. argument list
      class(t4), intent(in) :: pars
      real, intent(in)      :: x
      !.. function result
      real :: y

      y = pars%m_b * x + pars%m_a * x**2

      return

   end function t1p_t4

   subroutine set_b_t4(pars, b)

      class(t4), intent(inout) :: pars
      real, intent(in)         :: b

      pars%m_b = b

      return

   end subroutine set_b_t4

end module temp_m

program p

   use temp_m, only : t3, t4

   implicit none

   type(t3) :: foo
   type(t4) :: bar

   call foo%set_a( 1.0 )

   print *, " foo%t1p = ", foo%t1p( 2.0 )
   print *, " foo%t2p = ", foo%t2p( 3.0 )

   call bar%set_a( 2.0 )
   call bar%set_b( 3.0 )

   print *, " bar%t1p = ", bar%t1p( 1.0 )

   stop

end program p

  foo%t1p =  4.000000
  foo%t2p =  27.00000
  bar%t1p =  5.000000
Press any key to continue . . .

[rudi-gaelzer]

FortranFan wrote:

I'll now step back.   But note the crucial aspect is for you to perform a thorough object-oriented analysis and design (OOAD) of your "classes" and only then start to write code which follows the design; this while knowing there is quite a bit of flexibility available at your disposal - it's up to you to design your classes well; it is not about right or wrong but it is more about functionality, utility, efficiency, elegance, aesthetics. etc. - think of architecture which is usually mentioned in connection with OO design.  Look into UML diagrams, if you're not using them already - these can be of great help with development, implementation, maintenance, and communication of OO classes.

Hey, FortranFan, thanks for all your help.

I think now I cracked it. Your suggestions and tips helped a lot and, as I suspected, the solution is far simpler than I was thinking originally.  I attached a simple example, but I also managed to write another version where, using allocatable arrays, the code can handle an arbitrary number of particle species, each one with its own methods, and have all of them passed to the disp_eq routine that can handle all of them, as I wanted.

Just a last question.  Regarding OOAD and UML diagrams, is Rouson's book "Scientific Software Design" a good place to start?  I noticed that it has examples in F2003.

Again, thanks a lot.

[FortranFan]

rudi-gaelzer wrote:

..  Just a last question.  Regarding OOAD and UML diagrams, is Rouson's book "Scientific Software Design" a good place to start?  I noticed that it has examples in F2003. ..

Yes, and you can look at the references therein also.

Re: your attached example in the latest post, it still looks conceptually similar to that in the original post.  See message #2 above.  I'd view program tes_poly as being a "consumer" of particle class and from a OO design point-of-view, recommended practice is to create a demarcation between what a consumer can do and what the class does; a general OO view is that the consumer should only have to "instantiate" the class and invoke methods on the instantiated object, it should not directly access/assign to data nor do bindings, etc.  This goes into the whole line of thinking in terms of "information hiding" and the benefits it can provide in the long term.  If you subscribe to this notion, then you should avoid having to do "electron%pdf => zfn" and so forth on the caller side.  Just my 2 cents,