Decomposição de Cholesky com OpenMP

Eu tenho um projeto onde nós resolvemos o inverso de matrizes densas definidas positivas (sobre 3000×3000) grandes usando a decomposição de Cholesky . O projeto está em Java e usamos estão usando a biblioteca CERN Colt BLAS . A criação de perfil do código mostra que a decomposição de Cholesky é o gargalo.

Eu decidi tentar paralelizar a decomposição de Cholesky usando o OpenMP e usá-lo como uma DLL em Java (com JNA). Comecei com o código de decomposição Cholesky em C do Rosetta Code .

O que eu notei é que os valores em uma coluna, exceto o elemento diagonal, são independentes. Então decidi calcular os elementos diagonais em série e o resto dos valores da coluna em paralelo. Eu também troquei a ordem dos loops para que o loop interno fosse executado sobre as linhas e o loop externo sobre as colunas. A versão serial é um pouco mais lenta que a do RosettaCode, mas a versão paralela é seis vezes mais rápida que a versão RosettaCode no meu sistema de 4 núcleos (8 HT). Usando a DLL em Java acelera nossos resultados por seis vezes também. Aqui está o código:

double *cholesky(double *A, int n) { double *L = (double*)calloc(n * n, sizeof(double)); if (L == NULL) exit(EXIT_FAILURE); for (int j = 0; j <n; j++) { double s = 0; for (int k = 0; k < j; k++) { s += L[j * n + k] * L[j * n + k]; } L[j * n + j] = sqrt(A[j * n + j] - s); #pragma omp parallel for for (int i = j+1; i <n; i++) { double s = 0; for (int k = 0; k < j; k++) { s += L[i * n + k] * L[j * n + k]; } L[i * n + j] = (1.0 / L[j * n + j] * (A[i * n + j] - s)); } } return L; } 

Você pode encontrar o código completo para testar isso em http://coliru.stacked-crooked.com/a/6f5750c20d456da9

Inicialmente, pensei que o compartilhamento falso seria um problema quando os elementos restantes de uma coluna eram pequenos em comparação com o número de encadeamentos, mas esse não parece ser o caso. eu tentei

 #pragma omp parallel for schedule(static, 8) // a cache line is 8 doubles 

Eu não encontrei exemplos claros de como paralelizar a decomposição de Choleskey. Não sei se o que fiz é ideal. Por exemplo, funcionará bem em um sistema NUMA?

Talvez uma abordagem baseada em tarefas seja melhor em geral? Nos slides 7-9 em http://courses.engr.illinois.edu/cs554/fa2013/notes/07_cholesky.pdf, há um exemplo de decomposição em cholesky paralelo usando “tarefas refinadas”. Não está claro para mim como implementar isso ainda.

Eu tenho duas perguntas, específicas e gerais. Você tem alguma sugestão sobre como melhorar minha implementação da Decomposição de Cholesky com o OpenMP? Você pode sugerir uma implementação diferente da Decomposição de Cholesky com o OpenMP, por exemplo, com tarefas?

Edit: conforme solicitado aqui é a function AVX que eu usei para calcular s . Não ajudou

 double inner_sum_AVX(double *li, double *lj, int n) { __m256d s4; int i; double s; s4 = _mm256_set1_pd(0.0); for (i = 0; i < (n & (-4)); i+=4) { __m256d li4, lj4; li4 = _mm256_loadu_pd(&li[i]); lj4 = _mm256_loadu_pd(&lj[i]); s4 = _mm256_add_pd(_mm256_mul_pd(li4, lj4), s4); } double out[4]; _mm256_storeu_pd(out, s4); s = out[0] + out[1] + out[2] + out[3]; for(;i<n; i++) { s += li[i]*lj[i]; } return s; } 

    Eu consegui fazer o SIMD trabalhar com a decomposição de Cholesky. Eu fiz isso usando o ladrilho como eu usei antes na multiplicação de matrizes. A solução não foi trivial. Aqui estão os tempos para uma matriz de 5790×5790 no meu sistema Ivy Bridge de 4 núcleos / 8 HT (eff = GFLOPS / (pico GFLOPS)):

     double floating point peak GFLOPS 118.1 1 thread time 36.32 s, GFLOPS 1.78, eff 1.5% 8 threads time 7.99 s, GFLOPS 8.10, eff 6.9% 4 threads+AVX time 1.36 s, GFLOPS 47.64, eff 40.3% 4 threads MKL time 0.68 s, GFLOPS 95.14, eff 80.6% // from LAPACKE_dpotrf single floating point peak GFLOPS 236.2 1 thread time 33.88 s, GFLOPS 1.91, eff 0.8% 8 threads time 4.74 s, GFLOPS 13.64, eff 5.8% 4 threads+AVX time 0.78 s, GFLOPS 82.61, eff 35.0% 

    O novo método é 25 vezes mais rápido para o dobro e 40 vezes mais rápido para o single. A eficiência é de cerca de 35-40% do pico FLOPS agora. Com a multiplicação de matrizes, chego a 70% com o AVX no meu próprio código. Eu não sei o que esperar da decomposição de Cholesky. O algoritmo é parcialmente serial (ao calcular o bloco diagonal, chamado triangle no meu código abaixo) ao contrário da multiplicação de matrizes.

    Atualização: estou dentro de um fator para 2 dos MKL. Eu não sei se eu deveria estar orgulhoso disso ou envergonhado por isso, mas aparentemente meu código ainda pode ser melhorado significativamente. Eu encontrei uma tese de doutorado sobre isso, que mostra que o meu algoritmo de bloco é uma solução comum, então eu consegui reinventar a roda.

    Eu uso 32×32 tiles para double e 64×64 tiles para float. Eu também reordenar a memory para que cada bloco seja contíguo e seja sua transposição. Eu defini uma nova function de produção de matriz. A multiplicação de matrizes é definida como:

     C_i,j = A_i,k * B_k,j //sum over k 

    Eu percebi que no algoritmo de Cholesky há algo muito similar

     C_j,i = A_i,k * B_j,k //sum over k 

    Ao escrever a transposição dos tiles, eu pude usar quase exatamente a minha function otimizada para multiplicação de matrizes (eu só tive que mudar uma linha de código). Aqui está a function principal:

     reorder(tmp,B,n2,bs); for(int j=0; j 

    Aqui estão as outras funções. Eu tenho seis funções de produto para SSE2, AVX e FMA, cada uma com a versão dupla e flutuante. Eu só mostro o da AVX e dobro aqui:

     template  void triangle(Type *A, int n) { for (int j = 0; j < n; j++) { Type s = 0; for(int k=0; k void block(Type *A, Type *B, int n) { for (int j = 0; j  void reorder(Type *A, Type *B, int n, int bs) { int nb = n/bs; int stride = bs*bs; //printf("%d %d %d\n", bs, nb, stride); #pragma omp parallel for schedule(static) for(int i=0; i void reorder_inverse(Type *A, Type *B, int n, int bs) { int nb = n/bs; int stride = bs*bs; //printf("%d %d %d\n", bs, nb, stride); #pragma omp parallel for schedule(static) for(int i=0; i