Добрый день сообщество CF.

Я начал изучать потоки и столкнулся с проблемами. Недавно я нашёл алгоритм максимального потока минимальной стоимости.

Алгоритм

1) Для каждого ребра, соединяющего вершины u и v добавим обратное, соединяющее вершины v и u, с пропускной способностью равной нулю и стоимостью, равной минус стоимости прямого ребра.

2) Введем для каждой вершины u потенциал φu и расстояние от источника du. Пусть сначала они будут равны нулю.

3) Для каждого ребра e (в том числе и для обратного), соединяющего вершины u и v, введем вес weight(e) = cost(e) + φu — φv (именно для определения веса и нужны потенциалы — без них вес мог бы быть отрицательным для обратных ребер и нельзя было бы применить алгоритм Дейкстры), где cost(e) — стоимость ребра.

4) С помощью алгоритма Дейкстры находим расстояния от источника до всех вершин (в поиске участвуют только те ребра, поток по которым не максимальный). Если нет пути от источника до стока, значит максимальный поток найден и можно завершить алгоритм. Если есть — идем дальше.

5) Выставляем потенциалы для каждой вершины равными расстояниям от источника, которые нашли на предыдущем шаге. 6) Восстанавливаем кратчайший путь.

7) Находим максимальный дополняющий поток на этом пути.

8) Пропускаем найденный поток через кратчайший путь. Когда поток пропускается через ребро, то через соответствующее обратное (или прямое) ребро пропускается поток равный по модулю, но противоположный по значению. Таким образом, обратное ребро характеризует, на сколько можно уменьшить поток через прямое ребро, и какая от этого будет польза.

9) Итерация алгоритма закончена, идем на шаг 3.

Вот моя реализация:

