Ilya Tychinin (ejuo) wrote,
Ilya Tychinin
ejuo

Categories:
  • Music:

Умножение разреженных матриц

Как-то реализовывал параллельное умножение/сложение разреженных матриц. Побочным эффектом у меня вышел алгоритм, который совершенно не параллелится, но, мне кажется, тоже имеющий право на существование (не смотря на недостатки мне он очень глянется :) ).
Внутреннее представление матриц выглядит следующим образом:Внутреннее представление матриц выглядит следующим образом:



Имеется класс Matr:

class element
{
public:
int i,j;
float el;
};
class Matr
{
public:
element *matrix;
Matr(int s);
};


В котором осуществляется хранение элементов в форме массива, доступ к элементу происходит через отдельный класс. Нулевой элемент содержит размер разреженной матрицы.
Имеем исходную формулу

cij =

aikbkj

 

k

 

Как не трудно увидеть, умножение происходит по равным (k) колоночным и строчным индексам соответственно элементов aik и bkj. Будем выбирать такие элементы, при этом запоминая индекс, который «помнит» строчные(i) и колоночные(j) индексы. При этом он будет сохраняться в исходном состоянии, если (i) и (j) при следующих проходах цикла ему уже встречались (рассматриваемое представление хранит элементы в неупорядоченном(как есть) виде). Это отсечёт лишние проходы по циклу. Индекс представлен на основе односвязного списка.
 

 


int Size = (int)(Mat1.matrix[0].el);
int i,j;
int ii=0;
bool flag;
node *ch, *st, *ptr, *add;
ch = new node;
ch->i=0;
ch->j=0;
ch->ind=0;
ch->nextind=NULL;
st = ch;
for (i = 1;i<=Size;i++)
{
Rez.matrix[i].el = 0;
for (j = 1;j<=Size;j++)
{
if(Mat1.matrix[i].j == Mat2.matrix[j].i)
{
ch = st;
flag = false;
while(ch->nextind != NULL)
{

if(ch->i == Mat1.matrix[i].i && ch->j == Mat2.matrix[j].j)
{
flag = true;
ptr = ch;
}
ch=ch->nextind;
}
if(flag == false)
{
ii++;
add = new node;
add->i = Mat1.matrix[i].i;
add->j = Mat2.matrix[j].j;
add->ind = ii;
add->nextind = NULL;
ch->nextind = add;
}
ptr = ch->nextind;
Rez.matrix[ptr->ind].el=Rez.matrix[ptr->ind].el+Mat1.matrix[i].el*Mat2.matrix[j].el;
Rez.matrix[ptr->ind].i = Mat1.matrix[i].i;
Rez.matrix[ptr->ind].j = Mat2.matrix[j].j;
}
}
}
Проход по двойному циклу и через условие по не полному тройному даёт впечатляющие результаты. Если в матрице достаточно мало количество элементов, дающих в результате произведения не ноль, то время исполнения алгоритма не превосходит стандартную реализацию. Что можем пронаблюдать в следующей таблице (SPARC 64 1.08GHz, 32Gb, Solaris), время в сек:


размер

время

512x512

0,002211

768x768

0,004948

1024x1024

0,008891

1280x1280

0,013835

1536x1536

0,019734

1792x1792

0,026843

2048x2048

0,03511

2304x2304

0,044345

2560x2560

0,054737

2816x2816

0,066221

3072x3072

0,078796

 

 Для сравнения, таблица с классическим умножением плотных матриц в один поток:

 

размер

(k,i)*(k,j)

(i,k)*(k,j)

(i,k)*(j,k)

   

256x256

0,297322

0,179069

0,10361

 

512x512

2,44983

1,77681

0,813149

Статич.

768x768

8,12796

5,11412

3,52115

 

1024x1024

128,949

70,6321

10,0176

 

1280x1280

392,13

25,1306

18,8094

 

256x256

0,275139

0,248916

0,240391

 

512x512

2,33973

1,98403

1,91013

Динамич.

768x768

7,29593

6,7515

6,48599

 

1024x1024

24,9564

16,6434

15,3675

 

1280x1280

407,498

34,2095

30,3765



Формирование конечной матрицы происходит динамически, что даёт некоторые преимущества (проход только по индексу, который выбирает те элементы, которые дадут не нулевой результат), но и имеет недостатки.
Попытки распараллеливания не дают результатов, так как каждой нити(thread) требуются сведения о текущем индексе в конечной матрице. Что требует постоянной синхронизации и дополнительного времени.

 
Tags: c++, prog
Subscribe
  • Post a new comment

    Error

    Comments allowed for friends only

    Anonymous comments are disabled in this journal

    default userpic

    Your reply will be screened

    Your IP address will be recorded 

  • 0 comments