Maximum weighted pairing algorithm for complete graph

Question

The mathematical problem

Let there be 2n persons, and C(i,j) the "cost" of having i and j work together (the function C is quick to compute, in my case it is a given matrix, and is symmetric). The question is to find the arrangement of 2n pairs of persons that minimizes the sum of the costs of each pair.

This should be done in polynomial complexity in n, and implemented relatively easily in the Scilab language (input : cost matrix, output : pairings, for instance a n-by-2 matrix of indexes). I am aware that "relatively easily" is subject to interpretation...

Previous research

This problem is actually solved by the Blossom algorithm. See for instance this paper.

However, this (and its variants) looks like a nightmare to implement. My real problem is for n=20, so although brute force (= trying all possible pairings) is not OK (brute-forcing n=8 took an hour on my computer), pretty much anything better than brute force should do the trick; if I can avoid one week of coding at the cost of one hour of computation I'm in.

I was thinking along the lines of using the Hungarian/Munkres algorithm on a 2n-by-2n array filling the diagonal with +%inf and other elements by the symmetric cost matrix, then somehow selecting from the resulting permutation a relevant pairing, but I fail to find a reliable way to do this. (Note, the Hungarian algorithm is already coded for a separate section, so you may use it without cost to the "easy to implement" requirement.)

I hope that compared to the blossom-algorithm problem, the completeness of the graph allows for some shortcuts... (Edit: see DE's comment below, this is wrong for semi-obvious reasons)

Answer 1

I do not know Scilab I am afraid, but if you are willing to use Python it is very easy as the Networkx library provides support for this function:

import networkx as nx
import networkx.algorithms.matching as matching

def C(i,j):
    return i*j

n=40
G=nx.Graph()
for i in range(n):
    for j in range(n):
        G.add_edge(i,j,weight = -C(i,j))
M = matching.max_weight_matching(G,maxcardinality=True)
for i in M:
    print i,'with',M[i]

This code prints out the answer within a second.

The function C defines the cost of pairing i with j. Note that the weights are set to -C(i,j) in order to transform the max_weight_matching into a min_weight_matching algorithm.