Параллельные вычисления

#pragma comment( lib, "mpich.lib" )

#include "mpi.h"
#include "stdlib.h"
#include "stdio.h"
#include "math.h"

double function( double x ) {
    double f_x;
    f_x = x;
    return sin(f_x) + 1;
}

int main( int argc, char* argv[], int size, int rank ) {
    double a, b, dx, p;
    int n, i;
    long int my_n, r;
    double my_a, my_b, my_s = 0, my_ds = 0, s = 0, t1, t2;
   
    MPI_Status st;
    MPI_Init( &argc, &argv );

    MPI_Barrier( MPI_COMM_WORLD );
    t1 = MPI_Wtime();

    MPI_Comm_size( MPI_COMM_WORLD, &size ); // общее кол-во процессов
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );

    if( rank == 0 ) {                       // номер текущего процесса равен нулю
        FILE *fp;
        char s[100];
        if (( fp = fopen( "text.txt", "rt" ) ) == NULL ) {
            printf( "File not open!\n" );
            exit( 1 );
        }
        fgets( s, 99, fp );
        a = atof( s );
        fgets( s, 99, fp );
        b = atof( s );
        fgets( s, 99, fp );
        n = atoi( s );

        for( i = 1; i < size; i++ ) {       // отправляем значения длины отрезка и числа дробления во все процессы
            MPI_Send( &a, 1, MPI_DOUBLE, i, i,     MPI_COMM_WORLD );
            MPI_Send( &b, 1, MPI_DOUBLE, i, i * 2, MPI_COMM_WORLD );
            MPI_Send( &n, 1, MPI_INT,    i, i * 3, MPI_COMM_WORLD );
        }
    } else {                                // принимаем во всех процессах значения длины отрезка и числа дробления
            MPI_Recv( &a, 1, MPI_DOUBLE, 0, rank,     MPI_COMM_WORLD, &st );
            MPI_Recv( &b, 1, MPI_DOUBLE, 0, rank * 2, MPI_COMM_WORLD, &st );
            MPI_Recv( &n, 1, MPI_INT,    0, rank * 3, MPI_COMM_WORLD, &st );   
    }

    r = n % size;
    printf ( "size = %d, r = %d\n", size, r );

    dx = ( b - a ) / n;
    if( rank < r ) {                        // вычисляем значения концов отрезка и число дроблении
            my_n = n / size + 1;            // число отрезков dx в процессе
            my_a = a + rank * my_n * dx;    // начало рабочего отрезка в процессе
            my_b = my_a + dx * my_n;        // конец рабочего отрезка в процессе
    } else {
            my_n = n / size;
            my_a = a + r * dx * ( my_n + 1 ) + ( rank - r ) * my_n * dx;
            my_b = my_a + dx * my_n;
    }
   
    my_s = dx * ( function(my_a)/2 + function(my_b)/2 );
    for ( i = 1; i < my_n; i++ ) {
        my_s = my_s + dx* function(my_a+i*dx);// средний кусочек площади
    }
   
    if( rank > 0 ) {                        // отправляем вычисленное значение в нулевой процесс
            MPI_Send( &my_s, 1, MPI_DOUBLE, 0, rank, MPI_COMM_WORLD );
    } else {                                // получаем итоговый результат
        for( i = 1; i < size; i++ ) {
            MPI_Recv( &my_s, 1, MPI_DOUBLE, i, i, MPI_COMM_WORLD, &st );
            s = s + my_s;
        }
    }
   
    // погрешность
    p = ( 1.00 / ( 12.00 * n * n ) )* ( b - a )*( b - a )*( b - a );

    MPI_Barrier( MPI_COMM_WORLD );
    t2 = MPI_Wtime();
   
    if( rank == 0 ) {
        printf( "Значение интеграла функции на отрезке [%f, %f] равно %f. \nВремя работы %f\n", a, b, s, t2 - t1 );
        printf( "Погрешность вычислений: %f\n", p );
    }

    MPI_Finalize();
    return 0;
}