gistfile1.c

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <iostream>
#include "mpi.h"

unsigned long micros = 0;
float millis = 0.0;
clock_t start, end;

/*   run it with
 *   cd ~/Cpp/Architecture_Lab_3; mpicxx main.c -o lab_3; mpirun -n 1 ./lab_3
 */

#define MATRIX_SIZE 1000
#define MAX_MATRIX_SIZE_TO_SHOW_STATS 10
#define MAX_ABS_VALUE 10

double a[MATRIX_SIZE][MATRIX_SIZE];
double L[MATRIX_SIZE][MATRIX_SIZE];
double U[MATRIX_SIZE][MATRIX_SIZE];
int proc_count;

void generate_matrix()
{
    int i, j;
    for(i = 0; i < MATRIX_SIZE; ++i){
        for(j = 0; j < MATRIX_SIZE; ++j){
            a[i][j] = rand() % MAX_ABS_VALUE + 1;
        }

    }
}

void output_matrix(double a[MATRIX_SIZE][MATRIX_SIZE]){
    for(int i = 0; i < MATRIX_SIZE; ++i){
        for(int j = 0; j < MATRIX_SIZE; ++j){
            printf("%7.3f ", a[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}

void output_matrix_wolfram(double a[MATRIX_SIZE][MATRIX_SIZE]){
    printf("\n{");
    for(int i = 0; i < MATRIX_SIZE; ++i){
        printf("{");
        for(int j = 0; j < MATRIX_SIZE; ++j){
            printf("%f ", a[i][j]);
            if(j < MATRIX_SIZE - 1)
                printf(",");
        }
        printf("}");
        if(i < MATRIX_SIZE - 1)
            printf(",");
    }
    printf("}\n\n");
}

int get_number_of_rows_this_process_is_responsible_for(int proc_num) {
    int avg = (int)floor(1. * MATRIX_SIZE / proc_count);
    int remainder = MATRIX_SIZE - proc_count * avg;
    return remainder > proc_num ? avg + 1 : avg;
}

void transpose(){
    double temp;
    for(int i = 0; i < MATRIX_SIZE; ++i)
        for(int j = i; j < MATRIX_SIZE; ++j){
            temp = a[i][j];
            a[i][j] = a[j][i];
            a[j][i] = temp;
        }
}

void print_multiply_matrices(double A[MATRIX_SIZE][MATRIX_SIZE], double B[MATRIX_SIZE][MATRIX_SIZE]) {
    double C[MATRIX_SIZE][MATRIX_SIZE];
    int i, j, k;
    for (i = 0; i < MATRIX_SIZE; i++) {
        for (j = 0; j < MATRIX_SIZE; j++) {
            C[i][j] = 0;
            for (k=0; k < MATRIX_SIZE; k++)
                C[i][j] += A[i][k]*B[k][j];
        }
    }
    output_matrix(C);
}


int main(int argc, char *argv[])
{ 
    start = clock();

    int n = MATRIX_SIZE; //TODO: move s param

    int rank;
    int i, j, k;

    double * data_flattened = new double[MATRIX_SIZE * MATRIX_SIZE];

    //which process whould take care of nth line
    int * map = new int[n];
    int * displs = new int[n];
    int * counts = new int[n];

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &proc_count);

    //init matrix
    if(rank == 0){
        generate_matrix();
        if(MATRIX_SIZE <= MAX_MATRIX_SIZE_TO_SHOW_STATS){
            printf("\nA:\n");
            output_matrix(a);
            transpose();
        }
    }

    //distribute rows between processed
    for(i = 0; i < n;){
        int starting_position = 0;
        for(j = 0; j < proc_count; ++j){
            for(int ii = 0; ii < get_number_of_rows_this_process_is_responsible_for(j); ++ii){
                map[i++] = j;
            }
            counts[j] = MATRIX_SIZE * get_number_of_rows_this_process_is_responsible_for(j);
            displs[j] = starting_position;
            starting_position += MATRIX_SIZE * get_number_of_rows_this_process_is_responsible_for(j);
        }
    }

    //send to all other processes
    MPI_Bcast(a, n*n, MPI_DOUBLE, 0, MPI_COMM_WORLD);

    for(k = 0; k < MATRIX_SIZE; k++)
    {
        if(map[k] == rank){
            for(i = k + 1; i < MATRIX_SIZE; i++)
            {
                //if we're responsible for k, change it
                a[k][i] /= a[k][k];
            }
        }
        MPI_Bcast(&a[k][k + 1], n - k - 1, MPI_DOUBLE,  map[k], MPI_COMM_WORLD);

        for(i = k + 1; i < MATRIX_SIZE; i++)
        {
            if(map[i] == rank){
                for(j = k + 1; j < MATRIX_SIZE; j++) {
                    a[i][j] -= a[i][k] * a[k][j];
                }
            }
        }
    }

    // Printing the entries of the matrix
    i = 0;
    for(j = 0; j < proc_count; ++j){
        int row_count = get_number_of_rows_this_process_is_responsible_for(j);
        if(rank == j){
            MPI_Gatherv(&a[i], MATRIX_SIZE * row_count, MPI_DOUBLE,
                       data_flattened, counts, displs,
                        MPI_DOUBLE, 0, MPI_COMM_WORLD);
        }
        i += row_count;
    }

    MPI_Barrier(MPI_COMM_WORLD);

    if(rank == 0){
        if(MATRIX_SIZE <= MAX_MATRIX_SIZE_TO_SHOW_STATS){
            for(i = 0; i < MATRIX_SIZE; ++i){
                for(j = 0; j < MATRIX_SIZE; ++j){
                    if(i == j){
                        L[i][j] = 1;
                        U[i][j] = data_flattened[j * MATRIX_SIZE + i];
                    } else if(i > j){
                        L[i][j] = data_flattened[j * MATRIX_SIZE + i];
                        U[i][j] = 0;
                    } else {
                        L[i][j] = 0;
                        U[i][j] = data_flattened[j * MATRIX_SIZE + i];
                    }
                }
            }
            printf("\nL:\n");
            output_matrix(L);
            printf("\nU:\n");
            output_matrix(U);
            printf("\nLU:\n");
            print_multiply_matrices(L, U);
        }

        end = clock();
        micros = end - start;
        std::cout << "Потрачено времени: " << micros /1000 << " мсекунд" << std::endl;
    }

    MPI_Finalize();
    return 0;
}