Skip to content

Instantly share code, notes, and snippets.

@strrife
Created May 20, 2015 20:15
Show Gist options
  • Save strrife/181e6848e87367a5c66a to your computer and use it in GitHub Desktop.
Save strrife/181e6848e87367a5c66a to your computer and use it in GitHub Desktop.
#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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment