Symmetri c Умножение блочной матрицы - PullRequest
Новая
       38

Symmetri c Умножение блочной матрицы

0 голосов
/ 05 апреля 2020

Я пытаюсь умножить две симметричные матрицы c ( MATRIX_SIZE x MATRIX_SIZE ). Я хочу выполнить умножение блочной матрицы (разделить матрицу на несколько BLOCK_SIZE x BLOCK_SIZE матриц и умножить соответствующие блоки). Я написал некоторый код, но хочу улучшить его и хранить блоки выше основной диагонали, но у меня нет никаких идей. Ребята, можете ли вы помочь, если это возможно? 

#define IND(A, x, y) A[y*MATRIX_SIZE+x]
void block_mult2(double*& A, double*& B, double*& C){
int i, j, k, i0, j0, k0;
for (i = 0; i < MATRIX_SIZE; i += BLOCK_SIZE)
for (j = 0; j < MATRIX_SIZE; j += BLOCK_SIZE)
for (k = 0; k < MATRIX_SIZE; k += BLOCK_SIZE)
    for (i0 = i; i0 < min(BLOCK_SIZE+i, MATRIX_SIZE); i0++)
        for (j0 = j; j0 < min(BLOCK_SIZE+j, MATRIX_SIZE); j0++)
            for (k0 = k; k0 < min(BLOCK_SIZE+k, MATRIX_SIZE); k0++)
                IND(C, i0, j0) += IND(A, i0, k0) * IND(B, k0, j0);
}

Ответы [ 2 ]

0 голосов
/ 05 апреля 2020

Можете ли вы использовать существующие пакеты линейной алгебры? Если вы имеете дело с примитивными типами, например, double BLAS, вероятно, является наиболее оптимальным способом для go, но может иметь крутой кривой обучения. Для высокооптимизированной, но очень удобной библиотеки Eigen - один из моих любимых вариантов для таких задач в c ++.

Я очень рекомендую использовать существующий пакет линейной алгебры (даже не обязательно Я упомянул). Было бы легче выбросить sh ваших идей, поскольку фактическая реализация позаботилась о пакете. Не говоря уже о том, что такие пакеты существуют годами (несколько десятилетий в случае BLAS) и должны быть очень ОЧЕНЬ хороши в таких задачах. Если вы действительно не знаете, что делаете (имейте в виду очень очень конкретную задачу c с конкретными оптимизациями c, в которые вы можете кодировать) Я сомневаюсь, что вы можете оптимизировать так же легко, как и эти библиотеки, самостоятельно (если вообще) , Даже тогда, есть анализ затрат и выгод, чтобы рассмотреть: собираюсь ли я делать это сам по себе, сколько времени займет существующий хороший пакет?

Хотя я настоятельно рекомендую не делать это самостоятельно, если Вы обязательно должны сделать это сами, один вопрос, который неясен: все ли блоки одинакового размера? Кроме того, в какой форме хранятся матрицы, столбец или мажор строки? Предполагаемые блоки одинакового размера, и у вас есть основная форма строки, эскиз того, что вы могли бы сделать, это перебрать блоки и перевести умножение блочных блоков в обобщенную c функцию умножения матриц. Я сбрасываю double*& и прохожу только указатели double*. operator[] должен позаботиться о ссылке на правильное местоположение, но проверьте, что я правильно выполнил арифметику c внутри [], а также:

РЕДАКТИРОВАТЬ: Если A и B сохраняя только верхние блоки tri angular Я исправил код 

//Assuming all blocks are the same size
//Assuming matrix stored in row major form

#define NUMBER_OF_BLOCKS = MATRIX_SIZE/BLOCK_SIZE

void block_mult2(double* A, double* B, double* C){
  for(size_t i=0; i<NUMBER_OF_BLOCKS; i++)
    for(size_t j=0; j<NUMBER_OF_BLOCKS; j++)
      for(size_t k=0; k<NUMBER_OF_BLOCKS; k++)
        mult2(A[min(i,j)*BLOCK_SIZE*NUMBER_OF_BLOCKS + max(i,j)*BLOCK_SIZE],
              B[min(j,k)*BLOCK_SIZE*NUMBER_OF_BLOCKS + max(j,k)*BLOCK_SIZE],
              C[i*BLOCK_SIZE*NUMBER_OF_BLOCKS + k*BLOCK_SIZE]);
  return;
}

void mult2(double* A, double* B, double* C){
  for(size_t i=0; i<BLOCK_SIZE; i++)
    for(size_t j=0; j<BLOCK_SIZE; j++)
      for(size_t k=0; k<BLOCK_SIZE; k++)
        C[i*BLOCK_SIZE+k] = A[min(i,j)*BLOCK_SIZE+max(i,j)]*B[min(j,k)*BLOCK_SIZE+max(j,k)];
  return;
}

Не могу не подчеркнуть, насколько я рекомендую вам отбросить все это и потратить немного времени на изучение пакета линейной алгебры. Вы избавите себя от множества технических вопросов (например, только что возникших: правильно ли я сделал арифметику указателей c?), И вы могли бы использовать этот пакет для еще многих задач. Думаю, это пойдет на пользу вашей общей работе.

0 голосов
/ 05 апреля 2020 
for(int jj=0;jj<N;jj+= s){
    for(int kk=0;kk<N;kk+= s){
            for(int i=0;i<N;i++){
                    for(int j = jj; j<((jj+s)>N?N:(jj+s)); j++){
                            temp = 0;
                            for(int k = kk; k<((kk+s)>N?N:(kk+s)); k++){
                                    temp += a[i][k]*b[k][j];
                            }
                            c[i][j] += temp;
                    }
            }
     }
 }

Прошу прощения за этот фиктивный код, но вы можете считать N вашим BLOCK_SIZE

...