#include<iostream>  /// программа на случай неориентированного графа
#include<cstdio>
#include<cmath>
#include<vector>
#include<fstream>
using namespace std;
int p[103];   /// массив предков для восстановления кратчайших путей
int dist[103];  /// массив кратчайших путей из точки истока
int fi[103];   /// массив потециалов
bool flag[103];  
int last;  /// номер последнего элемента в куче
int n,m;   /// кол-во вершин и рёбер в графе
int start,finish; /// start - точка истока   finish - точка стока
struct kucha
{
    int v,c; /// v - номер вершины,   c - длина пути
};
kucha heap[100003]; /// куча для Дейкстры
struct rebro
{
    int v;     /// вершина в которую ведёт ребро
    int flow;  /// пропускная способность ребра
    int rflow; /// кол-во единиц потока уже проходящих через это ребро
    int cost;  /// стоимость ребра в случае использования, не зависит от кол-ва потока через него
};
vector < vector <rebro> > g;  /// способ хранения графа
void add(int x,int y,int flow,int cost) /// добавляем ребро
{
    rebro z;
    z.v=y;
    z.flow=flow;
    z.cost=cost;
    z.rflow=0;
    g[x].push_back(z);
    z.flow=0;
    z.cost=-cost;
    z.v=x;
    g[y].push_back(z);
}
int maxflow;  /// max поток
int mincost;  /// min цена
int main()
{
    int x,y,fl,cos;
    int i,j;
    scanf("%d %d",&n,&m);
    scanf("%d %d",&start,&finish);
    g.resize(n);
    start--;
    finish--;
    for(i=0;i<m;i++)  /// чтение графа
    {
        scanf("%d %d %d %d",&x,&y,&fl,&cos);  /// fl - пропускная способность ребра  cos - его стоимость
        x--;
        y--;
        add(x,y,fl,cos);
        add(y,x,fl,cos);
    }
    int yy;
    int addflow; /// добавочный поток
    while(true)
    {
        for(i=0;i<n;i++)
            dist[i]=1000000000;
        for(i=0;i<n;i++)
            flag[i]=false;
        for(i=0;i<2*m;i++)
            heap[i].c=1000000000;
        dist[start]=0;
        heap[1].v=start;
        heap[1].c=0;
        last=1;
        while(last>0) /// Дейкстра с потенциалами
        {
            x=heap[1].v;
            if(!flag[x])
            {
                flag[x]=true;
                for(i=0;i<g[x].size();i++)
                    if(g[x][i].rflow<g[x][i].flow && !flag[g[x][i].v])
                    {
                        if(dist[g[x][i].v]>dist[x]+g[x][i].cost+fi[x]-fi[g[x][i].v])
                        {
                            dist[g[x][i].v]=dist[x]+g[x][i].cost+fi[x]-fi[g[x][i].v];
                            p[g[x][i].v]=x;
                            yy=dist[g[x][i].v];
                            last++;
                            j=last;
                            while(j>1 && yy<heap[j>>1].c)
                            {
                                heap[j]=heap[j>>1];
                                j>>=1;
                            }
                            heap[j].c=yy;
                            heap[j].v=g[x][i].v;
                        }
                    }
            }
            j=1;
            while(j*2<last && (heap[last].c>heap[j*2].c || heap[last].c>heap[j*2+1].c))
            {
                if(heap[j*2].c<=heap[j*2+1].c)
                {
                    heap[j]=heap[j*2];
                    j*=2;
                }
                else
                {
                    heap[j]=heap[j*2+1];
                    j*=2;
                    j++;
                }
            }
            heap[j]=heap[last--];
        }
        if(dist[finish]==1000000000) /// если нет пути выходим из алгоритма
            break;
        for(i=0;i<n;i++) /// расстановка потенциалов
            fi[i]=dist[i];
        x=finish;
        addflow=1000000000;
        while(x!=start) /// восстановление кратчайшего пути и поиск добавочного потока
        {
            for(i=0;i<g[p[x]].size();i++)
                if(g[p[x]][i].v==x && g[p[x]][i].flow!=0)
                {
                    addflow=min(addflow,g[p[x]][i].flow-g[p[x]][i].rflow);
                    break;
                }
            x=p[x];
        }
        maxflow+=addflow;
        x=finish;
        while(x!=start) ///  сколько придётся заплатить за увеличение потока + увеличение потока
        {
            for(i=0;i<g[p[x]].size();i++)
                if(g[p[x]][i].v==x && g[p[x]][i].flow!=0)
                {
                    if(g[p[x]][i].rflow==0)
                        mincost+=g[p[x]][i].cost;
                    g[p[x]][i].rflow+=addflow;
                    for(j=0;j<g[x].size();j++)
                        if(g[x][j].v==p[x])
                        {
                            g[x][j].rflow-=addflow;
                            break;
                        }
                    break;
                }
            x=p[x];
        }
    }
    cout<<maxflow<<endl;
    cout<<mincost<<endl;
    return 0;
}  /// самый серьёзный тест на котором я проверял правильность
/*
5 6
1 5
1 3 3 5
1 4 1 1
1 2 4 10
2 4 5 1
3 4 2 2
4 5 3 1
*/

У меня остались вопросы:

По реализации:

1) Возможно удобнее хранить граф не vector < vector > g; , а rebro g[MAXN][MAXN]; ( матрица смежности ) + vector < vector > gg; (список смежности для скорости)?

2) Может можно избавиться от rflow в struct rebro?

С радостью выслушаю способы оптимизации моего кода.

По алгоритму:

1) Время работы O(n*m*log(n))?

2) Как быть если ограничен поток в вершине и задана её стомость?

3) Пригоден ли вообще такой алгоритм для случая, когда стоимоть ребра задана ввиде d денежных единиц на единицу потока и пропускная способность ребра >=1 ? Если да, то как это сделать?

Не могу найти задач на CF, которые могли бы проверить мою реализацию алгоритма на правильность. Пожалуйста, подскажите парочку таких на CF. Зааранее благодарю тех, кто поможет. Извиняюсь за большой текст.

Upd: Спасибо всем тем, кто ответил: CountZero,Miras321.

Особенно благодарен cmd,Perlik.

Благодаря их помощи:

1) заработал алгоритм 2) я понял, что решают задачи на потоки.