flow/general-weighted-matching.hpp

• View this file on GitHub
• Last update: 2026-07-22 11:58:01+07:00
• Include: #include "flow/general-weighted-matching.hpp"

Depends on

• template.hpp

Required by

• verify/yosupo/graph/assignment.cpp

Verified with

• verify/yosupo/graph/general_weighted_matching.test.cpp

Code

#pragma once
#include "template.hpp"

/**
 * Description: General weighted matching using blossom algorithm.
 * Time: O(|V|^3) \text{(fast in practice)}
 */

template<int N,class T>
struct GeneralWeightedMatching{
static const T INF=numeric_limits<T>::max();
struct Edge{
  int u,v;
  T w;
  Edge(){}
  Edge(int u,int v,T w):u(u),v(v),w(w){}
};
int n,n_x;
Edge g[N*2][N*2];
T lab[N*2];
int match[N*2],slack[N*2],st[N*2],pa[N*2];
int flo_from[N*2][N+1],S[N*2],vis[N*2];
vector<int> flo[N*2];
queue<int> q;
T e_delta(const Edge &e){return lab[e.u]+lab[e.v]-g[e.u][e.v].w*2;}
void update_slack(int u,int x){
  if(!slack[x]||e_delta(g[u][x])<e_delta(g[slack[x]][x]))slack[x]=u;
}
void set_slack(int x){
  slack[x]=0;
  for(int u=1;u<=n;u++)
    if(g[u][x].w>0&&st[u]!=x&&S[st[u]]==0)update_slack(u,x);
}
void q_push(int x){
  if(x<=n)q.push(x);
  else for(auto xs:flo[x])q_push(xs);
}
void set_st(int x,int b){
  st[x]=b;
  if(x>n)for(auto xs:flo[x])set_st(xs,b);
}
int get_pr(int b,int xr){
  int pr=find(ALL(flo[b]),xr)-flo[b].begin();
  if(pr&1){reverse(1+ALL(flo[b]));return SZ(flo[b])-pr;}
  return pr;
}
void set_match(int u,int v){
  Edge e=g[u][v];match[u]=g[u][v].v;if(u<=n)return;
  int xr=flo_from[u][e.u],pr=get_pr(u,xr);
  for(int i=0;i<pr;i++)set_match(flo[u][i],flo[u][i^1]);
  set_match(xr,v);rotate(flo[u].begin(),pr+ALL(flo[u]));
}
void augment(int u,int v){
  while(true){
    int xnv=st[match[u]];set_match(u,v);
    if(!xnv)return;
    set_match(xnv,st[pa[xnv]]);u=st[pa[xnv]];v=xnv;
  }
}
int get_lca(int u,int v){
  static int t=0;
  for(++t;u||v;swap(u,v)){
    if(u==0)continue;
    if(vis[u]==t)return u;
    vis[u]=t;u=st[match[u]];
    if(u)u=st[pa[u]];
  }
  return 0;
}
void add_blossom(int u,int lca,int v){
  int b=n+1;while(b<=n_x&&st[b])++b;
  if(b>n_x)++n_x; //new blossom;
  lab[b]=0;S[b]=0;match[b]=match[lca];flo[b]=vector<int>{lca};
  for(int x=u,y;x!=lca;x=st[pa[y]])
    flo[b].eb(x),flo[b].eb(y=st[match[x]]),q_push(y);
  reverse(1+ALL(flo[b]));
  for(int x=v,y;x!=lca;x=st[pa[y]])
    flo[b].eb(x),flo[b].eb(y=st[match[x]]),q_push(y);
  set_st(b,b);
  for(int x=1;x<=n_x;x++)g[b][x].w=g[x][b].w=0;
  for(int x=1;x<=n;x++)flo_from[b][x]=0;
  for(auto xs:flo[b]){
    for(int x=1;x<=n_x;x++)
      if(g[b][x].w==0||e_delta(g[xs][x])<e_delta(g[b][x]))
        g[b][x]=g[xs][x],g[x][b]=g[x][xs];
    for(int x=1;x<=n;x++)if(flo_from[xs][x])flo_from[b][x]=xs;
  }
  set_slack(b);
}
void expand_blossom(int b){
  for(auto xs:flo[b])set_st(xs,xs);
  int xr=flo_from[b][g[b][pa[b]].u],pr=get_pr(b,xr);
  for(int i=0;i<pr;i+=2){
    int xs=flo[b][i],xns=flo[b][i+1];
    pa[xs]=g[xns][xs].u;S[xs]=1;S[xns]=0;
    slack[xs]=0,set_slack(xns);q_push(xns);
  }
  S[xr]=1,pa[xr]=pa[b];
  for(int i=pr+1;i<SZ(flo[b]);i++)S[flo[b][i]]=-1,set_slack(flo[b][i]);
  st[b]=0;
}
bool on_found_edge(const Edge &e){
  int u=st[e.u],v=st[e.v];
  if(S[v]==-1){
    pa[v]=e.u;S[v]=1;
    int nu=st[match[v]];
    slack[v]=slack[nu]=S[nu]=0;q_push(nu);
  }else if(S[v]==0){
    int lca=get_lca(u,v);
    if(!lca)return augment(u,v),augment(v,u),true;
    add_blossom(u,lca,v);
  }
  return false;
}
bool matching(){
  for(int i=1;i<=n_x;i++)S[i]=-1;
  for(int i=1;i<=n_x;i++)slack[i]=-0;
  q=queue<int>();
  for(int x=1;x<=n_x;x++)if(st[x]==x&&!match[x])pa[x]=S[x]=0,q_push(x);
  if(q.empty())return false;
  while(true){
    while(!q.empty()){
      int u=q.front();q.pop();if(S[st[u]]==1)continue;
      for(int v=1;v<=n;v++)if(g[u][v].w>0&&st[u]!=st[v]){
        if(e_delta(g[u][v])==0){if(on_found_edge(g[u][v]))return true;}
        else update_slack(u,st[v]);
      }
    }
    T d=INF;
    for(int b=n+1;b<=n_x;b++)if(st[b]==b&&S[b]==1)d=min(d,lab[b]/2);
    for(int x=1;x<=n_x;x++)if(st[x]==x&&slack[x]){
      if(S[x]==-1)d=min(d,e_delta(g[slack[x]][x]));
      else if(S[x]==0)d=min(d,e_delta(g[slack[x]][x])/2);
    }
    for(int u=1;u<=n;u++){
      if(S[st[u]]==0){if(lab[u]<=d)return 0;lab[u]-=d;}
      else if(S[st[u]]==1)lab[u]+=d;
    }
    for(int b=n+1;b<=n_x;b++)if(st[b]==b){
      if(S[st[b]]==0)lab[b]+=d*2;
      else if(S[st[b]]==1)lab[b]-=d*2;
    }
    q=queue<int>();
    for(int x=1;x<=n_x;x++)
      if(st[x]==x&&slack[x]&&st[slack[x]]!=x&&e_delta(g[slack[x]][x])==0)
        if(on_found_edge(g[slack[x]][x]))return true;
    for(int b=n+1;b<=n_x;b++)
      if(st[b]==b&&S[b]==1&&lab[b]==0)expand_blossom(b);
  }
  return false;
}
pair<T,int> solve(){
  for(int i=1;i<=n;i++)match[i]=0;
  n_x=n;int n_matches=0;T w_max=0,tot_weight=0;
  for(int u=0;u<=n;u++)st[u]=u,flo[u].clear();
  for(int u=1;u<=n;u++)for(int v=1;v<=n;v++)
    flo_from[u][v]=u==v?u:0,w_max=max(w_max,g[u][v].w);
  for(int u=1;u<=n;u++)lab[u]=w_max;
  while(matching())++n_matches;
  for(int u=1;u<=n;u++)if(match[u]&&match[u]<u)tot_weight+=g[u][match[u]].w;
  return make_pair(tot_weight,n_matches);
}
void add_edge(int u,int v,T w){
  if(w<=g[u][v].w)return;
  g[u][v].w=g[v][u].w=w;
}
void init(int _n){
  n=_n;
  for(int u=1;u<=n;u++)for(int v=1;v<=n;v++)g[u][v]=Edge(u,v,T(0));
}

};

