MODULE Network_Matrix_Operations
      USE Precision
      USE Error_Handling
      USE System_Monitors
      IMPLICIT NONE
      PUBLIC::NodesExcessFlows,ArcsVoltages,Multiply_ADAt,Diagonal_ADAt
      PRIVATE
      
      CONTAINS
      
      
      SUBROUTINE NodesExcessFlows(heads_tails,arcs_mask,node_offset,&
      arcs_flows,excess_flows,initialize)
      INTEGER(KIND=i_wp),INTENT(IN)::node_offset
      INTEGER(KIND=i_wp),DIMENSION(:,:),INTENT(IN)::heads_tails
      LOGICAL(KIND=l_wp),DIMENSION(:),INTENT(IN),OPTIONAL::arcs_mask
      REAL(KIND=r_wp),DIMENSION(:),INTENT(IN)::arcs_flows
      REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT)::excess_flows
      LOGICAL,INTENT(IN),OPTIONAL::initialize
      
      INTEGER::head,tail,arc,node
      LOGICAL::first_time
      
      first_time=.TRUE.
      IF(PRESENT(initialize))first_time=initialize
      
      IF(first_time)excess_flows=0.0_r_wp
      IF(PRESENT(arcs_mask))THEN
      
      DO arc=INT(LBOUND(heads_tails,DIM=2),KIND=i_wp),INT(UBOUND(heads_tails,DIM=2),KIND=i_wp)
      IF(arcs_mask(arc))THEN
      head=heads_tails(1,arc)
      tail=heads_tails(2,arc)
      excess_flows(head)=excess_flows(head)+arcs_flows(arc)
      excess_flows(tail)=excess_flows(tail)-arcs_flows(arc)
      END IF
      END DO
      
      
      ELSE
      
      DO arc=INT(LBOUND(heads_tails,DIM=2),KIND=i_wp),INT(UBOUND(heads_tails,DIM=2),KIND=i_wp)
      head=heads_tails(1,arc)
      tail=heads_tails(2,arc)
      excess_flows(head)=excess_flows(head)+arcs_flows(arc)
      excess_flows(tail)=excess_flows(tail)-arcs_flows(arc)
      END DO
      
      
      END IF
      
      END SUBROUTINE NodesExcessFlows
      
      
      
      SUBROUTINE ArcsVoltages(heads_tails,node_offset,arcs_mask,&
      nodes_potentials,arcs_voltages)
      INTEGER(KIND=i_wp),INTENT(IN)::node_offset
      INTEGER(KIND=i_wp),DIMENSION(:,:),INTENT(IN)::heads_tails
      LOGICAL(KIND=l_wp),DIMENSION(:),INTENT(IN),OPTIONAL::arcs_mask
      REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(IN)::nodes_potentials
      REAL(KIND=r_wp),DIMENSION(:),INTENT(OUT)::arcs_voltages
      
      INTEGER::head,tail,arc,node
      
      IF(PRESENT(arcs_mask))THEN
      
      WHERE(arcs_mask)&
      arcs_voltages=nodes_potentials(heads_tails(1,:))-nodes_potentials(heads_tails(2,:))
      
      
      ELSE
      
      arcs_voltages=nodes_potentials(heads_tails(1,:))-nodes_potentials(heads_tails(2,:))
      
      
      END IF
      
      END SUBROUTINE ArcsVoltages
      
      
      
      SUBROUTINE Multiply_ADAt(heads_tails,node_offset,arcs_conductances,&
      arcs_mask,nodes_flows,nodes_potentials,initialize)
      INTEGER(KIND=i_wp),INTENT(IN)::node_offset
      INTEGER(KIND=i_wp),DIMENSION(:,:),INTENT(IN)::heads_tails
      REAL(KIND=r_wp),DIMENSION(:),INTENT(IN)::arcs_conductances
      LOGICAL(KIND=l_wp),DIMENSION(:),INTENT(IN),OPTIONAL::arcs_mask
      REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(IN)::nodes_potentials
      REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT)::nodes_flows
      LOGICAL,INTENT(IN),OPTIONAL::initialize
      
      INTEGER::head,tail,arc,node,n_arcs
      REAL(KIND=r_wp)::arc_flow
      LOGICAL::first_time
      
      first_time=.TRUE.
      IF(PRESENT(initialize))first_time=initialize
      
      IF(first_time)nodes_flows=0.0_r_wp
      IF(PRESENT(arcs_mask))THEN
      
      DO arc=INT(LBOUND(heads_tails,DIM=2),KIND=i_wp),INT(UBOUND(heads_tails,DIM=2),KIND=i_wp)
      IF(arcs_mask(arc))THEN
      head=heads_tails(1,arc)
      tail=heads_tails(2,arc)
      arc_flow=arcs_conductances(arc)*(nodes_potentials(head)-nodes_potentials(tail))
      nodes_flows(head)=nodes_flows(head)+arc_flow
      nodes_flows(tail)=nodes_flows(tail)-arc_flow
      END IF
      END DO
      
      
      ELSE
      
      DO arc=INT(LBOUND(heads_tails,DIM=2),KIND=i_wp),INT(UBOUND(heads_tails,DIM=2),KIND=i_wp)
      head=heads_tails(1,arc)
      tail=heads_tails(2,arc)
      arc_flow=arcs_conductances(arc)*(nodes_potentials(head)-nodes_potentials(tail))
      nodes_flows(head)=nodes_flows(head)+arc_flow
      nodes_flows(tail)=nodes_flows(tail)-arc_flow
      END DO
      
      
      END IF
      
      n_arcs=INT(SIZE(heads_tails,DIM=2),KIND=i_wp)
      n_reads=n_reads+7*n_arcs
      n_writes=n_writes+2*n_arcs
      n_flops=n_flops+4*n_arcs
      
      END SUBROUTINE Multiply_ADAt
      
      
      
      SUBROUTINE Diagonal_ADAt(heads_tails,node_offset,&
      arcs_conductances,nodes_conductances,nodes_resistances,&
      arcs_mask,initialize)
      INTEGER(KIND=i_wp),INTENT(IN)::node_offset
      INTEGER(KIND=i_wp),DIMENSION(:,:),INTENT(IN)::heads_tails
      REAL(KIND=r_wp),DIMENSION(:),INTENT(IN)::arcs_conductances
      REAL(KIND=r_wp),DIMENSION(-node_offset:),INTENT(OUT),OPTIONAL::&
      nodes_conductances,nodes_resistances
      LOGICAL(KIND=l_wp),DIMENSION(:),INTENT(IN),OPTIONAL::arcs_mask
      LOGICAL,INTENT(IN),OPTIONAL::initialize
      
      INTEGER::head,tail,arc,node,n_arcs
      LOGICAL::first_time
      
      first_time=.TRUE.
      IF(PRESENT(initialize))first_time=initialize
      
      IF(first_time)nodes_conductances=0.0_r_wp
      IF(PRESENT(arcs_mask))THEN
      
      IF(PRESENT(nodes_conductances))THEN
      
      DO arc=INT(LBOUND(heads_tails,DIM=2),KIND=i_wp),INT(UBOUND(heads_tails,DIM=2),KIND=i_wp)
      IF(arcs_mask(arc))THEN
      head=heads_tails(1,arc)
      tail=heads_tails(2,arc)
      nodes_conductances(head)=nodes_conductances(head)+arcs_conductances(arc)
      nodes_conductances(tail)=nodes_conductances(tail)+arcs_conductances(arc)
      END IF
      END DO
      ELSE IF(PRESENT(nodes_resistances))THEN
      
      DO arc=INT(LBOUND(heads_tails,DIM=2),KIND=i_wp),INT(UBOUND(heads_tails,DIM=2),KIND=i_wp)
      IF(arcs_mask(arc))THEN
      head=heads_tails(1,arc)
      tail=heads_tails(2,arc)
      nodes_resistances(head)=nodes_resistances(head)+arcs_conductances(arc)
      nodes_resistances(tail)=nodes_resistances(tail)+arcs_conductances(arc)
      END IF
      END DO
      nodes_resistances=1.0_r_wp/(nodes_resistances+EPSILON(1.0_r_wp))
      END IF
      ELSE
      
      IF(PRESENT(nodes_conductances))THEN
      
      DO arc=INT(LBOUND(heads_tails,DIM=2),KIND=i_wp),INT(UBOUND(heads_tails,DIM=2),KIND=i_wp)
      head=heads_tails(1,arc)
      tail=heads_tails(2,arc)
      nodes_conductances(head)=nodes_conductances(head)+arcs_conductances(arc)
      nodes_conductances(tail)=nodes_conductances(tail)+arcs_conductances(arc)
      END DO
      ELSE IF(PRESENT(nodes_resistances))THEN
      
      DO arc=INT(LBOUND(heads_tails,DIM=2),KIND=i_wp),INT(UBOUND(heads_tails,DIM=2),KIND=i_wp)
      head=heads_tails(1,arc)
      tail=heads_tails(2,arc)
      nodes_resistances(head)=nodes_resistances(head)+arcs_conductances(arc)
      nodes_resistances(tail)=nodes_resistances(tail)+arcs_conductances(arc)
      END DO
      nodes_resistances=1.0_r_wp/(nodes_resistances+EPSILON(1.0_r_wp))
      END IF
      END IF
      
      END SUBROUTINE Diagonal_ADAt
      
      
      
      END MODULE Network_Matrix_Operations