Duża macierz a czas

Witam, muszę zrobić dużą macierz - załóżmy 4000x4000 i ją przemnożyć. Zrobiłem na wskaźnikach, bo normalnie wywalało program. I teraz co? Mnoży mi się kilkanaście minut. Jest możliwość przyśpieszenia?

#include <iostream>

using namespace std;

int main(int argc, char** argv) 
{
	int m = 4000;
	
    float **tab  = new float * [m];
	float **tab2 = new float * [m];
	float **tab3 = new float * [m];
	
	for(int i = 0; i < m ; i++)
	{
		tab[i]  = new float[m];
		tab2[i] = new float[m];
		tab3[i] = new float[m];
	}
	
	for(int i = 0; i < m; i++)
	{
		for(int j = 0; j < m ; j++)
		{
			tab[i][j]  = 1;
			tab2[i][j] = 2;
			tab3[i][j] = 0;
		}
	}
	
	for(int i = 0; i < m; i++)
	{
		for(int j = 0; j < m; j++)
		{
			int s = 0;
			for(int k = 0; k < m; k++)
				s += tab[i][k] * tab2[k][j];
				tab3[i][j] = s;
		}
	}
	cout << tab3[m][m] << " ";
	
	system("pause");
	return 0;
}

Ułóż operacje tak, aby odczyty były przyjazne dla cache procesora. Zamiast używać kilku tysięcy tablic, użyj jednej, ciągłej. Jako widok możesz użyć mojego prostego szablonu:

template<typename T>
class simple_2d_matrix_view
{
	T* data_;
	size_t width_;
	size_t height_;
 
public:
 
	simple_2d_matrix_view(T* ptr, size_t h, size_t w):
		data_{ptr},
		width_{w},
		height_{h}
	{}
 
	size_t width() const { return width_; }
	size_t height() const { return height_; }
 
	T& operator()(size_t h, size_t w) {
		assert(w < width_);
		assert(h < height_);
		return data_[width_ * h + w];
	}
 
	T const& operator()(size_t h, size_t w) const {
		return const_cast<simple_2d_matrix_view&>(*this)(h, w);
	}
};

kq napisał(a):

Ułóż operacje tak, aby odczyty były przyjazne dla cache procesora. Zamiast używać kilku tysięcy tablic, użyj jednej, ciągłej. Jako widok możesz użyć:

Ale ja muszę użyć dwóch tablic, bo w jeden są wartości załóżmy 1, a w drugiej 2. I chcę to przemnożyć. Zminimalizowałem kod do takiej postaci:

int main(int argc, char** argv) 
{
	short m = 4000;
	
    float **tab  = new float * [m];
	float **tab2 = new float * [m];
	
	for(int i = 0; i < m; ++i)
	{
        tab[i]  = new float[m];
		tab2[i] = new float[m];
		for(int j = 0; j < m ; ++j)
		{
			tab[i][j]  = 1;
			tab2[i][j] = 2;
			tab[i][j] = tab[i][j] * tab2[i][j];
			cout << tab[i][j] << " ";
		}
		cout << endl;
	}
	system("pause");
	return 0;
}

A miej sobie te dwie tablice. Krzaqowi chodziło o to, żeby dwuwymiarową tablicę mieć w jednym bloku pamięci, a nie rozsypaną po RAMie.

template<typename T>
class simple_2d_matrix_view
{
    T* data_;
    size_t width_;
    size_t height_;
 
public:
 
    simple_2d_matrix_view(T* ptr, size_t h, size_t w):
        data_{ptr},
        width_{w},
        height_{h}
    {}
 
    size_t width() const { return width_; }
    size_t height() const { return height_; }
 
    T& operator()(size_t h, size_t w) {
        assert(w < width_);
        assert(h < height_);
        return data_[width_ * h + w];
    }
 
    T const& operator()(size_t h, size_t w) const {
        return const_cast<simple_2d_matrix_view&>(*this)(h, w);
    }
};

template <typename T>
simple_2d_matrix_view(T*, size_t, size_t) -> simple_2d_matrix_view<T>;

int main(int argc, char** argv) 
{
    short m = 4000;
 
    auto tab  = std::make_unique<float[]>(m*m);
    auto tab2  = std::make_unique<float[]>(m*m);
    
    auto tab_view = simple_2d_matrix_view(tab.get(), m, m);
    auto tab2_view = simple_2d_matrix_view(tab2.get(), m, m);
 
    for(int i = 0; i < m; ++i)
    {
        for(int j = 0; j < m ; ++j)
        {
            tab_view(i, j)  = 1;
            tab2_view(i, j) = 2;
            tab_view(i, j) = tab_view(i, j) * tab2_view(i, j);
            cout << tab_view(i, j) << " ";
        }
        cout << endl;
    }
    system("pause");
    return 0;
}

miałeś też wycieki pamięci, którymi teraz zajmą się unique_ptry. (może się nie kompilować, pisałem z palca)

EDIT
pozwoliłem sobie dopisać deduction guide'a

template<class T>
class Matrix
{
public:
      // rule of zero
      
      Matrix(int r, int c, T value = T())
          : data(r * c, value)
          , rowCount(r)
          , columnCount(c)
      {}

      T& operator(int r, int c)
      {
             return data[r + rowCount * c];
      }

      int RowCount() const { return rowCount; }
      int ColumnCount() const { return columnCount; }

private:
     std::vector<T> data;
     int rowCount = 0;
     int columnCount = 0;
};

Inne sztuczki przyspieszenia małym wysiłkiem

• std::ios::sync_with_stdio(false);
• zmienić każde endl na '\n'
• cin.tie(nullptr); - ale z głową, jeśli wyświetlasz prośby o dane to się nie pojawią bez poprzedniego flush.

Inna sztuczka żeby przyspieszyć obliczania to druga klasa: TransposedMatrix, żeby mnożyć liczby umieszczone w ciągłych fragmentach pamięci - zredukujesz to "cache miss" i da ekstra przyspieszenie. Albo metody mutiplyWithTransposed i ReadTransposed by uzyskać ten sam efekt.

Inne sztuczki to użycie wielowątkowości i by obliczenia prowadzić równolegle, ale to jest trudniejsze i może być już za wiele dla ciebie.

Tak duże macierze mnoży się blokowo. Przykładowo poczytaj o BLAS (Basic Linear Algebra Subprograms).

Idea jest taka żeby mnożyć małe macierze i je łączyć. Macierz dzielimy na małe blokowe macierze tak żeby mieściły się w kolejnych warstwach cache. Czyli L2, L1 i samych rejestrach procesora. Dodatkowo małe macierze pakujemy w tablicę w takiej kolejności żeby odczyt był liniowy.