MODULE TAUCS_Constants
      USE Precision
      IMPLICIT NONE
      PUBLIC
      
      INTEGER,PARAMETER::i_taucs=i_word,r_taucs=r_dp
      INTEGER,PARAMETER::TAUCS_CCSDIM=6
      
      END MODULE TAUCS_Constants
      
      
      MODULE CCS_Laplacians
      USE Precision
      USE TAUCS_Constants
      USE Error_Handling
      USE System_Monitors
      USE Network_Matrix_Operations
      USE Graph_Algorithms
      IMPLICIT NONE
      PUBLIC::CreateCCSLaplacian,ComputeCCSLaplacian
      
      PRIVATE
      
      CONTAINS
      
      
      SUBROUTINE CreateCCSLaplacian(indexing_offset,arc_offset,node_offset,&
      heads_tails,neighbours,my_neighbours,incident_arcs,&
      nodes_reordering,nodes_renumbering,n_edges)
      IMPLICIT NONE
      INTEGER(KIND=i_wp),INTENT(IN)::node_offset,arc_offset,indexing_offset
      INTEGER(KIND=i_wp),DIMENSION(:,-arc_offset:),INTENT(IN)::heads_tails
      INTEGER(KIND=i_taucs),DIMENSION(indexing_offset:),INTENT(OUT)::&
      neighbours,my_neighbours,incident_arcs
      INTEGER(KIND=i_wp),DIMENSION(-node_offset:),INTENT(IN),OPTIONAL::&
      nodes_reordering,nodes_renumbering
      INTEGER(KIND=i_wp),INTENT(OUT),OPTIONAL::n_edges
      
      INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::nodes_degrees
      INTEGER(KIND=i_wp)::n_nodes,n_special_nodes,n_arcs,n_special_arcs,&
      n_vertices
      INTEGER(KIND=i_wp)::index,node,arc,head,tail,special_node,special_arc,&
      degree,head_neighbour,tail_neighbour,&
      node_shift,arc_shift,arc_counter,node_counter
      INTEGER::alloc_status,status
      LOGICAL::graph_OK,reorder_nodes
      
      n_special_nodes=node_offset
      n_nodes=INT(SIZE(my_neighbours),KIND=i_wp)-2-n_special_nodes
      
      n_special_arcs=-INT(LBOUND(heads_tails,DIM=2),KIND=i_wp)
      n_arcs=INT(UBOUND(heads_tails,DIM=2),KIND=i_wp)
      n_vertices=n_nodes+n_special_nodes+1
      node_shift=n_special_nodes+indexing_offset
      arc_shift=n_special_arcs+indexing_offset
      
      reorder_nodes=(PRESENT(nodes_reordering).AND.PRESENT(nodes_renumbering))
      
      
      
      
      IF(.NOT.ALLOCATED(nodes_degrees))THEN
      ALLOCATE(nodes_degrees(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1_i_wp,&
      caller="CreateCCSLaplacian",alloc_status=alloc_status)
      END IF
      
      
      nodes_degrees=1
      DO arc=-n_special_arcs,n_arcs
      head=heads_tails(1,arc)
      tail=heads_tails(2,arc)
      
      IF(reorder_nodes)THEN
      head=nodes_renumbering(head)
      tail=nodes_renumbering(tail)
      END IF
      
      
      IF(head<tail)THEN
      nodes_degrees(head)=nodes_degrees(head)+1
      ELSE IF(tail<head)THEN
      nodes_degrees(tail)=nodes_degrees(tail)+1
      END IF
      END DO
      
      arc_counter=indexing_offset
      DO node=-n_special_nodes,n_nodes
      my_neighbours(node+node_shift)=arc_counter
      arc_counter=arc_counter+nodes_degrees(node)
      END DO
      my_neighbours(n_vertices+indexing_offset)=arc_counter
      IF(PRESENT(n_edges))n_edges=arc_counter-1
      
      
      DO arc=n_arcs,-n_special_arcs,-1
      
      head=heads_tails(1,arc)
      tail=heads_tails(2,arc)
      
      IF(reorder_nodes)THEN
      head=nodes_renumbering(head)
      tail=nodes_renumbering(tail)
      END IF
      
      IF(head<tail)THEN
      nodes_degrees(head)=nodes_degrees(head)-1
      head_neighbour=my_neighbours(head+node_shift)+nodes_degrees(head)
      neighbours(head_neighbour)=tail+node_shift
      incident_arcs(head_neighbour)=arc
      ELSE IF(tail<head)THEN
      nodes_degrees(tail)=nodes_degrees(tail)-1
      tail_neighbour=my_neighbours(tail+node_shift)+nodes_degrees(tail)
      neighbours(tail_neighbour)=head+node_shift
      incident_arcs(tail_neighbour)=arc
      END IF
      
      END DO
      
      
      DO node=-n_special_nodes,n_nodes
      neighbours(my_neighbours(node+node_shift))=node+node_shift
      incident_arcs(my_neighbours(node+node_shift))=0
      END DO
      
      
      IF(ALLOCATED(nodes_degrees))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(nodes_degrees),KIND=i_wp),mold=1_i_wp)
      DEALLOCATE(nodes_degrees,STAT=alloc_status)
      END IF
      
      END SUBROUTINE CreateCCSLaplacian
      
      
      
      SUBROUTINE ComputeCCSLaplacian(indexing_offset,arc_offset,&
      arcs_conductances,edges_conductances,&
      neighbours,my_neighbours,incident_arcs)
      IMPLICIT NONE
      INTEGER(KIND=i_wp),INTENT(IN)::indexing_offset,arc_offset
      INTEGER(KIND=i_taucs),DIMENSION(indexing_offset:),INTENT(IN)::&
      neighbours,my_neighbours,incident_arcs
      REAL(KIND=r_wp),DIMENSION(-arc_offset:),INTENT(IN)::arcs_conductances
      REAL(KIND=r_taucs),DIMENSION(indexing_offset:),INTENT(OUT)::edges_conductances
      
      INTEGER(KIND=i_wp)::n_vertices
      INTEGER(KIND=i_wp)::node,arc,vertex,edge,first,second
      
      n_vertices=INT(SIZE(my_neighbours),KIND=i_wp)-1
      
      DO vertex=indexing_offset,n_vertices+indexing_offset-1
      edges_conductances(my_neighbours(vertex))=0.0_r_taucs
      END DO
      
      DO vertex=indexing_offset,n_vertices+indexing_offset-1
      first=my_neighbours(vertex)
      DO edge=my_neighbours(vertex)+1,my_neighbours(vertex+1)-1
      arc=incident_arcs(edge)
      edges_conductances(edge)=-arcs_conductances(arc)
      second=my_neighbours(neighbours(edge))
      
      edges_conductances(first)=edges_conductances(first)+arcs_conductances(arc)
      edges_conductances(second)=edges_conductances(second)+arcs_conductances(arc)
      END DO
      END DO
      
      END SUBROUTINE ComputeCCSLaplacian
      
      
      
      END MODULE CCS_Laplacians
      
      
      
      MODULE TAUCS_Interface
      USE Precision
      USE TAUCS_Constants
      USE Error_Handling
      USE System_Monitors
      USE Vector_Operations
      USE Graph_Algorithms
      USE Network_Data_Types
      USE CCS_Laplacians
      IMPLICIT NONE
      PRIVATE
      PUBLIC::TAUCS_InitializePreconditioner,TAUCS_DestroyPreconditioner,&
      TAUCS_CreatePreconditioner,TAUCS_FreePreconditioner,TAUCS_ApplyPreconditioner
      PUBLIC::i_taucs,r_taucs
      
      
      ML_EXTERNAL::taucs_ccs_write_ijv,taucs_printf,&
      taucs_ccs_assemble,taucs_logfile,&
      taucs_ccs_free,taucs_supernodal_factor_free,&
      taucs_ccs_solve_llt,taucs_supernodal_solve_llt,&
      taucs_ccs_factor_llt,taucs_ccs_factor_llt_mf,taucs_ccs_factor_llt_ll,&
      taucs_ccs_order,taucs_free_order,taucs_ccs_permute_symmetrically,&
      taucs_amwb_preconditioner_create,taucs_supernodal_factor_free_numeric,&
      taucs_ccs_factor_llt_symbolic,taucs_ccs_factor_llt_numeric,&
      taucs_permute_vector_sp_to_sp,taucs_permute_vector_sp_to_dp,&
      taucs_permute_vector_dp_to_dp,taucs_permute_vector_dp_to_sp
      INTEGER::taucs_ccs_write_ijv,taucs_ccs_factor_llt_symbolic,taucs_ccs_factor_llt_numeric,&
      taucs_ccs_factor_llt,taucs_ccs_factor_llt_mf,taucs_ccs_factor_llt_ll,&
      taucs_ccs_solve_llt,taucs_supernodal_solve_llt,&
      taucs_amwb_preconditioner_create,taucs_ccs_permute_symmetrically
      INTERFACE taucs_permute_vector
      SUBROUTINE taucs_permute_vector_sp_to_sp(perm,x,perm_x,n)
      USE Precision
      INTEGER::perm,n
      REAL(KIND=r_sp)::x,perm_x
      END SUBROUTINE
      SUBROUTINE taucs_permute_vector_sp_to_dp(perm,x,perm_x,n)
      USE Precision
      INTEGER::perm,n
      REAL(KIND=r_sp)::x
      REAL(KIND=r_dp)::perm_x
      END SUBROUTINE
      SUBROUTINE taucs_permute_vector_dp_to_dp(perm,x,perm_x,n)
      USE Precision
      INTEGER::perm,n
      REAL(KIND=r_dp)::x,perm_x
      END SUBROUTINE
      SUBROUTINE taucs_permute_vector_dp_to_sp(perm,x,perm_x,n)
      USE Precision
      INTEGER::perm,n
      REAL(KIND=r_dp)::x
      REAL(KIND=r_sp)::perm_x
      END SUBROUTINE
      END INTERFACE taucs_permute_vector
      
      PUBLIC::TAUCS_CCS_Matrix,TAUCS_Network_Preconditioner
      INTEGER,PARAMETER,PUBLIC::TAUCS_LLt=1,TAUCS_MWB=2
      
      INTEGER,PARAMETER,PUBLIC::TAUCS_incomplete_fact=1,TAUCS_complete_fact_ll=2,TAUCS_complete_fact_mf=3
      
      
      TYPE TAUCS_CCS_Matrix
      PRIVATE
      INTEGER(KIND=i_taucs),DIMENSION(TAUCS_CCSDIM)::matrix
      INTEGER(KIND=i_taucs)::n_vertices=0,n_edges=0,n_eliminated=0
      INTEGER(KIND=i_taucs),DIMENSION(:),POINTER::col_ptr,row_ind,incident,eliminated
      
      REAL(KIND=r_taucs),DIMENSION(:),POINTER::val
      ENDTYPE
      
      TYPE TAUCS_Network_Preconditioner
      TYPE(Directed_Graph),POINTER::graph
      TYPE(TAUCS_CCS_Matrix)::ccs_matrix
      TYPE(Fill_Reducing_Ordering),POINTER::graph_ordering
      CHARACTER(LEN=25)::log_file=""
      INTEGER(KIND=i_wp),DIMENSION(:),POINTER::eliminated_nodes
      INTEGER::ordered_matrix=0,ccs_factor=0,ccs_mwb=0,perm=0,invperm=0
      
      INTEGER::precond_method=TAUCS_LLt,factorization=TAUCS_incomplete_fact
      REAL::subtree_ratio=0.0
      REAL(KIND=r_taucs)::droptol=0.0_r_taucs,subgraphs=1.0_r_taucs,&
      diagonal_conductance=1.0_r_taucs
      CHARACTER(LEN=10)::ordering_method="none"
      LOGICAL::reorder_graph=.FALSE.,allocated_eliminated=.FALSE.,&
      allocated_xb_buffers=.FALSE.
      REAL(KIND=r_wp),DIMENSION(:),POINTER::conductances=>NULL()
      REAL(KIND=r_taucs),DIMENSION(:),POINTER::flows=>NULL(),potentials=>NULL()
      
      REAL(KIND=r_wp),DIMENSION(:),POINTER::excess_flows=>NULL(),excess_potentials=>NULL()
      ENDTYPE
      
      CONTAINS
      
      
      SUBROUTINE TAUCS_InitializePreconditioner(taucs_preconditioner)
      IMPLICIT NONE
      TYPE(TAUCS_Network_Preconditioner),INTENT(INOUT),TARGET::taucs_preconditioner
      
      
      TYPE(Directed_Graph),POINTER::graph
      INTEGER::alloc_status,result
      INTEGER(KIND=i_taucs)::n_vertices,n_edges,second_node
      INTEGER(KIND=i_wp)::n_special_nodes,n_nodes,n_special_arcs,n_arcs
      LOGICAL::valid_graph
      
      graph=>taucs_preconditioner%graph
      valid_graph=.TRUE.
      IF(.NOT.ASSOCIATED(graph))THEN
      valid_graph=.FALSE.
      ELSE
      IF(.NOT.ASSOCIATED(graph%heads_tails))valid_graph=.FALSE.
      END IF
      IF(.NOT.ASSOCIATED(taucs_preconditioner%conductances))valid_graph=.FALSE.
      IF(.NOT.valid_graph)THEN
      CALL CriticalError(message="Incomplete or no graph passed to TAUCS library",&
      caller="TAUCS_InitializePreconditioner")
      RETURN
      END IF
      
      n_special_nodes=graph%n_special_nodes
      n_nodes=graph%n_nodes
      n_special_arcs=graph%n_special_arcs
      n_arcs=graph%n_arcs
      
      n_vertices=n_nodes+n_special_nodes+1
      n_edges=n_vertices+n_arcs+n_special_arcs+1
      
      
      taucs_preconditioner%reorder_graph=.FALSE.
      IF(ASSOCIATED(taucs_preconditioner%graph_ordering))THEN
      IF(ASSOCIATED(taucs_preconditioner%graph_ordering%nodes_reordering).AND.ASSOCIATED(taucs_preconditioner%graph_ordering%nodes&
     &_renumbering))THEN
      taucs_preconditioner%reorder_graph=.TRUE.
      END IF
      END IF
      IF(taucs_preconditioner%reorder_graph)taucs_preconditioner%ordering_method="none"
      
      
      IF((r_taucs/=r_wp).OR.(taucs_preconditioner%ordering_method/="none").OR.(taucs_preconditioner%reorder_graph))THEN
      WRITE(message_log_unit,*)"Allocating x and b buffers in TAUCS"
      taucs_preconditioner%allocated_xb_buffers=.TRUE.
      
      
      IF(.NOT.ASSOCIATED(taucs_preconditioner%flows))THEN
      ALLOCATE(taucs_preconditioner%flows(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_taucs,&
      caller="TAUCS_InitializePreconditioner",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.ASSOCIATED(taucs_preconditioner%potentials))THEN
      ALLOCATE(taucs_preconditioner%potentials(-n_special_nodes:n_nodes),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_nodes)-(-n_special_nodes)+1,mold=1.0_r_taucs,&
      caller="TAUCS_InitializePreconditioner",alloc_status=alloc_status)
      END IF
      ELSE
      CALL VectorPointerCopy(source=taucs_preconditioner%excess_potentials,target=taucs_preconditioner%potentials)
      CALL VectorPointerCopy(source=taucs_preconditioner%excess_flows,target=taucs_preconditioner%potentials)
      END IF
      
      
      IF(.NOT.ASSOCIATED(taucs_preconditioner%ccs_matrix%col_ptr))THEN
      ALLOCATE(taucs_preconditioner%ccs_matrix%col_ptr(0:n_vertices),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_vertices)-(0)+1,mold=1_i_taucs,&
      caller="TAUCS_InitializePreconditioner",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ASSOCIATED(taucs_preconditioner%ccs_matrix%row_ind))THEN
      ALLOCATE(taucs_preconditioner%ccs_matrix%row_ind(0:n_edges-1),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_edges-1)-(0)+1,mold=1_i_taucs,&
      caller="TAUCS_InitializePreconditioner",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ASSOCIATED(taucs_preconditioner%ccs_matrix%incident))THEN
      ALLOCATE(taucs_preconditioner%ccs_matrix%incident(0:n_edges-1),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_edges-1)-(0)+1,mold=1_i_taucs,&
      caller="TAUCS_InitializePreconditioner",alloc_status=alloc_status)
      END IF
      
      IF(.NOT.ASSOCIATED(taucs_preconditioner%ccs_matrix%val))THEN
      ALLOCATE(taucs_preconditioner%ccs_matrix%val(0:n_edges-1),STAT=alloc_status)
      CALL RecordAllocation(n_elements=(n_edges-1)-(0)+1,mold=1.0_r_taucs,&
      caller="TAUCS_InitializePreconditioner",alloc_status=alloc_status)
      END IF
      
      
      IF(.NOT.taucs_preconditioner%reorder_graph)THEN
      CALL CreateCCSLaplacian(indexing_offset=0,&
      arc_offset=n_special_arcs,node_offset=n_special_nodes,&
      heads_tails=graph%heads_tails,neighbours=taucs_preconditioner%ccs_matrix%row_ind,&
      my_neighbours=taucs_preconditioner%ccs_matrix%col_ptr,incident_arcs=taucs_preconditioner%ccs_matrix%incident,&
      n_edges=n_edges)
      ELSE
      CALL CreateCCSLaplacian(indexing_offset=0,&
      arc_offset=n_special_arcs,node_offset=n_special_nodes,&
      heads_tails=graph%heads_tails,neighbours=taucs_preconditioner%ccs_matrix%row_ind,&
      my_neighbours=taucs_preconditioner%ccs_matrix%col_ptr,incident_arcs=taucs_preconditioner%ccs_matrix%incident,&
      nodes_reordering=taucs_preconditioner%graph_ordering%nodes_reordering,nodes_renumbering=taucs_preconditioner%graph_ordering%&
     &nodes_renumbering,&
      n_edges=n_edges)
      END IF
      taucs_preconditioner%ccs_matrix%n_vertices=n_vertices
      taucs_preconditioner%ccs_matrix%n_edges=n_edges
      
      IF(taucs_preconditioner%log_file/="")THEN
      CALL taucs_logfile(TRIM(taucs_preconditioner%log_file)//CHAR(0))
      CALL taucs_printf("The Laplacian has %d non-zero entries \n"//CHAR(0),%VAL(INT(taucs_preconditioner%ccs_matrix%n_edges)))
      END IF
      
      
      IF(.NOT.ASSOCIATED(taucs_preconditioner%eliminated_nodes))THEN
      taucs_preconditioner%allocated_eliminated=.TRUE.
      ALLOCATE(taucs_preconditioner%eliminated_nodes(1:1))
      taucs_preconditioner%eliminated_nodes=-graph%n_special_nodes
      END IF
      taucs_preconditioner%ccs_matrix%n_eliminated=SIZE(taucs_preconditioner%eliminated_nodes)
      ALLOCATE(taucs_preconditioner%ccs_matrix%eliminated(1:taucs_preconditioner%ccs_matrix%n_eliminated))
      IF(.NOT.taucs_preconditioner%reorder_graph)THEN
      taucs_preconditioner%ccs_matrix%eliminated=taucs_preconditioner%eliminated_nodes+graph%n_special_nodes
      ELSE
      taucs_preconditioner%ccs_matrix%eliminated=taucs_preconditioner%graph_ordering%nodes_renumbering(taucs_preconditioner%elimin&
     &ated_nodes)+graph%n_special_nodes
      END IF
      
      
      CALL taucs_ccs_assemble(taucs_preconditioner%ccs_matrix%matrix,%VAL(n_vertices),&
      taucs_preconditioner%ccs_matrix%col_ptr(0),taucs_preconditioner%ccs_matrix%row_ind(0),taucs_preconditioner%ccs_matrix%val(0)&
     &)
      
      PreFactor:IF((taucs_preconditioner%precond_method==TAUCS_LLt))THEN
      
      
      FactorSymbolic:IF((taucs_preconditioner%factorization/=TAUCS_incomplete_fact).AND.(taucs_preconditioner%factorization/=TAUCS&
     &_complete_fact_ll))THEN
      
      
      ReorderMatrix:IF(taucs_preconditioner%ordering_method=="none")THEN
      
      taucs_preconditioner%ccs_factor=taucs_ccs_factor_llt_symbolic(taucs_preconditioner%ccs_matrix%matrix)
      
      ELSE ReorderMatrix
      
      
      CALL taucs_ccs_order(taucs_preconditioner%ccs_matrix%matrix,taucs_preconditioner%perm,taucs_preconditioner%invperm,TRIM(tauc&
     &s_preconditioner%ordering_method)//CHAR(0));
      
      IF((taucs_preconditioner%perm==0).OR.(taucs_preconditioner%invperm==0))CALL CriticalError(message="Ordering failed in TAUCS"&
     &,&
      caller="TAUCS_InitializePreconditioner")
      
      
      
      taucs_preconditioner%ordered_matrix=taucs_ccs_permute_symmetrically(taucs_preconditioner%ccs_matrix%matrix,%VAL(taucs_precon&
     &ditioner%perm),%VAL(taucs_preconditioner%invperm))
      IF(taucs_preconditioner%ordered_matrix==0)CALL CriticalError(message="Permuting failed in TAUCS",&
      caller="TAUCS_CreatePreconditioner")
      taucs_preconditioner%ccs_factor=taucs_ccs_factor_llt_symbolic(%VAL(taucs_preconditioner%ordered_matrix))
      
      CALL taucs_ccs_free(%VAL(taucs_preconditioner%ordered_matrix))
      taucs_preconditioner%ordered_matrix=0
      
      END IF ReorderMatrix
      
      ELSE FactorSymbolic
      
      IF(taucs_preconditioner%ordering_method/="none")THEN
      
      CALL taucs_ccs_order(taucs_preconditioner%ccs_matrix%matrix,taucs_preconditioner%perm,taucs_preconditioner%invperm,TRIM(tauc&
     &s_preconditioner%ordering_method)//CHAR(0));
      
      IF((taucs_preconditioner%perm==0).OR.(taucs_preconditioner%invperm==0))CALL CriticalError(message="Ordering failed in TAUCS"&
     &,&
      caller="TAUCS_InitializePreconditioner")
      END IF
      
      END IF FactorSymbolic
      
      END IF PreFactor
      
      END SUBROUTINE TAUCS_InitializePreconditioner
      
      
      
      SUBROUTINE TAUCS_DestroyPreconditioner(taucs_preconditioner)
      
      IMPLICIT NONE
      TYPE(TAUCS_Network_Preconditioner),INTENT(INOUT),TARGET::taucs_preconditioner
      
      INTEGER::alloc_status,result
      
      IF((taucs_preconditioner%perm/=0).AND.(taucs_preconditioner%invperm/=0))THEN
      CALL taucs_free_order(%VAL(taucs_preconditioner%perm),%VAL(taucs_preconditioner%invperm))
      taucs_preconditioner%perm=0;
      taucs_preconditioner%invperm=0
      END IF
      
      IF(taucs_preconditioner%ccs_factor/=0)THEN
      IF(taucs_preconditioner%factorization==TAUCS_incomplete_fact)THEN
      CALL taucs_ccs_free(%VAL(taucs_preconditioner%ccs_factor))
      ELSE
      CALL taucs_supernodal_factor_free(%VAL(taucs_preconditioner%ccs_factor))
      END IF
      taucs_preconditioner%ccs_factor=0
      END IF
      
      DEALLOCATE(taucs_preconditioner%ccs_matrix%eliminated,STAT=alloc_status)
      IF(taucs_preconditioner%allocated_eliminated)&
      DEALLOCATE(taucs_preconditioner%eliminated_nodes,STAT=alloc_status)
      
      IF(taucs_preconditioner%allocated_xb_buffers)THEN
      
      IF(ASSOCIATED(taucs_preconditioner%flows))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(taucs_preconditioner%flows),KIND=i_wp),mold=1.0_r_taucs)
      DEALLOCATE(taucs_preconditioner%flows,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(taucs_preconditioner%potentials))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(taucs_preconditioner%potentials),KIND=i_wp),mold=1.0_r_taucs)
      DEALLOCATE(taucs_preconditioner%potentials,STAT=alloc_status)
      END IF
      ELSE
      NULLIFY(taucs_preconditioner%flows)
      NULLIFY(taucs_preconditioner%potentials)
      END IF
      
      
      IF(ASSOCIATED(taucs_preconditioner%ccs_matrix%col_ptr))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(taucs_preconditioner%ccs_matrix%col_ptr),KIND=i_wp),mold=1_i_taucs)
      DEALLOCATE(taucs_preconditioner%ccs_matrix%col_ptr,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(taucs_preconditioner%ccs_matrix%row_ind))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(taucs_preconditioner%ccs_matrix%row_ind),KIND=i_wp),mold=1_i_taucs)
      DEALLOCATE(taucs_preconditioner%ccs_matrix%row_ind,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(taucs_preconditioner%ccs_matrix%incident))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(taucs_preconditioner%ccs_matrix%incident),KIND=i_wp),mold=1_i_taucs)
      DEALLOCATE(taucs_preconditioner%ccs_matrix%incident,STAT=alloc_status)
      END IF
      
      IF(ASSOCIATED(taucs_preconditioner%ccs_matrix%val))THEN
      CALL RecordAllocation(n_elements=-INT(SIZE(taucs_preconditioner%ccs_matrix%val),KIND=i_wp),mold=1.0_r_taucs)
      DEALLOCATE(taucs_preconditioner%ccs_matrix%val,STAT=alloc_status)
      END IF
      
      END SUBROUTINE TAUCS_DestroyPreconditioner
      
      
      
      SUBROUTINE TAUCS_CreatePreconditioner(taucs_preconditioner)
      USE Random_Numbers,ONLY:RandomUniform
      IMPLICIT NONE
      TYPE(TAUCS_Network_Preconditioner),INTENT(INOUT),TARGET::taucs_preconditioner
      
      TYPE(Directed_Graph),POINTER::graph
      INTEGER(KIND=i_wp)::arc,node
      INTEGER(KIND=i_taucs)::edge,vertex,index,seed
      INTEGER::alloc_status,result
      
      graph=>taucs_preconditioner%graph
      
      
      CALL ComputeCCSLaplacian(indexing_offset=0,arc_offset=graph%n_special_arcs,&
      arcs_conductances=taucs_preconditioner%conductances,edges_conductances=taucs_preconditioner%ccs_matrix%val,&
      neighbours=taucs_preconditioner%ccs_matrix%row_ind,my_neighbours=taucs_preconditioner%ccs_matrix%col_ptr,incident_arcs=taucs&
     &_preconditioner%ccs_matrix%incident)
      DO index=1,SIZE(taucs_preconditioner%ccs_matrix%eliminated)
      vertex=taucs_preconditioner%ccs_matrix%eliminated(index)
      
      taucs_preconditioner%ccs_matrix%val(taucs_preconditioner%ccs_matrix%col_ptr(vertex))=taucs_preconditioner%ccs_matrix%val(tau&
     &cs_preconditioner%ccs_matrix%col_ptr(vertex))+taucs_preconditioner%diagonal_conductance
      END DO
      
      Factorize:IF(taucs_preconditioner%precond_method==TAUCS_LLt)THEN
      
      ReorderMatrix:IF(taucs_preconditioner%ordering_method=="none")THEN
      
      SELECT CASE(taucs_preconditioner%factorization)
      CASE(TAUCS_incomplete_fact)
      taucs_preconditioner%ccs_factor=taucs_ccs_factor_llt(taucs_preconditioner%ccs_matrix%matrix,%VAL(taucs_preconditioner%dropto&
     &l),%VAL(0))
      CASE(TAUCS_complete_fact_ll)
      taucs_preconditioner%ccs_factor=taucs_ccs_factor_llt_ll(taucs_preconditioner%ccs_matrix%matrix)
      CASE DEFAULT
      result=taucs_ccs_factor_llt_numeric(taucs_preconditioner%ccs_matrix%matrix,%VAL(taucs_preconditioner%ccs_factor))
      IF(result/=0)THEN
      CALL CriticalError(message="Numerical factorization failed in TAUCS",&
      caller="TAUCS_CreatePreconditioner")
      END IF
      END SELECT
      
      ELSE ReorderMatrix
      
      
      taucs_preconditioner%ordered_matrix=taucs_ccs_permute_symmetrically(taucs_preconditioner%ccs_matrix%matrix,%VAL(taucs_precon&
     &ditioner%perm),%VAL(taucs_preconditioner%invperm))
      IF(taucs_preconditioner%ordered_matrix==0)CALL CriticalError(message="Permuting failed in TAUCS",&
      caller="TAUCS_CreatePreconditioner")
      
      SELECT CASE(taucs_preconditioner%factorization)
      CASE(TAUCS_incomplete_fact)
      taucs_preconditioner%ccs_factor=taucs_ccs_factor_llt(%VAL(taucs_preconditioner%ordered_matrix),%VAL(taucs_preconditioner%dro&
     &ptol),%VAL(0))
      CASE(TAUCS_complete_fact_ll)
      taucs_preconditioner%ccs_factor=taucs_ccs_factor_llt_ll(%VAL(taucs_preconditioner%ordered_matrix))
      CASE DEFAULT
      result=taucs_ccs_factor_llt_numeric(%VAL(taucs_preconditioner%ordered_matrix),%VAL(taucs_preconditioner%ccs_factor))
      IF(result/=0)THEN
      CALL CriticalError(message="Numerical factorization failed in TAUCS",&
      caller="TAUCS_CreatePreconditioner")
      END IF
      END SELECT
      
      CALL taucs_ccs_free(%VAL(taucs_preconditioner%ordered_matrix))
      taucs_preconditioner%ordered_matrix=0
      
      END IF ReorderMatrix
      
      ELSE
      
      
      taucs_preconditioner%subgraphs=MIN(taucs_preconditioner%graph%n_nodes-1_i_wp,MAX(1_i_wp,&
      INT(taucs_preconditioner%graph%n_nodes*taucs_preconditioner%subtree_ratio,i_wp)))
      CALL RandomUniform(seed)
      
      
      taucs_preconditioner%ccs_mwb=taucs_amwb_preconditioner_create(taucs_preconditioner%ccs_matrix%matrix,%VAL(seed),%VAL(taucs_p&
     &reconditioner%subgraphs))
      
      CALL taucs_ccs_order(%VAL(taucs_preconditioner%ccs_mwb),taucs_preconditioner%perm,taucs_preconditioner%invperm,TRIM(taucs_pr&
     &econditioner%ordering_method)//CHAR(0));
      
      IF((taucs_preconditioner%perm==0).OR.(taucs_preconditioner%invperm==0))CALL CriticalError(message="Ordering failed in TAUCS"&
     &,&
      caller="TAUCS_CreatePreconditioner")
      taucs_preconditioner%ordered_matrix=taucs_ccs_permute_symmetrically(%VAL(taucs_preconditioner%ccs_mwb),&
      %VAL(taucs_preconditioner%perm),%VAL(taucs_preconditioner%invperm))
      IF(taucs_preconditioner%ordered_matrix==0)CALL CriticalError(message="Permuting failed in TAUCS",&
      caller="TAUCS_CreatePreconditioner")
      CALL taucs_ccs_free(%VAL(taucs_preconditioner%ccs_mwb))
      taucs_preconditioner%ccs_mwb=0
      
      SELECT CASE(taucs_preconditioner%factorization)
      CASE(TAUCS_incomplete_fact)
      taucs_preconditioner%ccs_factor=taucs_ccs_factor_llt(%VAL(taucs_preconditioner%ordered_matrix),%VAL(taucs_preconditioner%dro&
     &ptol),%VAL(0))
      CASE(TAUCS_complete_fact_ll)
      taucs_preconditioner%ccs_factor=taucs_ccs_factor_llt_ll(%VAL(taucs_preconditioner%ordered_matrix))
      CASE DEFAULT
      taucs_preconditioner%ccs_factor=taucs_ccs_factor_llt_mf(%VAL(taucs_preconditioner%ordered_matrix))
      END SELECT
      
      CALL taucs_ccs_free(%VAL(taucs_preconditioner%ordered_matrix))
      taucs_preconditioner%ordered_matrix=0
      
      END IF Factorize
      
      IF(taucs_preconditioner%ccs_factor==0)THEN
      CALL CriticalError(message="Factorization failed in TAUCS",&
      caller="TAUCS_CreatePreconditioner")
      END IF
      
      END SUBROUTINE TAUCS_CreatePreconditioner
      
      
      
      SUBROUTINE TAUCS_FreePreconditioner(taucs_preconditioner)
      
      IMPLICIT NONE
      TYPE(TAUCS_Network_Preconditioner),INTENT(INOUT),TARGET::taucs_preconditioner
      
      INTEGER::alloc_status,result
      
      IF(taucs_preconditioner%precond_method/=TAUCS_LLt)THEN
      IF((taucs_preconditioner%perm/=0).AND.(taucs_preconditioner%invperm/=0))THEN
      CALL taucs_free_order(%VAL(taucs_preconditioner%perm),%VAL(taucs_preconditioner%invperm))
      taucs_preconditioner%perm=0;
      taucs_preconditioner%invperm=0
      END IF
      END IF
      
      IF(taucs_preconditioner%ccs_factor/=0)THEN
      IF(taucs_preconditioner%factorization==TAUCS_incomplete_fact)THEN
      CALL taucs_ccs_free(%VAL(taucs_preconditioner%ccs_factor))
      taucs_preconditioner%ccs_factor=0
      ELSE IF(taucs_preconditioner%factorization==TAUCS_complete_fact_ll)THEN
      CALL taucs_supernodal_factor_free(%VAL(taucs_preconditioner%ccs_factor))
      taucs_preconditioner%ccs_factor=0
      ELSE
      CALL taucs_supernodal_factor_free_numeric(%VAL(taucs_preconditioner%ccs_factor))
      END IF
      END IF
      
      IF(taucs_preconditioner%ccs_mwb/=0)THEN
      CALL taucs_ccs_free(%VAL(taucs_preconditioner%ccs_mwb))
      taucs_preconditioner%ccs_mwb=0
      END IF
      
      IF(taucs_preconditioner%ordered_matrix/=0)THEN
      CALL taucs_ccs_free(%VAL(taucs_preconditioner%ordered_matrix))
      taucs_preconditioner%ordered_matrix=0
      END IF
      
      END SUBROUTINE TAUCS_FreePreconditioner
      
      
      
      SUBROUTINE TAUCS_ApplyPreconditioner(taucs_preconditioner)
      IMPLICIT NONE
      TYPE(TAUCS_Network_Preconditioner),INTENT(INOUT),TARGET::taucs_preconditioner
      INTEGER(KIND=i_wp)::arc,node
      INTEGER(KIND=i_taucs)::edge,vertex,offset
      INTEGER::alloc_status,result
      
      offset=-taucs_preconditioner%graph%n_special_nodes
      
      IF(taucs_preconditioner%reorder_graph)THEN
      CALL VectorPermute(indexing_offset=offset,&
      source=taucs_preconditioner%excess_flows,target=taucs_preconditioner%flows,permutation=taucs_preconditioner%graph_ordering%n&
     &odes_reordering)
      ELSE IF(taucs_preconditioner%ordering_method/="none")THEN
      CALL taucs_permute_vector(perm=%VAL(taucs_preconditioner%perm),x=taucs_preconditioner%excess_flows(offset),&
      perm_x=taucs_preconditioner%flows(offset),n=%VAL(taucs_preconditioner%ccs_matrix%n_vertices))
      ELSE
      IF(taucs_preconditioner%allocated_xb_buffers)CALL VectorPointerCopy(source=taucs_preconditioner%excess_flows,target=taucs_pr&
     &econditioner%flows)
      END IF
      IF(taucs_preconditioner%factorization==TAUCS_incomplete_fact)THEN
      result=taucs_ccs_solve_llt(%VAL(taucs_preconditioner%ccs_factor),&
      taucs_preconditioner%potentials(offset),taucs_preconditioner%flows(offset))
      ELSE
      result=taucs_supernodal_solve_llt(%VAL(taucs_preconditioner%ccs_factor),&
      taucs_preconditioner%potentials(offset),taucs_preconditioner%flows(offset))
      END IF
      IF(result/=0)THEN
      CALL CriticalError(message="Triangular solve failed in TAUCS",&
      caller="TAUCS_ApplyPreconditioner")
      END IF
      
      IF(taucs_preconditioner%reorder_graph)THEN
      CALL VectorPermute(indexing_offset=offset,&
      source=taucs_preconditioner%potentials,target=taucs_preconditioner%excess_potentials,permutation=taucs_preconditioner%graph_&
     &ordering%nodes_renumbering)
      ELSE IF(taucs_preconditioner%ordering_method/="none")THEN
      CALL taucs_permute_vector(perm=%VAL(taucs_preconditioner%invperm),perm_x=taucs_preconditioner%excess_potentials(offset),&
      x=taucs_preconditioner%potentials(offset),n=%VAL(taucs_preconditioner%ccs_matrix%n_vertices))
      ELSE
      IF(taucs_preconditioner%allocated_xb_buffers)CALL VectorPointerCopy(source=taucs_preconditioner%potentials,target=taucs_prec&
     &onditioner%excess_potentials)
      END IF
      END SUBROUTINE TAUCS_ApplyPreconditioner
      
      
      
      END MODULE TAUCS_Interface