MODULE LSNNO_Interface USE Precision USE Error_Handling USE System_Monitors USE Network_Data_Structures IMPLICIT NONE PRIVATE PUBLIC::InitializeLSNNO,CreateLSNNO,DestroyLSNNO,CallLSNNO,CalculateLSNNOPotentials INTEGER,ALLOCATABLE,DIMENSION(:)::TO,FR INTEGER::NA,NN,NEL,LFUVAL,LIWK,LWK,IPDEVC,LFC,INFO,IFLAG,INFORM INTEGER,ALLOCATABLE,DIMENSION(:)::ELVAR,ELPTR,GPTR,HPTR,IWK,ELST,FCALC,XST INTEGER::IFC(3),IT(3) REAL(KIND=r_dp)::FX REAL(KIND=r_dp),ALLOCATABLE,DIMENSION(:)::FUVAL,WK,X CHARACTER,PUBLIC,SAVE::initial_flow='I' INTEGER,PUBLIC,SAVE::MAXIT=1000,WHAT=2,FREQ=1 REAL(KIND=r_dp),PUBLIC,SAVE::EPSF=1.0E-3_r_dp LOGICAL,PUBLIC,SAVE::PREC=.TRUE.,NIPST=.TRUE. INTEGER,PUBLIC,SAVE::ELST_low_cost=1,ELST_high_cost=1,& ELST_regular_cost=1 NAMELIST/LSNNOOptions/initial_flow,MAXIT,EPSF,& ELST_low_cost,ELST_high_cost,ELST_regular_cost,& PREC,NIPST,WHAT,FREQ INTERFACE SUBROUTINE LSNNO(NA,NN,NEL,ELVAR,ELPTR,ELST,FR,TO,X,& XST,FX,FUVAL,LFUVAL,GPTR,HPTR,INFORM,FCALC,& LFC,IFC,IT,MAXIT,EPSF,PREC,NIPST,IWK,& LIWK,WK,LWK,IPDEVC,WHAT,FREQ,INFO,IFLAG) USE Precision IMPLICIT NONE INTEGER,INTENT(IN)::NA,NN,NEL,LFUVAL,MAXIT,LIWK,LWK,& IPDEVC,WHAT,FREQ INTEGER,INTENT(OUT)::LFC,INFO,IFLAG INTEGER,INTENT(INOUT)::INFORM INTEGER,DIMENSION(NA),INTENT(IN)::FR,TO INTEGER,DIMENSION(NEL),INTENT(IN)::ELVAR INTEGER,DIMENSION(NEL+1),INTENT(IN)::ELPTR INTEGER,DIMENSION(NEL+1),INTENT(OUT)::GPTR,HPTR INTEGER,DIMENSION(3)::IFC,IT INTEGER,DIMENSION(LIWK)::IWK INTEGER,DIMENSION(NEL),INTENT(INOUT)::ELST,FCALC INTEGER,DIMENSION(NA),INTENT(INOUT)::XST REAL(KIND=r_dp),INTENT(IN)::EPSF REAL(KIND=r_dp),INTENT(OUT)::FX REAL(KIND=r_dp),DIMENSION(LFUVAL),INTENT(IN)::FUVAL REAL(KIND=r_dp),DIMENSION(LWK),INTENT(OUT)::WK REAL(KIND=r_dp),DIMENSION(NA),INTENT(INOUT)::X LOGICAL,INTENT(IN)::PREC,NIPST END SUBROUTINE END INTERFACE CONTAINS SUBROUTINE InitializeLSNNO() REWIND(UNIT=program_options_unit) READ(UNIT=program_options_unit,NML=LSNNOOptions,IOSTAT=error_status) IF(error_status/=0)& CALL CriticalError(message="NAMELIST LSNNOOptions was not read successfully"//& " from file "//TRIM(options_file),caller="InitializeLSNNO") WRITE(UNIT=message_log_unit,NML=LSNNOOptions) END SUBROUTINE InitializeLSNNO SUBROUTINE CreateLSNNO() IMPLICIT NONE INTEGER::arc,special_arc,node,alloc_status NA=n_arcs+n_special_arcs NN=n_nodes+n_special_nodes NEL=NA LFUVAL=3*NA ALLOCATE(TO(NA),STAT=alloc_status) CALL RecordAllocation(n_elements=NA,& mold=1,caller="CreateLSNNO",alloc_status=alloc_status) ALLOCATE(FR(NA),STAT=alloc_status) CALL RecordAllocation(n_elements=NA,& mold=1,caller="CreateLSNNO",alloc_status=alloc_status) ALLOCATE(X(NA),STAT=alloc_status) CALL RecordAllocation(n_elements=NA,mold=1.0_r_dp,& caller="CreateLSNNO",alloc_status=alloc_status) ALLOCATE(ELVAR(NEL),STAT=alloc_status) CALL RecordAllocation(n_elements=NEL,mold=1,& caller="CreateLSNNO",alloc_status=alloc_status) ALLOCATE(ELPTR(NEL+1),STAT=alloc_status) CALL RecordAllocation(n_elements=NEL+1,mold=1,& caller="CreateLSNNO",alloc_status=alloc_status) ALLOCATE(ELST(NEL),STAT=alloc_status) CALL RecordAllocation(n_elements=NEL,mold=1,& caller="CreateLSNNO",alloc_status=alloc_status) ALLOCATE(XST(NA),STAT=alloc_status) CALL RecordAllocation(n_elements=NA,mold=1,& caller="CreateLSNNO",alloc_status=alloc_status) ALLOCATE(GPTR(NEL+1),STAT=alloc_status) CALL RecordAllocation(n_elements=NEL+1,mold=1,& caller="CreateLSNNO",alloc_status=alloc_status) ALLOCATE(HPTR(NEL+1),STAT=alloc_status) CALL RecordAllocation(n_elements=NEL+1,mold=1,& caller="CreateLSNNO",alloc_status=alloc_status) ALLOCATE(FCALC(NEL),STAT=alloc_status) CALL RecordAllocation(n_elements=NEL,mold=1,& caller="CreateLSNNO",alloc_status=alloc_status) ALLOCATE(FUVAL(LFUVAL),STAT=alloc_status) CALL RecordAllocation(n_elements=LFUVAL,mold=1.0_r_dp,& caller="CreateLSNNO",alloc_status=alloc_status) END SUBROUTINE CreateLSNNO SUBROUTINE CallLSNNO(ElementalCosts) IMPLICIT NONE INTERFACE PURE SUBROUTINE ElementalCosts(flow,potential,resistance,cost,& arc_index,arguments_status,tolerance) USE Precision REAL(KIND=r_wp),INTENT(INOUT)::flow,potential,resistance,cost INTEGER(KIND=i_wp),INTENT(IN)::arc_index CHARACTER(LEN=4),INTENT(IN)::arguments_status REAL(KIND=r_wp),INTENT(IN)::tolerance END SUBROUTINE ElementalCosts END INTERFACE INTEGER::alloc_status INTEGER::element,arc,arc_index,node,special_arc,special_node,head,tail REAL(KIND=r_wp)::dummy IPDEVC=message_log_unit LWK=0 ALLOCATE(WK(LWK),STAT=alloc_status) LIWK=0 ALLOCATE(IWK(LIWK),STAT=alloc_status) FR(1:n_arcs)=heads_tails(1,1:n_arcs) TO(1:n_arcs)=heads_tails(2,1:n_arcs) DO special_arc=1,n_special_arcs head=heads_tails(1,-special_arc) tail=heads_tails(2,-special_arc) IF(head<0)THEN FR(n_arcs+special_arc)=n_nodes+ABS(head) ELSE FR(n_arcs+special_arc)=head END IF IF(tail<0)THEN TO(n_arcs+special_arc)=n_nodes+ABS(tail) ELSE TO(n_arcs+special_arc)=tail END IF END DO IF(initial_flow=='I'.OR.initial_flow=='i')THEN X=1.0_r_dp ELSE X(1:n_arcs)=arcs_flows(1:n_arcs) X(n_arcs+1:(n_arcs+n_special_arcs))=arcs_flows((-1):(-n_special_arcs):(-1)) END IF DO arc=1,NA ELVAR(arc)=arc ELPTR(arc)=arc END DO ELPTR(NA+1)=NA+1 ELST(1:n_arcs)=ELST_regular_cost DO special_arc=1,n_special_arcs SELECT CASE(special_arcs_status(-special_arc)) CASE(low_cost_arc) ELST(n_arcs+special_arc)=ELST_low_cost CASE(high_cost_arc) ELST(n_arcs+special_arc)=ELST_high_cost CASE DEFAULT ELST(n_arcs+special_arc)=ELST_regular_cost ENDSELECT END DO XST=-1 INFORM=0 ReverseLoop:DO CALL LSNNO(NA,NN,NEL,ELVAR,ELPTR,ELST,FR,TO,X,& XST,FX,FUVAL,LFUVAL,GPTR,HPTR,INFORM,FCALC,& LFC,IFC,IT,MAXIT,EPSF,PREC,NIPST,IWK,& LIWK,WK,LWK,IPDEVC,WHAT,FREQ,INFO,IFLAG) IF(IFLAG==1)THEN WRITE(message_log_unit,*)"Requested size of WK in LSNNO:",INFO DEALLOCATE(WK) IF(LWK>0)CALL RecordAllocation(n_elements=-LWK,mold=1.0_r_dp) LWK=INFO*10 ALLOCATE(WK(LWK),STAT=alloc_status) CALL RecordAllocation(n_elements=LWK,mold=1.0_r_dp,& caller="CallLSNNO",alloc_status=alloc_status) CYCLE ReverseLoop ELSE IF(IFLAG==2)THEN WRITE(message_log_unit,*)"Requested size of IWK in LSNNO:",INFO DEALLOCATE(IWK) IF(LIWK>0)CALL RecordAllocation(n_elements=-LIWK,mold=1) LIWK=INFO*10 ALLOCATE(IWK(LIWK),STAT=alloc_status) CALL RecordAllocation(n_elements=LIWK,mold=1,& caller="CallLSNNO",alloc_status=alloc_status) CYCLE ReverseLoop ELSE IF(IFLAG>=3)THEN WRITE(error_print_unit,*)"LSSNO's IFLAG=",IFLAG WRITE(error_log_unit,*)"LSSNO's IFLAG=",IFLAG CALL NonCriticalError("LSNNO aborted with IFLAG>=3",caller="CallLSNNO") EXIT ReverseLoop END IF IF(INFORM>0)THEN SELECT CASE(INFORM) CASE(1) DO element=1,LFC arc=FCALC(element) IF(arc>n_arcs)THEN arc_index=n_arcs-arc ELSE arc_index=arc END IF CALL ElementalCosts(X(arc),dummy,dummy,FUVAL(arc),& arc_index,"FDDC",EPSILON(1.0_r_wp)) END DO CASE(2) DO element=1,LFC arc=FCALC(element) IF(arc>n_arcs)THEN arc_index=n_arcs-arc ELSE arc_index=arc END IF CALL ElementalCosts(X(arc),FUVAL(NA+arc),dummy,dummy,& arc_index,"FCDD",EPSILON(1.0_r_wp)) END DO CASE(3) DO element=1,LFC arc=FCALC(element) IF(arc>n_arcs)THEN arc_index=n_arcs-arc ELSE arc_index=arc END IF CALL ElementalCosts(X(arc),dummy,FUVAL(2*NA+arc),dummy,& arc_index,"FDCD",EPSILON(1.0_r_wp)) END DO CASE(4) DO element=1,LFC arc=FCALC(element) IF(arc>n_arcs)THEN arc_index=n_arcs-arc ELSE arc_index=arc END IF CALL ElementalCosts(X(arc),FUVAL(NA+arc),FUVAL(2*NA+arc),dummy,& arc_index,"FCCD",EPSILON(1.0_r_wp)) END DO CASE DEFAULT CALL CriticalError("INFORM>4 in LSNNO",caller="CallLSNNO") ENDSELECT ELSE WRITE(message_print_unit,*)"LSNNO returned with success!!!" arcs_flows(1:n_arcs)=X(1:n_arcs) arcs_flows((-1):(-n_special_arcs):(-1))=X(n_arcs+1:(n_arcs+n_special_arcs)) arcs_flows(0)=0.0_r_dp EXIT ReverseLoop END IF END DO ReverseLoop END SUBROUTINE CallLSNNO SUBROUTINE DestroyLSNNO IMPLICIT NONE INTEGER::alloc_status DEALLOCATE(TO,STAT=alloc_status) CALL RecordAllocation(n_elements=-NA,& mold=1,caller="DestroyLSNNO",alloc_status=alloc_status) DEALLOCATE(FR,STAT=alloc_status) CALL RecordAllocation(n_elements=-NA,& mold=1,caller="DestroyLSNNO",alloc_status=alloc_status) DEALLOCATE(X,STAT=alloc_status) CALL RecordAllocation(n_elements=-NA,mold=1.0_r_dp,& caller="DestroyLSNNO",alloc_status=alloc_status) DEALLOCATE(ELVAR,STAT=alloc_status) CALL RecordAllocation(n_elements=-NEL,mold=1,& caller="DestroyLSNNO",alloc_status=alloc_status) DEALLOCATE(ELPTR,STAT=alloc_status) CALL RecordAllocation(n_elements=-(NEL+1),mold=1,& caller="DestroyLSNNO",alloc_status=alloc_status) DEALLOCATE(ELST,STAT=alloc_status) CALL RecordAllocation(n_elements=-NEL,mold=1,& caller="DestroyLSNNO",alloc_status=alloc_status) DEALLOCATE(XST,STAT=alloc_status) CALL RecordAllocation(n_elements=-NA,mold=1,& caller="DestroyLSNNO",alloc_status=alloc_status) DEALLOCATE(GPTR,STAT=alloc_status) CALL RecordAllocation(n_elements=-(NEL+1),mold=1,& caller="DestroyLSNNO",alloc_status=alloc_status) DEALLOCATE(HPTR,STAT=alloc_status) CALL RecordAllocation(n_elements=-(NEL+1),mold=1,& caller="DestroyLSNNO",alloc_status=alloc_status) DEALLOCATE(FCALC,STAT=alloc_status) CALL RecordAllocation(n_elements=-NEL,mold=1,& caller="DestroyLSNNO",alloc_status=alloc_status) DEALLOCATE(FUVAL,STAT=alloc_status) CALL RecordAllocation(n_elements=-LFUVAL,mold=1.0_r_dp,& caller="DestroyLSNNO",alloc_status=alloc_status) DEALLOCATE(WK,STAT=alloc_status) CALL RecordAllocation(n_elements=-LWK,mold=1.0_r_dp,& caller="DestroyLSNNO",alloc_status=alloc_status) DEALLOCATE(IWK,STAT=alloc_status) CALL RecordAllocation(n_elements=-LIWK,mold=1,& caller="DestroyLSNNO",alloc_status=alloc_status) END SUBROUTINE DestroyLSNNO SUBROUTINE CalculateLSNNOPotentials(ElementalCosts,potential_mismatch,flow_mismatch) USE Network_Matrix_Operations USE Graph_Algorithms USE Network_Spanning_Trees IMPLICIT NONE INTERFACE PURE SUBROUTINE ElementalCosts(flow,potential,resistance,cost,& arc_index,arguments_status,tolerance) USE Precision REAL(KIND=r_wp),INTENT(INOUT)::flow,potential,resistance,cost INTEGER(KIND=i_wp),INTENT(IN)::arc_index CHARACTER(LEN=4),INTENT(IN)::arguments_status REAL(KIND=r_wp),INTENT(IN)::tolerance END SUBROUTINE ElementalCosts END INTERFACE REAL(KIND=r_wp),DIMENSION(2),INTENT(OUT),OPTIONAL::& potential_mismatch,flow_mismatch REAL(KIND=r_wp),DIMENSION(:),ALLOCATABLE::temp_buffer INTEGER(KIND=i_wp)::node,arc INTEGER::alloc_status REAL(KIND=r_wp)::dummy DO arc=-n_special_arcs,n_arcs CALL ElementalCosts(flow=arcs_flows(arc),potential=arcs_voltages(arc),& resistance=dummy,cost=dummy,arc_index=arc,& arguments_status="FCDD",tolerance=EPSILON(1.0_r_wp)) END DO CALL CreateSpanningTree(tree_type=min_cost_tree,arcs_weights=ABS(arcs_voltages),& weights_distribution="Random") nodes_potentials(0:1)=0.0_r_wp CALL PropagateNodesPotentials(arc_offset=n_special_arcs,node_offset=n_special_nodes,& heads_tails=heads_tails,orientations=tree_arcs_orientations,& parents=tree_nodes_parents,level_ordering=tree_nodes_ordering,& arcs_potential_drops=arcs_voltages,nodes_potentials=nodes_potentials) IF(PRESENT(potential_mismatch))THEN ALLOCATE(temp_buffer(-n_special_arcs:n_arcs),STAT=alloc_status) CALL RecordAllocation(n_elements=n_special_arcs+n_arcs+1,mold=1.0_r_wp,& caller="CalculateLSNNOPotentials",alloc_status=alloc_status) CALL ArcsVoltages(heads_tails=heads_tails,node_offset=n_special_nodes,& nodes_potentials=nodes_potentials,arcs_voltages=temp_buffer) potential_mismatch=PercentMismatch(arcs_voltages,temp_buffer) CALL RecordAllocation(n_elements=-INT(SIZE(temp_buffer),KIND=i_wp),mold=1.0_r_wp) DEALLOCATE(temp_buffer) END IF IF(PRESENT(flow_mismatch))THEN ALLOCATE(temp_buffer(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=n_special_nodes+n_nodes+1,mold=1.0_r_wp,& caller="CalculateLSNNOPotentials",alloc_status=alloc_status) CALL NodesExcessFlows(heads_tails=heads_tails,node_offset=n_special_nodes,& arcs_flows=arcs_flows,excess_flows=temp_buffer) flow_mismatch=PercentMismatch(supplies_demands,temp_buffer) CALL RecordAllocation(n_elements=-INT(SIZE(temp_buffer),KIND=i_wp),mold=1.0_r_wp) DEALLOCATE(temp_buffer) END IF CONTAINS FUNCTION PercentMismatch(A,B)RESULT(norms) USE Precision REAL(KIND=r_wp),DIMENSION(:),INTENT(IN)::A,B REAL(KIND=r_wp),DIMENSION(2)::norms INTEGER(KIND=i_wp)::element,n_elements REAL(KIND=i_wp)::magnitude n_elements=INT(SIZE(A),KIND=i_wp) norms=0.0_r_wp magnitude=0.0_r_wp DO element=INT(LBOUND(A,DIM=1),KIND=i_wp),INT(UBOUND(A,DIM=1),KIND=i_wp) norms(1)=norms(1)+(A(element)-B(element))**2 norms(2)=MAX(norms(2),ABS(A(element)-B(element))) magnitude=magnitude+0.5_r_wp*(ABS(A(element))+ABS(B(element))) END DO magnitude=magnitude/REAL(n_elements,r_wp) norms(1)=SQRT(norms(1)/REAL(n_elements,r_wp))/magnitude norms(2)=norms(2)/magnitude END FUNCTION PercentMismatch END SUBROUTINE CalculateLSNNOPotentials END MODULE LSNNO_Interface FUNCTION RHS(node)RESULT(supply_demand) USE Precision USE Network_Data_Structures IMPLICIT NONE INTEGER,INTENT(IN)::node REAL(KIND=r_dp)::supply_demand IF(node>n_nodes)THEN supply_demand=supplies_demands(n_nodes-node) ELSE supply_demand=supplies_demands(node) END IF END FUNCTION RHS FUNCTION XLOWER(arc)RESULT(lb) USE Precision IMPLICIT NONE INTEGER,INTENT(IN)::arc REAL(KIND=r_dp)::lb lb=-1.0_r_dp/EPSILON(1.0_r_dp) END FUNCTION XLOWER FUNCTION XUPPER(arc)RESULT(ub) USE Precision IMPLICIT NONE INTEGER,INTENT(IN)::arc REAL(KIND=r_dp)::ub ub=1.0_r_dp/EPSILON(1.0_r_dp) END FUNCTION XUPPER SUBROUTINE RANGE(IEL,ELDIM,INTDIM,W1,W2,MODE) USE Precision USE Error_Handling IMPLICIT NONE INTEGER,INTENT(IN)::IEL,ELDIM,MODE INTEGER,INTENT(OUT)::INTDIM REAL(KIND=r_dp),DIMENSION(*),INTENT(INOUT)::W1,W2 IF(MODE==0)THEN INTDIM=ELDIM ELSE Call CriticalError("RANGE called by LSNNO with non-zero MODE",caller="RANGE") END IF END SUBROUTINE RANGE