ускорить вычисление адреса

а что оно должно делать?

процедура по заданным векторам b_row и b_col, а так же ряду других параметров, по некоторым формулам (математические выражения) заполняет матрицу matrix. В b_row и b_col координаты

 
x_{00},x_{01}, x_{10},x_{11}, ...x_{i0},x_{i1},...x_{n_row-1 0},x_{n_row-1 1}

Как-нибудь так:

#pragma omp parallel for
   for(i=0; i<n_row-1; ++i) {
       int ii = 2*i + b_row;
       xip0=*(ii);
       xip1=*(ii+1);
       xin0=*(ii+2);
       xin1=*(ii+3);
        for(int j=0; j<n_col-1; ++j){
           int jj = b_col + 2*j;
           xjp0=*(jj);
           xjp1=*(jj+1);
           xjn0=*(jj+2);
           xjn1=*(jj+3);

Думаю, что не поможет.

Если этих операций нет в mkl и подобных библиотеках, то единственный выход — заюзать sse и вручную всё оптимизировать.

Ну, у меня в подобных вычислениях до трехкратного выигрыша в скорости дает на четырех ядрах. На GPU было бы еще быстрее, если бы не узкое место - копирование больших объемов памяти между CPU и GPU.

параллелить я пока оставляю в стороне, да и в распоряжении есть компьютер с одним ядром.

не пойму почему вот это не проходит

int ii =  b_row;
for(i=0; i<n_row-1; ++i) {
       xip0=*(ii);
       xip1=*(ii+1);
       xin0=*(ii+2);
       xin1=*(ii+3);
       ii += 2;
...

ну и, надеюсь правильно считаю что подобный код быстрее кода самого первого поста.

Смотря что ты там делаешь. В для обычного O(n^2) алгоритма работы с матрицей (примеры: умножение матрицы на вектор, обращение верхнетреугольной матрицы ...) ускорение будет только если матрица не слишком маленькая (чтобы накладные расходы на создание потоков ничего не съели) и в то же время не слишком большая, чтобы _полностью_ влезать в кеш. Если во внутреннем цикле делается действительно какая-то сложная работа над элементами матрицы (вычисление тригонометрических функций и прочее), то ситуация будет на много лучше.

Потому что такую тривиальщину компилятор сам в состоянии выполнить. Можно вообще забить на это всё и писать везде просто и понятно a[i*n+j] ... и всё равно будет также по скорости.

Можно вообще забить на это всё и писать везде просто и понятно a[i*n+j] ... и всё равно будет также по скорости.

Это стоит проверить, у меня он такое место не соптимизировал. Но может я что-нить в ключах компилеру не указал.

Проверялось неоднократно — на линейной алгебре одни фиг выходит.

Это точно. Иногда даже выделение временной переменной для «ускорения» на самом деле тормозит вычисления. Сам только вчера сидел, мучился с оптимизацией openmp. А еще бывает, что параллельное вычисление по скорости вообще не дает выигрыша. Например, код

	GLfloat *normal = norms->data;
	for(y = 0; y < H; y++){
		ptr = &tr_norms[Wstride*(y+1)+1];
		for(x = 0; x < W; x++, normal+=3, ptr+=6){
			addVectors(normal, ptr);
			normalize(normal);
		}
	}

	int WW = W*3;
	#pragma omp parallel for private(y, ptr, normal)
	for(y = 0; y < H; y++){
		ptr = &tr_norms[Wstride*(y+1)+1];
		normal = &norms->data[y*WW];
		for(x = 0; x < W; x++, normal+=3, ptr+=6){
			addVectors(normal, ptr);
			normalize(normal);
		}
	}

Если в addVectors(normal, ptr); normalize(normal); делается что-то сложное, что приводит к сложности O(W*H*W), то надо курить про блочные алгоритмы (поможет книжка Деммеля). Если там делается что-то простое, что не приводит к увеличению сложности, то ничего сделать нельзя.

Простенькие арифметические функции. Но все-таки, теоретически параллелизация должна была бы ускорить вычисления, ан нет...

Ага, в теории на идеальном процессоре с бесконечным объемом кеша. Читай Деммеля, у него подробно описано почему так работает на реальном железе.

> Простенькие арифметические функции. Но все-таки, теоретически параллелизация должна была бы ускорить вычисления

С чего бы это? Параллеизация увеличивает нагрузку на подсистему памяти и добавляет накладных расходов на всякие синхронизации. Поэтому тривиальных задачах иногда бывает, что получается замедление в число нитей...

Деммель Дж. Вычислительная линейная алгебра. Теория и приложения

Эту?

да

пусть процессор с компилятором вычисляет адреса :) а программер тем временем должен писать ясный код

я бы делал примерно так:

...
// create temporary in stack (or in hip)
double **matr=(double **)alloca(sizeof(double *)*n_cols;
// place rows
for(i=0;i<n_rows;c++)
  matr[i]=matrix+n_cols;
// вот оно щастье :
...
matr[y][x]=15; // теперь нормально обращаемся к элементам
matr[3][3]=matr[1][1];  // и не ипём моск 

У тебя первая размерность matr неизвестна, поэтому делать matr[x][y] нельзя.

n_row и n_col вроде как переданны в аргументах. Насколько понял матрица получается n_row*n_col - всё известно.

Компилятору неизвестно, мужик.

а зачем это ему ? Вы С с Fortran` не путаете ?

математическая матрица задана как вектор (массив заданной размерности) --- это дано (пишу процедуру для numpy)

Как, по твоему, компилятор вычислит адрес matr[2][3], зная только то что mattr — это double **?

А что там про Фортран?

Как, по твоему, компилятор вычислит адрес matr[2][3]

хер знает, но код скомпилелся. Нет времени над ним думать и проверять valgrind-ом.

Сегфолтится, причём в предсказуемом месте :). Короче, код нуждается в доработке :). Ну а адрес оно, похоже, как base+x*y определило, но я не вручаюсь.

Оно пытается получить доступ к *(*(matr + x) + y).

Надо указывать размерность, ну.

double matr[n_rows][n_cols];

> void solve_matrix( const int n_row, const double* b_row, const int n_col, const double* b_col, ..., double* matrix)

inline void solve_matrix( const int n_row, const double* b_row, const int n_col, const double* b_col, ..., double* matrix)

я думал что это дело применимо в случае

int simply_funct(params){ return exspression; }

А что там про Фортран?

насколько помню в нём положение элемента вычисляется именно как ТС делает, то есть на основе размерностей, а вот в С немного совсем по другому.

вобщем почитайте K&R про указатели и массивы в С.

то что я написал, содержит ошибки (даже различимые на глаз), повнимательнее надо. Общая идея - делаете вектор double **matr нужного размера (по строкам или по столбцам), инициализуете его адресами соответсвующих строк(или столбцов) из плоского предствления matrix - и с этого момента работаете с естественным для С представлением (через указатель+индекс). И компилер сумеет оптимизировать эти обращения вместе с кодом который с ними работает, и себя не парите лишними вычислениями и кеш процессору не забивается ненужной арифметикой.

//inline
uw_t
stec_div_by_x(uw_t *aR, const uw_t *aX, const uw_t x, const sz_t N){
  const uw_t c_x = (STEC_MAX_UW - (STEC_MAX_UW % x)) / x; // max(rcx) == STEC_MAX_UW * (x-1)/x < STEC_MAX_UW
  const uw_t d_x = (1+(STEC_MAX_UW % x));// max(d_x) + 1 <= x
  sz_t i  = N;
  uw_t r = 0, rem, rcx = 0;

  assert(x <= g_L); // (x-1)*d_x + (x-1) <= (x-1)*x + (x-1) = (x-1)*(x+1) = x*x - 1, g_L*g_L - 1 < STEC_MAX_UW
                    //  --> we need x <= g_L (g_L = STEC_MAX_UW >> sizeof(uw_t)*4;)
  //aR += N; aX += N;
  do{
    i--;// aR--; aX--;
    /* uw_t rcx = r * c_x; */
    *(aR + i)  = *(aX + i) / x;
    rem        = *(aX + i) % x;
    *(aR + i) += rcx; //
    // STEC_MAX_UW/x +  (x-1)*(STEC_MAX_UW /x) <= STEC_MAX_UW*(1+x-1)/x == STEC_MAX_UW  <= STEC_MAX_UW
    r *= d_x; // r*d_x + rem < STEC_MAX_UW (see assertion above)
    r += rem;
    if(r >= x){
      rem = r / x;
      r   = r % x;
      *(aR + i) += rem; // STEC_MAX_UW + (x-1) - worst case (more precise calculations show the worst case is STEC_MAX_UW + 1)
      if(*(aR + i) < rem)
         stec_inc((aR + i) + 1, 1, N - i - 1); // next optimisation is to replace stec_inc with (*(aR + i + 1))++; but some proof is necessary
    }
    rcx = r * c_x;
  }while(i);

  return r;
}

Екземпл. Работает примерно в 4 раза медленнее без инлайна.

И работать будет медленнее

Теперь понял, я не до конца въехал что происходит.

Почитать нужно Вам, уважаемый. А ошибку свою так и не увидели. Сие печально.

> что здесь не так?

double ptr_j = b_col;

vs

double* ptr_j = b_col;

Это же не рабочий код?

И да,

for(i=0; i<n_row-1; ++i)

for(j=0; j<n_col-1; ++j)

должны быть без "-1", если, конечно, ты всё правильно делаешь.

Это же нормальное действие ptr_j=x_j? Можно как-то в программе указать компилятору (ключи отключающие предупреждения передавать компилятору не хочу), что здесь не нужно предупреждать (если конечно код нормальный :-) )?

void matrix_maker_for_test( const int n_row, const double* x_row, const int n_col, const double* x_col, ... , double* mtx) 
{
  double ti0, ti1, ... tj0, tj1 , ... ;
  for(i=0; i<n_row; ++i) {
       ti0 = *(x_row++);
       ti1 = *(x_row++);
       ...
       ptr_j = x_col;
       for(j=0; j<n_col-1; ++j){
       tj0 = *(ptr_j++);
       tj1 = *(ptr_j++);
       tj2 = *(ptr_j++);
...

компилятор на строку ptr_j = x_col выдает предупреждение

 warning: assignment discards qualifiers from pointer target type