MODULE SSCNO_Interface USE Precision USE Error_Handling USE System_Monitors USE Random_Numbers USE Initialization_Termination USE Network_Data_Structures USE Network_Spanning_Trees USE Network_Data_Types USE Network_Matrix_Operations USE Conjugate_Gradient USE Dual_Network_Solvers USE Dual_Line_Minimizers USE Dual_Newton_SSCNO USE TAUCS_Interface USE TAUCS_Constants USE SCOTCH_Interface USE SCOTCH_Constants USE CHACO_Interface USE CHACO_Constants IMPLICIT NONE PRIVATE PUBLIC::InitializeSSCNO,CreateSSCNO,DestroySSCNO,ProfileSSCNO,CallSSCNO TYPE(Dual_Network_Problem),POINTER,PUBLIC::SSCNO_problem=>NULL() REAL(KIND=r_wp),DIMENSION(:),ALLOCATABLE,TARGET,PUBLIC::& arcs_conductances,arcs_resistances,arcs_excess_voltages,arcs_excess_flows REAL(KIND=r_wp),DIMENSION(:),ALLOCATABLE,TARGET,PUBLIC::& nodes_excess_potentials,nodes_excess_flows CHARACTER,PUBLIC,SAVE::& SSCNO_lagrangian='DTN',SSCNO_linear_solver='Direct',SSCNO_preconditioner='Diagonal',& SSCNO_ls_method='Newton',SSCNO_factorization='Multifrontal',& SSCNO_ordering='SCOTCH',SSCNO_CHACO_method='KL',SSCNO_ST_mapping='Default' INTEGER,PUBLIC,SAVE::SSCNO_max_iterations=10,SSCNO_max_cg_iterations=-1,& SSCNO_max_ls_iterations=10,PCG_log_unit=0 INTEGER(KIND=i_wp),PUBLIC,SAVE::SSCNO_ST_height=0,SSCNO_ST_degree=2 REAL(KIND=r_wp),DIMENSION(4),PUBLIC,SAVE::SSCNO_tolerances=-1.0_r_wp REAL(KIND=r_wp),PUBLIC,SAVE::SSCNO_tolerance=1.0E-3_r_wp,& SSCNO_TAUCS_droptol=0.0_r_wp,SSCNO_TAUCS_subtree=0.0_r_wp,& typical_voltage=1.0_r_wp,typical_flow=1.0_r_wp,typical_resistance=1.0_r_wp,& low_conductance=-1.0_r_wp,high_conductance=-1.0_r_wp REAL(KIND=r_wp),DIMENSION(2),PUBLIC,SAVE::& SSCNO_forcing=(/0.1_r_wp,1.0_r_wp/),SSCNO_ls_forcing=(/0.1_r_wp,1.0_r_wp/) CHARACTER(LEN=10),PUBLIC,SAVE::SSCNO_TAUCS_ordering="metis" CHARACTER(LEN=SCOTCH_MAX_LEN),PUBLIC,SAVE::& SSCNO_SCOTCH_mapping=& "b{strat=m{low=(gx|h{pass=10}),asc=f{bal=0.1,type=b}x,rat=0.9,type=h,vert=50}}",& SSCNO_SCOTCH_ordering=& "n{sep=/(vert>100)?m{low=h,asc=f{bal=0.1},vert=25};,ole=h{halo=y,cmin=32,cmax=100,frat=0.08},ose=g}" CHARACTER(LEN=25),PUBLIC,SAVE::SSCNO_TAUCS_log_file="" LOGICAL,SAVE,PUBLIC::debug_SSCNO=.FALSE.,profile_SSCNO=.TRUE.,& SSCNO_warm_starts=.FALSE. LOGICAL,SAVE::allocated_problem=.FALSE. NAMELIST/SSCNOOptions/SSCNO_lagrangian,SSCNO_linear_solver,SSCNO_ls_method,& SSCNO_factorization,SSCNO_ordering,SSCNO_SCOTCH_mapping,SSCNO_SCOTCH_ordering,& SSCNO_max_iterations,SSCNO_max_cg_iterations,SSCNO_max_ls_iterations,& SSCNO_tolerance,SSCNO_tolerances,SSCNO_forcing,SSCNO_ls_forcing,& PCG_log_unit,SSCNO_preconditioner,SSCNO_ST_height,SSCNO_ST_degree,SSCNO_ST_mapping,& SSCNO_TAUCS_droptol,SSCNO_TAUCS_subtree,SSCNO_TAUCS_ordering,SSCNO_CHACO_method,& typical_voltage,typical_flow,typical_resistance,low_conductance,high_conductance,& debug_SSCNO,profile_SSCNO,SSCNO_warm_starts,SSCNO_TAUCS_log_file CONTAINS SUBROUTINE InitializeSSCNO() REWIND(UNIT=program_options_unit) READ(UNIT=program_options_unit,NML=SSCNOOptions,IOSTAT=error_status) IF(error_status/=0)& CALL CriticalError(message="NAMELIST SSCNOOptions was not read successfully"//& " from file "//TRIM(options_file),caller="InitializeSSCNO") WRITE(UNIT=message_log_unit,NML=SSCNOOptions) WRITE(UNIT=message_log_unit,FMT=*) END SUBROUTINE InitializeSSCNO SUBROUTINE CreateSSCNO() IMPLICIT NONE TYPE(Dual_Network_Problem),POINTER::dual_problem TYPE(Directed_Graph),POINTER::graph TYPE(Spanning_Tree),POINTER::tree TYPE(Network_SC_Cost),POINTER::cost_function TYPE(Network_Arrays),POINTER::arrays TYPE(Dual_Network_Preconditioner),POINTER::preconditioner TYPE(Dual_Network_Factorization),POINTER::factorization INTEGER::alloc_status,precond_method,mapping_method,global_method,local_method IF(.NOT.ASSOCIATED(SSCNO_problem))THEN allocated_problem=.TRUE. ALLOCATE(SSCNO_problem,STAT=alloc_status) END IF dual_problem=>SSCNO_problem IF(.NOT.ALLOCATED(arcs_conductances))THEN ALLOCATE(arcs_conductances(-n_special_arcs:n_arcs),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1.0_r_wp,& caller="InitializeSSCNO",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(arcs_excess_voltages))THEN ALLOCATE(arcs_excess_voltages(-n_special_arcs:n_arcs),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1.0_r_wp,& caller="InitializeSSCNO",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(nodes_excess_flows))THEN ALLOCATE(nodes_excess_flows(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_wp,& caller="InitializeSSCNO",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(nodes_excess_potentials))THEN ALLOCATE(nodes_excess_potentials(-n_special_nodes:n_nodes),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_wp,& caller="InitializeSSCNO",alloc_status=alloc_status) END IF IF(.NOT.ALLOCATED(arcs_excess_flows))THEN ALLOCATE(arcs_excess_flows(-n_special_arcs:n_arcs),STAT=alloc_status) CALL RecordAllocation(n_elements=(n_arcs)-(-n_special_arcs)+1,mold=1.0_r_wp,& caller="InitializeSSCNO",alloc_status=alloc_status) END IF dual_problem%PIN=1 graph=>dual_problem%problem%graph cost_function=>dual_problem%problem%cost_function arrays=>dual_problem%problem%arrays graph%heads_tails=>heads_tails graph%n_special_nodes=n_special_nodes graph%n_special_arcs=n_special_arcs graph%n_nodes=n_nodes graph%n_arcs=n_arcs dual_problem%n_fixed_potentials=n_fixed_potential_nodes arrays%supplies_demands=>supplies_demands arrays%nodes_potentials=>nodes_potentials arrays%nodes_excess_potentials=>nodes_excess_potentials arrays%nodes_excess_flows=>nodes_excess_flows arrays%arcs_conductances=>arcs_conductances arrays%arcs_voltages=>arcs_voltages arrays%arcs_flows=>arcs_flows arrays%arcs_excess_voltages=>arcs_excess_voltages arrays%arcs_excess_flows=>arcs_excess_flows arrays%nodes_coordinates=>nodes_coords SELECT CASE(SSCNO_lagrangian) CASE('D','d','N','n') dual_problem%method=SSCNO_TDN CASE DEFAULT dual_problem%method=SSCNO_TSQP ENDSELECT SELECT CASE(SSCNO_linear_solver) CASE('I','i','C','c','P','p') dual_problem%system%method=dual_network_pcg CASE DEFAULT dual_problem%system%method=dual_network_LLt ENDSELECT UsePCG:IF(dual_problem%system%method==dual_network_pcg)THEN WRITE(*,*)"Using PCG as a linear solver" preconditioner=>dual_problem%system%preconditioner dual_problem%system%allocated_solver=.TRUE. ALLOCATE(dual_problem%system%system%solver) dual_problem%system%system%solver%log_unit=PCG_log_unit SELECT CASE(SSCNO_preconditioner) CASE('N','n') dual_problem%system%preconditioner%method=no_precond CASE('B','b') precond_method=MST_BHBt_precond CASE('R','r') precond_method=MST_LDLt_precond CASE('S','s') precond_method=ST_precond CASE('P','p') precond_method=ST_MST_precond CASE('V','v') precond_method=TAUCS_MWB_precond CASE('T','t') precond_method=TAUCS_LLt_precond CASE DEFAULT precond_method=diagonal_precond ENDSELECT preconditioner%method=precond_method PartitioningMethod:IF((precond_method==ST_precond).OR.& (precond_method==ST_MST_precond))THEN preconditioner%allocated_support_tree=.TRUE. ALLOCATE(preconditioner%support_tree) preconditioner%support_tree%height=SSCNO_ST_height preconditioner%support_tree%degree=SSCNO_ST_degree SELECT CASE(SSCNO_ST_mapping) CASE('M','m') mapping_method=Matching_ordering CASE('S','s') mapping_method=SCOTCH_ordering CASE('C','c') mapping_method=CHACO_ordering CASE DEFAULT mapping_method=Default_ordering END SELECT preconditioner%support_tree%mapping_method=mapping_method IF(mapping_method==CHACO_ordering)THEN ALLOCATE(preconditioner%support_tree%chaco_mapping) SELECT CASE(SSCNO_CHACO_method) CASE('Multilevel') global_method=CHACO_multilevel local_method=CHACO_local_KL CASE('Inertial') IF(.NOT.ASSOCIATED(arrays%nodes_coordinates))& arrays%nodes_coordinates=>nodes_coords global_method=CHACO_inertial local_method=CHACO_local_none CASE('Spectral') global_method=CHACO_spectral local_method=CHACO_local_none CASE DEFAULT global_method=CHACO_linear local_method=CHACO_local_KL END SELECT preconditioner%support_tree%chaco_mapping%strategy%global_method=global_method preconditioner%support_tree%chaco_mapping%strategy%local_method=local_method ELSE IF(mapping_method==SCOTCH_ordering)THEN ALLOCATE(preconditioner%support_tree%scotch_mapping) preconditioner%support_tree%scotch_mapping%scotch_mapping%& mapping_string=TRIM(SSCNO_SCOTCH_mapping) END IF END IF PartitioningMethod IF((precond_method==TAUCS_LLt_precond).OR.& (precond_method==TAUCS_MWB_precond))THEN preconditioner%allocated_taucs=.TRUE. ALLOCATE(preconditioner%taucs_preconditioner) preconditioner%taucs_preconditioner%ordering_method=TRIM(SSCNO_TAUCS_ordering) preconditioner%taucs_preconditioner%log_file=TRIM(SSCNO_TAUCS_log_file) END IF ELSE UsePCG WRITE(*,*)"Using Cholesky factorization as a linear solver" factorization=>dual_problem%system%factorization factorization%allocated_factorization=.TRUE. ALLOCATE(factorization%taucs_factorization) factorization%taucs_factorization%log_file=TRIM(SSCNO_TAUCS_log_file) SELECT CASE(SSCNO_factorization) CASE('I','i') factorization%factorization=TAUCS_fact_llt CASE('L','l') factorization%factorization=TAUCS_fact_ll CASE DEFAULT factorization%factorization=TAUCS_fact_mf END SELECT SELECT CASE(SSCNO_ordering) CASE('N','n') factorization%ordering=no_fact_ordering CASE('T','t') factorization%ordering=TAUCS_fact_ordering factorization%taucs_factorization%ordering_method=TRIM(SSCNO_TAUCS_ordering) CASE DEFAULT factorization%allocated_scotch_ordering=.TRUE. ALLOCATE(factorization%scotch_ordering) factorization%scotch_ordering%scotch_ordering%& ordering_string=TRIM(SSCNO_SCOTCH_ordering) factorization%taucs_factorization%ordering_method="none" END SELECT END IF UsePCG IF(profile_SSCNO)THEN dual_problem%timers=SSCNO_Timers(5,6,7,8,9) dual_problem%system%timers=Dual_System_Timers(CG_Timers(15,16,17,18),& 11 ,12,13,14,31,32,33) dual_problem%line_search%timers=Dual_Search_Timers(21,22,23) END IF debug_dual_search=debug_SSCNO dual_problem%log_SSCNO=debug_SSCNO CALL InitializeDualSSCNO(dual_problem=dual_problem) END SUBROUTINE CreateSSCNO SUBROUTINE DestroySSCNO IMPLICIT NONE TYPE(Network_SC_Cost),POINTER::cost_function INTEGER::alloc_status CALL DestroyDualSSCNO(dual_problem=SSCNO_problem) IF(ALLOCATED(arcs_conductances))THEN CALL RecordAllocation(n_elements=-INT(SIZE(arcs_conductances),KIND=i_wp),mold=1.0_r_wp) DEALLOCATE(arcs_conductances,STAT=alloc_status) END IF IF(ALLOCATED(arcs_excess_voltages))THEN CALL RecordAllocation(n_elements=-INT(SIZE(arcs_excess_voltages),KIND=i_wp),mold=1.0_r_wp) DEALLOCATE(arcs_excess_voltages,STAT=alloc_status) END IF IF(ALLOCATED(nodes_excess_flows))THEN CALL RecordAllocation(n_elements=-INT(SIZE(nodes_excess_flows),KIND=i_wp),mold=1.0_r_wp) DEALLOCATE(nodes_excess_flows,STAT=alloc_status) END IF IF(ALLOCATED(nodes_excess_potentials))THEN CALL RecordAllocation(n_elements=-INT(SIZE(nodes_excess_potentials),KIND=i_wp),mold=1.0_r_wp) DEALLOCATE(nodes_excess_potentials,STAT=alloc_status) END IF IF(ALLOCATED(arcs_excess_flows))THEN CALL RecordAllocation(n_elements=-INT(SIZE(arcs_excess_flows),KIND=i_wp),mold=1.0_r_wp) DEALLOCATE(arcs_excess_flows,STAT=alloc_status) END IF IF(allocated_problem)THEN DEALLOCATE(SSCNO_problem,STAT=alloc_status) END IF END SUBROUTINE DestroySSCNO SUBROUTINE ProfileSSCNO() INTEGER::linear_solver SELECT CASE(SSCNO_linear_solver) CASE('I','i','C','c','P','p') linear_solver=dual_network_pcg CASE DEFAULT linear_solver=dual_network_LLt ENDSELECT IF(profile_SSCNO)THEN WRITE(message_print_unit,*)"Profiling timing results:" WRITE(message_print_unit,*)"______________________________________________:" WRITE(message_print_unit,*)"SSCNO initialization :",& ReadTimer(5) IF(linear_solver==dual_network_pcg)THEN WRITE(message_print_unit,*)"---> Preconditioner initialization :",& ReadTimer(11) ELSE WRITE(message_print_unit,*)"---> Fill-reducing reordering in SCOTCH :",& ReadTimer(31) END IF WRITE(message_print_unit,*)"TCGN iterations :" WRITE(message_print_unit,*)"--->Array updates :",& ReadTimer(6) WRITE(message_print_unit,*)"--->Conductance calculation :",& ReadTimer(7) WRITE(message_print_unit,*)"--->Line Search :",& ReadTimer(9) WRITE(message_print_unit,*)"------>Root bracketing :",& ReadTimer(22) WRITE(message_print_unit,*)"------>Root finding solver :",& ReadTimer(23) WRITE(message_print_unit,*)"--------->Cost/derivatives evaluation :",& ReadTimer(21) WRITE(message_print_unit,*)"--->Solving Newton's linear system :",& ReadTimer(8) IF(linear_solver==dual_network_pcg)THEN WRITE(message_print_unit,*)"------>Preconditioner construction :" WRITE(message_print_unit,*)"--------->Combinatorial reoptimization :",& ReadTimer(12) WRITE(message_print_unit,*)"--------->Creation (factorization) :",& ReadTimer(13) WRITE(message_print_unit,*)"------>PCG iteration :",& ReadTimer(14) WRITE(message_print_unit,*)"--------->Preconditioning :",& ReadTimer(15) WRITE(message_print_unit,*)"--------->Dot products :",& ReadTimer(16) WRITE(message_print_unit,*)"--------->Vector updates :",& ReadTimer(17) WRITE(message_print_unit,*)"--------->Matrix-vector products :",& ReadTimer(18) ELSE WRITE(message_print_unit,*)"------>LLT numerical factorization :",& ReadTimer(32) WRITE(message_print_unit,*)"------>LLT triangular solve :",& ReadTimer(33) END IF WRITE(message_print_unit,*)"______________________________________________:" END IF END SUBROUTINE ProfileSSCNO SUBROUTINE CallSSCNO(ElementalCosts,initialize_guess) IMPLICIT NONE INTERFACE SUBROUTINE ElementalCosts(cost_function,arguments_status,tolerance,& arcs_indices,arcs_flows,arcs_voltages,arcs_resistances,arcs_costs) USE Precision USE Network_Data_Types TYPE(Network_SC_Cost),INTENT(INOUT),OPTIONAL::cost_function CHARACTER(LEN=4),INTENT(IN)::arguments_status REAL(KIND=r_wp),INTENT(IN)::tolerance INTEGER(KIND=i_wp),DIMENSION(2),INTENT(IN)::arcs_indices REAL(KIND=r_wp),DIMENSION(arcs_indices(1):),INTENT(INOUT),OPTIONAL::& arcs_flows,arcs_voltages,arcs_resistances,arcs_costs END SUBROUTINE ElementalCosts END INTERFACE LOGICAL,INTENT(IN),OPTIONAL::initialize_guess TYPE(Dual_Network_Problem),POINTER::dual_problem TYPE(Network_SC_Cost),POINTER::cost_function TYPE(TAUCS_Network_Preconditioner),POINTER::taucs_preconditioner INTEGER::alloc_status,linear_solver,precond_method LOGICAL::random_guess random_guess=.FALSE. IF(PRESENT(initialize_guess))random_guess=initialize_guess dual_problem=>SSCNO_problem cost_function=>dual_problem%problem%cost_function linear_solver=dual_problem%system%method dual_problem%max_iterations=SSCNO_max_iterations dual_problem%error_tolerances=SSCNO_Error_Tolerances(SSCNO_tolerance,& SSCNO_tolerances(1),SSCNO_tolerances(2),SSCNO_tolerances(3),SSCNO_tolerances(4)) dual_problem%forcing_prefactor=SSCNO_forcing(1) dual_problem%forcing_exponent=SSCNO_forcing(2) dual_problem%ls_tolerance_prefactor=SSCNO_ls_forcing(1) dual_problem%ls_tolerance_exponent=SSCNO_ls_forcing(2) dual_problem%line_search%root_finder%max_iterations=SSCNO_max_ls_iterations dual_problem%line_search%root_finder%method=SSCNO_ls_method dual_problem%warm_starts=SSCNO_warm_starts cost_function%typical_voltage=typical_voltage cost_function%typical_flow=typical_flow cost_function%typical_resistance=typical_resistance cost_function%low_conductance=low_conductance cost_function%high_conductance=high_conductance UsePCG:IF(linear_solver==dual_network_pcg)THEN dual_problem%system%system%solver%max_iterations=SSCNO_max_cg_iterations precond_method=dual_problem%system%preconditioner%method IF((precond_method==TAUCS_LLt_precond).OR.& (precond_method==TAUCS_MWB_precond))THEN taucs_preconditioner=>dual_problem%system%preconditioner%taucs_preconditioner taucs_preconditioner%droptol=SSCNO_TAUCS_droptol taucs_preconditioner%subtree_ratio=SSCNO_TAUCS_subtree END IF END IF UsePCG IF(random_guess)THEN InitialGuess:IF(dual_problem%method==SSCNO_TDN)THEN CALL RandomUniform(nodes_potentials(-(n_special_nodes-n_fixed_potential_nodes):n_nodes),& range=(/-typical_voltage,typical_voltage/)) nodes_potentials(0)=0.0_r_wp arcs_flows=0.0_r_wp ELSE InitialGuess CALL RandomUniform(arcs_flows,& range=(/-typical_flow,typical_flow/)) arcs_flows(0)=0.0_r_wp nodes_potentials(-(n_special_nodes-n_fixed_potential_nodes):n_nodes)=0.0_r_wp END IF InitialGuess END IF IF(n_fixed_potential_nodes>0)THEN nodes_potentials(-n_special_nodes:-(n_special_nodes-n_fixed_potential_nodes))=& supplies_demands(-n_special_nodes:-(n_special_nodes-n_fixed_potential_nodes)-1) supplies_demands(-n_special_nodes:-(n_special_nodes-n_fixed_potential_nodes)-1)=0.0_r_wp END IF CALL SolveDualSSCNO_TCGN(dual_problem=dual_problem,ElementalCosts=ElementalCosts) END SUBROUTINE CallSSCNO END MODULE SSCNO_Interface