This is a corrected version for something I missed and some typesetting
issues with the code (that was Outlooks fault) of the message I just sent:
Hi,
Here is my attempt to use F2K syntax to give (in my view) a robust and nice
way of implementing a reverse communication interface. But I am not really
sure that what I am proposing here will work with the F2K standard, in
particular because I am not sure how type-bound procedures are passed as
dummy-procedure arguments, and after reading some of section 12.4 in the F2K
draft-proposal, I did not find an answer to my questions. In particular, I
suspect the routine Multiply below should in fact be made a polymorphic
pointer procedure, but I am not sure how and why...Experts, please look over
this when you have time.
I will use my previous Conjugate Gradient code as an example, with the main
solver routine being called CG, placed in a module Conjugate_Gradient. In
CG, which solves the system Ax=b, what is required from the user is to
multiply a given vector Vx with some matrix A and store the result in Vy.
The extents of all vectors are (n_lb : n_ub), which is an argument to CG. We
will agree that the user will provide a routine Multiply to CG that will do
this. This whole argument is about how to interface Mulitply.
First issue is how to declare and provide Vx and Vy to the user. They could
be passed on to Multiply as assumed-shape array arguments. I dislike this
approach for reasons that become obvious when doing parallel calculations.
In my case for example, to perform this multiplication it is required that
additional, *shadow* space (shadow edges) be allocated on Vx and Vy, where
non-local variables will be copied to avoid, say, write conflicts on a
shared-memory machine, or to hoist communication out of the loops on a
distributed-memory machine. The point is, that the user needs to worry about
how exactly Vx and Vy are allocated. I see two approaches:
1. Make Vx and Vy public array pointers in Conjugate_Gradient, and let the
user associate them with the correct portions of the arrays he needs to
operate on (the non-shadow region). This is rather elegant, but there is a
problem--Vx and Vy may be used in intensive loops in CG, so that aliasing
may be an optimizaiton hinderance. To resolve this, I will use the trick I
mentioned yesterday of making CG a wrapper around another routine which will
be passed, as done in the code below. This has some dangers though as
Multiply has effects on arguments it is not using, and some optimizers may
trip over this, though I believe the standard allows for that so long as
there is no PURE attribute on Multiply?
2. Make Vx and Vy allocatable arrays that the user should allocate according
to his needs and make them public--less elegant, but avoids some issues and
a wrapping process.
Here is the code for CG which makes Multiply have no arguments whatsoever
(and so relies on its effect on global data). The first usage of OOP
features comes when one considers what happens if one needs to use several
CG's at the same time (say a CG as a preconditioner for another CG). Then
only two arrays Vx and Vy are not enough, so we need a class type solution,
in which each CG routine gets an associated Vx and Vy arguments.
MODULE Conjugate_Gradient
PRIVATE
...
TYPE, PUBLIC, EXTENSIBLE :: CG_Solver
...
REAL, DIMENSION(:), POINTER, SAVE :: Vx, Vy
CONTAINS
PROCEDURE, PASS_OBJ :: CG
END TYPE
...
CONTAINS
...
! This is a wrapper routine for the real solver
SUBROUTINE CG(CG_Solver, Multiply, rhs_vector, solution_vector, ...)
CLASS(CG_Solver), INTENT(INOUT) :: CG_Solver ! PASSED IMPLICTLY
INTERFACE ! This may be *wrong*.
! A pointer procedure may be what I am thinking of
SUBROUTINE Multiply() ! No arguments whatsoever
END SUBROUTINE Multiply
END INTERFACE
REAL, DIMENSION(:), INTENT(IN) :: rhs_vector ! Vector b
REAL, DIMENSION(:), INTENT(INOUT) :: solution_vector
...
! To avoid aliasing, I call CG_internal with Vx and Vy passed on
CALL CG_internal(Multiply, rhs_vector, solution_vector,...,&
CG_Solver%Vx, CG_Solver%Vy)
END SUBROUTINE CG
! This routine has a classical interface and is array-based, so it is
efficient
! and can be borrowed from older code
SUBROUTINE CG_internal(Multiply, rhs_vector, solution_vector,...,Vx, Vy)
INTERFACE
SUBROUTINE Multiply() ! No arguments whatsoever
END SUBROUTINE Multiply
END INTERFACE
REAL, DIMENSION(:), INTENT(IN) :: rhs_vector ! Vector b
REAL, DIMENSION(:), INTENT(INOUT) :: &
solution_vector ! Vector x--contains initial guess on entry
REAL, DIMENSION(:), INTENT(INOUT) :: Vx , Vy ! No aliasing
....
Iterate: DO ! Start iteration in CG
...
IF(converged) EXIT Iterate
...
! This has side effects on Vx and Vy???
! Effectively Vy=A*Vx
CALL Multiply()
...
solution_vector=...Vy
Vx=....
...
END DO Iterate
...
END SUBROUTINE CG_internal
...
END MODULE Conjugate_Gradient
I hope you are still with me, cause there is more!!! Now we want to actually
write the Multiply routine. I will take a very simplified example from my
own work, namely, where the muliplication is an operation associated with a
graph, which is characterized by the number of nodes n, the number of arcs
m, an array of heads and tails heads_tails(2,m) and an array of costs for
the arcs arcs_costs(m). I make the arrays pointers so that they can be
associated with arrays that the user creates in his own ways someplace else,
and use the same trick as above to avoid aliasing optimization problems. One
of the operations related to this graph will be related to our
multiplication, namely CalculatePotentials:
MODULE Directed_Graphs
PRIVATE
...
TYPE, PUBLIC, EXTENSIBLE :: Graph
INTEGER :: n, m, dn ! Number of nodes, arcs and shadow nodes needed
INTEGER, DIMENSION(:,:), POINTER, SAVE :: heads_tails
REAL, DIMENSION(:), POINTER, SAVE :: arcs_costs
...
CONTAINS
PROCEDURE, PASS_OBJ :: CalculatePotentials
...
END TYPE
...
CONTAINS
...
SUBROUTINE CalculatePotentials(Graph, nodes_potentials, arcs_flows,...)
CLASS(Graph), INTENT(IN) :: Graph ! Passed implicitly
REAL, DIMENSION(:), INTENT(OUT) :: nodes_potentials ! This is Vy
REAL, DIMENSION(:), INTENT(IN) :: nodes_flows ! This is Vx
...
! This is an efficient array-based routine
CALL CalculatePotentials_Internal(Graph%n,Graph%m,Graph%dn, &
Graph%heads_tails, Graph%arcs_costs, &
nodes_potentials, nodes_flows)
...
END SUBROUTINE CalculatePotentials
SUBROUTINE CalculatePotentials_Internal(n,m,dn, &
heads_tails, arcs_costs, &
nodes_potentials, nodes_flows)
INTEGER, INTENT(IN) :: n,m,dn
INTEGER, DIMENSION(:,:), INTENT(IN) :: heads_tails ! From Graph
REAL, DIMENSION(:), INTENT(IN) :: arcs_costs ! From Graph
REAL, DIMENSION(:), INTENT(OUT) :: nodes_potentials ! This is Vy
REAL, DIMENSION(:), INTENT(IN) :: nodes_flows ! This is Vx
...
! All intensive and important calculations go here
...
END SUBROUTINE CalculatePotentials_Internal
...
END MODULE Directed_Graphs
OK, are you still with me? Now comes the module that combines these two into
a graph (network) optimization problem, where CG is used to solve some
linear system related to a directed graph:
MODULE Network_Optimization
USE Directed_Graphs
USE Conjugate_Gradient
PRIVATE
...
! Do I need POINTER attributes here somewhere (on the classes)?
TYPE, PUBLIC, EXTENSIBLE :: Network_Problem
CLASS(CG_Solver) :: My_CG_Solver ! The solver
CLASS(Graph) :: My_Graph ! The Graph associated with this problem
! These arrays are Vx and Vy but with the shadows on them as well
REAL, DIMENSION(:), POINTER, SAVE :: nodes_potentials ! This is Vy
REAL, DIMENSION(:), POINTER, SAVE :: nodes_flows ! This is Vx
...
CONTAINS
PROCEDURE, PASS_OBJ :: Multiply
...
END TYPE
...
CONTAINS
...
SUBROUTINE Multiply(Problem)
! The only argument Problem is passed implicitly--exactly what I
want
CLASS(Network_Problem), INTENT(INOUT) :: Problem
....
INTEGER :: n
n=Problem%My_Graph%n
! The CG Vx and Vy should not contain the shadow in them!
Problem%My_CG_Solver%Vy => Problem%nodes_potentials(1:n)
Problem%My_CG_Solver%Vx => Problem%nodes_flows(1:n)
CALL Problem%My_Graph%CalculatePotentials &
(nodes_potentials=Problem%nodes_potentials,
nodes_flows=Problem%nodes_flows,...)
END SUBROUTINE Multiply
...
END MODULE Network_Optimization
Now I am going to ask you for help. Can the routine Network_Problem%Multiply
be passed to the routine CG, where it was declared:
INTERFACE
SUBROUTINE Multiply() ! No arguments whatsoever
END SUBROUTINE Multiply
END INTERFACE
i.e. how will this routine know that one of its arguments is to be passed
implictly. In C++, this is obvious because the procedure Multiply is a
*part* of the class, but as I read the F2K proposal, I am not sure if this
is the case?
Now comes an example caller, here a main program, to illustrate the usage of
this example package:
PROGRAM GiveItAShot
USE Network_Optimization
...
INTEGER, PARAMETER :: n=..., m=... ! Number of nodes and arcs
INTEGER, PARAMETER :: dn=... ! The *shadow* size needed in Multiply
INTEGER, DIMENSION(2,m), TARGET :: heads_tails ! For Graph
REAL, DIMENSION(m), TARGET :: arcs_costs ! For Graph
REAL, DIMENSION(n+dn), TARGET :: &
nodes_excess_potentials, nodes_excess_flows ! Vx and Vy
REAL, DIMENSION(n) :: supplies_demands, nodes_potentials ! For b and x
CLASS(Network_Problem) :: Sample_Problem ! The center of attention
...
! Assign values to heads_tails and arcs_costs, as well as
! supplies_demands and a guess for nodes_potentials
Sample_Problem%My_Graph%heads_tails=>heads_tails
Sample_Problem%My_Graph%arcs_costs=>arcs_costs
! These pointers should not contain the shadow in them
Sample_Problem%nodes_flows=>nodes_excess_flows
Sample_Problem%nodes_potentials=>nodes_excess_potentials
! My nightmare--CAN Sample_Problem%Multiply be passed here?
! Would a pointer procedure help fix the problem?
CALL Sample_Problem%CG(Multiply=Sample_Problem%Multiply, &
rhs_vector=supplies_demands, solution_vector=nodes_potentials ,
...)
END PROGRAM GiveItAShot
So that is my on-the-fly attempt. It took me a while to get through this,
and I am not even sure if this will work. If it does, I like it, and would
appreaciate any input from all of you that got this far.
Thank you very much,
Aleksandar
__________________________________
Aleksandar Donev
http://www.pa.msu.edu/~donev/
[log in to unmask]
(517) 432-6770
Department of Physics and Astronomy
Michigan State University
East Lansing, MI 48824-1116
_____________________________________________
|
