Hola tengo que hacer una multiplicación de matrices en C usando MPI(message parsing interface), donde tengo que hacer A*B=C, donde con la matriz A tengo que hacer un broadcast y con B hacer un scatter para repartirla de manera local y crear una submatriz C_local donde guardará los resultados locales y luego hacer un gather para juntar los resultados en C. Bien el producto se tiene que hacer con cblas_dgem. El caso es que me da un error de segmentation fault y no entiendo porqué. El trabajo es a repartir entre 4 nodos.
#include <mpi.h>
#include <stdlib.h>
#include <stdio.h>
#include <cblas.h>
#include <math.h> // Para fabs
// Función para generar matrices aleatorias
void generar_matriz(double* matriz, int filas, int columnas) {
for (int i = 0; i < filas * columnas; i++) {
matriz[i] = (double)rand() / RAND_MAX;
}
}
// Función secuencial para multiplicar matrices de forma clásica
void multiplicar_matrices_secuencial(double* A, double* B, double* C, int n) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
C[i * n + j] = 0.0;
for (int k = 0; k < n; k++) {
C[i * n + j] += A[i * n + k] * B[k * n + j];
}
}
}
}
// Función para comparar dos matrices y verificar si son iguales
int comparar_matrices(double* C1, double* C2, int n) {
double tolerancia = 1e-6; // Diferencia permitida por errores numéricos
for (int i = 0; i < n * n; i++) {
if (fabs(C1[i] - C2[i]) > tolerancia) {
return 0; // Las matrices no son iguales
}
}
return 1; // Las matrices son iguales
}
int main(int argc, char** argv) {
int n = atoi(argv[1]); // Toma el tamaño n de la línea de comandos
int rank, size;
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
double *A, *B, *B_local, *C_local, *C_paralelo, *C_secuencial;
if (rank == 0) {
// P0 genera las matrices A y B
A = (double*)malloc(n * n * sizeof(double));
B = (double*)malloc(n * n * sizeof(double));
generar_matriz(A, n, n);
generar_matriz(B, n, n);
// Para almacenar el resultado final de la versión paralela
C_paralelo = (double*)malloc(n * n * sizeof(double));
// Para almacenar el resultado secuencial
C_secuencial = (double*)malloc(n * n * sizeof(double));
// Multiplicación secuencial clásica para comparar
multiplicar_matrices_secuencial(A, B, C_secuencial, n);
}
B_local = (double*)malloc(n * (n / 4) * sizeof(double));
C_local = (double*)malloc(n * (n / 4) * sizeof(double));
// Distribución de la matriz B en bloques y replicación de A
MPI_Scatter(B, n * (n / 4), MPI_DOUBLE, B_local, n * (n / 4), MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(A, n * n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
// Calculo local de C_local usando dgemm de BLAS
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
n, n / 4, n,
1.0, A, n, B_local, n,
0.0, C_local, n);
// Reunión de los resultados
MPI_Gather(C_local, n * (n / 4), MPI_DOUBLE, C_paralelo, n * (n / 4), MPI_DOUBLE, 0, MPI_COMM_WORLD);
if (rank == 0) {
// Comparar el resultado paralelo con el secuencial
if (comparar_matrices(C_secuencial, C_paralelo, n)) {
printf("Las matrices calculadas por la versión secuencial y paralela son iguales.\n");
} else {
printf("Las matrices calculadas por la versión secuencial y paralela son DIFERENTES.\n");
}
}
// Liberar la memoria
if (rank == 0) {
free(A);
free(B);
free(C_paralelo);
free(C_secuencial);
}
free(B_local);
free(C_local);
MPI_Finalize();
return 0;
}