#line 2 "template.hpp"
#include<bits/stdc++.h>

using namespace std;

#define pb push_back
#define eb emplace_back
#define mp make_pair
#define mt make_tuple
#define fi first
#define se second

#define ALL(a) a.begin(),a.end()
#define RALL(a) a.rbegin(),a.rend()
#define SORT(a) sort(ALL(a))
#define RSORT(a) sort(RALL(a))
#define REV(a) reverse(ALL(a))
#define UNI(a) a.erase(unique(ALL(a)),a.end())
#define SZ(a) (int)(a.size())
#define LB(a,x) (int)(lower_bound(ALL(a),x)-a.begin())
#define UB(a,x) (int)(upper_bound(ALL(a),x)-a.begin())
#define MIN(a) *min_element(ALL(a))
#define MAX(a) *max_element(ALL(a))

using ll = long long;
using db = long double;
using i128 = __int128_t;
using u32 = uint32_t;
using u64 = uint64_t;

const int INF=INT_MAX/2;
const ll LINF=LLONG_MAX/4;
const db DINF=numeric_limits<db>::infinity();
const int MOD=998244353;
const int MOD2=1000000007;
const db EPS=1e-9;
const db PI=acos(db(-1));

template<class T>
using PQ = priority_queue<T,vector<T>,greater<T>>;

#define vv(T,a,n,...) vector<vector<T>> a(n,vector<T>(__VA_ARGS__))
#define vvv(T,a,n,m,...) vector<vector<vector<T>>> a(n,vector<vector<T>>(m,vector<T>(__VA_ARGS__)))
#define vvvv(T,a,n,m,k,...) vector<vector<vector<vector<T>>>> a(n,vector<vector<vector<T>>>(m,vector<vector<T>>(k,vector<T>(__VA_ARGS__))))

template<class T,class U>
bool chmin(T &a,U b){return b<a?a=b,1:0;}
template<class T,class U>
bool chmax(T &a,U b){return a<b?a=b,1:0;}
template<class T,class U>
T SUM(const U &a){return accumulate(ALL(a),T{});}

mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
mt19937_64 rng64(chrono::steady_clock::now().time_since_epoch().count());
#line 3 "flow/general-weighted-matching.hpp"

/**
 * Description: General weighted matching using blossom algorithm.
 * Time: O(|V|^3) \text{(fast in practice)}
 */

template<int N,class T>
struct GeneralWeightedMatching{
static const T INF=numeric_limits<T>::max();
struct Edge{
  int u,v;
  T w;
  Edge(){}
  Edge(int u,int v,T w):u(u),v(v),w(w){}
};
int n,n_x;
Edge g[N*2][N*2];
T lab[N*2];
int match[N*2],slack[N*2],st[N*2],pa[N*2];
int flo_from[N*2][N+1],S[N*2],vis[N*2];
vector<int> flo[N*2];
queue<int> q;
T e_delta(const Edge &e){return lab[e.u]+lab[e.v]-g[e.u][e.v].w*2;}
void update_slack(int u,int x){
  if(!slack[x]||e_delta(g[u][x])<e_delta(g[slack[x]][x]))slack[x]=u;
}
void set_slack(int x){
  slack[x]=0;
  for(int u=1;u<=n;u++)
    if(g[u][x].w>0&&st[u]!=x&&S[st[u]]==0)update_slack(u,x);
}
void q_push(int x){
  if(x<=n)q.push(x);
  else for(auto xs:flo[x])q_push(xs);
}
void set_st(int x,int b){
  st[x]=b;
  if(x>n)for(auto xs:flo[x])set_st(xs,b);
}
int get_pr(int b,int xr){
  int pr=find(ALL(flo[b]),xr)-flo[b].begin();
  if(pr&1){reverse(1+ALL(flo[b]));return SZ(flo[b])-pr;}
  return pr;
}
void set_match(int u,int v){
  Edge e=g[u][v];match[u]=g[u][v].v;if(u<=n)return;
  int xr=flo_from[u][e.u],pr=get_pr(u,xr);
  for(int i=0;i<pr;i++)set_match(flo[u][i],flo[u][i^1]);
  set_match(xr,v);rotate(flo[u].begin(),pr+ALL(flo[u]));
}
void augment(int u,int v){
  while(true){
    int xnv=st[match[u]];set_match(u,v);
    if(!xnv)return;
    set_match(xnv,st[pa[xnv]]);u=st[pa[xnv]];v=xnv;
  }
}
int get_lca(int u,int v){
  static int t=0;
  for(++t;u||v;swap(u,v)){
    if(u==0)continue;
    if(vis[u]==t)return u;
    vis[u]=t;u=st[match[u]];
    if(u)u=st[pa[u]];
  }
  return 0;
}
void add_blossom(int u,int lca,int v){
  int b=n+1;while(b<=n_x&&st[b])++b;
  if(b>n_x)++n_x; //new blossom;
  lab[b]=0;S[b]=0;match[b]=match[lca];flo[b]=vector<int>{lca};
  for(int x=u,y;x!=lca;x=st[pa[y]])
    flo[b].eb(x),flo[b].eb(y=st[match[x]]),q_push(y);
  reverse(1+ALL(flo[b]));
  for(int x=v,y;x!=lca;x=st[pa[y]])
    flo[b].eb(x),flo[b].eb(y=st[match[x]]),q_push(y);
  set_st(b,b);
  for(int x=1;x<=n_x;x++)g[b][x].w=g[x][b].w=0;
  for(int x=1;x<=n;x++)flo_from[b][x]=0;
  for(auto xs:flo[b]){
    for(int x=1;x<=n_x;x++)
      if(g[b][x].w==0||e_delta(g[xs][x])<e_delta(g[b][x]))
        g[b][x]=g[xs][x],g[x][b]=g[x][xs];
    for(int x=1;x<=n;x++)if(flo_from[xs][x])flo_from[b][x]=xs;
  }
  set_slack(b);
}
void expand_blossom(int b){
  for(auto xs:flo[b])set_st(xs,xs);
  int xr=flo_from[b][g[b][pa[b]].u],pr=get_pr(b,xr);
  for(int i=0;i<pr;i+=2){
    int xs=flo[b][i],xns=flo[b][i+1];
    pa[xs]=g[xns][xs].u;S[xs]=1;S[xns]=0;
    slack[xs]=0,set_slack(xns);q_push(xns);
  }
  S[xr]=1,pa[xr]=pa[b];
  for(int i=pr+1;i<SZ(flo[b]);i++)S[flo[b][i]]=-1,set_slack(flo[b][i]);
  st[b]=0;
}
bool on_found_edge(const Edge &e){
  int u=st[e.u],v=st[e.v];
  if(S[v]==-1){
    pa[v]=e.u;S[v]=1;
    int nu=st[match[v]];
    slack[v]=slack[nu]=S[nu]=0;q_push(nu);
  }else if(S[v]==0){
    int lca=get_lca(u,v);
    if(!lca)return augment(u,v),augment(v,u),true;
    add_blossom(u,lca,v);
  }
  return false;
}
bool matching(){
  for(int i=1;i<=n_x;i++)S[i]=-1;
  for(int i=1;i<=n_x;i++)slack[i]=-0;
  q=queue<int>();
  for(int x=1;x<=n_x;x++)if(st[x]==x&&!match[x])pa[x]=S[x]=0,q_push(x);
  if(q.empty())return false;
  while(true){
    while(!q.empty()){
      int u=q.front();q.pop();if(S[st[u]]==1)continue;
      for(int v=1;v<=n;v++)if(g[u][v].w>0&&st[u]!=st[v]){
        if(e_delta(g[u][v])==0){if(on_found_edge(g[u][v]))return true;}
        else update_slack(u,st[v]);
      }
    }
    T d=INF;
    for(int b=n+1;b<=n_x;b++)if(st[b]==b&&S[b]==1)d=min(d,lab[b]/2);
    for(int x=1;x<=n_x;x++)if(st[x]==x&&slack[x]){
      if(S[x]==-1)d=min(d,e_delta(g[slack[x]][x]));
      else if(S[x]==0)d=min(d,e_delta(g[slack[x]][x])/2);
    }
    for(int u=1;u<=n;u++){
      if(S[st[u]]==0){if(lab[u]<=d)return 0;lab[u]-=d;}
      else if(S[st[u]]==1)lab[u]+=d;
    }
    for(int b=n+1;b<=n_x;b++)if(st[b]==b){
      if(S[st[b]]==0)lab[b]+=d*2;
      else if(S[st[b]]==1)lab[b]-=d*2;
    }
    q=queue<int>();
    for(int x=1;x<=n_x;x++)
      if(st[x]==x&&slack[x]&&st[slack[x]]!=x&&e_delta(g[slack[x]][x])==0)
        if(on_found_edge(g[slack[x]][x]))return true;
    for(int b=n+1;b<=n_x;b++)
      if(st[b]==b&&S[b]==1&&lab[b]==0)expand_blossom(b);
  }
  return false;
}
pair<T,int> solve(){
  for(int i=1;i<=n;i++)match[i]=0;
  n_x=n;int n_matches=0;T w_max=0,tot_weight=0;
  for(int u=0;u<=n;u++)st[u]=u,flo[u].clear();
  for(int u=1;u<=n;u++)for(int v=1;v<=n;v++)
    flo_from[u][v]=u==v?u:0,w_max=max(w_max,g[u][v].w);
  for(int u=1;u<=n;u++)lab[u]=w_max;
  while(matching())++n_matches;
  for(int u=1;u<=n;u++)if(match[u]&&match[u]<u)tot_weight+=g[u][match[u]].w;
  return make_pair(tot_weight,n_matches);
}
void add_edge(int u,int v,T w){
  if(w<=g[u][v].w)return;
  g[u][v].w=g[v][u].w=w;
}
void init(int _n){
  n=_n;
  for(int u=1;u<=n;u++)for(int v=1;v<=n;v++)g[u][v]=Edge(u,v,T(0));
}

};

Back to top page