[C - Openmp] Calcolo determinante in modo parallelo è troppo lento
Ultimamente ho iniziato a guardare la libreria openmp per la programmazione multithreading in C. Dopo gli esempi banali, ho provato a parallelizzare in modo semplice un mio programma per il calcolo del determinante di una matrice, aggiungendo la parallelizzazione del calcolo dei determinanti dei minori:
Il programma fornisce il risultato corretto, ma provandolo anche con un singolo thread (compilandolo senza fornigli l'opzione -fopenmp), ho notato che va molto più velocemente in questo caso che nel caso parallelo, ho misurato ad esempio i seguenti tempi;
Qualcuno sa trovare l'errore o darmi una motivazione di questo cosi vistoso rallentamento?
#include "stdio.h"
#include "math.h"
// SubMatrix Without Raw I And Column J
void submatrix (double matrix[], double m_tmp[], int length, int I, int J);
// Matrix Determinant
double matrix_det (double matrix[], int length);
// Main
int main () {
int LENGTH = 12;
double matrix[LENGTH*LENGTH];
int i = 0;
// Initializing The Matrix
for (i = 0; i < LENGTH*LENGTH; i++) {
matrix[i] = i*i;
}
printf ("\n");
// Determinant Calculation
double det = matrix_det (matrix, LENGTH);
printf ("Det = %lf\n", det);
return 0;
}
// Matrix Determinant
/*
* Calculate The Determinant Of A Matrix
*
*/
double matrix_det (double matrix[], int length) {
double det = 0.0;
if (length > 0) {
// Length == 1, Return The Only Number
if (length == 1) {
det = matrix[0];
}
// Length == 2, Apply Definition
if (length == 2) {
det = matrix[0]*matrix[3] - matrix[1]*matrix[2];
}
// Length > 2, Laplace Rule
if (length > 2) {
// TMP Matrix To Hold The Current Submatrix
double m_tmp[(length-1)*(length-1)];
// Foor Loop Variable
int z = 0;
// TMP Matrix Value
double m_value = 0.0;
// Single Minor Computation
double res = 0.0;
#pragma omp parallel for private(z,m_value,res) reduction(+:det)
for (z = 0; z < length; z++) {
// Submatrix Calculation
submatrix (matrix, m_tmp, length, 0, z);
// m_value = pow (-1, z)*matrix[z]
if ( z%2 == 0 ) {
m_value = matrix[z];
} else {
m_value = -matrix[z];
}
// Determinant Calculation
// Separated From det To Minimize Time In Critical Section
res = m_value*matrix_det(m_tmp, (length - 1));
// Adding Partial Result To Total Determinant
det += res;
}
}
}
return det;
}
// SubMatrix Without Row I And Column J
void submatrix (double matrix[], double m_tmp[], int length, int I, int J) {
// Virtual Row Index
int i = 0;
// Virtual Column Index
int j = 0;
// TMP Matrix Current Position
int a = 0;
// Matrix Current Position
int b = 0;
for (i = 0; i < length; i++) {
if (i != I) {
for (j = 0; j < length; j++) {
if (j != J) {
m_tmp[a] = matrix[b];
a++; // Next matrix_tmp position
b++; // Next matrix position
} else { // Next matrix position
b++;
}
}
} else {
b += length; // Next matrix position, skipping the whole row
}
}
}
Il programma fornisce il risultato corretto, ma provandolo anche con un singolo thread (compilandolo senza fornigli l'opzione -fopenmp), ho notato che va molto più velocemente in questo caso che nel caso parallelo, ho misurato ad esempio i seguenti tempi;
Single Thread - MATRIX DIMENSION = 12 real 0m12.517s user 0m12.500s sys 0m0.000s Multithread (4 Threads) - MATRIX DIMENSION = 12 real 0m20.521s user 0m52.420s sys 0m25.340s Single Thread - MATRIX DIMENSION = 13 real 2m41.931s user 2m41.728s sys 0m0.004s Multithread (4 Threads) - MATRIX DIMENSION = 13 real 5m6.844s user 11m18.524s sys 5m26.280s
Qualcuno sa trovare l'errore o darmi una motivazione di questo cosi vistoso rallentamento?
Risposte
A prima vista direi che il problema potrebbe essere la copia delle sotto-matrici.
"apatriarca":
A prima vista direi che il problema potrebbe essere la copia delle sotto-matrici.
Se ti riferisci al fatto che la matrice temporanea è una variabile condivisa tra tutti i thread e non privata per ognuno, ho provato ad inserire la definizione direttamente nel ciclo for (non si può definire come private un vettore in openmp), ma non cambia nulla, da sempre gli stessi tempi.
L'algoritmo scelto non è particolarmente efficiente e neanche particolarmente cache-friendly. Inoltre usa moltissima memoria.
Un metodo più veloce è l'uso della fattorizzazione LU. Oggi pomeriggio ho buttato giù questo. Non sono un esperto di OpenMP e sono un po' fuori allenamento con la decomposizione. Quindi non sono sicuro che sia corretta e con la matrice che hai scelto mi da problemi per grosse dimensioni. Sul mio pc calcola il determinante per una 2300x2300 in 6s contro i 16s senza OpenMP. In 65s calcola il determinante dell'identità di dimensione 5000x5000 (in debug). In release va ovviamente anche meglio (lo riduce di circa 1/3). Comunque penso si possa fare molto meglio.
Un metodo più veloce è l'uso della fattorizzazione LU. Oggi pomeriggio ho buttato giù questo. Non sono un esperto di OpenMP e sono un po' fuori allenamento con la decomposizione. Quindi non sono sicuro che sia corretta e con la matrice che hai scelto mi da problemi per grosse dimensioni. Sul mio pc calcola il determinante per una 2300x2300 in 6s contro i 16s senza OpenMP. In 65s calcola il determinante dell'identità di dimensione 5000x5000 (in debug). In release va ovviamente anche meglio (lo riduce di circa 1/3). Comunque penso si possa fare molto meglio.
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <omp.h>
double my_abs(double const a)
{
return (a < 0.)? -a : a;
}
double matrix_det(double const m[], uint32_t const n);
double factorize (double m[], uint32_t const n);
// Funzioni per inizializzare la matrice
void set_matrix(double m[], uint32_t const n);
void set2identity(double m[], uint32_t const n);
int main()
{
uint32_t const MatDim = 2300;
uint32_t const MatSize = MatDim * MatDim;
double * matrix = calloc(MatSize, sizeof(double));
if(matrix)
{
set_matrix(matrix, MatDim);
//set2identity(matrix, MatDim);
double const det = matrix_det(matrix, MatDim);
printf("Det = %f\n", det);
free(matrix);
}
else
{
perror("calloc() :");
}
}
void set_matrix(double m[], uint32_t const n)
{
uint32_t const n2 = n*n;
#pragma omp parallel for
for (uint32_t i = 0; i < n2; ++i)
m[i] = (i*i + 0.)/(n+i);
}
void set2identity(double m[], uint32_t const n)
{
uint32_t const n2 = n*n;
#pragma omp parallel for
for (uint32_t i = 0; i < n2; i+=(n+1))
m[i] = 1.;
}
double matrix_det(double const m[], uint32_t const n)
{
if(n>3)
{
uint32_t const n2 = n*n;
double * matrix = calloc(n2, sizeof(double));
if(matrix)
{
// copia della matrice per non sovrascriverla
#pragma omp parallel for
for (uint32_t i = 0; i < n2; ++i)
matrix[i] = m[i];
double det = factorize(matrix, n);
if(det != 0)
{
#pragma omp parallel for reduction(*:det)
for (uint32_t i = 0; i < n; ++i)
{
uint32_t const ind = i * n + i;
double value = matrix[ ind ];
det *= value;
}
}
free(matrix);
return det;
}
else
{
perror("calloc() :");
return 0.;
}
}
else if(n==3)
{
return ( m[0]*(m[4] * m[8] - m[5] * m[7])
- m[1]*(m[3] * m[8] - m[5] * m[6])
+ m[2]*(m[3] * m[7] - m[4] * m[6]) );
}
else if(n == 2)
{
return (m[0] * m[3] - m[1] * m[2]);
}
else if(n == 1)
{
return m[0];
}
else return 0.;
}
double factorize (double m[], uint32_t const n)
{
double parity = 1.;
uint32_t const n2 = n*n;
for( uint32_t k = 0; k < n-1; ++k )
{
// ricerca massimo
double max = my_abs( m[k*n + k] );
uint32_t ip = k;
for(uint32_t i = k+1; i < n; ++i)
{
double const t = my_abs( m[i*n + k] );
if(t > max)
{
max = t;
ip = i;
}
}
if(max == 0.)
{
return 0.;
}
if(ip != k)
{
parity *= -1;
#pragma omp parallel for
for(uint32_t i = 0; i < n; ++i)
{
double const t = m[k*n + i];
m[k*n + i] = m[ip*n + i];
m[ip*n + i] = t;
}
}
double const pe = m[k * n + k];
#pragma omp parallel for
for(uint32_t i = k+1; i < n; ++i)
m[i*n + k] /= pe;
#pragma omp parallel for
for(uint32_t i = k+1; i < n; ++i)
{
double const mik = m[i*n + k];
for(uint32_t j = k+1; j<n; ++j)
{
m[i*n + j] -= mik * m[k*n + j];
}
}
}
return parity;
}
"vict85":
L'algoritmo scelto non è particolarmente efficiente e neanche particolarmente cache-friendly. Inoltre usa moltissima memoria.
Un metodo più veloce è l'uso della fattorizzazione LU. Oggi pomeriggio ho buttato giù questo. Non sono un esperto di OpenMP e sono un po' fuori allenamento con la decomposizione. Quindi non sono sicuro che sia corretta e con la matrice che hai scelto mi da problemi per grosse dimensioni. Sul mio pc calcola il determinante per una 2300x2300 in 6s contro i 16s senza OpenMP. In 65s calcola il determinante dell'identità di dimensione 5000x5000 (in debug). In release va ovviamente anche meglio (lo riduce di circa 1/3). Comunque penso si possa fare molto meglio.
Grazie per la risposta[nota]Quando ho tempo guardo bene la tua implementazione pareallela della decomposizione LU[/nota].
So che l'algoritmo scelto non è dei migliori, ma dato che lo avevo già implementato tempo fa in modo sequenziale, ho pensato che fosse un ottimo punto di partenza. Consideralo un toy problem del tipo parallelizza questo algoritmo già fatto (nota anche che non è presente praticamente nessuna ottimizzazione o trucchetto strano, diciamo che è un'implementazione out of the book).
Deve esserci qualche errore nella logica della gestione della parallelizzazione, perchè l'algoritmo è lo stesso e l'unica differenza è che il calcolo dei determinanti dei minori avviene in parallelo e non in modo sequenziale.
La questione della copia delle sottomatrici penalizza entrambi i due modi di esecuzione (sarebbe più rapida una riduzione in loco), quindi non dovrebbe essere questo il problema.
Quell'algoritmo non è parallelizzabile, specialmente usando OpenMP. Insomma puoi parallelizzarlo usando pthread, MPI o TBB, ma non saprei dire con che risultati. Insomma dovresti farlo gestendo manualmente sincronizzazione e memoria di ogni thread. È sbagliata l'idea secondo cui prima fai le cose in modo che funzionino e poi parallelizzi, devi direttamente pensare il problema in termini paralleli (seppur possa essere comodo scrivere prima una versione sequenziale adatta ad essere parallelizzata).
Detto questo l'unico modo che mi viene in mente per renderlo parallelizzabile è girare l'algoritmo e usare un po' di memoria aggiuntiva. In questo modo si può rendere l'algoritmo più parallelizzabile (nel ciclo più interno) e di complessità \(\displaystyle O(2^n) \) invece che \(\displaystyle O(n!) \), ma usando \(\displaystyle O(n!) \) memoria.
Consideriamo infatti i minori \(\displaystyle 2\times 2 \) che stanno nelle ultime due righe. I determinanti di questi minori sono calcolati tantissime volte. L'idea consiste quindi nel calcolare questi determinanti una volta sola (ogni calcolo è indipendente cosa che rende le cose ottimali). A questo punto puoi calcolare i determinanti dei minori \(\displaystyle 3\times 3 \) usando i valori appena trovati. Ci sono dei particolari di questo sistema a cui bisogna fare un po' di attenzione per non moltiplicare troppo la complessità (devi poter trovare velocemente l'indice dell'elemento che stai cercando). In ogni caso anche questi determinanti sono indipendenti e quindi si può usare OpenMP per parallelizzarli. A questo punto i minori \(\displaystyle 2\times 2 \) non ti sono più necessari quindi puoi riscriverci sopra. E così via fino a che non calcoli il determinante.
Detto questo l'unico modo che mi viene in mente per renderlo parallelizzabile è girare l'algoritmo e usare un po' di memoria aggiuntiva. In questo modo si può rendere l'algoritmo più parallelizzabile (nel ciclo più interno) e di complessità \(\displaystyle O(2^n) \) invece che \(\displaystyle O(n!) \), ma usando \(\displaystyle O(n!) \) memoria.
Consideriamo infatti i minori \(\displaystyle 2\times 2 \) che stanno nelle ultime due righe. I determinanti di questi minori sono calcolati tantissime volte. L'idea consiste quindi nel calcolare questi determinanti una volta sola (ogni calcolo è indipendente cosa che rende le cose ottimali). A questo punto puoi calcolare i determinanti dei minori \(\displaystyle 3\times 3 \) usando i valori appena trovati. Ci sono dei particolari di questo sistema a cui bisogna fare un po' di attenzione per non moltiplicare troppo la complessità (devi poter trovare velocemente l'indice dell'elemento che stai cercando). In ogni caso anche questi determinanti sono indipendenti e quindi si può usare OpenMP per parallelizzarli. A questo punto i minori \(\displaystyle 2\times 2 \) non ti sono più necessari quindi puoi riscriverci sopra. E così via fino a che non calcoli il determinante.
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo