MMVSS optimisation

#include <iostream>
#include <iomanip>
#include <string>
#include <cmath>

using namespace std;

struct matrix_elem
{
    double value;
    int i;
    int j;
};

template<typename T>
T ** alloc_matrix(int nrows, int ncols)
{
    T ** temp_matrix = new T*[nrows];
    for (int i = 0; i < nrows; ++i)
    {
        temp_matrix[i] = new T[ncols];
        for (int j = 0; j < ncols; ++j)
        {
            temp_matrix[i][j] = 0;
        }
    }

    return temp_matrix;
}

template<typename T>
void delete_matrix(T ** matrix, int nrows)
{
    for (int i = 0; i < nrows; ++i)
    {
        delete[] matrix[i];
    }
    delete[] matrix;
}

template<typename T>
T read_value()
{
    T result;
    if (cin >> result)
        return result;

    cout << "wrong value" << endl;
    return read_value<T>();
}

template<typename T>
void read_matrix(T ** matrix, int nrows, int ncols)
{
    if (ncols != 1)
    {
        for (int i = 0; i < nrows; ++i)
        {
            for (int j = 0; j < ncols; ++j)
            {
                matrix[i][j] = read_value<T>();
            }
        }
    }
    else
    {
        for (int i = 0; i < nrows; i++)
        {
            *matrix[i] = read_value<T>();
        }
    }
}

template<typename T>
void print_matrix(T ** matrix, int nrows, int ncols)
{
    for (int i = 0; i < nrows; ++i)
    {
        for (int j = 0; j < ncols; ++j)
        {
            cout << setw(9) << matrix[i][j] << " ";
        }
        cout << endl;
    }
}

void calculate_DEL(double ** DEL,
                    double ** A,
                    double ** b0,
                    int nrows,
                    int ncols,
                    int L,
                    int dc,
                    int Optimal_i,
                    int Optimal_j)
{
    double Bm = 0.0;
    double Am = 0.0;

    for (int i = 0; i < nrows; ++i)
    {
        for (int j = 0; j < ncols; ++j)
        {
            Bm = b0[i][j];
            Am = A[i][j];

            if (Optimal_i == i && Optimal_j == j)
                Bm += dc;

            DEL[i][j] = L / (Bm - Am);
        }
    }
}

void calculate_dl(double ** dl,
                    double ** DEL,
                    int ** R,
                    int nrows,
                    int ncols)
{
    double sum = 0.0;

    for (int i = 1; i <= nrows; ++i)
    {
        for (int j = 1; j <= ncols; ++j)
        {
            sum = 0.0;

            if (i == j)
            {
                sum = DEL[i - 1][j - 1];
            }
            else
            {
                int vert = i;
                int next_vert = j;

                while (vert != j) {
                    next_vert = R[vert - 1][j - 1];

                    sum += DEL[vert - 1][next_vert - 1];

                    vert = next_vert;
                }
            }

            dl[i - 1][j - 1] = sum;
        }
    }
}

double calculate_sum(double ** dl,
                    int nrows,
                    int ncols,
                    double Topt)
{
    double sum = 0.0;

    for (int i = 0; i < nrows; ++i)
    {
        for (int j = 0; j < ncols; ++j)
        {
            sum += (dl[i][j] - Topt) * (dl[i][j] - Topt);
        }
    }

    return sum;
}

matrix_elem find_min_O(double ** O,
                                int nrows,
                                int ncols,
                                double pred_min_value)
{
    matrix_elem new_min = {pred_min_value, 0, 0};

    for (int i = 0; i < nrows; ++i)
    {
        for (int j = 0; j < ncols; ++j)
        {
            if (O[i][j] != 0)
            {
                double temp_min = O[i][j];

                if (temp_min < pred_min_value)
                {
                    new_min.value = temp_min;
                    new_min.i = i;
                    new_min.j = j;
                }
            }
        }
    }
    return new_min;
}

void calculate_optimal(double ** O,
                        double ** dl,
                        double ** b0,
                        double ** A,
                        double ** B,
                        int ** R,
                        double ** DEL,
                        int nrows,
                        int ncols,
                        double Topt,
                        int dc,
                        int L)
{
    for (int i = 0; i < nrows; ++i)
    {
        for (int j = 0; j < ncols; ++j)
        {
            if (b0[i][j] != 0)
            {
                calculate_DEL(DEL, A, b0, nrows, ncols, L, dc, i, j);
            }
            calculate_dl(dl, DEL, R, nrows, ncols);

            O[i][j] = calculate_sum(dl, nrows, ncols, Topt);
        }
    }
}

void start_optimisation(double ** O,
                        double ** dl,
                        double ** b0,
                        double ** A,
                        double ** B,
                        double ** DEL,
                        int ** R,
                        int nrows,
                        int ncols,
                        double Topt,
                        int L,
                        int dc)
{
    double O_pred = 10E16;
    matrix_elem min = { 10E16, 0, 0 };

    int iteration = 0;

    do
    {
        cout << "iteration: " << ++iteration << endl;
        O_pred = min.value;

        calculate_optimal(O, dl, b0,
                            A, B, R,
                            DEL, nrows, ncols,
                            Topt, dc, L);

        min = { O[0][0], 0, 0 };

        min = find_min_O(O, nrows, ncols, min.value);

        b0[min.i][min.j] += dc;

    } while (min.value < O_pred);
}

int main()
{
    double T0 = 0.1;    //seconds
    int L = 200 * 8;        //byte
    int a0 = 85600;     //kbit/s
    double q = 98.0;    //probability

    int nverts = 20;    //number of graph verts
    double y0 = 0.1;    //Erlang

    double Topt = T0 / 2;   //optimal delay
    int dc = 10000;         //step

    int ** R = alloc_matrix<int>(nverts, nverts);
    double ** A = alloc_matrix<double>(nverts, nverts);
    double ** B = alloc_matrix<double>(nverts, nverts);
    double ** b0 = alloc_matrix<double>(nverts, nverts);
    double ** DEL = alloc_matrix<double>(nverts, nverts);
    double ** dl = alloc_matrix<double>(nverts, nverts);
    double ** O = alloc_matrix<double>(nverts, nverts);

    cout << "R" << endl;
    read_matrix(R, nverts, nverts);
    cout << "A" << endl;
    read_matrix(A, nverts, nverts);
    cout << "B" << endl;
    read_matrix(B, nverts, nverts);

    for (int i = 0; i < nverts; ++i)
    {
        for (int j = 0; j < nverts; ++j)
        {
            b0[i][j] = B[i][j];
        }
    }

    cout << "START OPTIMISATION" << endl;

    start_optimisation(O, dl, b0,
                        A, B, DEL,
                        R, nverts, nverts,
                        Topt, L, dc);

    cout << "B0" << endl;
    print_matrix<double>(b0, nverts, nverts);
    cout << "DEL" << endl;
    print_matrix<double>(DEL, nverts, nverts);
    cout << "DL" << endl;
    print_matrix<double>(dl, nverts, nverts);
    cout << "O" << endl;
    print_matrix<double>(O, nverts, nverts);

    cout << "press any key to exit..." << endl;
    string exit;
    cin >> exit;

    delete_matrix(R, nverts);
    delete_matrix(A, nverts);
    delete_matrix(B, nverts);
    delete_matrix(b0, nverts);
    delete_matrix(DEL, nverts);
    delete_matrix(dl, nverts);
    delete_matrix(O, nverts);
    return 0;
}