@z This file was created by Aleksandar Donev as part of the Network Optimization project. Feel free to use any portion of it and contact me at donev@pa.msu.edu @x \Title{Near-Neighbour Lattice Networks} \author{Aleksandar Donev} \date{February 2000} @*0 Module |Lattice_Geometry|. This module will set up a lattice with a simple (one-site) basis and incoorporate appropriate boundary conditions. Although it is meant to be primarily used for hypercube lattices, with minor to no modifications the source code is applicable to any $n$-dimensional simple lattice, such as a triangular lattice. The lattice is returned as a general network, i.e. a heads and tails array |heads_tails| (from |Network_Data_Structures|), with appropriate coordinate information for the nodes in the array |nodes_coords|, as well as status identifiers telling what kind of arc/node each one is in the arrays |arcs_status| and |nodes_status|. The routine |InitializeLattice| reads the parameter file and should be called once in the beginning, while the very bulky and complex routine |CreateLatticeNetwork| can be called many times to create lattices with different parameters (such as dilution or sizes of the lattice), so long as it is always called subsequently after a call to |DestroyLatticeNetwork|, which deallocates the arrays that are created in |CreateLatticeNetwork|. I will explain some of the important publicly available module data as I go along. For now, I want to say that after the call to |CreateLatticeNetwork| each arc and node will have a status which is related to their place in the geometry. The possible kinds are defined as one-byte integers below, such as |interior_arc| or |plate_arc| or |injection_node| etc. The arrays |arc_status| and |nodes_status|] are one-byte integer arc and nodal arrays which can contain an integer from 0-255 giving the status (kind) indicator for each node and arc. The module |Lattice_Geometry| defines the status of arcs and nodes based on their geometrical position and role. The module |Network_Geometry| will use this module and base the creation of the arrays in |Network_Data_Structures| upon these geometrical statuses. Basically, a dummy node is labeled as |empty_node|. A node in the interior of the lattice is a |interior_node|, while a node that is on a plate boundary is a |plate_node|. A node can also be connected to $S$ or $T$, as in |to_source_node| or |to_sink_node|, or it can be a |source_node| or a |sink_node|. It may also be an node where flow is injected (|injection_node|), or where flow is taken out (|outflow_node|). The number of injection and outflow nodes is recorded in |n_injection_nodes| and |n_outflow_nodes| because flow must be conserved so these will be needed when assigning supplies/demands later on. {\em Only sources and sinks are treated as special nodes} since it is likely that these will have fixed-potential or other special boundary conditions applied to them. Similarly, an arc can be a |dummy_arc| or an interior |interior_arc|, it can be a wrap-around arc |periodic_bc_arc| (these would not be plotted usually, for example). On the other hand, they can be arcs connecting the source/sink plates or nodes, |plate_arc|, |source_arc| or |sink_arc|. {\em Plate, source and sink arcs are treated as special arcs}, because it is expected these should be zero cost arcs later on. {\em Higher numbers for the status indicators superseed smaller ones}, that is, if a node is both a plate and an injection node, it will be marked as |injection_node|. @ @f dim _dim @f source _source @a MODULE Lattice_Geometry @; USE Precision @; USE Error_Handling @; USE System_Monitors @; USE Random_Numbers @; USE Graph_Algorithms @; USE Network_Data_Structures @; USE Network_Graphics @; IMPLICIT NONE @; PUBLIC :: InitializeLattice, CreateLatticeNetwork, DestroyLatticeNetwork @; PRIVATE @; INTEGER, PARAMETER, PUBLIC :: n_directions = NDIR @; // Number of arcs originating from each lattice site (lattice coordination number). INTEGER, DIMENSION (n_dim), PUBLIC, SAVE :: lengths @; // The sizes of the lattice box CHARACTER, PUBLIC, SAVE :: component_extraction='C' @; // Either |"Cluster"| or |"Backbone"| INTEGER, DIMENSION (n_dim), PUBLIC, SAVE :: & source_layer_thickness=1, sink_layer_thickness=1, & injection_layer_thickness=1, outflow_layer_thickness=1 @; INTEGER, DIMENSION (n_dim, n_directions), PARAMETER :: displacements= & RESHAPE(SOURCE=(/NO_DISPLACEMENTS/), SHAPE=(/n_dim, n_directions/)) @; INTEGER, PARAMETER, PUBLIC :: bc_periodic=1, bc_injection=2, bc_plate=3, bc_source_sink=4 @; LOGICAL, DIMENSION(4,n_dim), PUBLIC, SAVE :: boundaries=.FALSE. @; // Boundary conditions, for each dimension specifyin order whether the boundary is: {\em periodic, injection, plate} or {\em source-sink ($ST$)} boundary. INTEGER(KIND=i_wp), DIMENSION (n_dim), PUBLIC :: sources, sinks @; INTEGER(KIND=i_wp), PUBLIC, SAVE :: n_injection_nodes, n_outflow_nodes @; // Counters REAL, PUBLIC, SAVE :: node_dilution=-1.0, arc_dilution=-1.0 @; // Dilutions for nodes and arcs INTEGER(KIND=i_byte), DIMENSION(:), ALLOCATABLE, PUBLIC :: arcs_status, nodes_status @; // Geometric status of each arc and node. I keep these in bytes to save memory ! INTEGER(KIND=i_byte), PARAMETER, PUBLIC :: empty_node=0, interior_node=1, & plate_node=2, to_source_node=3, to_sink_node=4, & injection_node=5, outflow_node=6, source_node=7, sink_node=8 @; // Possible kinds of nodes for |nodes_status| INTEGER(KIND=i_byte), PARAMETER, PUBLIC :: empty_arc=0, interior_arc=1, periodic_bc_arc=2, & plate_arc=3, source_arc=4, sink_arc=5 @; // Possible kinds of arcs for |arc_status| NAMELIST /LatticeOptions/ lengths, component_extraction, source_layer_thickness, sink_layer_thickness, injection_layer_thickness, outflow_layer_thickness, node_dilution, arc_dilution, boundaries @; // Geometry parameters in options file. CONTAINS @; @@; @@; @@; END MODULE Lattice_Geometry @; @*1 Setting up the Geometry. The main job of the routine |InitializeLattice| is to read the namelist {\tt LatticeOptions} from the options file. The main parameters that can and should be set via this file are: \begin{description} \item[|lengths|] which is an array of size |n_dim| containing the lengths of the lattice box in number of sites \item[|component_extraction|] tells how to extract the actual network from the randomly diluted lattice (shich may be disconnected, for example). Either the largest connected cluster (connected component) can be extracted, |"Cluster"|, or the largest backbone (multi-cycle) if |"Backbone"|. This backbone cluster does not have dangling ends and is made up entirely of cycles. \item[|node_dilution| and |arc_dilution|] which are default reals and give the desired dilutions for the nodes and/or arcs in the graph. For example, a |node_dilution=0.75| means that $25$ ($100-75$) percent of the nodes will be removed in the initial phase (during the creation of the lattice). More nodes will likely be diluted when a single connected-component is extracted from the graph. {\em Negative dilutions imply no dilution at all.} \item[|boundaries|] which is a logical array of size |(4,n_dim)| where there is one column for each dimension and each column contains four indicators that say, in order, whether the boundary along the dimension is: {\em periodic, injection, plate} or {\em source-sink} boundary (these are given by |bc_periodic=1, bc_injection=2, bc_plate=3, bc_source_sink=4| above). \\ So a first column of |F,T,T,F| in |boundaries| says that the boundary is free (no periodic BC's), that the nodes on the boundaries |x=1| will be injection sites and connected in a plate with special arcs, while the nodes with |x=lengths[1]| will be outflow sites and on a plate, and that there will be no source-sink pair in the |x|-direction. \item[|sources| and |sinks|] are two pointer-like arrays that store the locations of the sources and sinks in each dimension. These will be needed when setting up the supplies and demands for the special nodes, since one source/sink can only belong to one dimension. Notice that an injection/outflow node can belong to more than one dimension (such as corner or edge nodes), and thereofore I do not attempt to distinguish these by dimension, but rather count them all together in |n_injection_nodes| and |n_outflow_nodes|. \item[|source_layer_thickness|, |sink_layer_thickness|, |injection_layer_thickness| and |outflow_layer_thickness|] give the thickness (in terms of lattice coordinates) of the layers to which the source and sink are connected and the thicknesses of the layers of injection/outflow nodes, {\em for each dimension} (just use 0 or 1 for dimensions which do not have an $ST$ or injection/outflow boundary). Note that these can be zero, which implies that there is no sink/source or injection/outflow sites on that side of the lattice box. This allows a great variety of geometries to be simulated with this program. Notice that plate boundaries always have thickness 1 since I see no point in having a thicker plate. \end{description} Two important geometry parameters have been left out as |PARAMETER|s to be inputed via FWEB's macros, and can not be changed at run time. These are the number of dimensions |n_dim| and the type of lattice (square nearest-neighbor, second nearest-neighbor, eight-neighbor, triangular, etc.), which is determined by the |n_directions| |n_dim|-dimensional vectors stored as rows in the array |displacements|. These displacements give the displacements (in lattice coordinates) of the neighbors of each node. For example, in a nearest neighbour hypercubic lattice, these are just the unit vectors $(1,0,...)$, $(0,1,...)$, etc., but for other lattices these may be more complicated. The choice to make these fixed parameters was in order to facilitate compiler optimization of the short loops they create and avoid a lot of small allocatable arrays. @ @= SUBROUTINE InitializeLattice () @; IMPLICIT NONE @; REWIND(UNIT=program_options_unit) @; // Rewind the options file READ(UNIT=program_options_unit, NML=LatticeOptions, IOSTAT=error_status) @; IF(error_status!=0) @~ & CALL CriticalError(message="NAMELIST LatticeOptions was not read successfully"\/ & " from file "\/TRIM(options_file), caller="InitializeLattice") @; WRITE(UNIT=message_log_unit, NML=LatticeOptions) @; END SUBROUTINE InitializeLattice @; @*1 Creating the Lattice. To create a lattice as a network, we need to assign the heads-and-tails arrays for the network, along with the correct number of nodes |n_nodes| and arcs |n_arcs|, and in our case assign coordinates to the nodes in |nodes_coordinates|. One has to choose a default numbering for the nodes, and in my case this is done by numbering the points (sites) of the lattice in the $x$ direction first, then in $y$, etc. The following two macros do the conversion from coordinates to site number and vice-versa: @ @m _CoordsToPoint(point_coords, point_num, dim) @; point_num = 1 @; DO dim = 1, n_dim @; point_num = point_num + (point_coords[dim]-1) * lengths_product[dim] @; END DO @; @m _PointToCoords(point_num, point_coords, dim, point) @; point = point_num-1 @; DO dim = n_dim, 1, -1 @; point_coords [dim] = point/lengths_product[dim] + 1 @; // Integer coefficient point = point - (point_coords[dim]-1) * lengths_product[dim] @; // Integer remainder--|MOD| operation inlined END DO @; @*2 Procedure |CreateLatticeNetwork|. Here is the routine |CreateLatticeNetwork|. I can not comment too much on it since it is rather involved and coded according to my own standards. I also partially had future parallelization with HPF or OpenMP in mind, but this is a different issue. The routine will first allocate enough space (meaning full lattice size) for all needed arrays, then it will calculate the connectivity (heads-tails array) for the lattice and the coordinates of the nodes, all the while carefully considering the boundary conditions and marking each node and arc with the correct status. The number of nodes and arcs after this phase will correspond closely to the desired dilution. Then, the routine calls a connected-components-labeling algorithm from |Graph_Algorithms| to find and label the largest connected-component in the graph. Finally, it will extract this component, throwing away all arcs and nodes which do not belong to it (suitably renumbering them). If there is no dilution then the rather costly extraction step is skipped. @ @= SUBROUTINE CreateLatticeNetwork () @; IMPLICIT NONE @; INTEGER(KIND=i_wp), DIMENSION(:), ALLOCATABLE :: points_count, nodes_labels, nodes_count @; // Temporaries /* We need some temporary buffers for reallocating some of the arrays later on: */ INTEGER(KIND=i_wp), DIMENSION(:,:), ALLOCATABLE :: heads_tails_buffer @; REAL, DIMENSION(:,:), ALLOCATABLE :: coords_buffer @; INTEGER(KIND=i_byte), DIMENSION(:), ALLOCATABLE :: status_buffer @; INTEGER(KIND=i_wp) :: n_points, n_links, n_special_points, n_special_links, & n_selected_nodes, n_selected_arcs, n_selected_special_arcs @; // Counters of sites and edges INTEGER(KIND=i_wp) :: link, point, node, arc, neighbor_point, node_counter, arc_counter, & point_counter, point_temp, selected_label, selected_size @; INTEGER :: dim, direction, dim_temp @; // Temporaries and loop counters INTEGER(KIND=i_byte) :: arc_status, node_status @; // Temporaries INTEGER(KIND=i_wp), DIMENSION (n_dim) :: lengths_product INTEGER, DIMENSION (n_dim) :: neighbor_coords, point_coords @; // Temporaries REAL :: dice, real_temp @; // Random dice for dilution LOGICAL :: dilute_nodes, dilute_arcs, extraction_success, & node_on_plate, special_arc, renumber_nodes @; // Indicators @@; @@; @@; END SUBROUTINE CreateLatticeNetwork @; @*2 Initial Allocation or Network Arrays. The following piece of code finds the full (maximal) size of the lattice as well as the maximal number of special arcs and nodes, and then allocates the network arrays with these full sizes. It would have been nice if all reallocations were avoided in the case of a fully connected network. But counting the number of arcs that will be present under the plethora of boundary conditionsand general lattices I allow turned out to be a formiddable task, especially because of the arcs that cross the box boundaries (which can, but do not have to be, wrapped around in a periodic fashion). The following piece of code was placed below and caused a lot of debugging trouble because of double-counting of certain edge and vertex arcs: @ @= DO dim=1, n_dim @; // Check boundaries for crossing arcs PeriodicArcs: IF(.NOT.boundaries[bc_periodic,dim]) THEN @; // Plate boundary DO direction=1, n_directions @; IF(displacements[dim, direction]!=0) THEN @; // An arc that crosses the box! n_links=n_links-n_points/lengths[dim] @; // These are absent END IF @; END DO @; END IF PeriodicArcs @; END DO @; @ Here is the code which works and is simple, but does require reallocation in certain cases: @= lengths_product[1]=1 @; // $lengths\_product_k=\prod_{i=0}^k {lengths_i}$ DO dim = 2, n_dim @; // |lengths_product| is needed when numbering the nodes lengths_product[dim] = lengths_product[dim-1]*lengths[dim-1] @; END DO @; n_points = PRODUCT (lengths) @; // Total number of sites in the lattice n_links = n_points * n_directions @; // Total number of connections sources=0 ;@~ sinks=0 ;@; // The node numbers for the sources and sinks--needed later on /* Boundary conditions have to be scanned for possible plate or source/sink boundary conditions, which require special arcs and nodes:*/ n_special_points=0 ;@~ n_special_links=0 ;@; DO dim=1, n_dim @; // Check boundaries for all dimensions PlateArcs: IF(boundaries[bc_plate,dim]) THEN @; // Plate boundary DO direction=1, n_directions @; IF(displacements[dim, direction]==0) THEN @; // An arc in plate! n_special_links = n_special_links + 2*n_points/lengths[dim] @; // Zero-cost arcs n_links=n_links-2*n_points/lengths[dim] @; // {\em Possible} special arcs END IF @; END DO @; END IF PlateArcs @; SourcesSinks: IF(boundaries[bc_source_sink,dim]) THEN @; // Add a source IF(source_layer_thickness[dim] > 0) THEN @; n_special_links = n_special_links + & source_layer_thickness[dim]* n_points / lengths[dim] @; // Arcs from $S$ n_special_points=n_special_points + 1 @; // $S$ is a special nodes END IF @; IF(sink_layer_thickness[dim] > 0) THEN @; // Add a sink n_special_links = n_special_links + & sink_layer_thickness[dim]* n_points / lengths[dim] @; // Arcs from $S$ n_special_points=n_special_points + 1 @; // $S$ is a special nodes END IF @; END IF SourcesSinks @; END DO @; /* Now we need to allocate all the arrays with their maximal possible sizes: */ ALLOCATE(nodes_coords(n_dim,-n_special_points : n_points), STAT=error_status) @; CALL RecordAllocation(n_elements=n_dim*(n_points+n_special_points+1), mold=1.0, & caller="CreateLatticeNetwork", alloc_status=error_status) @; nodes_coords(1:n_dim,0)=0.0 @; // Dummy node ALLOCATE(nodes_status(-n_special_points : n_points), STAT=error_status) @; CALL RecordAllocation(n_elements=(n_points+n_special_points+1), mold=1_i_byte, & caller="CreateLatticeNetwork", alloc_status=error_status) @; ALLOCATE(heads_tails(2,-n_special_links : n_links), STAT=error_status) @; CALL RecordAllocation(n_elements=2*(n_links+n_special_links+1), mold=1_i_wp, & caller="CreateLatticeNetwork", alloc_status=error_status) @; heads_tails(1:2,0)=0 @; // Initialize zeroth arc to be loop to the dummy node ALLOCATE(arcs_status(-n_special_links : n_links), STAT=error_status) @; CALL RecordAllocation(n_elements=(n_links+n_special_links+1), mold=1_i_byte, & caller="CreateLatticeNetwork", alloc_status=error_status) @; nodes_status=empty_node ;@~ arcs_status=empty_arc @; // Initialize dummy ones @*2 Lattice Geometry Calculation. The connectivity of the network is calculated in an involved way because we allow for both node and arc dilution. The major loop is over nodes, and for each node I calculate its neighbors and assign the proper heads-and-tails and nodal coordinates. This piece of code is a monster when it comes to parallelization due to the constant use of indirections: @ @= dilute_nodes= ((node_dilution<1.0) .AND. (node_dilution>0.0)) @; // Dilution is in $(0,1)$ dilute_arcs= ((arc_dilution<1.0) .AND. (arc_dilution>0.0)) @; // Dilution is in $(0,1)$ /* When nodes are diluted, their numbering changes from the default because some of them will be absent (diluted). We need to make a mask and so a prefix counting operation to find the new numbering of the nodes, and store this for later use: */ IF(dilute_nodes) THEN @; // |points_count=COUNT_PREFIX(points_mask)| in essence ALLOCATE(points_count(n_points), STAT=error_status) @; CALL RecordAllocation(n_elements=n_points, mold=1_i_wp, & caller="CreateLatticeNetwork", alloc_status=error_status) @; points_count=0 @; // Zero means a diluted--empty node point_counter=0 @; // Just a counter for |points_count| DiluteNodes: DO point=1, n_points @; // See which nodes are diluted node_on_plate=.FALSE. @; _PointToCoords(point, point_coords, dim_temp, point_temp) @; // Get the coordinates PlateNodes: DO dim=1, n_dim @; // Is this a plate node? IF(boundaries[bc_plate,dim]) THEN @; // Is this point on a plate? IF(point_coords[dim]==1 .OR. point_coords[dim]==lengths[dim]) & node_on_plate=.TRUE. @; // Plate nodes are not diluted END IF @; END DO PlateNodes @; IF(.NOT.node_on_plate) THEN @; // Dilute the node if dice says so CALL RandomUniform(dice) @; // In $[0,1]$ IF (dice>node_dilution) @~ CYCLE DiluteNodes @; // Throw away the node END IF @; point_counter=point_counter+1 @; points_count[point] = point_counter @; END DO DiluteNodes @; END IF @; n_nodes=0 ;@~ n_arcs=0 ;@~ n_special_nodes=0 ;@~ n_special_arcs=0 ;@; // Initialize SourceSinkNodes: DO dim=1, n_dim @; // Check the boundary conditions for special nodes IF(boundaries[bc_source_sink,dim]) THEN @; //Initialize the source and sink IF(source_layer_thickness[dim] > 0) THEN @; // Add $S$ n_special_nodes=n_special_nodes+1 @; // $S$ node sources[dim]=-n_special_nodes @; // Save this location for later use nodes_status[-n_special_nodes]=source_node @; // This is a source node nodes_coords[1:n_dim,-n_special_nodes]=0.5*@E REAL(lengths+1) @; // An attempt to plot $ST$ nodes nicely. nodes_coords[dim,-n_special_nodes]=1.0-0.1 * @E REAL(lengths[dim]) @; // Place the sink or source 10 percent to the side of the boundary END IF @; IF(sink_layer_thickness[dim] > 0) THEN @; // Add $T$ n_special_nodes=n_special_nodes+1 @; // $T$ node sinks[dim]=-n_special_nodes @; // Store for later use nodes_status[-n_special_nodes]=sink_node @; // This is a sink node nodes_coords[1:n_dim,-n_special_nodes]=0.5*@E REAL(lengths+1) @; nodes_coords[dim,-n_special_nodes]=1.1*@E REAL(lengths[dim]) @; // For plotting END IF @; END IF @; END DO SourceSinkNodes @; /* Now I start the major loop over nodes and add the arcs to the neighbors one by one: */ ForAllPoints: DO point=1, n_points @; // |FORALL| masked with |points_mask| node_status=interior_node @; IF(dilute_nodes) THEN @; // Check dice mask for dilution IF(points_count[point]==0) @~ CYCLE ForAllPoints @; // This site is diluted n_nodes=points_count[point] @; // The new node number ELSE @; // Go on n_nodes=point @; // No need to renumber nodes END IF @; _PointToCoords(point, point_coords, dim_temp, point_temp) @; // Get the coordinates nodes_coords[1:n_dim,n_nodes]=@E REAL(point_coords) @; // Conversion--could be more complicated if the topology is more complicated, such as a triangular lattice PlateOrInjectionNodes: DO dim=1, n_dim @; // Check if this is a plate or injection node IF(boundaries[bc_plate,dim]) THEN @; // Check for plate node IF(point_coords[dim]==1 .OR. point_coords[dim]==lengths[dim]) THEN @; node_status=MAX(node_status, plate_node) @; // Plate node END IF @; END IF @; IF(boundaries[bc_injection,dim]) THEN @; // Is this node an injection site? IF(point_coords[dim] <= injection_layer_thickness[dim] ) THEN @; node_status=MAX(node_status, injection_node) @; // Injection site ELSE IF (point_coords[dim] > lengths[dim]-outflow_layer_thickness[dim]) THEN @; node_status=MAX(node_status, outflow_node) @; // Outflow site END IF @; END IF @; IF(boundaries[bc_source_sink,dim]) THEN @; // Check for connections to $ST$ nodes IF(point_coords[dim] <= source_layer_thickness[dim] ) THEN @; // Connect to $S$--here all source arcs start at $S$ and all sink arcs end at $T$ by convention node_status=MAX(node_status, to_source_node) @; n_special_arcs=n_special_arcs+1 @; arcs_status[-n_special_arcs]=MAX(arcs_status[-n_special_arcs],source_arc) @; IF(dilute_nodes) THEN @; heads_tails[1:2,-n_special_arcs]=(/sources[dim],points_count[point]/) @; ELSE @; heads_tails[1:2,-n_special_arcs]=(/sources[dim],point/) @; END IF @; ELSE IF(point_coords[dim] > lengths[dim]-sink_layer_thickness[dim]) THEN @; // Connect to $T$ node_status=MAX(node_status, to_sink_node) @; n_special_arcs=n_special_arcs+1 @; arcs_status[-n_special_arcs]=MAX(arcs_status[-n_special_arcs],sink_arc) @; IF(dilute_nodes) THEN @; heads_tails[1:2,-n_special_arcs]=(/points_count[point],sinks[dim]/) @; ELSE @; heads_tails[1:2,-n_special_arcs]=(/point,sinks[dim]/) @; END IF @; END IF @; END IF @; END DO PlateOrInjectionNodes @; nodes_status[n_nodes]=node_status @; // Kind of node /* Now I start adding the arcs that originate at this node (have the node as their head) one by one by looking at all the directions specified in |displacements|: */ ForAllDirections: DO direction = 1, n_directions @; // Different directions arc_status=interior_arc @; special_arc=.FALSE. @; // An ordinary arc neighbor_coords=point_coords+displacements[:,direction] @; // The neighbour coordinates Plates: DO dim=1, n_dim @; // Check for special arcs BoundaryPlate: IF(boundaries[bc_plate,dim]) THEN @; // On the plate? IF(((point_coords[dim]==lengths[dim]) .AND. (neighbor_coords[dim]==lengths[dim])) & @\ .OR. ((point_coords[dim]==1) .AND. (neighbor_coords[dim]==1))) THEN @; special_arc=.TRUE. @; // A plate arc arc_status=MAX(plate_arc, arc_status) @; // I will not dilute plate arcs END IF @; END IF BoundaryPlate @; END DO Plates @; IF(.NOT.special_arc.AND.dilute_arcs) THEN @; // Throw dilution dice for arcs CALL RandomUniform(dice) @; // In $[0,1]$ IF(dice>arc_dilution) @~ CYCLE ForAllDirections @; // This arc is diluted END IF @; // Go on Periodics: DO dim=1, n_dim @; // Check if crossing a periodic boundary BoundaryCrossing: IF((neighbor_coords[dim]>lengths[dim]) & .OR. (neighbor_coords[dim]<1) ) THEN @; // Outside of square IF(boundaries[bc_periodic,dim]) THEN @; IF(neighbor_coords[dim]>lengths[dim]) & @\ neighbor_coords[dim]=MOD(neighbor_coords[dim],lengths[dim]) @; // Wrap around the lattice in this dimension IF(neighbor_coords[dim]<1) & @\ neighbor_coords[dim]=lengths[dim]+MOD(neighbor_coords[dim],lengths[dim]) @; // Wrap around arc_status=MAX(periodic_bc_arc, arc_status) @; // A wrap-around arc ELSE @; CYCLE ForAllDirections @; // This arc is not allowed END IF @; END IF BoundaryCrossing @; END DO Periodics @; _CoordsToPoint(neighbor_coords, neighbor_point, dim_temp) @; IF(dilute_nodes) THEN @; // Nodes are diluted--there is danger of empty tail! Also, the nodes were renumbered in |points_count| IF(points_count[neighbor_point]==0) @~ CYCLE ForAllDirections @; // The tail is empty IF(.NOT.special_arc) THEN @; // A regular arc n_arcs=n_arcs+1 @; heads_tails[1:2,n_arcs]=(/points_count[point], points_count[neighbor_point]/) @; arcs_status[n_arcs]=arc_status @; ELSE @; // A special arc n_special_arcs=n_special_arcs+1 @; heads_tails[1:2,-n_special_arcs]=(/points_count[point], points_count[neighbor_point]/) @; arcs_status[-n_special_arcs]=arc_status @; END IF @; ELSE @; // There is no node dilution, danger of empty tail or node renumbering IF(.NOT.special_arc) THEN @; n_arcs=n_arcs+1 @; heads_tails[1:2,n_arcs]=(/point, neighbor_point/) @; arcs_status[n_arcs]=arc_status @; ELSE @; n_special_arcs=n_special_arcs+1 @; heads_tails[1:2,-n_special_arcs]=(/point, neighbor_point/) @; arcs_status[-n_special_arcs]=arc_status @; END IF @; END IF @; END DO ForAllDirections @; END DO ForAllPoints @; IF(dilute_nodes) THEN @; // Free the prefix count array CALL RecordAllocation(n_elements=-_SIZE(points_count, i_wp), & mold=1_i_wp, caller="CreateLatticeNetwork") @; DEALLOCATE(points_count) @; END IF @; @*2 Connected-Component Extraction. The following code section calls the routine |ConnectedComponents| to find the largest connected component and extracts it. This is skipped if there was no dilution. Then the code recreates (reallocates) the lattice network with just the largest connected component extracted from the original lattice. Evem when there is no dilution for certain combinations of boundary conditions the network arrays will need to be reallocated: @ @= n_injection_nodes=0 @; // We have to count these because of the constraint $\sum{b_i}=0$ n_outflow_nodes=0 @; ExtractSelectedComponent: IF(dilute_nodes .OR. dilute_arcs) THEN @; // Extract a single connected component WRITE(message_log_unit,*) "Starting connected-component extraction with:" WRITE(message_log_unit,*) "n_arcs=", n_arcs, " n_nodes=", n_nodes, & " n_special_nodes=", n_special_nodes, " n_special_arcs=", n_special_arcs @; ALLOCATE(nodes_mask(-n_special_nodes:n_nodes), STAT=error_status) @; CALL RecordAllocation(n_elements=n_special_nodes+n_nodes+1, mold=.TRUE._l_wp, & caller="CreateNetwork", alloc_status=error_status) @; ALLOCATE(arcs_mask(-n_special_arcs:n_arcs)) @; CALL RecordAllocation(n_elements=n_special_arcs+n_arcs+1, mold=.TRUE._l_wp, & caller="CreateNetwork", alloc_status=error_status) @; ALLOCATE(nodes_labels(-n_special_nodes:n_nodes), STAT=error_status) @; // The connected-components labels for the nodes CALL RecordAllocation(n_elements=n_special_nodes+n_nodes+1, mold=1_i_wp, & caller="CreateLatticeNetwork", alloc_status=error_status) @; nodes_labels(0)=0 @; // Dummy node SELECT CASE(component_extraction) @; CASE('b','B') @; CALL BackboneCycles(heads_tails=heads_tails(:,-n_special_arcs:n_arcs), & node_offset=n_special_nodes, supernodes=nodes_labels, & largest_supernode=selected_label, & largest_supernode_size=selected_size) @; // Label the backbone cycles WRITE(message_log_unit,*) "The largest backbone component has: ", & selected_size, " nodes" @; CASE DEFAULT @; CALL ConnectedComponents(heads_tails=heads_tails(:,-n_special_arcs:n_arcs), & node_offset=n_special_nodes, labels=nodes_labels, & largest_component_label=selected_label, & largest_component_size=selected_size) @; // Label the connected components WRITE(message_log_unit,*) "The largest connected component has: ", & selected_size, " nodes" @; ENDSELECT @; n_selected_nodes=0 @; // The number of nodes in the extracted component. All special nodes and arcs must be in the selected connected component, so I don't count them explicitly n_selected_arcs=0 @; // The number of arcs n_selected_special_arcs=0 @; // Number of special arcs in the extracted component extraction_success=.TRUE. @; /* Now we mask as |.TRUE.| all the nodes and arcs that belong to the largest component. I also check whether this componet contains all the special nodes (sources and sinks and plate nodes in it), and whether there is at least one injection and outflow site in it (if these were requested in |boundaries|). If not, this means the lattice is too diluted and the program is terminated with a critical error: */ MaskSelectedNodes: DO node=-n_special_nodes, n_nodes @; // Select only the nodes in the major component SelectNodes: IF(nodes_labels[node]==selected_label) THEN @; // Select this node nodes_mask[node] = .TRUE._l_wp @; IF(node>0) @~ n_selected_nodes=n_selected_nodes+1 @; IF(nodes_status[node]==injection_node) THEN @; // Count the injection/outflow nodes n_injection_nodes=n_injection_nodes+1 @; ELSE IF(nodes_status[node]==outflow_node) THEN @; // We want to count these n_outflow_nodes=n_outflow_nodes+1 @; END IF @; ELSE @; // If this was a special or plate node, then abort the mission IF( (node<0) .OR. (nodes_status[node]==plate_node) ) THEN @; extraction_success =.FALSE. @; EXIT MaskSelectedNodes @; // This is a critical error! END IF @; nodes_mask[node]=.FALSE._l_wp @; END IF SelectNodes @; END DO MaskSelectedNodes @; CheckConnectedness: IF(.NOT.extraction_success) THEN @; // The graph is too diluted CALL CriticalError(message="The sources/sinks and/or plate nodes were not in the largest connected component", caller="CreateLatticeNetwork") @; ELSE @; DO dim=1, n_dim @; IF(boundaries[bc_injection,dim]) THEN @; // There were injection boundaries IF((injection_layer_thickness[dim]>0) .AND. (n_injection_nodes==0)) & CALL CriticalError(message="No injection sites in the largest component", caller="CreateLatticeNetwork") @; IF((outflow_layer_thickness[dim]>0) .AND. (n_outflow_nodes==0)) & CALL CriticalError(message="No outflow sites in the largest component", caller="CreateLatticeNetwork") @; END IF @; END DO @; @%% ELSE IF( @E REAL(ABS(n_injection_nodes-n_outflow_nodes)) > & @%% 0.1*@E REAL(n_injection_nodes+n_outflow_nodes) ) THEN @; // Imbalance @%% CALL Warning(message="An imbalance in the number of injection and outflow sites larger than 20 percent", caller="CreateLatticeNetwork") @; @%% WRITE(warning_print_unit,*) "n_injection_nodes=", n_injection_nodes, & @%% "n_outflow_nodes=", n_outflow_nodes @; @%% WRITE(warning_log_unit,*) "n_injection_nodes=", n_injection_nodes, & @%% "n_outflow_nodes=", n_outflow_nodes @; END IF CheckConnectedness @; CALL RecordAllocation(n_elements=-_SIZE(nodes_labels, i_wp), & mold=1_i_wp, caller="CreateLatticeNetwork") @; DEALLOCATE(nodes_labels) @; // We don't need these any more MaskSelectedSpecialArcs: DO arc=-n_special_arcs, -1 @; // Mask empty arcs IF(nodes_mask(heads_tails[1,arc]).AND.nodes_mask(heads_tails[2,arc])) THEN @; n_selected_special_arcs=n_selected_special_arcs+1 @; arcs_mask[arc]=.TRUE._l_wp @; // Arc is in the right component ELSE @; arcs_mask[arc]=.FALSE._l_wp @; // Throw away this arc END IF @; END DO MaskSelectedSpecialArcs @; MaskSelectedArcs: DO arc=1, n_arcs @; // Mask empty arcs IF(nodes_mask(heads_tails[1,arc]).AND.nodes_mask(heads_tails[2,arc])) THEN @; n_selected_arcs=n_selected_arcs+1 @; arcs_mask[arc]=.TRUE._l_wp @; // Arc is in the right component ELSE @; arcs_mask[arc]=.FALSE._l_wp @; // Throw away this arc END IF @; END DO MaskSelectedArcs @; /* Now I check whether reallocation, packing and node renumbering are really neccessary and then perform them if they are. This section of code is really too complicated for what it is worth, since the cost of the reallocation is not that great: */ renumber_nodes=.FALSE. @; // Are nodes renumbered? IF((n_selected_nodes!=n_points).OR.(n_special_nodes!=n_special_points)) THEN @; // Pack WRITE(message_log_unit,*) "Reallocating node arrays..." @@; IF(n_selected_nodes!=n_nodes) THEN @; // Repacking/renumbering of nodes is needed renumber_nodes=.TRUE. @; WRITE(message_log_unit,*) "Repacking node arrays..." ALLOCATE(nodes_count(-n_special_nodes:n_nodes), STAT=error_status) @; CALL RecordAllocation(n_elements=n_special_nodes+n_nodes+1, mold=1_i_wp, & caller="CreateLatticeNetwork", alloc_status=error_status) @; @@; @@; ELSE @; // Rather unlikely! @@; END IF @; @@; END IF @; CALL RecordAllocation(n_elements=-_SIZE(nodes_mask, i_wp), & mold=.TRUE._l_wp, caller="CreateLatticeNetwork") @; DEALLOCATE(nodes_mask) @; IF((n_selected_arcs!=n_links).OR.(n_selected_special_arcs!=n_special_links)) THEN @; // We need to repackage the arc arrays into smaller arrays WRITE(message_log_unit,*) "Reallocating arc arrays..." @@; IF((n_selected_arcs!=n_arcs).OR.(n_selected_special_arcs!=n_special_arcs)) THEN @; WRITE(message_log_unit,*) "Repacking arc arrays..." IF(renumber_nodes) THEN @; @@; ELSE @; @@; END IF @; ELSE @; // Very unlikely! IF(renumber_nodes) THEN @; @@; ELSE @; @@; END IF @; END IF @; @@; END IF @; IF(renumber_nodes) THEN @; // Deallocate the prefix-count array CALL RecordAllocation(n_elements=-_SIZE(nodes_count, i_wp), & mold=1_i_wp, caller="CreateLatticeNetwork") @; DEALLOCATE(nodes_count) @; // We don't need these any more END IF @; CALL RecordAllocation(n_elements=-_SIZE(arcs_mask, i_wp), & mold=.TRUE._l_wp, caller="CreateLatticeNetwork") @; DEALLOCATE(arcs_mask) @; n_nodes=n_selected_nodes @; // Update $n$ n_arcs=n_selected_arcs @; // Update $m$ n_special_arcs=n_selected_special_arcs @; ELSE @; // If there is no dilution there is no need for connected-components extraction n_selected_nodes=n_nodes @; // All nodes are connected n_selected_arcs=n_arcs @; // All arcs are in the selected component n_selected_special_arcs=n_special_arcs @; IF((n_nodes!=n_points).OR.(n_special_nodes!=n_special_points)) THEN @; // Nodes WRITE(message_log_unit,*) "Reallocating node arrays..." @@; @@; @@; END IF @; IF((n_arcs!=n_links).OR.(n_special_arcs!=n_special_links)) THEN @; // Reallocate for arcs WRITE(message_log_unit,*) "Reallocating arc arrays..." @@; @@; @@; END IF @; IF( ANY(boundaries[bc_injection,:]) ) THEN @; // Count the number of in/out sites DO node=1, n_nodes @; // Explicit counting to avoid double-counting IF(nodes_status[node]==injection_node) THEN @; // Count the injection/outflow nodes n_injection_nodes=n_injection_nodes+1 @; ELSE IF(nodes_status[node]==outflow_node) THEN @; // We want to count these n_outflow_nodes=n_outflow_nodes+1 @; END IF @; END DO @; END IF @; END IF ExtractSelectedComponent @; IF( (n_injection_nodes!=0) .OR. (n_outflow_nodes!=0) ) & WRITE(message_log_unit,*) "There are ", n_injection_nodes," injection sites, and ", & n_outflow_nodes, " outflow sites in the network." @*2 Reallocating and Packing the Network Arrays. Now that we know which nodes and arcs to keep and which to throw away, we need to actually pack the network arrays and reallocate them with the correct sizes. Since Fortran does not allow direct reallocation, this piece here is rather messy and long, and temporary reallocation buffers must be used at the cost of one extra memory transfer. Care must also be taken because nodes are again renumbered because some will be thrown away, and the head and tails arrays must be updated to reflect this. In HPF, the routine |PACK| or an analog will be used to do this, at a rather large cost of memory movement and lots of indirections. @ This code will allocate the copy buffers required for packing the arc-arrays: @= ALLOCATE(heads_tails_buffer(2,-n_selected_special_arcs:n_selected_arcs), & STAT=error_status) @; CALL RecordAllocation(n_elements=2*(n_selected_arcs+n_selected_special_arcs+1), & mold=1_i_wp, caller="CreateLatticeNetwork", alloc_status=error_status) @; ALLOCATE(status_buffer(-n_selected_special_arcs:n_selected_arcs), STAT=error_status) @; CALL RecordAllocation(n_elements=n_selected_arcs+n_selected_special_arcs+1, mold=1_i_byte, & caller="CreateLatticeNetwork", alloc_status=error_status) @; heads_tails_buffer[1:2,0]=0 @; // Null element status_buffer[0]=empty_arc @; // Dummy arc @ Once the buffers have been assigned values (see macros below), the network arrays need to be reallocated and the buffers copied into the new arrays: @= CALL RecordAllocation(n_elements=-_SIZE(heads_tails, i_wp), & mold=1_i_wp, caller="CreateLatticeNetwork") @; DEALLOCATE(heads_tails) @; // Start the re-allocation process ALLOCATE(heads_tails(2,-n_selected_special_arcs:n_selected_arcs), & STAT=error_status) @; // New size CALL RecordAllocation(n_elements=2*(n_selected_arcs+n_selected_special_arcs+1), & mold=1_i_wp, caller="CreateLatticeNetwork", alloc_status=error_status) @; heads_tails=heads_tails_buffer @; // Copy the buffer back into the array CALL RecordAllocation(n_elements=-_SIZE(heads_tails_buffer, i_wp), & mold=1_i_wp, caller="CreateLatticeNetwork") @; DEALLOCATE(heads_tails_buffer) @; CALL RecordAllocation(n_elements=-_SIZE(arcs_status, i_wp), & mold=1_i_byte, caller="CreateLatticeNetwork") @; DEALLOCATE(arcs_status) @; // Start the re-allocation process ALLOCATE(arcs_status(-n_selected_special_arcs:n_selected_arcs), & STAT=error_status) @; // New size CALL RecordAllocation(n_elements=n_selected_arcs+n_selected_special_arcs+1, mold=1_i_byte, & caller="CreateLatticeNetwork", alloc_status=error_status) @; arcs_status=status_buffer @; // Copy the buffer CALL RecordAllocation(n_elements=-_SIZE(status_buffer, i_wp), & mold=1_i_byte, caller="CreateLatticeNetwork") @; DEALLOCATE(status_buffer) @; @ When nodes are packed, they are renumbered, according to |nodes_labels=COUNT_PREFIX(nodes_mask)|. Here is this operation inlined: @= nodes_count[0]=0 @; DO node=-n_special_nodes, -1 @; nodes_count[node]=node @; // Since they are all present END DO @; node_counter=0 @; DO node=1, n_nodes @; // |COUNT_PREFIX| IF(nodes_mask[node]) THEN @; // Count this node node_counter=node_counter+1 @; nodes_count[node]=node_counter @; END IF @; // No need to count the rest of them END DO @; @ How the arc or node arrays are packed into the buffer depends on whether there is dilution or not and whether nodes were renumbered (which has to be reflected in the new heads and tails arrays). Here is the piece of code that will assign values to the buffers when there is dilution @= arc_counter=0 @; DO arc=-1, -n_special_arcs, -1 @; // Sometimes we need to pack these as well IF(arcs_mask[arc]) THEN @; // Pack this arc arc_counter=arc_counter-1 @; heads_tails_buffer[1:2,arc_counter]=nodes_count(heads_tails[1:2,arc]) @; // Renumber status_buffer[arc_counter]=arcs_status[arc] @; END IF @; END DO @; arc_counter=0 @; PackArcsRenumbered: DO arc=1, n_arcs IF(arcs_mask[arc]) THEN @; // Pack this arc arc_counter=arc_counter+1 @; heads_tails_buffer[1:2,arc_counter]=nodes_count(heads_tails[1:2,arc]) @; status_buffer[arc_counter]=arcs_status[arc] @; END IF @; END DO PackArcsRenumbered @; @ When nodes are not renumbered packing is a bit easier: @= arc_counter=0 @; DO arc=-1, -n_special_arcs, -1 @; // Sometimes we need to pack these as well IF(arcs_mask[arc]) THEN @; // Pack this arc arc_counter=arc_counter-1 @; heads_tails_buffer[1:2,arc_counter]=heads_tails[1:2,arc] @; // Renumber status_buffer[arc_counter]=arcs_status[arc] @; END IF @; END DO @; arc_counter=0 @; PackArcs: DO arc=1, n_arcs IF(arcs_mask[arc]) THEN @; // Pack this arc arc_counter=arc_counter+1 @; heads_tails_buffer[1:2,arc_counter]=heads_tails[1:2,arc] @; status_buffer[arc_counter]=arcs_status[arc] @; END IF @; END DO PackArcs @; @ When there is no dilution, there is no need to perform a connected-components extraction and to pack the arrays, but simply reallocate them by copying the old contents to the buffers, deallocating, again allocating and recopying the buffer back into the arrays: @= DO arc=-n_selected_special_arcs, n_selected_arcs @; heads_tails_buffer[1:2,arc]=nodes_count(heads_tails[1:2,arc]) @; END DO @; status_buffer=arcs_status[-n_selected_special_arcs:n_selected_arcs] @; @ @= heads_tails_buffer=heads_tails[1:2,-n_selected_special_arcs:n_selected_arcs] @; status_buffer=arcs_status[-n_selected_special_arcs:n_selected_arcs] @; @ Packing the nodes is very similar to packing the arcs, only there is no need to update any arrays or renumber anything, but just pack them (renumbering of the nodes was taken care of above). Here are the corresponding pieces of code: @ @= ALLOCATE(coords_buffer(n_dim,-n_special_nodes:n_selected_nodes), STAT=error_status) @; CALL RecordAllocation(n_elements=n_dim*(n_special_nodes+n_selected_nodes+1), mold=1.0, & caller="CreateLatticeNetwork", alloc_status=error_status) @; ALLOCATE(status_buffer(-n_special_nodes:n_selected_nodes), STAT=error_status) @; CALL RecordAllocation(n_elements=n_special_nodes+n_selected_nodes+1, mold=1_i_byte, & caller="CreateLatticeNetwork", alloc_status=error_status) @; coords_buffer[1:n_dim,0]=0.0 @; // Null element status_buffer[0]=empty_node @; // Dummy node @ @= CALL RecordAllocation(n_elements=-_SIZE(nodes_coords, i_wp), & mold=1.0, caller="CreateLatticeNetwork") @; DEALLOCATE(nodes_coords) @; // Start the re-allocation process ALLOCATE(nodes_coords(n_dim,-n_special_nodes:n_selected_nodes), STAT=error_status) @; CALL RecordAllocation(n_elements=n_dim*(n_special_nodes+n_selected_nodes+1), mold=1.0, & caller="CreateLatticeNetwork", alloc_status=error_status) @; nodes_coords=coords_buffer @; // Copy buffer CALL RecordAllocation(n_elements=-_SIZE(coords_buffer, i_wp), & mold=1.0, caller="CreateLatticeNetwork") @; DEALLOCATE(coords_buffer) @; CALL RecordAllocation(n_elements=-_SIZE(nodes_status, i_wp), & mold=1_i_byte, caller="CreateLatticeNetwork") @; DEALLOCATE(nodes_status) @; // Start the re-allocation process ALLOCATE(nodes_status(-n_special_nodes:n_selected_nodes), STAT=error_status) @; CALL RecordAllocation(n_elements=n_special_nodes+n_selected_nodes+1, mold=1_i_byte, & caller="CreateLatticeNetwork", alloc_status=error_status) @; nodes_status=status_buffer @; // Copy buffer CALL RecordAllocation(n_elements=-_SIZE(status_buffer, i_wp), & mold=1_i_byte, caller="CreateLatticeNetwork") @; DEALLOCATE(status_buffer) @; @ @= coords_buffer[1:n_dim,:-1]=nodes_coords[1:n_dim,:-1] // Just a copy for special nodes status_buffer[:-1]=nodes_status[:-1] @; // No need to |PACK| for special nodes node_counter=0 @; PackNodes: DO node=1, n_nodes IF(nodes_mask[node]) THEN @; // Pack this node node_counter=node_counter+1 @; coords_buffer[1:n_dim,node_counter]=nodes_coords[1:n_dim,node] @; status_buffer[node_counter]=nodes_status[node] @; END IF @; END DO PackNodes @; @ @= coords_buffer=nodes_coords[1:n_dim,-n_special_nodes:n_selected_nodes] @; status_buffer=nodes_status[-n_special_nodes:n_selected_nodes] @; @*1 Destroying the Lattice. The routine |DestroyNetworkLattice| should be called to deallocate the arrays that were allocated in |CreateNetworkLattice|: @ @= SUBROUTINE DestroyLatticeNetwork () @; IMPLICIT NONE @; CALL RecordAllocation(n_elements=-_SIZE(nodes_coords, i_wp), & mold=1.0, caller="DestroyLatticeNetwork") @; DEALLOCATE(nodes_coords) @; // Allocate enough space CALL RecordAllocation(n_elements=-_SIZE(nodes_status, i_wp), & mold=1_i_byte, caller="DestroyLatticeNetwork") @; DEALLOCATE(nodes_status) @; // Allocate enough space CALL RecordAllocation(n_elements=-_SIZE(heads_tails, i_wp), & mold=1_i_wp, caller="DestroyLatticeNetwork") @; DEALLOCATE(heads_tails) @; CALL RecordAllocation(n_elements=-_SIZE(arcs_status, i_wp), & mold=1_i_byte, caller="DestroyLatticeNetwork") @; DEALLOCATE(arcs_status) @; END SUBROUTINE DestroyLatticeNetwork @; @%%