Jacobi Algorithm

#include <mpi.h>
#include <iostream>
#include "math.h"
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <ctime>

using namespace std;

/*
int read_data(double** a) {
ifstream f_inp("matrix.txt");
int n = 0;
f_inp >> n;
for (int j = 0; j < n; j++){
for (int k = 0; k < n; k++){
f_inp >> a[j][k];
}
}
f_inp.close();
return n;
}
*/


int main(int argc, char *argv[]) {

    int rank;
    int size;
    int rc;

    if (rc = MPI_Init(&argc, &argv))
    {
        cout << "ошибка запуска, выполнение остановлено " << endl;
        MPI_Abort(MPI_COMM_WORLD, rc);
    }

    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    if (rank == 0)
        cout << "number of processes: " << size << endl;

    if (size < 2){
        MPI_Abort(MPI_COMM_WORLD, MPI_ERR_OTHER);
    }

    for (int count = 0; count < 10; count++){
        int n = 30 + 30 * count;
        double **a;
        double *d;

        d = new double[n];
        a = new double*[n];

        clock_t t1, t2;

        if (rank == 0){

            for (int i = 0; i < n; i++)
                a[i] = new double[n];

            for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                a[i][j] = rand() % 101 - 50;
            /*
            int num = read_data(a);
            if (num == 0){
            cout << "read_data error" << endl;
            system("pause");
            return -1;
            }
            */

            t1 = clock();
        }
//////////////////////////////////////////////////////////
//          Jacobi
//////////////////////////////////////////////////////////

        double *b;
        double *z;
        b = new double[n];
        z = new double[n];

        if (rank == 0)
        for (int i = 0; i < n; ++i){
            z[i] = 0.;
            b[i] = d[i] = a[i][i];
        }

        for (int i = 0; i < 50; ++i)
        {
            double sm = 0.;
            double tresh;
            int p, q;
            bool BR = false;

            if (rank == 0){
                for (p = 0; p < n - 1; ++p)
                {
                    for (q = p + 1; q < n; ++q)
                        sm += fabs(a[p][q]);
                }
                if (sm == 0)
                    BR = true;
            }

            MPI_Bcast(&BR, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);
            if (BR)
                break;

            if (rank == 0)
                tresh = i < 3 ? 0.2 * sm / (n*n) : 0.;

            MPI_Bcast(&tresh, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

            for (p = 0; p < n - 1; ++p)
            {
                for (q = p + 1; q < n; ++q)
                {
                    double apqdqdp[3];
                    if (rank == 0){
                        apqdqdp[0] = a[p][q];
                        apqdqdp[1] = d[q];
                        apqdqdp[2] = d[p];
                    }
                    MPI_Bcast(apqdqdp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);

                    double g = 1e12 * fabs(apqdqdp[0]);
                    if (i >= 3 && fabs(apqdqdp[1]) > g && fabs(apqdqdp[1]) > g){
                        apqdqdp[0] = 0.;
                        if (rank == 0)
                            a[p][q] = 0.;
                    }
                    else
                    if (fabs(apqdqdp[0]) > tresh)
                    {
                        const double theta = 0.5 * (apqdqdp[1] - apqdqdp[2]) / apqdqdp[0];
                        double t = 1. / (fabs(theta) + sqrt(1. + theta*theta));
                        if (theta < 0)
                            t = -t;
                        const double c = 1. / sqrt(1. + t*t);
                        const double s = t * c;
                        const double tau = s / (1. + c);
                        const double h = t * apqdqdp[0];
                        if (rank == 0){
                            z[p] -= h;
                            z[q] += h;
                            d[p] -= h;
                            d[q] += h;
                            a[p][q] = 0.;
                        }

                        if (rank == 0){
                            double gh[2];
                            double ajpq[2];

                            for (int j = 0; j < p; j++){
                                gh[0] = a[j][p];
                                gh[1] = a[j][q];
                                MPI_Send(gh, 2, MPI_DOUBLE, j % (size - 1) + 1, j, MPI_COMM_WORLD);
                            }
                            for (int j = p + 1; j < q; j++){
                                gh[0] = a[p][j];
                                gh[1] = a[j][q];
                                MPI_Send(gh, 2, MPI_DOUBLE, j % (size - 1) + 1, j, MPI_COMM_WORLD);
                            }
                            for (int j = q + 1; j < n; j++){
                                gh[0] = a[p][j];
                                gh[1] = a[q][j];
                                MPI_Send(gh, 2, MPI_DOUBLE, j % (size - 1) + 1, j, MPI_COMM_WORLD);
                            }

                            MPI_Status status;
                            for (int k = 0; k < n-2; k++){
                                MPI_Recv(ajpq, 2, MPI_DOUBLE, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
                                int jtag = status.MPI_TAG;
                                if (jtag < p){
                                    a[jtag][p] = ajpq[0];
                                    a[jtag][q] = ajpq[1];
                                }
                                else if (jtag > p && jtag < q){
                                    a[p][jtag] = ajpq[0];
                                    a[jtag][q] = ajpq[1];
                                }
                                else if (jtag > q && jtag < n){
                                    a[p][jtag] = ajpq[0];
                                    a[q][jtag] = ajpq[1];
                                }
                            }
                        }
                        else
                        {
                            double gh[2], ajpq[2];

                            for (int j = rank - 1; j < n; j += size - 1)
                            {
                                if (j != p && j != q){
                                    MPI_Recv(gh, 2, MPI_DOUBLE, 0, j, MPI_COMM_WORLD, NULL);
                                    ajpq[0] = gh[0] - s * (gh[1] + gh[0] * tau);
                                    ajpq[1] = gh[1] + s * (gh[0] - gh[1] * tau);
                                    MPI_Send(ajpq, 2, MPI_DOUBLE, 0, j, MPI_COMM_WORLD);
                                }
                            }
                        }
                    }
                }
            }


            if (rank == 0)
                for (int p = 0; p < n; ++p){
                    d[p] = (b[p] += z[p]);
                    z[p] = 0.;
                }
        }

        free(z);
        free(b);

//////////////////////////////////////////////////////////////////////

        if (rank == 0){
            t2 = clock();

            double seconds = ((double)t2 - t1) / ((double)CLOCKS_PER_SEC);

            FILE* f_outp = fopen("outp.txt", "a");
            fprintf(f_outp, "* %f\n", seconds);
            cout << '*' << seconds << endl;
            fclose(f_outp);

            /* вывод результата: */
            //for (int i = 0; i < num; i++)
            //cout << d[i] << endl;
        }
        MPI_Barrier(MPI_COMM_WORLD);
    }
    MPI_Finalize();
}