PROGRAM Saturation USE Precision USE Error_Handling USE System_Monitors USE Initialization_Termination USE Network_Data_Structures USE Lattice_Geometry USE Network_Geometry USE Simple_Graphics USE Network_Graphics USE Network_Spanning_Trees USE Graph_Algorithms USE Graph_Algorithms_F77 USE jj_i_Cost_Parameters USE jj_i_Cost_Functions USE Lattice_Network_Optimization USE Conjugate_Gradient USE Dual_Network_Solvers USE Dual_Line_Minimizers USE Dual_Newton_SSCNO USE SSCNO_Interface IMPLICIT NONE INTEGER(KIND=i_wp),DIMENSION(:,:),ALLOCATABLE::special_heads_tails INTEGER(KIND=i_byte),DIMENSION(:),ALLOCATABLE::nodes_cuts REAL(KIND=r_wp),DIMENSION(:),ALLOCATABLE::arcs_cuts_capacities,cuts_capacities INTEGER(KIND=i_byte),DIMENSION(:),ALLOCATABLE::arcs_cuts INTEGER(KIND=i_byte)::head_cut,tail_cut,n_cuts=1 INTEGER::arc,node,sample,point,instance CHARACTER(LEN=60)::title,plot_title(3)="" CHARACTER(LEN=50)::file_name="",file_type="CONS",file_extension="dummy" REAL(KIND=r_wp)::applied_voltages_sum,min_arc_value,max_arc_value REAL(KIND=r_wp),DIMENSION(2)::current_interval=(/0.0_r_wp,1.0_r_wp/),& voltage_interval=(/0.0_r_wp,1.0_r_wp/) INTEGER::n_samples=1,n_points=2,n_instances=1 INTEGER::alloc_status,source_type REAL(KIND=r_wp),DIMENSION(4)::axis REAL(KIND=r_wp),DIMENSION(:),ALLOCATABLE::total_currents,total_voltages,& averaged_currents,averaged_voltages REAL(KIND=r_wp),DIMENSION(:,:),ALLOCATABLE::applied_voltages,injected_currents LOGICAL::plot_samples=.FALSE.,plot_potentials=.FALSE.,& change_current=.FALSE.,change_voltage=.TRUE.,save_plot WRITE(*,*)"Enter change_current and change_voltage (T or F)" READ(*,*)change_current,change_voltage IF(change_current)THEN WRITE(*,*)"Enter the interval for the total injected current:" READ(*,*)current_interval END IF IF(change_voltage)THEN WRITE(*,*)"Enter the interval for the total applied voltage:" READ(*,*)voltage_interval END IF WRITE(*,*)"Enter n_points, n_cuts" READ(*,*)n_points,n_cuts WRITE(*,*)"Plot every sample (T or F):" READ(*,*)plot_samples IF(plot_samples)THEN WRITE(*,*)"Plot the nodes potentials as well?:" READ(*,*)plot_potentials END IF WRITE(*,*)"Enter file output type (CONS, XWIN, PNG, POST, PSCL)" READ(*,*)file_type SELECT CASE(file_type) CASE("POST","PSCL") file_extension="ps" CASE("PNG") file_extension="png" CASE DEFAULT file_extension="dummy" ENDSELECT CALL StartProgram ALLOCATE(total_currents(n_points),total_voltages(n_points)) ALLOCATE(averaged_currents(n_points),averaged_voltages(n_points)) ALLOCATE(applied_voltages(n_points,n_samples),injected_currents(n_points,n_samples)) CALL InitializeLatticeNetworkProblem CALL InitializeSSCNO CALL ResetTimer(100) CALL StartTimer(100) DO point=1,n_points total_currents(point)=current_interval(1)+& REAL(point-1)*(current_interval(2)-current_interval(1))/REAL(n_points-1) total_voltages(point)=voltage_interval(1)+& REAL(point-1)*(voltage_interval(2)-voltage_interval(1))/REAL(n_points-1) END DO injected_currents=0.0_r_wp applied_voltages=0.0_r_wp SourceType:DO source_type=1,2 IF(source_type==1.AND.(.NOT.change_current))CYCLE SourceType IF(source_type==2.AND.(.NOT.change_voltage))CYCLE SourceType SampleLoop:DO sample=1,n_samples WRITE(message_print_unit,*)"->Doing sample #",sample IF(source_type==1)THEN sources_status(1)='R' sinks_status(1)='R' ELSE IF(source_type==2)THEN sources_status(1)='F' sinks_status(1)='F' END IF CALL CreateLatticeNetworkProblem(create_tree=.TRUE.) CALL CreateSSCNO IF(.NOT.ALLOCATED(nodes_cuts))THEN ALLOCATE(nodes_cuts(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_byte,& caller="Breakdown",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(arcs_cuts))THEN ALLOCATE(arcs_cuts(-n_special_arcs:n_arcs),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1_i_byte,& caller="Breakdown",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(arcs_cuts_capacities))THEN ALLOCATE(arcs_cuts_capacities(-n_special_arcs:n_arcs),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1.0_r_wp,& caller="Breakdown",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(cuts_capacities))THEN ALLOCATE(cuts_capacities(1:n_cuts),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_cuts)-(1)+1,mold=1.0_r_wp,& caller="Breakdown",alloc_status=alloc_status) END IF ALLOCATE(special_heads_tails(1:2,-n_special_arcs:(-1)),STAT=alloc_status) CALL RecordAllocation(n_elements=INT(SIZE(special_heads_tails),KIND=i_wp),mold=1_i_wp,& caller="Saturation",alloc_status=alloc_status) InstanceLoop:DO instance=1,n_instances CALL AssignCostParameters() special_heads_tails=heads_tails(1:2,-n_special_arcs:(-1)) heads_tails(1:2,(-1):-n_special_arcs:(-1))=& SSCNO_problem%problem%special_subgraph%heads_tails(1:2,1:n_special_arcs) CALL LabelingMinCut(node_offset=n_special_nodes,arc_offset=n_special_arcs,& heads_tails=heads_tails,arcs_capacities=arcs_cost_parameters(1,1:n_arcs),& source_sink=(/sources(1),sinks(1)/),nodes_cuts=nodes_cuts,arcs_cuts=arcs_cuts(1:n_arcs),& n_cuts=n_cuts,cuts_capacities=cuts_capacities) heads_tails(1:2,-n_special_arcs:(-1))=special_heads_tails WHERE(cuts_capacities<0.0_r_wp)cuts_capacities=0.0_r_wp nodes_cuts(1:n_nodes)=nodes_cuts(1:n_nodes)-MINVAL(nodes_cuts(1:n_nodes))+1 WRITE(*,*)"Cuts capacities:",cuts_capacities DO arc=1,n_arcs IF(arcs_cuts(arc)>0)THEN arcs_cuts_capacities(arc)=cuts_capacities(arcs_cuts(arc)) arcs_cuts(arc)=n_cuts-arcs_cuts(arc)+1 ELSE arcs_cuts_capacities(arc)=0.0_r_wp END IF END DO arcs_cuts_capacities(-n_special_arcs:0)=0.0_r_wp arcs_mask=.TRUE._l_wp arcs_mask(:0)=.FALSE._l_wp DO arc=1,n_arcs arcs_mask(arc)=(arcs_mask(arc).AND.(arcs_status(arc)/=periodic_bc_arc)) END DO nodes_mask=.TRUE._l_wp nodes_mask(:0)=.FALSE._l_wp CALL InitNetworkGraphics(file="Screen",file_type="CONS",& x_label="Current flow \Huge{$\rightarrow$}",& page_size=(/5000,5000/),label_format="5N0",& color_table="RAIN",colorbar_position="Horizontal",& axis_labels_format=(/"NONE","NONE","NONE","NONE"/)) CALL PlotNetwork2D(heads_tails=heads_tails(:,1:n_arcs),& node_offset=-1,node_coords=nodes_coords(:,1:n_nodes),& node_mask=nodes_mask(1:n_nodes),arc_mask=arcs_mask(1:n_arcs),& arc_values=REAL(arcs_cuts(1:n_arcs),r_wp),& node_values=REAL(nodes_cuts(1:n_nodes),r_wp),& resize_nodes=.TRUE.,color_nodes=.TRUE.,& resize_arcs=.TRUE.,color_arcs=.TRUE.,& node_size_range=(/-HUGE(1.0_r_wp)/50,HUGE(1.0_r_wp)/20/),& node_colorbar_format="5F0",& arc_colorbar_format="5F0",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) nodes_mask(0:n_nodes)=.FALSE._l_wp nodes_mask(-n_special_nodes:-1)=.TRUE._l_wp arcs_mask(-1:-n_special_arcs:(-1))=.TRUE. DO arc=-n_special_arcs,-1 arcs_mask(arc)=(arcs_mask(arc).AND.(arcs_status(arc)/=periodic_bc_arc)) END DO CALL PlotNetwork2D(heads_tails=SSCNO_problem%problem%special_subgraph%heads_tails,& node_offset=n_special_nodes,node_coords=nodes_coords,& node_mask=nodes_mask,arc_mask=arcs_mask(-1:-n_special_arcs:(-1)),& arc_values=REAL(arcs_cuts(-1:-n_special_arcs:(-1)),r_wp),& resize_arcs=.TRUE.,color_arcs=.FALSE.,& arc_size_range=(/0.0_r_wp,HUGE(1.0_r_wp)/10/),& node_colorbar_format="5E1",arc_colorbar_format="5E1",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) CALL EndNetworkGraphics() arcs_mask=.TRUE._l_wp arcs_mask(:0)=.FALSE._l_wp DO arc=1,n_arcs arcs_mask(arc)=(arcs_mask(arc).AND.(arcs_status(arc)/=periodic_bc_arc)) END DO nodes_mask=.TRUE._l_wp nodes_mask(:0)=.FALSE._l_wp CALL InitNetworkGraphics(file="Screen",file_type="CONS",& x_label="Current flow \Huge{$\rightarrow$}",& page_size=(/5000,5000/),label_format="5N0",& color_table="RAIN",colorbar_position="Horizontal",& axis_labels_format=(/"NONE","NONE","NONE","NONE"/)) CALL PlotNetwork2D(heads_tails=heads_tails(:,1:n_arcs),& node_offset=-1,node_coords=nodes_coords(:,1:n_nodes),& node_mask=nodes_mask(1:n_nodes),arc_mask=arcs_mask(1:n_arcs),& arc_values=REAL(arcs_cuts_capacities(1:n_arcs),r_wp),& node_values=REAL(nodes_cuts(1:n_nodes),r_wp),& resize_nodes=.TRUE.,color_nodes=.TRUE.,& resize_arcs=.TRUE.,color_arcs=.TRUE.,& node_size_range=(/-HUGE(1.0_r_wp)/50,HUGE(1.0_r_wp)/20/),& node_colorbar_format="5F0",& arc_colorbar_format="5F1",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) nodes_mask(0:n_nodes)=.FALSE._l_wp nodes_mask(-n_special_nodes:-1)=.TRUE._l_wp arcs_mask(-1:-n_special_arcs:(-1))=.TRUE. DO arc=-n_special_arcs,-1 arcs_mask(arc)=(arcs_mask(arc).AND.(arcs_status(arc)/=periodic_bc_arc)) END DO CALL PlotNetwork2D(heads_tails=SSCNO_problem%problem%special_subgraph%heads_tails,& node_offset=n_special_nodes,node_coords=nodes_coords,& node_mask=nodes_mask,arc_mask=arcs_mask(-1:-n_special_arcs:(-1)),& arc_values=REAL(arcs_cuts_capacities(-1:-n_special_arcs:(-1)),r_wp),& resize_arcs=.TRUE.,color_arcs=.FALSE.,& arc_size_range=(/0.0_r_wp,HUGE(1.0_r_wp)/10/),& node_colorbar_format="5E1",arc_colorbar_format="5E1",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) CALL EndNetworkGraphics() IF(file_type/="CONS".AND.file_type/="XWIN")THEN arcs_mask=.TRUE._l_wp arcs_mask(:0)=.FALSE._l_wp DO arc=1,n_arcs arcs_mask(arc)=(arcs_mask(arc).AND.(arcs_status(arc)/=periodic_bc_arc)) END DO nodes_mask=.TRUE._l_wp nodes_mask(:0)=.FALSE._l_wp CALL InitNetworkGraphics(file="MC_MF_cuts."//TRIM(file_extension),file_type=TRIM(file_type),& x_label="Current flow \Huge{$\rightarrow$}",& page_size=(/5000,5000/),label_format="5N0",& color_table="RAIN",colorbar_position="Horizontal",& axis_labels_format=(/"NONE","NONE","NONE","NONE"/)) CALL PlotNetwork2D(heads_tails=heads_tails(:,1:n_arcs),& node_offset=-1,node_coords=nodes_coords(:,1:n_nodes),& node_mask=nodes_mask(1:n_nodes),arc_mask=arcs_mask(1:n_arcs),& arc_values=REAL(arcs_cuts(1:n_arcs),r_wp),& node_values=REAL(nodes_cuts(1:n_nodes),r_wp),& resize_nodes=.TRUE.,color_nodes=.TRUE.,& resize_arcs=.TRUE.,color_arcs=.TRUE.,& node_size_range=(/-HUGE(1.0_r_wp)/50,HUGE(1.0_r_wp)/20/),& node_colorbar_format="5F0",& arc_colorbar_format="5F0",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) nodes_mask(0:n_nodes)=.FALSE._l_wp nodes_mask(-n_special_nodes:-1)=.TRUE._l_wp arcs_mask(-1:-n_special_arcs:(-1))=.TRUE. DO arc=-n_special_arcs,-1 arcs_mask(arc)=(arcs_mask(arc).AND.(arcs_status(arc)/=periodic_bc_arc)) END DO CALL PlotNetwork2D(heads_tails=SSCNO_problem%problem%special_subgraph%heads_tails,& node_offset=n_special_nodes,node_coords=nodes_coords,& node_mask=nodes_mask,arc_mask=arcs_mask(-1:-n_special_arcs:(-1)),& arc_values=REAL(arcs_cuts(-1:-n_special_arcs:(-1)),r_wp),& resize_arcs=.TRUE.,color_arcs=.FALSE.,& arc_size_range=(/0.0_r_wp,HUGE(1.0_r_wp)/10/),& node_colorbar_format="5E1",arc_colorbar_format="5E1",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) CALL EndNetworkGraphics() arcs_mask=.TRUE._l_wp arcs_mask(:0)=.FALSE._l_wp DO arc=1,n_arcs arcs_mask(arc)=(arcs_mask(arc).AND.(arcs_status(arc)/=periodic_bc_arc)) END DO nodes_mask=.TRUE._l_wp nodes_mask(:0)=.FALSE._l_wp CALL InitNetworkGraphics(file="MC_MF_capacities."//TRIM(file_extension),file_type=TRIM(file_type),& x_label="Current flow \Huge{$\rightarrow$}",& page_size=(/5000,5000/),label_format="5N0",& color_table="RAIN",colorbar_position="Horizontal",& axis_labels_format=(/"NONE","NONE","NONE","NONE"/)) CALL PlotNetwork2D(heads_tails=heads_tails(:,1:n_arcs),& node_offset=-1,node_coords=nodes_coords(:,1:n_nodes),& node_mask=nodes_mask(1:n_nodes),arc_mask=arcs_mask(1:n_arcs),& arc_values=REAL(arcs_cuts_capacities(1:n_arcs),r_wp),& node_values=REAL(nodes_cuts(1:n_nodes),r_wp),& resize_nodes=.TRUE.,color_nodes=.TRUE.,& resize_arcs=.TRUE.,color_arcs=.TRUE.,& node_size_range=(/-HUGE(1.0_r_wp)/50,HUGE(1.0_r_wp)/20/),& node_colorbar_format="5F0",& arc_colorbar_format="5F1",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) nodes_mask(0:n_nodes)=.FALSE._l_wp nodes_mask(-n_special_nodes:-1)=.TRUE._l_wp arcs_mask(-1:-n_special_arcs:(-1))=.TRUE. DO arc=-n_special_arcs,-1 arcs_mask(arc)=(arcs_mask(arc).AND.(arcs_status(arc)/=periodic_bc_arc)) END DO CALL PlotNetwork2D(heads_tails=SSCNO_problem%problem%special_subgraph%heads_tails,& node_offset=n_special_nodes,node_coords=nodes_coords,& node_mask=nodes_mask,arc_mask=arcs_mask(-1:-n_special_arcs:(-1)),& arc_values=REAL(arcs_cuts_capacities(-1:-n_special_arcs:(-1)),r_wp),& resize_arcs=.TRUE.,color_arcs=.FALSE.,& arc_size_range=(/0.0_r_wp,HUGE(1.0_r_wp)/10/),& node_colorbar_format="5E1",arc_colorbar_format="5E1",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) CALL EndNetworkGraphics() END IF PointLoop:DO point=1,n_points WRITE(message_print_unit,*)"---->Doing instance, point #",instance,point IF(source_type==1)THEN sources_flow_capacity(1)=total_currents(point) sinks_flow_capacity(1)=-total_currents(point) ELSE IF(source_type==2)THEN sources_flow_capacity(1)=0.5_r_wp*total_voltages(point) sinks_flow_capacity(1)=-0.5_r_wp*total_voltages(point) END IF CALL AssignSuppliesDemands() CALL CallSSCNO(ElementalCosts=jj_iElementalCosts,& initialize_guess=(point==1)) nodes_potentials=nodes_potentials-MINVAL(nodes_potentials) IF(source_type==1)THEN applied_voltages(point,sample)=applied_voltages(point,sample)+& (nodes_potentials(sources(1))-nodes_potentials(sinks(1))) ELSE IF(source_type==2)THEN injected_currents(point,sample)=injected_currents(point,sample)+& 0 .5_r_wp*(supplies_demands(sources(1))-supplies_demands(sinks(1))) END IF IF(plot_samples)THEN IF(source_type==2)THEN arcs_mask=.TRUE. arcs_mask(:0)=.FALSE._l_wp DO arc=1,n_arcs arcs_mask(arc)=(arcs_mask(arc).AND.(arcs_status(arc)/=periodic_bc_arc)) END DO CALL InitNetworkGraphics(file="Dummy",file_type="CONS",& x_label="Current flow \Huge{$\rightarrow$}",& page_size=(/5000,5000/),label_format="5N0",& color_table="RAIN",colorbar_position="Horizontal",& axis_labels_format=(/"NONE","NONE","NONE","NONE"/)) min_arc_value=MINVAL(ABS(arcs_voltages(1:n_arcs))) max_arc_value=MAXVAL(ABS(arcs_voltages(1:n_arcs))) nodes_mask=.FALSE. CALL PlotNetwork2D(heads_tails=heads_tails(:,1:n_arcs),& node_offset=-1,node_coords=nodes_coords(:,1:n_nodes),& node_mask=nodes_mask(1:n_nodes),arc_mask=arcs_mask(1:n_arcs),& arc_values=ABS(arcs_voltages(1:n_arcs)),& resize_arcs=.TRUE.,color_arcs=.TRUE.,& arc_colorbar_format="5E1",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) IF(plot_potentials)THEN nodes_mask=.TRUE. nodes_mask(:0)=.FALSE._l_wp arcs_mask=.FALSE._l_wp CALL PlotNetwork2D(heads_tails=heads_tails(:,1:n_arcs),& node_offset=-1,node_coords=nodes_coords(:,1:n_nodes),& node_mask=nodes_mask(1:n_nodes),arc_mask=arcs_mask(1:n_arcs),& node_values=& nodes_potentials(1:n_nodes),& resize_nodes=.TRUE.,color_nodes=.TRUE.,& node_size_range=(/-HUGE(1.0_r_wp)/50,HUGE(1.0_r_wp)/20/),& node_colorbar_format="5E1",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) END IF nodes_mask(0:n_nodes)=.FALSE._l_wp nodes_mask(-n_special_nodes:-1)=.TRUE._l_wp arcs_mask(-1:-n_special_arcs:(-1))=.TRUE. DO arc=-n_special_arcs,-1 arcs_mask(arc)=(arcs_mask(arc).AND.(arcs_status(arc)/=periodic_bc_arc)) END DO CALL PlotNetwork2D(heads_tails=SSCNO_problem%problem%special_subgraph%heads_tails,& node_offset=n_special_nodes,node_coords=nodes_coords,& node_mask=nodes_mask,arc_mask=arcs_mask(-1:-n_special_arcs:(-1)),& arc_values=arcs_cuts_capacities(-1:-n_special_arcs:(-1)),& node_values=& nodes_potentials,& resize_nodes=.TRUE.,color_nodes=.FALSE.,& resize_arcs=.TRUE.,color_arcs=.FALSE.,& arc_size_range=(/0.0_r_wp,HUGE(1.0_r_wp)/10/),& node_size_range=(/-HUGE(1.0_r_wp)/10,HUGE(1.0_r_wp)/100/),& node_colorbar_format="5E1",arc_colorbar_format="5E1",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) CALL EndNetworkGraphics() IF(file_type/="CONS".AND.file_type/="XWIN")THEN WRITE(*,*)"Save this plot?" READ(*,*)save_plot IF(save_plot)THEN WRITE(*,*) WRITE(file_name,"(E10.2)")total_currents(point) file_name="V_"//TRIM(ADJUSTL(file_name))//"."//TRIM(file_extension) arcs_mask=.TRUE. arcs_mask(:0)=.FALSE._l_wp DO arc=1,n_arcs arcs_mask(arc)=(arcs_mask(arc).AND.(arcs_status(arc)/=periodic_bc_arc)) END DO CALL InitNetworkGraphics(file=TRIM(file_name),file_type=TRIM(file_type),& x_label="Current flow \Huge{$\rightarrow$}",& page_size=(/5000,5000/),label_format="5N0",& color_table="RAIN",colorbar_position="Horizontal",& axis_labels_format=(/"NONE","NONE","NONE","NONE"/)) min_arc_value=MINVAL(ABS(arcs_voltages(1:n_arcs))) max_arc_value=MAXVAL(ABS(arcs_voltages(1:n_arcs))) nodes_mask=.FALSE. CALL PlotNetwork2D(heads_tails=heads_tails(:,1:n_arcs),& node_offset=-1,node_coords=nodes_coords(:,1:n_nodes),& node_mask=nodes_mask(1:n_nodes),arc_mask=arcs_mask(1:n_arcs),& arc_values=ABS(arcs_voltages(1:n_arcs)),& resize_arcs=.TRUE.,color_arcs=.TRUE.,& arc_colorbar_format="5E1",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) IF(plot_potentials)THEN nodes_mask=.TRUE. nodes_mask(:0)=.FALSE._l_wp arcs_mask=.FALSE._l_wp CALL PlotNetwork2D(heads_tails=heads_tails(:,1:n_arcs),& node_offset=-1,node_coords=nodes_coords(:,1:n_nodes),& node_mask=nodes_mask(1:n_nodes),arc_mask=arcs_mask(1:n_arcs),& node_values=& nodes_potentials(1:n_nodes),& resize_nodes=.TRUE.,color_nodes=.TRUE.,& node_size_range=(/-HUGE(1.0_r_wp)/50,HUGE(1.0_r_wp)/20/),& node_colorbar_format="5E1",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) END IF nodes_mask(0:n_nodes)=.FALSE._l_wp nodes_mask(-n_special_nodes:-1)=.TRUE._l_wp arcs_mask(-1:-n_special_arcs:(-1))=.TRUE. DO arc=-n_special_arcs,-1 arcs_mask(arc)=(arcs_mask(arc).AND.(arcs_status(arc)/=periodic_bc_arc)) END DO CALL PlotNetwork2D(heads_tails=SSCNO_problem%problem%special_subgraph%heads_tails,& node_offset=n_special_nodes,node_coords=nodes_coords,& node_mask=nodes_mask,arc_mask=arcs_mask(-1:-n_special_arcs:(-1)),& arc_values=arcs_cuts_capacities(-1:-n_special_arcs:(-1)),& node_values=& nodes_potentials,& resize_nodes=.TRUE.,color_nodes=.FALSE.,& resize_arcs=.TRUE.,color_arcs=.FALSE.,& arc_size_range=(/0.0_r_wp,HUGE(1.0_r_wp)/10/),& node_size_range=(/-HUGE(1.0_r_wp)/10,HUGE(1.0_r_wp)/100/),& node_colorbar_format="5E1",arc_colorbar_format="5E1",& vector_type=0,axis=(/0.0,0.0,REAL(lengths+1)/)) CALL EndNetworkGraphics() END IF END IF ELSE END IF ELSE END IF END DO PointLoop END DO InstanceLoop applied_voltages(point,sample)=applied_voltages(point,sample)/n_instances injected_currents(point,sample)=injected_currents(point,sample)/n_instances Call DestroySSCNO CALL DestroyLatticeNetworkProblem(destroy_tree=.TRUE.) IF(ALLOCATED(nodes_cuts))THEN CALL RecordAllocation(n_elements=-INT(SIZE(nodes_cuts),KIND=i_wp),mold=1_i_byte) DEALLOCATE(nodes_cuts,STAT=alloc_status) END IF IF(ALLOCATED(arcs_cuts))THEN CALL RecordAllocation(n_elements=-INT(SIZE(arcs_cuts),KIND=i_wp),mold=1_i_byte) DEALLOCATE(arcs_cuts,STAT=alloc_status) END IF IF(ALLOCATED(arcs_cuts_capacities))THEN CALL RecordAllocation(n_elements=-INT(SIZE(arcs_cuts_capacities),KIND=i_wp),mold=1.0_r_wp) DEALLOCATE(arcs_cuts_capacities,STAT=alloc_status) END IF IF(ALLOCATED(cuts_capacities))THEN CALL RecordAllocation(n_elements=-INT(SIZE(cuts_capacities),KIND=i_wp),mold=1.0_r_wp) DEALLOCATE(cuts_capacities,STAT=alloc_status) END IF CALL RecordAllocation(n_elements=-INT(SIZE(special_heads_tails),KIND=i_wp),mold=1_i_wp,& caller="Saturation") DEALLOCATE(special_heads_tails,STAT=alloc_status) END DO SampleLoop END DO SourceType CALL StopTimer(100) WRITE(*,*)"Solving",n_points*n_samples*n_instances,& " NO problems took a total of (s) :",ReadTimer(100) WRITE(*,*)"___________________________________________" CALL ProfileSSCNO WRITE(*,*)"___________________________________________" CALL InitGraphics(file="V_I."//TRIM(file_extension),file_type=TRIM(file_type),& n_plots=n_samples+1,tick_labels="9F1",legend_position="UL",& plot_title=(/"IV curve for a 2D network"/),& x_label="$I$",y_label="$V$") IF(change_voltage.AND.change_current)THEN axis=(/MIN(MINVAL(injected_currents),MINVAL(total_currents)),& MAX(MAXVAL(injected_currents),MAXVAL(total_currents)),& MIN(MINVAL(applied_voltages),MINVAL(total_voltages)),& MAX(MAXVAL(applied_voltages),MAXVAL(total_voltages))/) ELSE IF(change_voltage)THEN axis=(/MINVAL(injected_currents),MAXVAL(injected_currents),& MINVAL(total_voltages),MAXVAL(total_voltages)/) ELSE IF(change_current)THEN axis=(/MINVAL(total_currents),MAXVAL(total_currents),& MINVAL(applied_voltages),MAXVAL(applied_voltages)/) END IF IF(change_current)THEN averaged_voltages=SUM(applied_voltages,DIM=2)/n_samples CALL Plot2D(x=total_currents,y=averaged_voltages,plot_spec="L.B",axis=axis) DO sample=1,n_samples CALL Plot2D(x=total_currents,y=applied_voltages(:,sample),plot_spec="SCR",axis=axis) END DO END IF IF(change_voltage)THEN averaged_currents=SUM(injected_currents,DIM=2)/n_samples CALL Plot2D(y=total_voltages,x=averaged_currents,plot_spec="L:R",axis=axis) DO sample=1,n_samples CALL Plot2D(y=total_voltages,x=injected_currents(:,sample),plot_spec="SSB",axis=axis) END DO END IF CALL EndGraphics() CALL EndProgram DEALLOCATE(total_currents,total_voltages,averaged_currents,averaged_voltages) DEALLOCATE(applied_voltages,injected_currents) END PROGRAM Saturation