2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему
 
 Алгоритм Штрассена
Сообщение19.04.2014, 08:56 


20/10/12
235
Добрый день, уважаемые участники форума!
Вопрос совершенно ироничный, но у меня никак не получается заставить работать алгоритм Штрассена эффективнее перемножения матриц в лоб.
код: [ скачать ] [ спрятать ]
Используется синтаксис C++
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "header.h"

#define SIZE 512
//for square matrixes!
int main(void)
{
  double *test_a = malloc(SIZE*SIZE*sizeof(double)),*test_b = malloc(SIZE*SIZE*sizeof(double)),
  *result = malloc(SIZE*SIZE*sizeof(double)), dt;
  srand((unsigned)time(NULL));
   
  RandomInit(test_a, SIZE);
  //PrintMatrix(test_a, SIZE);
  RandomInit(test_b, SIZE);
  //PrintMatrix(test_b, SIZE);
  dt = clock();
  DefMatrixMultiply(test_a, test_b, result, SIZE);
  dt = clock()-dt;
  printf("Time is %lf:\n", dt);
  //PrintMatrix(result, SIZE);
  dt = clock();
  RecursiveMatrixMultiply(test_a, test_b, result, SIZE);
  dt = clock()-dt;
  printf("Time is %lf:\n", dt);
  //PrintMatrix(result, SIZE);
  return 0;
}

void DefMatrixMultiply(double *a, double *b, double *c, int size) //O(n^3)
{
  int i, j, k;
  for(i = 0; i < size; ++i)
   for(j = 0; j < size; ++j)
     for(c[size*i + j] = 0, k = 0; k < size; ++k)
       c[size*i + j] += a[size*i + k]*b[size*k + j];
}

void RecursiveMatrixMultiply(double *a, double *b, double *c, int size)
{
  if(size == 1)
    c[0] = a[0]*b[0];
  else
  {
    int i, j;
    double *S[10], *P[7], *sub_a[4], *sub_b[4];

    for(i = 0; i < 10; ++i)
    {
      if((S[i] = (double *)malloc(sizeof(double) * size * size/4))==NULL)
        printf("Access violation: Memory denied...\n");
      if(i < 4)
      {
        if((sub_a[i] = (double *)malloc(sizeof(double) * size * size/4))==NULL)
          printf("Access violation: Memory denied...\n");
        if((sub_b[i] = (double *)malloc(sizeof(double) * size * size/4))==NULL)
          printf("Access violation: Memory denied...\n");
      }
      if(i < 7)
        if((P[i] = (double *)malloc(sizeof(double) * size * size/4))==NULL)
          printf("Access violation: Memory denied...\n");
    }

    for(i = 0; i < size; ++i)
      for(j = 0; j < size; ++j)
        if(i < size/2 && j < size/2)
        {
          sub_a[0][(size/2)*i + j] = a[size*i + j];
          sub_b[0][(size/2)*i + j] = b[size*i + j];
        }
        else if(i < size/2 && j >= size/2)
        {
          sub_a[1][(size/2)*(i - 1) + j] = a[size*i + j];
          sub_b[1][(size/2)*(i - 1) + j] = b[size*i + j];
        }
        else if(i >= size/2 && j < size/2)
        {
          sub_a[2][(size/2)*(i - size/2) + j] = a[size*i + j];
          sub_b[2][(size/2)*(i - size/2) + j] = b[size*i + j];
        }
        else if(i >= size/2 && j >= size/2)
        {
          sub_a[3][(size/2)*(i-size/2-1) + j] = a[size*i + j];
          sub_b[3][(size/2)*(i-size/2-1) + j] = b[size*i + j];
        }

      for(i = 0; i < size/2; ++i)//parallel
       for(j = 0; j < size/2; ++j)
       {
          S[0][size/2*i+j] = sub_b[1][size/2*i+j] - sub_b[3][size/2*i+j];
          S[1][size/2*i+j] = sub_a[0][size/2*i+j] + sub_a[1][size/2*i+j];
          S[2][size/2*i+j] = sub_a[2][size/2*i+j] + sub_a[3][size/2*i+j];
          S[3][size/2*i+j] = sub_b[2][size/2*i+j] - sub_b[0][size/2*i+j];
          S[4][size/2*i+j] = sub_a[0][size/2*i+j] + sub_a[3][size/2*i+j];
          S[5][size/2*i+j] = sub_b[0][size/2*i+j] + sub_b[3][size/2*i+j];
          S[6][size/2*i+j] = sub_a[1][size/2*i+j] - sub_a[3][size/2*i+j];
          S[7][size/2*i+j] = sub_b[2][size/2*i+j] + sub_b[3][size/2*i+j];
          S[8][size/2*i+j] = sub_a[0][size/2*i+j] - sub_a[2][size/2*i+j];
          S[9][size/2*i+j] = sub_b[0][size/2*i+j] + sub_b[1][size/2*i+j];
       }
      RecursiveMatrixMultiply(sub_a[0], S[0]    , P[0], size/2);
      RecursiveMatrixMultiply(S[1]    , sub_b[3], P[1], size/2);
      RecursiveMatrixMultiply(S[2]    , sub_b[0], P[2], size/2);
      RecursiveMatrixMultiply(sub_a[3], S[3]    , P[3], size/2);
      RecursiveMatrixMultiply(S[4]    , S[5]    , P[4], size/2);
      RecursiveMatrixMultiply(S[6]    , S[7]    , P[5], size/2);
      RecursiveMatrixMultiply(S[8]    , S[9]    , P[6], size/2);
      for(i = 0; i < size; ++i)
        for(j = 0; j < size; ++j)
        {
         if(i < size/2 && j < size/2)
           c[size*i + j] = P[4][(size/2)*i + j]+P[3][(size/2)*i + j]-P[1][(size/2)*i + j]+P[5][(size/2)*i + j];
         else if(i < size/2 && j >= size/2)
           c[size*i + j] = P[0][(size/2)*(i - 1) + j] + P[1][(size/2)*(i - 1) + j];
         else if(i >= size/2 && j < size/2)
           c[size*i + j] = P[2][(size/2)*(i - size/2) + j] + P[3][(size/2)*(i - size/2) + j];
         else if(i >= size/2 && j >= size/2)
           c[size*i + j] = P[4][(size/2)*(i-size/2-1) + j] + P[0][(size/2)*(i-size/2-1) + j] - P[2][(size/2)*(i-size/2-1) + j] - P[6][(size/2)*(i-size/2-1) + j];
       }
    }
}
 

То есть он работает, но медленнее, а как только размер матрицы доходит до 512 и больше(где он и должен взять свое) переполняется память, которую я в этом алгоритме не умею очищать.
Пожалуйста, посоветуйте, как улучшить реализацию и очищать здесь память.
в header.h прячется:
код: [ скачать ] [ спрятать ]
Используется синтаксис C++
void DefMatrixMultiply(double *a, double *b, double *c, int size);
void RecursiveMatrixMultiply(double *a, double *b, double *c, int size);

void PrintMatrix(double *a, int size)
{
  int i, j;
  for(i = 0; i < size; ++i)
  {
    for(j = 0; j < size; ++j)
      printf("%.3lf ", a[size*i + j]);
    printf("\n");
  }
  printf("\n");
}

void RandomInit(double *a, int size)
{
  int i, j;
  for(i = 0; i < size; ++i)
    for(j = 0; j < size; ++j)
       a[size*i + j] = 0.5*(double)rand()/RAND_MAX;
}

 

 Профиль  
                  
 
 Re: Алгоритм Штрассена
Сообщение19.04.2014, 12:18 
Заслуженный участник
Аватара пользователя


06/10/08
6422
Простейшая оптимизация: не обязательно вычислять эти матрицы S все сразу, можно вычислить нужные для одного маленького произведения, перемножить и использовать ту же память для следующего. sub_a и sub_b вообще не нужны, кроме тех, которые используются в маленьких произведениях (при желании даже без них можно обойтись, выбирая нужные элементы из исходных матриц, но это уже сложнее).

 Профиль  
                  
 
 Re: Алгоритм Штрассена
Сообщение19.04.2014, 15:38 
Заслуженный участник


04/05/09
4582
Алгоритм Штрассена выгоден только для матриц, больших некоторого размера (порядка сотни), поэтому для умножения подматриц меньшего размера надо использовать обычный алгоритм, а не использовать рекурсию до конца, как у вас.
Ещё надо учесть, что непоследовательное обращение к памяти (по столбцам в одном из операндов обычного умножения) заметно медленнее, если матрица не влезает в кеш процессора (8-64 Кб). Поэтому для более эффективной работы с памятью стоит сначала транспонировать вторую матрицу, и модифицировать алгоритмы с учётом этого.

 Профиль  
                  
 
 Re: Алгоритм Штрассена
Сообщение19.04.2014, 16:43 


20/10/12
235
venco, спасибо, в принципе я добился, чего хотел
Xaositect, спасибо и вам, ваши советы мне предстоит понять

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 4 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group